Page MenuHomeHEPForge
Files F8308730

No OneTemporary
Actions

View Options

Index: trunk/examples/GenFitDs2KKpi.cc
===================================================================
--- trunk/examples/GenFitDs2KKpi.cc (revision 415)
+++ trunk/examples/GenFitDs2KKpi.cc (revision 416)
@@ -1,169 +1,169 @@
#include <cstdlib>
#include <iostream>
#include <vector>
#include "TFile.h"
#include "TH2.h"
#include "TString.h"
#include "TTree.h"
#include "LauSimpleFitModel.hh"
#include "LauBkgndDPModel.hh"
#include "LauDaughters.hh"
#include "LauEffModel.hh"
#include "LauIsobarDynamics.hh"
#include "LauMagPhaseCoeffSet.hh"
#include "LauVetoes.hh"
void usage( std::ostream& out, const TString& progName )
{
out<<"Usage:\n";
out<<progName<<" gen [nExpt = 1] [firstExpt = 0]\n";
out<<"or\n";
out<<progName<<" fit <iFit> [nExpt = 1] [firstExpt = 0]"<<std::endl;
}
int main( int argc, char** argv )
{
// Process command-line arguments
// Usage:
// ./GenFitDs2KKpi gen [nExpt = 1] [firstExpt = 0]
// or
// ./GenFitDs2KKpi fit <iFit> [nExpt = 1] [firstExpt = 0]
if ( argc < 2 ) {
usage( std::cerr, argv[0] );
return EXIT_FAILURE;
}
TString command = argv[1];
command.ToLower();
Int_t iFit(0);
Int_t nExpt(1);
Int_t firstExpt(0);
if ( command == "gen" ) {
if ( argc > 2 ) {
nExpt = atoi( argv[2] );
if ( argc > 3 ) {
firstExpt = atoi( argv[3] );
}
}
} else if ( command == "fit" ) {
if ( argc < 3 ) {
usage( std::cerr, argv[0] );
return EXIT_FAILURE;
}
iFit = atoi( argv[2] );
if ( argc > 3 ) {
nExpt = atoi( argv[3] );
if ( argc > 4 ) {
firstExpt = atoi( argv[4] );
}
}
} else {
usage( std::cerr, argv[0] );
return EXIT_FAILURE;
}
// If you want to use square DP histograms for efficiency,
// backgrounds or you just want the square DP co-ordinates
// stored in the toy MC ntuple then set this to kTRUE
Bool_t squareDP = kFALSE;
// This defines the DP => decay is D_s+ -> pi+ K+ K-
// Particle 1 = pi+
// Particle 2 = K+
// Particle 3 = K-
// The DP is defined in terms of m13Sq and m23Sq
LauDaughters* daughters = new LauDaughters("D_s+", "pi+", "K+", "K-", squareDP);
// Define the efficiency model (defaults to unity everywhere)
LauVetoes* vetoes = new LauVetoes();
LauEffModel* effModel = new LauEffModel(daughters, vetoes);
// Create the isobar model
LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel);
- LauAbsResonance* reson(0);
- reson = sigModel->addResonance("phi(1020)", 1, LauAbsResonance::RelBW); // resPairAmpInt = 1 => resonance mass is m23.
- reson = sigModel->addResonance("K*0(892)", 2, LauAbsResonance::RelBW);
- reson = sigModel->addResonance("NonReson", 0, LauAbsResonance::FlatNR);
+ //LauAbsResonance* reson(0);
+ /*reson =*/ sigModel->addResonance("phi(1020)", 1, LauAbsResonance::RelBW); // resPairAmpInt = 1 => resonance mass is m23.
+ /*reson =*/ sigModel->addResonance("K*0(892)", 2, LauAbsResonance::RelBW);
+ /*reson =*/ sigModel->addResonance("NonReson", 0, LauAbsResonance::FlatNR);
// Reset the maximum signal DP ASq value
// This will be automatically adjusted to avoid bias or extreme
// inefficiency if you get the value wrong but best to set this by
// hand once you've found the right value through some trial and error.
sigModel->setASqMaxValue(0.27);
// Create the fit model
LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel);
// Create the complex coefficients for the isobar model
// Here we're using the magnitude and phase form:
// c_j = a_j exp(i*delta_j)
std::vector<LauAbsCoeffSet*> coeffset;
coeffset.push_back( new LauMagPhaseCoeffSet("phi(1020)", 1.00, 0.00, kTRUE, kTRUE) );
coeffset.push_back( new LauMagPhaseCoeffSet("K*0(892)", 1.00, 0.00, kFALSE, kFALSE) );
coeffset.push_back( new LauMagPhaseCoeffSet("NonReson", 1.00, 0.00, kFALSE, kFALSE) );
for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
fitModel->setAmpCoeffSet(*iter);
}
// Set the signal yield and define whether it is fixed or floated
const Double_t nSigEvents = 10000.0;
LauParameter * signalEvents = new LauParameter("signalEvents",nSigEvents,-1.0*nSigEvents,2.0*nSigEvents,kFALSE);
fitModel->setNSigEvents(signalEvents);
// Set the number of experiments to generate or fit and which
// experiment to start with
fitModel->setNExpts( nExpt, firstExpt );
// Configure various fit options
// Switch on/off calculation of asymmetric errors.
fitModel->useAsymmFitErrors(kFALSE);
// Randomise initial fit values for the signal mode
fitModel->useRandomInitFitPars(kTRUE);
const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 );
// Switch on/off Poissonian smearing of total number of events
fitModel->doPoissonSmearing(haveBkgnds);
// Switch on/off Extended ML Fit option
fitModel->doEMLFit(haveBkgnds);
// Generate toy from the fitted parameters
//TString fitToyFileName("fitToyMC_Ds2KKpi_");
//fitToyFileName += iFit;
//fitToyFileName += ".root";
//fitModel->compareFitData(100, fitToyFileName);
// Write out per-event likelihoods and sWeights
//TString splotFileName("splot_Ds2KKpi_");
//splotFileName += iFit;
//splotFileName += ".root";
//fitModel->writeSPlotData(splotFileName, "splot", kFALSE);
// Set the names of the files to read/write
TString dataFile("gen-Ds2KKpi.root");
TString treeName("genResults");
TString rootFileName("");
TString tableFileName("");
if (command == "fit") {
rootFileName = "fitDs2KKpi_"; rootFileName += iFit;
rootFileName += "_expt_"; rootFileName += firstExpt;
rootFileName += "-"; rootFileName += (firstExpt+nExpt-1);
rootFileName += ".root";
tableFileName = "fitDs2KKpiResults_"; tableFileName += iFit;
} else {
rootFileName = "dummy.root";
tableFileName = "genDs2KKpiResults";
}
// Execute the generation/fit
fitModel->run( command, dataFile, treeName, rootFileName, tableFileName );
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:59 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022865
Default Alt Text
(5 KB)

Event Timeline