FairRoot/PandaRoot
macro/softrig/quickfsimana.C
Go to the documentation of this file.
1 // **********************************************************************************************
2 // Macro for running fast simulation + Software Trigger + simple analysis (PndSimpleCombinerTask)
3 // **********************************************************************************************
4 
5 // The parameters are
6 // -------------------
7 // quickfsimana.C( <pref>, <decfile>, <mom>, <decay>, [nevt], [res], [parms] )
8 // <pref> : output file names prefix
9 // <decfile> : EvtGen decfile; DPM/FTF/BOX uses DPM/FTF generator (inelastic mode) or box generator instead
10 // <mom> : pbar momentum; negative values are interpreted as -E_cm
11 // [decay] : the decay pattern to be reconstructed, e.g. 'phi -> K+ K-; D_s+ -> phi pi+'; emtpy string: only fast sim will be run
12 // [nevt] : number of events; default = 1000
13 // [res] : resonance/particle type for BOX generator (ignored when running DPM); default = 'pbarpSystem0'
14 // [parms] : parameters for the analysis, e.g. 'mwin=0.4:mwin(phi)=0.1:emin=0.1:pmin=0.1:qamc'
15 // [runST] : if 'true' runs Software Trigger (default: false)
16 // [runnum] : integer run number (default: 0)
17 // [mode] : arbitrary mode number (default: 0)
18 // -------------------
19 
20 
21 int quickfsimana(TString Prefix="", TString Decfile="", Float_t Mom=0., TString anadecay="",
22  Int_t nEvents = 1000, TString Resonance="pbarpSystem0", TString anaparms="", bool runST=false, int run=0 , int mode=0)
23 {
24  if (Prefix=="" || Decfile=="" || Mom==0. )
25  {
26  cout << "USAGE:\n";
27  cout << "quickfsimana.C( <pref>, <decfile>, <mom>, <decay>, [nevt], [res], [parms] )\n\n";
28  cout << " <pref> : output file names prefix\n";
29  cout << " <decfile> : EvtGen decfile; DPM/FTF/BOX uses DPM/FTF generator (inelastic mode) or box generator instead\n";
30  cout << " <mom> : pbar momentum; negative values are interpreted as -E_cm\n";
31  cout << " [decay] : the decay pattern to be reconstructed, e.g. 'phi -> K+ K-; D_s+ -> phi pi+'; emtpy string: only fast sim will be run\n";
32  cout << " [nevt] : number of events; default = 1000\n";
33  cout << " [res] : resonance/particle type for BOX generator (ignored when running DPM); default = 'pbarpSystem0'\n";
34  cout << " [parms] : parameters for the analysis, e.g. 'mwin=0.4:mwin(phi)=0.1:emin=0.1:pmin=0.1:qamc'\n";
35  cout << " [runST] : if 'true' runs Software Trigger (default: false)\n";
36  cout << " [runnum] : integer run number (default: 0)\n";
37  cout << " [mode] : arbitrary mode number (default: 0)\n\n";
38  return;
39  }
40 
41  // only run fast sim without analysis
42  bool simonly = (anadecay=="");
43 
44  // for submission to queue all blanks in decay string were replaced by '§'; now we replace again the other way around
45  anadecay.ReplaceAll("§," "); // if Mom<0, interprete as -E_cm double mp = 0.938272; if (Mom<0) { double X = (Mom*Mom-2*mp*mp)/(2*mp); Mom = sqrt(X*X-mp*mp); } // Allow shortcut for resonance if (Resonance=="pbp") Resonance = "pbarpSystem"; if (Resonance=="pbp0") Resonance = "pbarpSystem0"; anadecay.ReplaceAll("pbp","pbarpSystem"); anadecay.ReplaceAll("pbp0","pbarpSystem0"); // Prevent generator from throwing a lot of warnings TLorentzVector fIni(0,0,Mom,mp+sqrt(Mom*Mom+mp*mp)); TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",fIni.M(),kFALSE,0.1,0, "",88888,0); TDatabasePDG::Instance()->AddParticle("pbarpSystem0","pbarpSystem0",fIni.M(),kFALSE,0.1,0, "",88880,0); TDatabasePDG::Instance()->AddParticle("Z(3900)+","Z+",3.900,kFALSE,0.03,0, "",90000); TDatabasePDG::Instance()->AddParticle("Z(3900)-","Z-",3.900,kFALSE,0.03,0, "",-90000); //----- Switches for Fast Simulation Options ------------------------------ Bool_t enableSplitoff = true; // create e.-m. and hadronic split offs Bool_t mergeNeutrals = true; // merge neutrals (for merged pi0s) Bool_t electronBrems = true; // bremsstrahlung loss for electrons Bool_t useEventFilter = false; // enable Fast Sim event filter. *** Needs configuration (see below) *** Bool_t usePndEventFilter = false; // enable Panda event filter. *** Needs configuration (see below) *** //----- Presist simulation output ------------------------------ Bool_t persist = simonly; // if analysis is running, fsim output not needed //-----General settings----------------------------------------------- TString BaseDir = gSystem->Getenv("VMCWORKDIR"); TString splitpars = BaseDir+"/fsim/splitpars.dat"; gRandom->SetSeed(); //-----User Settings:------------------------------------------------- TString OutputFile = TString::Format("%s_%d_ana.root",Prefix.Data(), run); if (simonly) OutputFile = TString::Format("%s_%d_fsim.root",Prefix.Data(), run); gDebug = 0; // choose your event generator Bool_t UseEvtGenDirect = kTRUE; Bool_t UseFtf = kFALSE; Bool_t UseDpm = kFALSE; Bool_t UseBoxGenerator = kFALSE; // use DPM generator; default: inelastic @ pbarmom = mom if (Decfile.BeginsWith("DPM")) { UseEvtGenDirect = kFALSE; UseDpm = kTRUE; } // use FTF generator; if (Decfile.BeginsWith("FTF")) { UseEvtGenDirect = kFALSE; UseFtf = kTRUE; } // use BOX generator; default: single mu-, 0<tht<180, 0<phi<360, 0.1<p<mom if (Decfile=="BOX") { UseEvtGenDirect = kFALSE; UseBoxGenerator = kTRUE; } usePndEventFilter=UseDpm; Double_t MomMin = 0.1; // minimum momentum for box generator Double_t MomMax = Mom; // maximum " " // Start a stop watch TStopwatch timer; timer.Start(); // Create the Simulation run manager // -------------------------------- FairRunSim *fRun = new FairRunSim(); fRun->SetOutputFile(OutputFile.Data()); fRun->SetGenerateRunInfo(kFALSE); if (!simonly) fRun->SetUserConfig(BaseDir+"/tutorials/analysis/g3ConfigNoMC.C"); FairLogger::GetLogger()->SetLogToFile(kFALSE); // Create and Set Event Generator // ------------------------------- FairFilteredPrimaryGenerator* primGen = new FairFilteredPrimaryGenerator(); //FairPrimaryGenerator* primGen = new FairPrimaryGenerator(); fRun->SetGenerator(primGen); fRun->SetName("TGeant3"); if(UseBoxGenerator) { // Box Generator int Pdgcode = TDatabasePDG::Instance()->GetParticle(Resonance)->PdgCode(); FairBoxGenerator* boxGen = new FairBoxGenerator(Pdgcode, 1); // 211 = pion; 1 = multipl. boxGen->SetPRange(MomMin,MomMax); // GeV/c boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree] boxGen->SetThetaRange(0., 180.); // Polar angle in lab system range [degree] boxGen->SetXYZ(0., 0., 0.); //cm primGen->AddGenerator(boxGen); } if(UseDpm) { int mode = 0; if (Decfile=="DPM1") mode = 1; if (Decfile=="DPM2") mode = 2; PndDpmDirect *Dpm= new PndDpmDirect(Mom,mode); // 0 = inelastic, 1 = inelastic & elastic, 2 = elastic // since fastsim doesn't have a transport, let all long-living resonances decay by the generator Dpm->SetUnstable(111); // pi0 Dpm->SetUnstable(310); // K_S0 Dpm->SetUnstable(311); // K0 Dpm->SetUnstable(-311); // K0bar Dpm->SetUnstable(3122); // Lambda0 Dpm->SetUnstable(-3122); // anti-Lambda0 Dpm->SetUnstable(221); // eta*/ Dpm->SetUnstable(3222); // Sigma+ Dpm->SetUnstable(-3222); // anti-Sigma- Dpm->SetUnstable(3334); // Omega- primGen->AddGenerator(Dpm); } if(UseFtf) { bool noelastic = true; if (Decfile=="FTF1") noelastic=false; PndFtfDirect *Ftf = new PndFtfDirect("anti_proton", "G4_H", 1, "ftfp", Mom, 0, noelastic); primGen->AddGenerator(Ftf); } if(UseEvtGenDirect) { PndEvtGenDirect *EvtGen = new PndEvtGenDirect(Resonance, Decfile.Data(), Mom); EvtGen->SetStoreTree(kTRUE); primGen->AddGenerator(EvtGen); } // ------------- switch off the transport of particles primGen->DoTracking(kFALSE); //---------------------Create and Set the Field(s)---------- PndMultiField *fField= new PndMultiField("AUTO"); fRun->SetField(fField); // ********************************** // *** BEGIN Fast Simulation Task *** // ********************************** PndFastSim* fastSim = new PndFastSim(persist); // increasing verbosity increases the amount of console output (mainly for debugging) fastSim->SetVerbosity(0); // set PANDA event filters //----------------------------- if (usePndEventFilter) { // primGen->SetFilterMaxTries(100000); // for testing small number, for real produrction set usually to 9999999 or something very big // FairEvtFilterOnSingleParticleCounts* lambfilt= new FairEvtFilterOnSingleParticleCounts("PdgFilter"); // lambfilt->AndMaxPdgCodes(0, 3122, -3122); // filter out Lambda0 // lambfilt->AndMaxPdgCodes(0, 3222, -3222); // filter out Sigma+ // primGen->AndFilter(lambfilt); } // set event filters //----------------------------- if (useEventFilter) { // fastSim->SetMultFilter(type, min, max); requires min <= mult <= max // available types are: '+', '-', 'gam', 'pi0' (2 gam, 0.135 +- 0.03 GeV), 'eta' (2 gam, 0.547 +- 0.04 GeV), 'ks' (pi+ pi-, 0.497 +- 0.04 GeV) // fastSim->SetInvMassFilter(comb, m_min, m_max, mult); requires at least mult combined candidates with m_min < m < m_max // comb is a TString describing the combinatoric // - particle codes are: e+ e- mu+ mu- pi+ pi- k+ k- p+ p- gam pi0 ks eta, separated with blank // Examples: // - ("k+ k-", 0.98, 1.1, 2) : forms K+ K- candidate and requires >=2 in the given window // - ("ks k+ pi- cc", 2.8, 3.2,1 ) : forms ks k+ pi- / ks k- pi+ cands and req. at least one in window } // enable the merging of neutrals if they have similar direction fastSim->MergeNeutralClusters(mergeNeutrals); // enable bremsstahlung loss for electrons fastSim->EnableElectronBremsstrahlung(electronBrems); //enable the producting of parametrized neutral (hadronic) split offs // generate electro-magnetic / hadronic split offs in the EMC? switch off when running w/o EMC if (enableSplitoff) fastSim->EnableSplitoffs(splitpars.Data()); // simulates covariance matrix for fitting fastSim->SetUseFlatCov(true); // ----------------------------------------------------------------------------------- // Now add the detector components // ----------------------------------------------------------------------------------- // Full Panda Tracking: STT MVD GEM FTS (scanning over theta range from 0° to 160°) fastSim->AddDetector("ScFts", "thtMin=0. thtMax=5. ptmin=0.0 pmin=0.5 pRes=0.05 thtRes=0.002 phiRes=0.002 efficiency=0.80"); fastSim->AddDetector("ScMvdGem", "thtMin=5. thtMax=7.8 ptmin=0.1 pmin=0.0 pRes=0.03 thtRes=0.001 phiRes=0.001 efficiency=0.60"); fastSim->AddDetector("ScSttMvdGem", "thtMin=7.8 thtMax=20.9 ptmin=0.1 pmin=0.0 pRes=0.02 thtRes=0.001 phiRes=0.001 efficiency=0.85"); fastSim->AddDetector("ScSttMvd", "thtMin=20.9 thtMax=145. ptmin=0.1 pmin=0.0 pRes=0.02 thtRes=0.001 phiRes=0.001 efficiency=0.85"); fastSim->AddDetector("ScSttAlone", "thtMin=145. thtMax=159.5 ptmin=0.1 pmin=0.0 pRes=0.04 thtRes=0.006 phiRes=0.007 efficiency=0.25"); // Vertexing fastSim->AddDetector("ScVtxNoMvd", "thtMin=0. thtMax=5. ptmin=0.0 vtxRes=0.05 efficiency=1."); // efficiency=1: all tracks found in trackers will get a vertex information fastSim->AddDetector("ScVtxMvd", "thtMin=5. thtMax=145. ptmin=0.1 vtxRes=0.005 efficiency=1."); // efficiency=1: all tracks found in trackers will get a vertex information // EM Calorimeters w/ default parameters fastSim->AddDetector("EmcFS", "thtMin=0.05 thtMax=10.0 aPar=0.013 bPar=0.0283 Emin=0.01 dist=8.2"); fastSim->AddDetector("EmcFwCap", "thtMin=10.0 thtMax=22.0 Emin=0.01 dist=2.5"); fastSim->AddDetector("EmcBarrel", "thtMin=22.0 thtMax=142.0 Emin=0.01 barrelRadius=0.5"); fastSim->AddDetector("EmcBwCap", "thtMin=142.0 thtMax=160.0 Emin=0.01 dist=0.7"); // PID Target Spectrometer: STT, MVD, MDT, EMC, Barrel DIRC, Disc DIRC fastSim->AddDetector("MvdPid", "thtMin=5. thtMax=133.6 ptmin=0.1 dEdxResMulti=1. efficiency=1."); // Note: Bethe-Bloch-Landau-Gauss Prametrization from 2008 fastSim->AddDetector("SttPid", "thtMin=7.8 thtMax=159.5 ptmin=0.1 dEdxRes=0.15 efficiency=1."); // Note: dEdX parametrization from 2008 fastSim->AddDetector("ScMdtPidBarrel", "thtMin=10.0 thtMax=130.0 pmin=0.5 efficiency=0.95 misId=0.01"); fastSim->AddDetector("ScEmcPidFwCap", "thtMin=10.0 thtMax=22.0 ptmin=0.0 pmin=0.0 efficiency=1.0"); fastSim->AddDetector("ScEmcPidBarrel", "thtMin=22.0 thtMax=142.0 ptmin=0.2 pmin=0.0 efficiency=1.0"); fastSim->AddDetector("ScEmcPidBwCap", "thtMin=142.0 thtMax=160.0 ptmin=0.0 pmin=0.0 efficiency=1.0"); fastSim->AddDetector("DrcDisc", "thtMin=5.0 thtMax=22.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075"); fastSim->AddDetector("DrcBarrel", "thtMin=22.0 thtMax=140.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075"); // PID Fwd spectrometer PID: RICH, FwdMUO and EMC FS fastSim->AddDetector("ScEmcPidFS", "thtMin=0.5 thtMax=10.0 ptmin=0.0 pmin=0.5 efficiency=1.0"); fastSim->AddDetector("ScMdtPidForward", "thtMin=0.0 thtMax=10.0 pmin=0.5 efficiency=0.95 misId=0.01"); fastSim->AddDetector("Rich", "angleXMax=5.0 angleYMax=10.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075"); fRun->AddTask(fastSim); // ******************************** // *** END Fast Simulation Task *** // ******************************** // ***************************** // *** BEGIN SoftTriggerTask *** // ***************************** if (runST) { // this file contains the trigger line definitions TString triggercfg = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/triggerlines_fsim.cfg"; PndSoftTriggerTask *stTask = new PndSoftTriggerTask(Mom, mode, run, triggercfg); stTask->SetFastSimDefaults(); fRun->AddTask(stTask); } // *************************** // *** END SoftTriggerTask *** // *************************** // *********************************** // *** BEGIN PndSimpleCombinerTask *** // *********************************** if (!simonly) { PndSimpleCombinerTask *scTask = new PndSimpleCombinerTask(anadecay, anaparms+":algo=PidChargedProbability",Mom, run, mode); scTask->SetPidAlgo("PidChargedProbability"); fRun->AddTask(scTask); } // ********************************* // *** END PndSimpleCombinerTask *** // ********************************* //------------------------- Initialize the RUN ----------------- fRun->Init(); //------------------------- Run the Simulation ----------------- fRun->Run(nEvents); //------------------------- Write Filter Info to File ----------- if (usePndEventFilter) primGen->WriteEvtFilterStatsToRootFile(); //------------------------Print some info and exit---------------- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime); return 0; } "," ");
46 
47  // if Mom<0, interprete as -E_cm
48  double mp = 0.938272;
49  if (Mom<0)
50  {
51  double X = (Mom*Mom-2*mp*mp)/(2*mp);
52  Mom = sqrt(X*X-mp*mp);
53  }
54 
55  // Allow shortcut for resonance
56  if (Resonance=="pbp") Resonance = "pbarpSystem";
57  if (Resonance=="pbp0") Resonance = "pbarpSystem0";
58 
59  anadecay.ReplaceAll("pbp","pbarpSystem");
60  anadecay.ReplaceAll("pbp0","pbarpSystem0");
61 
62  // Prevent generator from throwing a lot of warnings
63  TLorentzVector fIni(0,0,Mom,mp+sqrt(Mom*Mom+mp*mp));
64  TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",fIni.M(),kFALSE,0.1,0, "",88888,0);
65  TDatabasePDG::Instance()->AddParticle("pbarpSystem0","pbarpSystem0",fIni.M(),kFALSE,0.1,0, "",88880,0);
66  TDatabasePDG::Instance()->AddParticle("Z(3900)+","Z+",3.900,kFALSE,0.03,0, "",90000);
67  TDatabasePDG::Instance()->AddParticle("Z(3900)-","Z-",3.900,kFALSE,0.03,0, "",-90000);
68 
69  //----- Switches for Fast Simulation Options ------------------------------
70  Bool_t enableSplitoff = true; // create e.-m. and hadronic split offs
71  Bool_t mergeNeutrals = true; // merge neutrals (for merged pi0s)
72  Bool_t electronBrems = true; // bremsstrahlung loss for electrons
73  Bool_t useEventFilter = false; // enable Fast Sim event filter. *** Needs configuration (see below) ***
74  Bool_t usePndEventFilter = false; // enable Panda event filter. *** Needs configuration (see below) ***
75 
76  //----- Presist simulation output ------------------------------
77  Bool_t persist = simonly; // if analysis is running, fsim output not needed
78 
79  //-----General settings-----------------------------------------------
80  TString BaseDir = gSystem->Getenv("VMCWORKDIR");
81  TString splitpars = BaseDir+"/fsim/splitpars.dat";
82  gRandom->SetSeed();
83 
84  //-----User Settings:-------------------------------------------------
85  TString OutputFile = TString::Format("%s_%d_ana.root",Prefix.Data(), run);
86  if (simonly) OutputFile = TString::Format("%s_%d_fsim.root",Prefix.Data(), run);
87 
88  gDebug = 0;
89 
90  // choose your event generator
91  Bool_t UseEvtGenDirect = kTRUE;
92  Bool_t UseFtf = kFALSE;
93  Bool_t UseDpm = kFALSE;
94  Bool_t UseBoxGenerator = kFALSE;
95 
96  // use DPM generator; default: inelastic @ pbarmom = mom
97  if (Decfile.BeginsWith("DPM"))
98  {
99  UseEvtGenDirect = kFALSE;
100  UseDpm = kTRUE;
101  }
102 
103  // use FTF generator;
104  if (Decfile.BeginsWith("FTF"))
105  {
106  UseEvtGenDirect = kFALSE;
107  UseFtf = kTRUE;
108  }
109 
110  // use BOX generator; default: single mu-, 0<tht<180, 0<phi<360, 0.1<p<mom
111  if (Decfile=="BOX")
112  {
113  UseEvtGenDirect = kFALSE;
114  UseBoxGenerator = kTRUE;
115  }
116 
117  usePndEventFilter=UseDpm;
118 
119  Double_t MomMin = 0.1; // minimum momentum for box generator
120  Double_t MomMax = Mom; // maximum " "
121 
122  // Start a stop watch
123  TStopwatch timer;
124  timer.Start();
125 
126  // Create the Simulation run manager
127  // --------------------------------
128  FairRunSim *fRun = new FairRunSim();
129  fRun->SetOutputFile(OutputFile.Data());
130  fRun->SetGenerateRunInfo(kFALSE);
131  if (!simonly) fRun->SetUserConfig(BaseDir+"/tutorials/analysis/g3ConfigNoMC.C");
132 
133  FairLogger::GetLogger()->SetLogToFile(kFALSE);
134 
135 
136  // Create and Set Event Generator
137  // -------------------------------
139  //FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
140  fRun->SetGenerator(primGen);
141  fRun->SetName("TGeant3");
142 
143  if(UseBoxGenerator)
144  { // Box Generator
145  int Pdgcode = TDatabasePDG::Instance()->GetParticle(Resonance)->PdgCode();
146  FairBoxGenerator* boxGen = new FairBoxGenerator(Pdgcode, 1); // 211 = pion; 1 = multipl.
147  boxGen->SetPRange(MomMin,MomMax); // GeV/c
148  boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree]
149  boxGen->SetThetaRange(0., 180.); // Polar angle in lab system range [degree]
150  boxGen->SetXYZ(0., 0., 0.); //cm
151  primGen->AddGenerator(boxGen);
152  }
153 
154  if(UseDpm)
155  {
156  int mode = 0;
157  if (Decfile=="DPM1") mode = 1;
158  if (Decfile=="DPM2") mode = 2;
159 
160  PndDpmDirect *Dpm= new PndDpmDirect(Mom,mode); // 0 = inelastic, 1 = inelastic & elastic, 2 = elastic
161  // since fastsim doesn't have a transport, let all long-living resonances decay by the generator
162  Dpm->SetUnstable(111); // pi0
163  Dpm->SetUnstable(310); // K_S0
164  Dpm->SetUnstable(311); // K0
165  Dpm->SetUnstable(-311); // K0bar
166  Dpm->SetUnstable(3122); // Lambda0
167  Dpm->SetUnstable(-3122); // anti-Lambda0
168  Dpm->SetUnstable(221); // eta*/
169  Dpm->SetUnstable(3222); // Sigma+
170  Dpm->SetUnstable(-3222); // anti-Sigma-
171  Dpm->SetUnstable(3334); // Omega-
172  primGen->AddGenerator(Dpm);
173  }
174 
175  if(UseFtf)
176  {
177  bool noelastic = true;
178  if (Decfile=="FTF1") noelastic=false;
179  PndFtfDirect *Ftf = new PndFtfDirect("anti_proton", "G4_H", 1, "ftfp", Mom, 0, noelastic);
180  primGen->AddGenerator(Ftf);
181  }
182 
183  if(UseEvtGenDirect)
184  {
185  PndEvtGenDirect *EvtGen = new PndEvtGenDirect(Resonance, Decfile.Data(), Mom);
186  EvtGen->SetStoreTree(kTRUE);
187  primGen->AddGenerator(EvtGen);
188  }
189 
190  // ------------- switch off the transport of particles
191  primGen->DoTracking(kFALSE);
192 
193  //---------------------Create and Set the Field(s)----------
194  PndMultiField *fField= new PndMultiField("AUTO");
195  fRun->SetField(fField);
196 
197 
198 
199  // **********************************
200  // *** BEGIN Fast Simulation Task ***
201  // **********************************
202 
203  PndFastSim* fastSim = new PndFastSim(persist);
204 
205  // increasing verbosity increases the amount of console output (mainly for debugging)
206  fastSim->SetVerbosity(0);
207 
208  // set PANDA event filters
209  //-----------------------------
210  if (usePndEventFilter)
211  {
212  // primGen->SetFilterMaxTries(100000); // for testing small number, for real produrction set usually to 9999999 or something very big
213  // FairEvtFilterOnSingleParticleCounts* lambfilt= new FairEvtFilterOnSingleParticleCounts("PdgFilter");
214  // lambfilt->AndMaxPdgCodes(0, 3122, -3122); // filter out Lambda0
215  // lambfilt->AndMaxPdgCodes(0, 3222, -3222); // filter out Sigma+
216  // primGen->AndFilter(lambfilt);
217  }
218 
219  // set event filters
220  //-----------------------------
221  if (useEventFilter)
222  {
223  // fastSim->SetMultFilter(type, min, max); requires min <= mult <= max
224  // available types are: '+', '-', 'gam', 'pi0' (2 gam, 0.135 +- 0.03 GeV), 'eta' (2 gam, 0.547 +- 0.04 GeV), 'ks' (pi+ pi-, 0.497 +- 0.04 GeV)
225 
226  // fastSim->SetInvMassFilter(comb, m_min, m_max, mult); requires at least mult combined candidates with m_min < m < m_max
227  // comb is a TString describing the combinatoric
228  // - particle codes are: e+ e- mu+ mu- pi+ pi- k+ k- p+ p- gam pi0 ks eta, separated with blank
229 
230  // Examples:
231  // - ("k+ k-", 0.98, 1.1, 2) : forms K+ K- candidate and requires >=2 in the given window
232  // - ("ks k+ pi- cc", 2.8, 3.2,1 ) : forms ks k+ pi- / ks k- pi+ cands and req. at least one in window
233  }
234 
235  // enable the merging of neutrals if they have similar direction
236  fastSim->MergeNeutralClusters(mergeNeutrals);
237 
238  // enable bremsstahlung loss for electrons
239  fastSim->EnableElectronBremsstrahlung(electronBrems);
240 
241  //enable the producting of parametrized neutral (hadronic) split offs
242  // generate electro-magnetic / hadronic split offs in the EMC? switch off when running w/o EMC
243  if (enableSplitoff) fastSim->EnableSplitoffs(splitpars.Data());
244 
245  // simulates covariance matrix for fitting
246  fastSim->SetUseFlatCov(true);
247 
248  // -----------------------------------------------------------------------------------
249  // Now add the detector components
250  // -----------------------------------------------------------------------------------
251 
252  // Full Panda Tracking: STT MVD GEM FTS (scanning over theta range from 0° to 160°)
253  fastSim->AddDetector("ScFts", "thtMin=0. thtMax=5. ptmin=0.0 pmin=0.5 pRes=0.05 thtRes=0.002 phiRes=0.002 efficiency=0.80");
254  fastSim->AddDetector("ScMvdGem", "thtMin=5. thtMax=7.8 ptmin=0.1 pmin=0.0 pRes=0.03 thtRes=0.001 phiRes=0.001 efficiency=0.60");
255  fastSim->AddDetector("ScSttMvdGem", "thtMin=7.8 thtMax=20.9 ptmin=0.1 pmin=0.0 pRes=0.02 thtRes=0.001 phiRes=0.001 efficiency=0.85");
256  fastSim->AddDetector("ScSttMvd", "thtMin=20.9 thtMax=145. ptmin=0.1 pmin=0.0 pRes=0.02 thtRes=0.001 phiRes=0.001 efficiency=0.85");
257  fastSim->AddDetector("ScSttAlone", "thtMin=145. thtMax=159.5 ptmin=0.1 pmin=0.0 pRes=0.04 thtRes=0.006 phiRes=0.007 efficiency=0.25");
258 
259  // Vertexing
260  fastSim->AddDetector("ScVtxNoMvd", "thtMin=0. thtMax=5. ptmin=0.0 vtxRes=0.05 efficiency=1."); // efficiency=1: all tracks found in trackers will get a vertex information
261  fastSim->AddDetector("ScVtxMvd", "thtMin=5. thtMax=145. ptmin=0.1 vtxRes=0.005 efficiency=1."); // efficiency=1: all tracks found in trackers will get a vertex information
262 
263  // EM Calorimeters w/ default parameters
264  fastSim->AddDetector("EmcFS", "thtMin=0.05 thtMax=10.0 aPar=0.013 bPar=0.0283 Emin=0.01 dist=8.2");
265  fastSim->AddDetector("EmcFwCap", "thtMin=10.0 thtMax=22.0 Emin=0.01 dist=2.5");
266  fastSim->AddDetector("EmcBarrel", "thtMin=22.0 thtMax=142.0 Emin=0.01 barrelRadius=0.5");
267  fastSim->AddDetector("EmcBwCap", "thtMin=142.0 thtMax=160.0 Emin=0.01 dist=0.7");
268 
269  // PID Target Spectrometer: STT, MVD, MDT, EMC, Barrel DIRC, Disc DIRC
270  fastSim->AddDetector("MvdPid", "thtMin=5. thtMax=133.6 ptmin=0.1 dEdxResMulti=1. efficiency=1."); // Note: Bethe-Bloch-Landau-Gauss Prametrization from 2008
271  fastSim->AddDetector("SttPid", "thtMin=7.8 thtMax=159.5 ptmin=0.1 dEdxRes=0.15 efficiency=1."); // Note: dEdX parametrization from 2008
272  fastSim->AddDetector("ScMdtPidBarrel", "thtMin=10.0 thtMax=130.0 pmin=0.5 efficiency=0.95 misId=0.01");
273  fastSim->AddDetector("ScEmcPidFwCap", "thtMin=10.0 thtMax=22.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
274  fastSim->AddDetector("ScEmcPidBarrel", "thtMin=22.0 thtMax=142.0 ptmin=0.2 pmin=0.0 efficiency=1.0");
275  fastSim->AddDetector("ScEmcPidBwCap", "thtMin=142.0 thtMax=160.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
276  fastSim->AddDetector("DrcDisc", "thtMin=5.0 thtMax=22.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
277  fastSim->AddDetector("DrcBarrel", "thtMin=22.0 thtMax=140.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
278 
279  // PID Fwd spectrometer PID: RICH, FwdMUO and EMC FS
280  fastSim->AddDetector("ScEmcPidFS", "thtMin=0.5 thtMax=10.0 ptmin=0.0 pmin=0.5 efficiency=1.0");
281  fastSim->AddDetector("ScMdtPidForward", "thtMin=0.0 thtMax=10.0 pmin=0.5 efficiency=0.95 misId=0.01");
282  fastSim->AddDetector("Rich", "angleXMax=5.0 angleYMax=10.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
283 
284 
285  fRun->AddTask(fastSim);
286 
287  // ********************************
288  // *** END Fast Simulation Task ***
289  // ********************************
290 
291 
292  // *****************************
293  // *** BEGIN SoftTriggerTask ***
294  // *****************************
295 
296  if (runST)
297  {
298  // this file contains the trigger line definitions
299  TString triggercfg = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/triggerlines_fsim.cfg";
300 
301  PndSoftTriggerTask *stTask = new PndSoftTriggerTask(Mom, mode, run, triggercfg);
302  stTask->SetFastSimDefaults();
303  fRun->AddTask(stTask);
304  }
305 
306  // ***************************
307  // *** END SoftTriggerTask ***
308  // ***************************
309 
310 
311  // ***********************************
312  // *** BEGIN PndSimpleCombinerTask ***
313  // ***********************************
314 
315  if (!simonly)
316  {
317  PndSimpleCombinerTask *scTask = new PndSimpleCombinerTask(anadecay, anaparms+":algo=PidChargedProbability",Mom, run, mode);
318  scTask->SetPidAlgo("PidChargedProbability");
319  fRun->AddTask(scTask);
320  }
321  // *********************************
322  // *** END PndSimpleCombinerTask ***
323  // *********************************
324 
325 
326  //------------------------- Initialize the RUN -----------------
327  fRun->Init();
328 
329  //------------------------- Run the Simulation -----------------
330  fRun->Run(nEvents);
331 
332  //------------------------- Write Filter Info to File -----------
333  if (usePndEventFilter) primGen->WriteEvtFilterStatsToRootFile();
334 
335  //------------------------Print some info and exit----------------
336  timer.Stop();
337  Double_t rtime = timer.RealTime();
338  Double_t ctime = timer.CpuTime();
339  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
340  return 0;
341 }
342 
PndMultiField * fField
Definition: sim_emc_apd.C:97
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
Int_t run
Definition: autocutx.C:47
void WriteEvtFilterStatsToRootFile(TFile *outputFile=NULL)
Writes all relevant event filter information to the output root file.
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
bool AddDetector(std::string name, std::string params="")
Definition: PndFastSim.cxx:313
void SetPidAlgo(TString algo)
static const double mp
Definition: mzparameters.h:11
void SetUseFlatCov(bool v=true)
Definition: PndFastSim.h:64
FairPrimaryGenerator * primGen
Definition: sim_emc_apd.C:81
FairRunAna * fRun
Definition: hit_dirc.C:58
Int_t mode
Definition: autocutx.C:47
void SetUnstable(int pdg)
Primary generator with added event filtering capabilities.
bool MergeNeutralClusters(bool merge=true, double par=0.389)
Definition: PndFastSim.h:61
void SetStoreTree(Bool_t store=true)
Double_t
void EnableElectronBremsstrahlung(bool brems=true)
Definition: PndFastSim.h:62
Int_t nEvents
Definition: hit_dirc.C:11
TStopwatch timer
Definition: hit_dirc.C:51
gDebug
Definition: sim_emc_apd.C:6
bool EnableSplitoffs(std::string fname="splitpars.dat")
Definition: PndFastSim.cxx:224
Double_t ctime
Definition: hit_dirc.C:114
FairBoxGenerator * boxGen
Definition: sim_emc_apd.C:85
double X
Definition: anaLmdDigi.C:68
void SetVerbosity(int vb)
Definition: PndFastSim.h:59
void quickfsimana(TString Prefix="", TString Decfile="", Float_t Mom=0., TString anadecay="", Int_t nEvents=1000, TString anaparms="", bool runST=false, int run=0, int runmode=0)
TLorentzVector fIni
Definition: full_core_ntp.C:29
Double_t rtime
Definition: hit_dirc.C:113