FairRoot/PandaRoot
production/dayone-2017/fastsim/prod_fsim.C
Go to the documentation of this file.
1 // Macro for running fast simulation
3 // to print usage information:
4 // root -l -b -q prod_fsim.C
5 //
6 // to run with specific setting:
7 // root -l -b -q 'prod_fsim.C("EvtD0D0b", 100, "D0toKpi.dec:pbarpSystem0", 12.)'
8 // root -l -b -q 'prod_fsim.C("DpmInel", 100, "DPM", 12.)'
9 // root -l -b -q 'prod_fsim.C("Box1Kp", 100, "BOX:type(321,1):p(0.1,10):tht(22,140):phi(0,360)",1.)'
10 
11 void getRange(TString par, double &min, double &max);
12 TString getInitialResonance(TString &fEvtGenFile);
13 int SplitString(TString s, TString delim, std::vector<TString> &toks);
14 
15 int prod_fsim(TString prefix="", Int_t nEvents = 100, TString inputGen="", Float_t pbeam = 0. , TString simopt="")
16 {
17  if (prefix=="" || inputGen=="" || pbeam==0.)
18  {
19  cout << "USAGE:\n";
20  cout << "prod_fsim.C( <pref>, <nevt>, <gen>, <pbeam>, <simopt> )\n\n";
21  cout << " <pref> : output file names prefix\n";
22  cout << " <nevt> : number of events\n";
23  cout << " <gen> : generator input: EvtGen decfile; DPM/FTF/BOX uses DPM/FTF generator (inelastic mode) or BOX generator instead\n";
24  cout << " DPM settings: DPM = inelastic only, DPM1 = inel. + elastic, DPM2 = elastic only\n";
25  cout << " FTF settings: FTF = inel. + elastic, FTF1 = inelastic only\n";
26  cout << " BOX settings: type[pdgcode,mult] and optional ranges 'p/tht/phi[min,max]' separated with colon; example: 'BOX:type[211,1]:p[1,5]:tht[45]:phi[90,210]'\n";
27  cout << " <pbeam> : pbar momentum (for BOX generator it still controls the magnetic field) \n";
28  cout << " <simopt> : option like e.g. \"SetupB:Filter2:EMC[22-120]\"\n\n";
29  cout << "Example 1 : root -l -b -q 'prod_fsim.C(\"EvtD0D0b\", 100, \"D0toKpi.dec:pbarpSystem0\", 12.)'\n";
30  cout << "Example 2 : root -l -b -q 'prod_fsim.C(\"DpmInel\", 100, \"DPM\", 12.)'\n";
31  cout << "Example 3 : root -l -b -q 'prod_fsim.C(\"SingleK\", 100, \"BOX:type[321,1]:p[0.1,10]:tht[22,140]:phi[0,360]\", 12.)'\n\n";
32 
33  return 0;
34  }
35 
36  // ----- BEGIN EVALUATE simopt string -----------
37  // set the detector day1 setup (further reduction of phase 1 setup)
38  // defaults
39  int setupmode = 0; // Setup A or B
40 
41  double thtmin = 22.; // min polar angle of EMC barrel
42  double thtmax = 142.; // max polar angle of EMC barrel
43 
44  double phimin = 0.; // min phi angle of EMC barrel (adds opposite site +180 automatically)
45  double phimax = 180.; // max phi angle of EMC barrel (adds opposite site +180 automatically)
46 
47  int filtermode = -1; // event filter (-1 = no filter)
48 
49  int rndseed = 0; // random seed
50 
51  simopt.ToLower();
52  std::vector<TString> opts;
53  SplitString(simopt,":",opts);
54 
55  for (auto op:opts)
56  {
57  // Setup A: Full with 75% GEM, 35% FTS, no EMC Barrel Xtals, no TOF Fwd, no MVD Pixels
58  // Setup B: Full with 75% GEM, 35% FTS, no EMC Barrel Xtals, no TOF Fwd, no MVD Pixels, no Muon, no Shashlyk
59  if (op=="setupb") setupmode=1;
60 
61  // tht range
62  if (op.BeginsWith("emc"))
63  {
64  op.ReplaceAll("emc[","");
65  op.ReplaceAll("]","");
66 
67  thtmin = ((TString)op(0,op.Index("-"))).Atof();
68  thtmax = ((TString)op(op.Index("-")+1,1000)).Atof();
69  }
70 
71  // phi range
72  if (op.BeginsWith("slc"))
73  {
74  op.ReplaceAll("slc[","");
75  op.ReplaceAll("]","");
76 
77  phimin = ((TString)op(0,op.Index("-"))).Atof();
78  phimax = ((TString)op(op.Index("-")+1,1000)).Atof();
79  }
80 
81  // event filter
82  if (op.BeginsWith("filter"))
83  {
84  op.ReplaceAll("filter","");
85  filtermode = op.Atoi();
86  }
87 
88  // random seed
89  if (op.BeginsWith("seed"))
90  {
91  op.ReplaceAll("seed","");
92  rndseed = op.Atoi();
93  }
94  }
95 
96  gRandom->SetSeed(rndseed);
97 
98 
99  cout <<"------------------------------"<<endl;
100  cout <<"Setup : "<<(setupmode==0?"A":"B")<<endl;
101  cout <<"Filter : "<<filtermode<<endl;
102  cout <<"EMC Barrel tht : "<<thtmin<<" < tht < "<<thtmax<<endl;
103  cout <<"EMC Barrel phi : "<<phimin<<" < phi < "<<phimax<<" OR "<<phimin+180.<<" phi "<<phimax+180<<endl;
104  cout <<"Rndm seed : "<<gRandom->GetSeed()<<endl;
105  cout <<"------------------------------"<<endl;
106 
107  // ----- END EVALUATE simopt string -----------
108 
109  // persist fast sim output?
110  bool persist = true;
111 
112  // if pbeam<0, interprete as -E_cm
113  double mp = 0.938272;
114  if (pbeam<0)
115  {
116  double X = (pbeam*pbeam-2*mp*mp)/(2*mp);
117  pbeam = sqrt(X*X-mp*mp);
118  }
119  double Ecm = sqrt(pow(sqrt(pbeam*pbeam + mp*mp) + mp,2) - pbeam*pbeam);
120 
121 
122  // Prevent generator from throwing a lot of warnings
123  TLorentzVector fIni(0,0,pbeam,mp+sqrt(pbeam*pbeam+mp*mp));
124  TDatabasePDG *pdg = TDatabasePDG::Instance();
125  pdg->AddParticle("pbarpSystem","pbarpSystem",fIni.M(),kFALSE,0.1,0, "",88888,0);
126  pdg->AddParticle("pbarpSystem0","pbarpSystem0",fIni.M(),kFALSE,0.1,0, "",88880,0);
127  pdg->AddParticle("eta_c1","eta_c1",4.3,kFALSE,0.02,0,"",999441);
128 
129 
130  //----- Switches for Simulation Options ------------------------------
131  Bool_t enableSplitoff = true; // create e.-m. and hadronic split offs
132  Bool_t mergeNeutrals = true; // merge neutrals (for merged pi0s)
133  Bool_t electronBrems = true; // bremsstrahlung loss for electrons
134 
135  //----- Switches for Event Filter Options ------------------------------
136  Bool_t usePndEventFilter = (filtermode >= 0); // enable Panda event filter. *** Needs configuration (see below) ***
137 
138 
139  //-----General settings-----------------------------------------------
140  TString BaseDir = gSystem->Getenv("VMCWORKDIR");
141  TString splitpars = BaseDir+"/fsim/splitpars.dat";
142 
143  //-----User Settings:-------------------------------------------------
144  TString OutputFile = prefix+"_fsim.root";
145  gDebug = 0;
146 
147 
148  // Start a stop watch
149  TStopwatch fTimer;
150  fTimer.Start();
151 
152  // --------------------------------
153  // Create the Simulation run manager
154  // --------------------------------
155  FairRunSim *fRun = new FairRunSim();
156  fRun->SetOutputFile(OutputFile.Data());
157  fRun->SetGenerateRunInfo(kFALSE);
158  //fRun->SetUserConfig(BaseDir+"/macro/prod/g3ConfigNoMC.C");
159 
160  FairLogger::GetLogger()->SetLogToFile(kFALSE);
161 
162 
163  // -------------------------------
164  // Create and Set Event Generator
165  // -------------------------------
166  FairFilteredPrimaryGenerator* primGen = new FairFilteredPrimaryGenerator();
167  if (!usePndEventFilter) primGen->SetVerbose(0);
168  fRun->SetGenerator(primGen);
169  fRun->SetName("TGeant3");
170 
171 
172  // ------------------------------------------------------
173  // Determine event generator according to inputGen
174  // ------------------------------------------------------
175 
176  // ------------------------------------------------------
177  // use DPM generator; default: inelastic @ pbarmom = mom
178  // ------------------------------------------------------
179  if (inputGen.BeginsWith("DPM",TString::kIgnoreCase) && !inputGen.EndsWith(".dec",TString::kIgnoreCase))
180  {
181  inputGen.ToLower();
182  int mode = 0;
183  if (inputGen=="dpm1") mode = 1;
184  if (inputGen=="dpm2") mode = 2;
185 
186  PndDpmDirect *Dpm= new PndDpmDirect(pbeam,mode); // 0 = inelastic, 1 = inelastic & elastic, 2 = elastic
187 
188  // since fastsim doesn't have a transport, let all long-living resonances decay by the generator
189  Dpm->SetUnstable(111); // pi0
190  Dpm->SetUnstable(310); // K_S0
191  Dpm->SetUnstable(311); // K0
192  Dpm->SetUnstable(-311); // K0bar
193  Dpm->SetUnstable(3122); // Lambda0
194  Dpm->SetUnstable(-3122); // anti-Lambda0
195  Dpm->SetUnstable(221); // eta*/
196  Dpm->SetUnstable(3222); // Sigma+
197  Dpm->SetUnstable(-3222); // anti-Sigma-
198  Dpm->SetUnstable(3334); // Omega-
199  primGen->AddGenerator(Dpm);
200  }
201 
202  // ------------------------------------------------------
203  // use FTF generator;
204  // ------------------------------------------------------
205  else if (inputGen.BeginsWith("FTF",TString::kIgnoreCase) && !inputGen.EndsWith(".dec",TString::kIgnoreCase))
206  {
207  inputGen.ToLower();
208  PndFtfDirect *Ftf = new PndFtfDirect("anti_proton", "G4_H", 1, "ftfp", pbeam, 0, (inputGen=="ftf1") );
209  primGen->AddGenerator(Ftf);
210  }
211 
212  // ------------------------------------------------------
213  // use BOX generator; defaults
214  // ------------------------------------------------------
215  else if (inputGen.BeginsWith("BOX",TString::kIgnoreCase) && !inputGen.EndsWith(".dec",TString::kIgnoreCase))
216  {
217  // Set Box generator defaults
218  Double_t BoxMomMin = 0.05; // minimum momentum for box generator
219  Double_t BoxMomMax = pbeam; // maximum " "
220  Double_t BoxThtMin = 0. ; // minimum theta for box generator
221  Double_t BoxThtMax = 180.; // maximum " "
222  Double_t BoxPhiMin = -180. ; // minimum phi for box generator
223  Double_t BoxPhiMax = 180.; // maximum " "
224  Bool_t BoxCosTht = false; // isotropic in cos(theta) instead theta
225 
226  Int_t BoxType = 13; // default particle muon
227  Int_t BoxMult = 1; // default particle multiplicity
228  Double_t type=0,mult=0; // ref. parameters for range function
229 
230  inputGen.ToLower();
231 
232  // Parse configuratio string
233  if (inputGen!="box")
234  {
235  inputGen.ReplaceAll("box","");
236  inputGen.ReplaceAll(" ","");
237  inputGen += ":";
238 
239  while (inputGen.Contains(":"))
240  {
241  TString curpar = inputGen(0,inputGen.Index(":"));
242  inputGen = inputGen(inputGen.Index(":")+1,1000);
243  curpar.ReplaceAll("[","("); curpar.ReplaceAll("]",")");
244 
245  if (curpar.BeginsWith("type(")) {getRange(curpar,type,mult); BoxType = (Int_t)type; BoxMult = (Int_t)mult; }
246  if (curpar.BeginsWith("p(")) getRange(curpar,BoxMomMin,BoxMomMax);
247  if (curpar.BeginsWith("tht(")) getRange(curpar,BoxThtMin,BoxThtMax);
248  if (curpar.BeginsWith("ctht(")) {getRange(curpar,BoxThtMin,BoxThtMax); BoxCosTht=true;}
249  if (curpar.BeginsWith("phi(")) getRange(curpar,BoxPhiMin,BoxPhiMax);
250  }
251  }
252 
253  cout <<"BOX generator range: type["<<BoxType<<","<<BoxMult<<"] p["<<BoxMomMin<<","<<BoxMomMax
254  <<"] tht["<<BoxThtMin<<","<<BoxThtMax<<"]"<<(BoxCosTht?"*":"")<<" phi["<<BoxPhiMin<<","<<BoxPhiMax<<"]"<<endl;
255 
256 
257  PndBoxGenerator* boxGen = new PndBoxGenerator(BoxType, BoxMult);
258  boxGen->SetDebug(0);
259 
260  boxGen->SetPRange(BoxMomMin,BoxMomMax); // GeV/c
261  boxGen->SetPhiRange(BoxPhiMin, BoxPhiMax); // Azimuth angle range [degree]
262  boxGen->SetThetaRange(BoxThtMin, BoxThtMax); // Polar angle in lab system range [degree]
263 
264  if (BoxCosTht) boxGen->SetCosTheta();
265 
266  boxGen->SetXYZ(0., 0., 0.); //cm
267  primGen->AddGenerator(boxGen);
268  }
269 
270  // ------------------------------------------------------
271  // EvtGen Generator
272  // ------------------------------------------------------
273  else
274  {
275  TString Resonance = getInitialResonance(inputGen);
276  Resonance.ReplaceAll("pbp","pbarpSystem");
277 
278  PndEvtGenDirect *EvtGen = new PndEvtGenDirect(Resonance, inputGen.Data(), pbeam);
279  EvtGen->SetStoreTree(kTRUE);
280  primGen->AddGenerator(EvtGen);
281  }
282 
283  // ------------- switch off the transport of particles
284  primGen->DoTracking(kFALSE);
285 
286  //---------------------Create and Set the Field(s)----------
287  PndMultiField *fField= new PndMultiField("AUTO");
288  fRun->SetField(fField);
289 
290  // Setup the Fast Simulation Task
291  //-----------------------------
292  PndFastSim* fastSim = new PndFastSim(persist);
293 
294  // increasing verbosity increases the amount of console output (mainly for debugging)
295  fastSim->SetVerbosity(0);
296 
297 
298  //-----------------------------
299  // set PANDA event filters
300  //-----------------------------
301 
302  if (usePndEventFilter)
303  {
304  // *** Example Configuration of Event Filter
305  cout <<"Using FairEventFilter"<<endl;
306  primGen->SetFilterMaxTries(100000);
307 
308  switch (filtermode)
309  {
310  case 0: // ******* J/psi -> e+ e- filter
311  {
312  // require 4 charged tracks
313  FairEvtFilterOnSingleParticleCounts* chrgFilter = new FairEvtFilterOnSingleParticleCounts("chrgFilter");
314  chrgFilter->AndMinCharge(4, FairEvtFilter::kCharged);
315  primGen->AndFilter(chrgFilter);
316 
317  // require 1 ee combination in the mass range 2.6 < m(ee) < 3.5 GeV
318  PndEvtFilterOnInvMassCounts* eeInv= new PndEvtFilterOnInvMassCounts("eeInvMFilter");
319  eeInv->SetPdgCodesToCombine( 11, -11);
320  eeInv->SetMinMaxInvMass( 2.6, 3.6 );
321  eeInv->SetMinMaxCounts(1,10000);
322  primGen->AndFilter(eeInv); //add filter to fFilterList
323  }
324  break;
325 
326  case 1: // ******* J/psi -> mu+ mu- filter
327  {
328  // require 4 charged tracks
329  FairEvtFilterOnSingleParticleCounts* chrgFilter = new FairEvtFilterOnSingleParticleCounts("chrgFilter");
330  chrgFilter->AndMinCharge(4, FairEvtFilter::kCharged);
331  primGen->AndFilter(chrgFilter);
332 
333  // require 1 ee combination in the mass range 2.7 < m(ee) < 3.5 GeV
334  PndEvtFilterOnInvMassCounts* mmInv= new PndEvtFilterOnInvMassCounts("mmInvMFilter");
335  mmInv->SetPdgCodesToCombine( 13, -13);
336  mmInv->SetMinMaxInvMass( 2.6, 3.6 );
337  mmInv->SetMinMaxCounts(1,10000);
338  primGen->AndFilter(mmInv); //add filter to fFilterList
339  }
340  break;
341 
342  case 2: // ******* phi -> K+ K- filter (two cand)
343  {
344  // require 4 charged tracks
345  FairEvtFilterOnSingleParticleCounts* chrgFilter = new FairEvtFilterOnSingleParticleCounts("chrgFilter");
346  chrgFilter->AndMinCharge(4, FairEvtFilter::kCharged);
347  primGen->AndFilter(chrgFilter);
348 
349  // require 2 KK combination in the mass range 0.9 < m(KK) < 1.3 GeV
350  PndEvtFilterOnInvMassCounts* kkInv= new PndEvtFilterOnInvMassCounts("kkInvMFilter");
351  kkInv->SetPdgCodesToCombine( 321, -321);
352  kkInv->SetMinMaxInvMass( 0.9, 1.2 );
353  kkInv->SetMinMaxCounts(2,10000);
354  primGen->AndFilter(kkInv); //add filter to fFilterList
355  }
356  break;
357 
358  case 3: // ******* pbar p -> e+ e- filter
359  {
360  // require ee combination in the mass range |m - Ecm| < 0.5
361  PndEvtFilterOnInvMassCounts* eeInv= new PndEvtFilterOnInvMassCounts("eeInvMFilter");
362  eeInv->SetPdgCodesToCombine( 11, -11);
363  eeInv->SetMinMaxInvMass( Ecm-0.5, Ecm+0.5 );
364  eeInv->SetMinMaxCounts(1,10000);
365  primGen->AndFilter(eeInv); //add filter to fFilterList
366  }
367  break;
368 
369  case 4: // ******* pbar p -> eta_c1 eta, eta_c1->chi_c 2pi0, chi_c -> J/psi gamma, => J/psi -> e+ e- filter
370  {
371  // require 1 ee combination in the mass range 2.6 < m(ee) < 3.5 GeV
372  PndEvtFilterOnInvMassCounts* eeInv= new PndEvtFilterOnInvMassCounts("eeInvMFilter");
373  eeInv->SetPdgCodesToCombine( 11, -11);
374  eeInv->SetMinMaxInvMass( 2.6, 3.6);
375  eeInv->SetMinMaxCounts(1,10000);
376  primGen->AndFilter(eeInv); //add filter to fFilterList
377  }
378  break;
379  default: break;
380  }
381  }
382 
383  // enable the merging of neutrals if they have similar direction
384  //-----------------------------
385  fastSim->MergeNeutralClusters(mergeNeutrals);
386 
387  // enable bremsstahlung loss for electrons
388  //-----------------------------
389  fastSim->EnableElectronBremsstrahlung(electronBrems);
390 
391  //enable the producting of parametrized neutral (hadronic) split offs
392  // generate electro-magnetic / hadronic split offs in the EMC? switch off when running w/o EMC
393  if (enableSplitoff) fastSim->EnableSplitoffs(splitpars.Data());
394 
395  fastSim->SetUseFlatCov(true);
396 
397 
398 
399  // -----------------------------------------------------------------------------------
400  // -----------------------------------------------------------------------------------
401  // ********* BEGIN Fast Simulation Configuration ********
402  // -----------------------------------------------------------------------------------
403  // -----------------------------------------------------------------------------------
404 
405 
406  // Tracking: Set up in parts of theta coverage. All modelled by PndFsmSimpleTracker.
407  // Mind: Numbers on resolution (pRes,thtRes,phiRes) and efficiency are guessed
408  // -----------------------------------------------------------------------------------
409  // - (Full Panda Tracking: STT MVD GEM FTS)
410 
411  // **** full tracking
412 // 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");
413 
414 // 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");
415 // 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");
416 // 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");
417 // 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");
418 
419  // **** reduced tracking
420  fastSim->AddDetector("ScMvdGem", "thtMin=5. thtMax=7.8 ptmin=0.1 pmin=0.0 pRes=0.04 thtRes=0.003 phiRes=0.003 efficiency=0.55");
421  fastSim->AddDetector("ScSttMvdGem", "thtMin=7.8 thtMax=20.9 ptmin=0.1 pmin=0.0 pRes=0.025 thtRes=0.002 phiRes=0.002 efficiency=0.82");
422  fastSim->AddDetector("ScSttMvd", "thtMin=20.9 thtMax=145. ptmin=0.1 pmin=0.0 pRes=0.03 thtRes=0.003 phiRes=0.003 efficiency=0.80");
423  fastSim->AddDetector("ScSttAlone", "thtMin=145. thtMax=159.5 ptmin=0.1 pmin=0.0 pRes=0.05 thtRes=0.006 phiRes=0.007 efficiency=0.25");
424 
425  // -----------------------------------------------------------------------------------
426  // Vertexing
427  // -----------------------------------------------------------------------------------
428  // **** full vertexing
429 // 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
430 // 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
431 
432  // **** reduced vertexing
433  fastSim->AddDetector("ScVtxNoMvd", "thtMin=5. thtMax=145. ptmin=0.1 vtxRes=0.02 efficiency=1."); // efficiency=1: all tracks found in trackers will get a vertex information
434 
435  // -----------------------------------------------------------------------------------
436  // EM Calorimeters w/ default parameters
437  // -----------------------------------------------------------------------------------
438  // (don't have to be set, just to list the available parameters
439  // -----------------------------------------------------------------------------------
440 
441  fastSim->AddDetector("EmcFwCap", "thtMin=10.0 thtMax=22.0 Emin=0.01 dist=2.5");
442  fastSim->AddDetector("EmcBwCap", "thtMin=142.0 thtMax=160.0 Emin=0.01 dist=0.7");
443  if (setupmode==0)
444  fastSim->AddDetector("EmcFS", "thtMin=0.05 thtMax=10.0 aPar=0.013 bPar=0.0283 Emin=0.01 dist=8.2");
445 
446  // configure polar angle range of EMC
447  if (thtmin<thtmax)
448  {
449  fastSim->AddDetector("EmcBarrel1", Form("thtMin=%.1f thtMax=%.1f phiMin=%.1f phiMax=%.1f Emin=0.01 barrelRadius=0.5", thtmin, thtmax, phimin, phimax));
450  fastSim->AddDetector("EmcBarrel2", Form("thtMin=%.1f thtMax=%.1f phiMin=%.1f phiMax=%.1f Emin=0.01 barrelRadius=0.5", thtmin, thtmax, phimin-180., phimax-180.));
451  }
452 
453  // -----------------------------------------------------------------------------------
454  // PID
455  // -----------------------------------------------------------------------------------
456 
457  // Cherenkovs
458  fastSim->AddDetector("DrcBarrel","thtMin=22.0 thtMax=140.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
459  //fastSim->AddDetector("DrcDisc","thtMin=5.0 thtMax=22.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
460  //fastSim->AddDetector("Rich","angleXMax=5.0 angleYMax=10.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
461 
462  // TOF
463  fastSim->AddDetector("Tof","thtMin=22.0 thtMax=140.0 dp=0.01");
464 
465  // Trackers with dE/dx
466  //Note: A dEdX parametrization from 2008
467  fastSim->AddDetector("SttPid","thtMin=7.8 thtMax=159.5 ptmin=0.1 dEdxRes=0.15 efficiency=1.");
468  //Note: A Bethe-Bloch-Landau-Gauss Prametrization from 2008
469  fastSim->AddDetector("MvdPid","thtMin=5. thtMax=133.6 ptmin=0.1 dEdxResMulti=1. efficiency=1.");
470 
471  // Muon counters
472  if (setupmode==0)
473  {
474  fastSim->AddDetector("ScMdtPidBarrel", "thtMin=10.0 thtMax=130.0 pmin=0.5 efficiency=0.95 misId=0.01");
475  fastSim->AddDetector("ScMdtPidForward","thtMin=0.0 thtMax=10.0 pmin=0.5 efficiency=0.95 misId=0.01");
476  }
477 
478  // EMCs
479  fastSim->AddDetector("ScEmcPidFwCap", "thtMin=10.0 thtMax=22.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
480  fastSim->AddDetector("ScEmcPidBwCap", "thtMin=142.0 thtMax=160.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
481 
482  if (setupmode==0)
483  fastSim->AddDetector("ScEmcPidFS", "thtMin=0.5 thtMax=10.0 ptmin=0.0 pmin=0.5 efficiency=1.0");
484 
485  // configure polar angle range of EMC
486  if (thtmin<thtmax)
487  {
488  fastSim->AddDetector("ScEmcPidBarrel1", Form("thtMin=%.1f thtMax=%.1f phiMin=%.1f phiMax=%.1f ptmin=0.2 pmin=0.0 efficiency=1.0", thtmin, thtmax, phimin, phimax));
489  fastSim->AddDetector("ScEmcPidBarrel2", Form("thtMin=%.1f thtMax=%.1f phiMin=%.1f phiMax=%.1f ptmin=0.2 pmin=0.0 efficiency=1.0", thtmin, thtmax, phimin-180., phimax-180.));
490  }
491 
492  // -----------------------------------------------------------------------------------
493  // ********* END Fast Simulation Configuration ********
494  // -----------------------------------------------------------------------------------
495 
496 
497 
498  fRun->AddTask(fastSim);
499 
500  //------------------------- Initialize the RUN -----------------
501  fRun->Init();
502 
503  //------------------------- Run the Simulation -----------------
504  fRun->Run(nEvents);
505 
506  //------------------------- Write Filter Info to File -----------
507  if (usePndEventFilter) primGen->WriteEvtFilterStatsToRootFile();
508 
509  //------------------------Print some info and exit----------------
510  fTimer.Stop();
511  FairSystemInfo sysInfo;
512  Float_t maxMemory=sysInfo.GetMaxMemory();
513  Double_t rtime = fTimer.RealTime();
514  Double_t ctime = fTimer.CpuTime();
515 
516  Float_t cpuUsage=ctime/rtime;
517 
518  cout << endl;
519  cout << "[INFO ] Macro call : prod_fsim.C(\""<<prefix<<"\", "<<nEvents<<", \""<<inputGen<<"\", "<<pbeam<<")" <<endl;
520  cout << "[INFO ] Generated Events : " <<primGen->GetNumberOfGeneratedEvents()<<endl;
521  cout << "[INFO ] Output file : " << OutputFile << endl;
522  cout << "[INFO ] Real time : " << rtime << " s, CPU time " << ctime << "s" << endl;
523  cout << "[INFO ] CPU usage : " << cpuUsage*100. << "%" << endl;
524  cout << "[INFO ] Max Memory : " << maxMemory << " MB" << endl;
525  cout << "[INFO ] Macro finished successfully." << endl<<endl;
526  return 0;
527 }
528 
529 // ------ Helper functions -------
530 
531 void getRange(TString par, double &min, double &max)
532 {
533  par.ReplaceAll(" ","");
534  par = par(par.Index("(")+1, par.Length()-par.Index("(")-2);
535 
536  TString smin=par, smax=par;
537 
538  if (par.Contains(","))
539  {
540  smin = par(0,par.Index(","));
541  smax = par(par.Index(",")+1,1000);
542  }
543 
544  min = smin.Atof();
545  max = smax.Atof();
546 
547  //if (min>max) {double tmp=min;min=max;max=tmp;}
548 }
549 
550 TString getInitialResonance(TString &fEvtGenFile)
551 {
552 
553  TString IniRes="";
554 
555  if (fEvtGenFile.Contains(":")) // is the initial resonance provide as <decfile>.dec:iniRes ?
556  {
557  IniRes = fEvtGenFile(fEvtGenFile.Index(":")+1,1000);
558  fEvtGenFile = fEvtGenFile(0,fEvtGenFile.Index(":"));
559  }
560 
561  if (IniRes=="") // we need to search the decay file
562  {
563  std::ifstream fs(fEvtGenFile.Data());
564  char line[250];
565 
566  while (fs)
567  {
568  fs.getline(line,249);
569  TString s(line);
570  s.ReplaceAll("\r","");
571  if (IniRes=="" && s.Contains("Decay "))
572  {
573  if (s.Contains("#")) s=s(0,s.Index("#"));
574  s.ReplaceAll("Decay ","");
575  s.ReplaceAll(" ","");
576  IniRes = s;
577  }
578  }
579  fs.close();
580  }
581 
582  return IniRes;
583 }
584 // -------------------------------------------
585 // Splits a TString into vector of strings; separation character contained in delim
586 int SplitString(TString s, TString delim, std::vector<TString> &toks)
587 {
588  toks.clear();
589  TObjArray *tok = s.Tokenize(delim);
590  int N = tok->GetEntries();
591  for (int i=0;i<N;++i)
592  {
593  TString token = (((TObjString*)tok->At(i))->String()).Strip(TString::kBoth);
594  toks.push_back(token);
595  }
596  return toks.size();
597 }
FairFilteredPrimaryGenerator::GetNumberOfGeneratedEvents
Int_t GetNumberOfGeneratedEvents()
returns the total (accepted + rejected) number of events generated by the event generators. If no event filters are used this number is equal to the number of simulated events.
Definition: FairFilteredPrimaryGenerator.h:136
fField
PndMultiField * fField
Definition: sim_emc_apd.C:97
PndEvtFilterOnInvMassCounts::SetMinMaxCounts
Bool_t SetMinMaxCounts(Int_t min, Int_t max)
Definition: PndEvtFilterOnInvMassCounts.cxx:29
i
Int_t i
Definition: run_full.C:25
FairFilteredPrimaryGenerator::WriteEvtFilterStatsToRootFile
void WriteEvtFilterStatsToRootFile(TFile *outputFile=NULL)
Writes all relevant event filter information to the output root file.
Definition: FairFilteredPrimaryGenerator.h:159
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
FairFilteredPrimaryGenerator::AndFilter
void AndFilter(FairEvtFilter *filter)
Register a non-veto event filter using a logical AND to connect with previously defined non-veto even...
Definition: FairFilteredPrimaryGenerator.h:72
PndFastSim::AddDetector
bool AddDetector(std::string name, std::string params="")
Definition: PndFastSim.cxx:313
s
TLorentzVector s
Definition: Pnd2DStar.C:50
par
Double_t par[3]
Definition: anaLmdCluster.C:146
getInitialResonance
TString getInitialResonance(TString &fEvtGenFile)
Definition: production/dayone-2017/fastsim/prod_fsim.C:550
FairFilteredPrimaryGenerator::SetFilterMaxTries
void SetFilterMaxTries(Int_t maxTries=99999)
Define the maximum number of times that this object should try to find an event which suits all event...
Definition: FairFilteredPrimaryGenerator.h:120
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
getRange
void getRange(TString par, double &min, double &max)
Definition: production/dayone-2017/fastsim/prod_fsim.C:531
mp
static const double mp
Definition: mzparameters.h:11
PndFastSim::SetUseFlatCov
void SetUseFlatCov(bool v=true)
Definition: PndFastSim.h:64
primGen
FairPrimaryGenerator * primGen
Definition: sim_emc_apd.C:81
PndEvtFilterOnInvMassCounts::SetMinMaxInvMass
Bool_t SetMinMaxInvMass(Double_t min, Double_t max)
Definition: PndEvtFilterOnInvMassCounts.cxx:59
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
mode
Int_t mode
Definition: autocutx.C:47
PndDpmDirect::SetUnstable
void SetUnstable(int pdg)
Definition: PndDpmDirect.cxx:134
FairEvtFilterOnSingleParticleCounts::AndMinCharge
Bool_t AndMinCharge(Int_t min, ChargeState charge)
Definition: FairEvtFilterOnSingleParticleCounts.h:130
FairFilteredPrimaryGenerator
Primary generator with added event filtering capabilities.
Definition: FairFilteredPrimaryGenerator.h:52
PndFastSim::MergeNeutralClusters
bool MergeNeutralClusters(bool merge=true, double par=0.389)
Definition: PndFastSim.h:61
PndEvtGenDirect::SetStoreTree
void SetStoreTree(Bool_t store=true)
Definition: PndEvtGenDirect.h:47
prod_fsim
int prod_fsim(TString prefix="", Int_t nEvents=100, TString inputGen="", Float_t pbeam=0., TString simopt="")
Definition: production/dayone-2017/fastsim/prod_fsim.C:15
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndFastSim::EnableElectronBremsstrahlung
void EnableElectronBremsstrahlung(bool brems=true)
Definition: PndFastSim.h:62
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
gDebug
gDebug
Definition: sim_emc_apd.C:6
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
PndFastSim::EnableSplitoffs
bool EnableSplitoffs(std::string fname="splitpars.dat")
Definition: PndFastSim.cxx:224
SplitString
int SplitString(TString s, TString delim, StrVec &toks)
Definition: invexp.C:43
ctime
Double_t ctime
Definition: hit_dirc.C:114
boxGen
FairBoxGenerator * boxGen
Definition: sim_emc_apd.C:85
X
double X
Definition: anaLmdDigi.C:68
PndFastSim::SetVerbosity
void SetVerbosity(int vb)
Definition: PndFastSim.h:59
PndEvtFilterOnInvMassCounts::SetPdgCodesToCombine
Bool_t SetPdgCodesToCombine(Int_t pdgCode1, Int_t pdgCode2, Int_t pdgCode3=kInvalidPdgCode, Int_t pdgCode4=kInvalidPdgCode, Int_t pdgCode5=kInvalidPdgCode)
Definition: PndEvtFilterOnInvMassCounts.cxx:119
fIni
TLorentzVector fIni
Definition: full_core_ntp.C:29
rtime
Double_t rtime
Definition: hit_dirc.C:113
mult
Double_t mult
Definition: GammaSpectraAnalysis.C:45
FairFilteredPrimaryGenerator::SetVerbose
void SetVerbose(Int_t verbose=12)
Set the level of commenting output.
Definition: FairFilteredPrimaryGenerator.h:177
Generated on Wed Apr 3 2019 08:44:04 for FairRoot/PandaRoot by   doxygen 1.8.5