FairRoot/PandaRoot
production/scripts/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 
14 int prod_fsim(TString prefix="", Int_t nEvents = 100, TString inputGen="", Float_t pbeam = 0. )
15 {
16  if (prefix=="" || inputGen=="" || pbeam==0.)
17  {
18  cout << "USAGE:\n";
19  cout << "prod_fsim.C( <pref>, <nevt>, <gen>, <pbeam> )\n\n";
20  cout << " <pref> : output file names prefix\n";
21  cout << " <nevt> : number of events\n";
22  cout << " <gen> : generator input: EvtGen decfile; DPM/FTF/BOX uses DPM/FTF generator (inelastic mode) or BOX generator instead\n";
23  cout << " DPM settings: DPM = inelastic only, DPM1 = inel. + elastic, DPM2 = elastic only\n";
24  cout << " FTF settings: FTF = inel. + elastic, FTF1 = inelastic only\n";
25  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";
26  cout << " <pbeam> : pbar momentum (for BOX generator it still controls the magnetic field) \n\n";
27  cout << "Example 1 : root -l -b -q 'prod_fsim.C(\"EvtD0D0b\", 100, \"D0toKpi.dec:pbarpSystem0\", 12.)'\n";
28  cout << "Example 2 : root -l -b -q 'prod_fsim.C(\"DpmInel\", 100, \"DPM\", 12.)'\n";
29  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";
30 
31  return 0;
32  }
33 
34  // persist fast sim output?
35  bool persist = true;
36 
37  // if pbeam<0, interprete as -E_cm
38  double mp = 0.938272;
39  if (pbeam<0)
40  {
41  double X = (pbeam*pbeam-2*mp*mp)/(2*mp);
42  pbeam = sqrt(X*X-mp*mp);
43  }
44 
45  // Prevent generator from throwing a lot of warnings
46  TLorentzVector fIni(0,0,pbeam,mp+sqrt(pbeam*pbeam+mp*mp));
47  TDatabasePDG *pdg = TDatabasePDG::Instance();
48  pdg->AddParticle("pbarpSystem","pbarpSystem",fIni.M(),kFALSE,0.1,0, "",88888,0);
49  pdg->AddParticle("pbarpSystem0","pbarpSystem0",fIni.M(),kFALSE,0.1,0, "",88880,0);
50 
51 
52  //----- Switches for Simulation Options ------------------------------
53  Bool_t enableSplitoff = true; // create e.-m. and hadronic split offs
54  Bool_t mergeNeutrals = true; // merge neutrals (for merged pi0s)
55  Bool_t electronBrems = true; // bremsstrahlung loss for electrons
56 
57  //----- Switches for Event Filter Options ------------------------------
58  Bool_t usePndEventFilter = false; // enable Panda event filter. *** Needs configuration (see below) ***
59 
60 
61  //-----General settings-----------------------------------------------
62  TString BaseDir = gSystem->Getenv("VMCWORKDIR");
63  TString splitpars = BaseDir+"/fastsim/splitpars.dat";
64  gRandom->SetSeed();
65 
66  //-----User Settings:-------------------------------------------------
67  TString OutputFile = prefix+"_fsim.root";
68  gDebug = 0;
69 
70 
71  // Start a stop watch
72  TStopwatch fTimer;
73  fTimer.Start();
74 
75  // --------------------------------
76  // Create the Simulation run manager
77  // --------------------------------
78  FairRunSim *fRun = new FairRunSim();
79  fRun->SetOutputFile(OutputFile.Data());
80  fRun->SetGenerateRunInfo(kFALSE);
81  //fRun->SetUserConfig(BaseDir+"/macro/production/scripts/g3ConfigNoMC.C");
82 
83  FairLogger::GetLogger()->SetLogToFile(kFALSE);
84 
85 
86  // -------------------------------
87  // Create and Set Event Generator
88  // -------------------------------
90  primGen->SetEventPrintFrequency(100);
91  if (!usePndEventFilter) primGen->SetVerbose(0);
92  fRun->SetGenerator(primGen);
93  fRun->SetName("TGeant3");
94 
95  //-----------------------------
96  // set PANDA event filters
97  //-----------------------------
98 
99  if (usePndEventFilter)
100  {
101  // *** Example Configuration of Event Filter
102  cout <<"Using PndEventFilter"<<endl;
103  primGen->SetVerbose(1);
104  primGen->AddFilter("(t+-;4..) && M(e+ e-; m[3.1,0.6])"); //require 4 tracks and at least one e+e- candidate in mass window [2.8,3.4]
105  }
106 
107  // ------------------------------------------------------
108  // Determine event generator according to inputGen
109  // ------------------------------------------------------
110 
111  // ------------------------------------------------------
112  // use DPM generator; default: inelastic @ pbarmom = mom
113  // ------------------------------------------------------
114  if (inputGen.BeginsWith("DPM",TString::kIgnoreCase) && !inputGen.EndsWith(".dec",TString::kIgnoreCase))
115  {
116  inputGen.ToLower();
117  int mode = 0;
118  if (inputGen=="dpm1") mode = 1;
119  if (inputGen=="dpm2") mode = 2;
120 
121  PndDpmDirect *Dpm= new PndDpmDirect(pbeam,mode); // 0 = inelastic, 1 = inelastic & elastic, 2 = elastic
122 
123  // since fastsim doesn't have a transport, let all long-living resonances decay by the generator
124  Dpm->SetUnstable(111); // pi0
125  Dpm->SetUnstable(310); // K_S0
126  Dpm->SetUnstable(311); // K0
127  Dpm->SetUnstable(-311); // K0bar
128  Dpm->SetUnstable(3122); // Lambda0
129  Dpm->SetUnstable(-3122); // anti-Lambda0
130  Dpm->SetUnstable(221); // eta*/
131  Dpm->SetUnstable(3222); // Sigma+
132  Dpm->SetUnstable(-3222); // anti-Sigma-
133  Dpm->SetUnstable(3334); // Omega-
134  primGen->AddGenerator(Dpm);
135  }
136 
137  // ------------------------------------------------------
138  // use FTF generator;
139  // ------------------------------------------------------
140  else if (inputGen.BeginsWith("FTF",TString::kIgnoreCase) && !inputGen.EndsWith(".dec",TString::kIgnoreCase))
141  {
142  inputGen.ToLower();
143  PndFtfDirect *Ftf = new PndFtfDirect("anti_proton", "G4_H", 1, "ftfp", pbeam, 0, (inputGen=="ftf1") );
144  primGen->AddGenerator(Ftf);
145  }
146 
147  // ------------------------------------------------------
148  // use BOX generator; defaults
149  // ------------------------------------------------------
150  else if (inputGen.BeginsWith("BOX",TString::kIgnoreCase) && !inputGen.EndsWith(".dec",TString::kIgnoreCase))
151  {
152  // Set Box generator defaults
153  Double_t BoxMomMin = 0.05; // minimum momentum for box generator
154  Double_t BoxMomMax = pbeam; // maximum " "
155  Double_t BoxThtMin = 0. ; // minimum theta for box generator
156  Double_t BoxThtMax = 180.; // maximum " "
157  Double_t BoxPhiMin = 0. ; // minimum phi for box generator
158  Double_t BoxPhiMax = 360.; // maximum " "
159  Bool_t BoxCosTht = false; // isotropic in cos(theta) instead theta
160 
161  Int_t BoxType = 13; // default particle muon
162  Int_t BoxMult = 1; // default particle multiplicity
163  Double_t type=0,mult=0; // ref. parameters for range function
164 
165  inputGen.ToLower();
166 
167  // Parse configuratio string
168  if (inputGen!="box")
169  {
170  inputGen.ReplaceAll("box","");
171  inputGen.ReplaceAll(" ","");
172  inputGen += ":";
173 
174  while (inputGen.Contains(":"))
175  {
176  TString curpar = inputGen(0,inputGen.Index(":"));
177  inputGen = inputGen(inputGen.Index(":")+1,1000);
178  curpar.ReplaceAll("[","("); curpar.ReplaceAll("]",")");
179 
180  if (curpar.BeginsWith("type(")) {getRange(curpar,type,mult); BoxType = (Int_t)type; BoxMult = (Int_t)mult; }
181  if (curpar.BeginsWith("p(")) getRange(curpar,BoxMomMin,BoxMomMax);
182  if (curpar.BeginsWith("tht(")) getRange(curpar,BoxThtMin,BoxThtMax);
183  if (curpar.BeginsWith("ctht(")) {getRange(curpar,BoxThtMin,BoxThtMax); BoxCosTht=true;}
184  if (curpar.BeginsWith("phi(")) getRange(curpar,BoxPhiMin,BoxPhiMax);
185  }
186  }
187 
188  cout <<"BOX generator range: type["<<BoxType<<","<<BoxMult<<"] p["<<BoxMomMin<<","<<BoxMomMax
189  <<"] tht["<<BoxThtMin<<","<<BoxThtMax<<"]"<<(BoxCosTht?"*":"")<<" phi["<<BoxPhiMin<<","<<BoxPhiMax<<"]"<<endl;
190 
191 
192  PndBoxGenerator* boxGen = new PndBoxGenerator(BoxType, BoxMult);
193  boxGen->SetDebug(0);
194 
195  boxGen->SetPRange(BoxMomMin,BoxMomMax); // GeV/c
196  boxGen->SetPhiRange(BoxPhiMin, BoxPhiMax); // Azimuth angle range [degree]
197  boxGen->SetThetaRange(BoxThtMin, BoxThtMax); // Polar angle in lab system range [degree]
198 
199  if (BoxCosTht) boxGen->SetCosTheta();
200 
201  boxGen->SetXYZ(0., 0., 0.); //cm
202  primGen->AddGenerator(boxGen);
203  }
204 
205  // ------------------------------------------------------
206  // EvtGen Generator
207  // ------------------------------------------------------
208  else
209  {
210  TString Resonance = getInitialResonance(inputGen);
211  Resonance.ReplaceAll("pbp","pbarpSystem");
212 
213  PndEvtGenDirect *EvtGen = new PndEvtGenDirect(Resonance, inputGen.Data(), pbeam);
214  EvtGen->SetStoreTree(kTRUE);
215  primGen->AddGenerator(EvtGen);
216  }
217 
218  // ------------- switch off the transport of particles
219  primGen->DoTracking(kFALSE);
220 
221  //---------------------Create and Set the Field(s)----------
222  PndMultiField *fField= new PndMultiField("AUTO");
223  fRun->SetField(fField);
224 
225  // Setup the Fast Simulation Task
226  //-----------------------------
227  PndFastSim* fastSim = new PndFastSim(persist);
228 
229  // increasing verbosity increases the amount of console output (mainly for debugging)
230  fastSim->SetVerbosity(0);
231 
232 
233  // enable the merging of neutrals if they have similar direction
234  //-----------------------------
235  fastSim->MergeNeutralClusters(mergeNeutrals);
236 
237  // enable bremsstahlung loss for electrons
238  //-----------------------------
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  fastSim->SetUseFlatCov(true);
246 
247 
248 
249  // -----------------------------------------------------------------------------------
250  // -----------------------------------------------------------------------------------
251  // ********* BEGIN Fast Simulation Configuration ********
252  // -----------------------------------------------------------------------------------
253  // -----------------------------------------------------------------------------------
254 
255 
256  // Tracking: Set up in parts of theta coverage. All modelled by PndFsmSimpleTracker.
257  // Mind: Numbers on resolution (pRes,thtRes,phiRes) and efficiency are guessed
258  // -----------------------------------------------------------------------------------
259  // - (Full Panda Tracking: STT MVD GEM FTS)
260 
261  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");
262  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");
263  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");
264  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");
265  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");
266 
267  // -----------------------------------------------------------------------------------
268  // Vertexing
269  // -----------------------------------------------------------------------------------
270  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
271  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
272 
273  // -----------------------------------------------------------------------------------
274  // EM Calorimeters w/ default parameters
275  // -----------------------------------------------------------------------------------
276  // (don't have to be set, just to list the available parameters
277  // -----------------------------------------------------------------------------------
278 
279  fastSim->AddDetector("EmcFwCap", "thtMin=10.0 thtMax=22.0 Emin=0.01 dist=2.5");
280  fastSim->AddDetector("EmcBwCap", "thtMin=142.0 thtMax=160.0 Emin=0.01 dist=0.7");
281  fastSim->AddDetector("EmcBarrel","thtMin=22.0 thtMax=142.0 Emin=0.01 barrelRadius=0.5");
282  fastSim->AddDetector("EmcFS", "thtMin=0.05 thtMax=10.0 aPar=0.013 bPar=0.0283 Emin=0.01 dist=8.2");
283 
284  // -----------------------------------------------------------------------------------
285  // PID
286  // -----------------------------------------------------------------------------------
287 
288  // Cherenkovs
289  fastSim->AddDetector("DrcBarrel","thtMin=22.0 thtMax=140.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
290  fastSim->AddDetector("DrcDisc","thtMin=5.0 thtMax=22.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
291  fastSim->AddDetector("Rich","angleXMax=5.0 angleYMax=10.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
292 
293  // Trackers with dE/dc
294  //Note: A dEdX parametrization from 2008
295  fastSim->AddDetector("SttPid","thtMin=7.8 thtMax=159.5 ptmin=0.1 dEdxRes=0.15 efficiency=1.");
296  //Note: A Bethe-Bloch-Landau-Gauss Prametrization from 2008
297  fastSim->AddDetector("MvdPid","thtMin=5. thtMax=133.6 ptmin=0.1 dEdxResMulti=1. efficiency=1.");
298 
299  // Muon counters
300  fastSim->AddDetector("ScMdtPidBarrel", "thtMin=10.0 thtMax=130.0 pmin=0.5 efficiency=0.95 misId=0.01");
301  fastSim->AddDetector("ScMdtPidForward","thtMin=0.0 thtMax=10.0 pmin=0.5 efficiency=0.95 misId=0.01");
302 
303  // EMCs
304  fastSim->AddDetector("ScEmcPidFwCap", "thtMin=10.0 thtMax=22.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
305  fastSim->AddDetector("ScEmcPidBwCap", "thtMin=142.0 thtMax=160.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
306  fastSim->AddDetector("ScEmcPidBarrel", "thtMin=22.0 thtMax=142.0 ptmin=0.2 pmin=0.0 efficiency=1.0");
307  fastSim->AddDetector("ScEmcPidFS", "thtMin=0.5 thtMax=10.0 ptmin=0.0 pmin=0.5 efficiency=1.0");
308 
309  // -----------------------------------------------------------------------------------
310  // ********* END Fast Simulation Configuration ********
311  // -----------------------------------------------------------------------------------
312 
313 
314 
315  fRun->AddTask(fastSim);
316 
317  //------------------------- Initialize the RUN -----------------
318  fRun->Init();
319 
320  //------------------------- Run the Simulation -----------------
321  fRun->Run(nEvents);
322 
323  //------------------------- Write Filter Info to File -----------
324  if (usePndEventFilter) primGen->WriteEvtFilterStatsToRootFile();
325 
326  //------------------------Print some info and exit----------------
327  fTimer.Stop();
328  FairSystemInfo sysInfo;
329  Float_t maxMemory=sysInfo.GetMaxMemory();
330  Double_t rtime = fTimer.RealTime();
331  Double_t ctime = fTimer.CpuTime();
332 
333  Float_t cpuUsage=ctime/rtime;
334 
335  cout << endl;
336  cout << "[INFO ] Macro call : prod_fsim.C(\""<<prefix<<"\", "<<nEvents<<", \""<<inputGen<<"\", "<<pbeam<<")" <<endl;
337  cout << "[INFO ] Generated Events : " <<primGen->GetNumberOfGeneratedEvents()<<endl;
338  cout << "[INFO ] Output file : " << OutputFile << endl;
339  cout << "[INFO ] Real time : " << rtime << " s, CPU time " << ctime << "s" << endl;
340  cout << "[INFO ] CPU usage : " << cpuUsage*100. << "%" << endl;
341  cout << "[INFO ] Max Memory : " << maxMemory << " MB" << endl;
342  cout << "[INFO ] Macro finished successfully." << endl<<endl;
343  return 0;
344 }
345 
346 // ------ Helper functions -------
347 
348 void getRange(TString par, double &min, double &max)
349 {
350  par.ReplaceAll(" ","");
351  par = par(par.Index("(")+1, par.Length()-par.Index("(")-2);
352 
353  TString smin=par, smax=par;
354 
355  if (par.Contains(","))
356  {
357  smin = par(0,par.Index(","));
358  smax = par(par.Index(",")+1,1000);
359  }
360 
361  min = smin.Atof();
362  max = smax.Atof();
363 
364  //if (min>max) {double tmp=min;min=max;max=tmp;}
365 }
366 
368 {
369 
370  TString IniRes="";
371 
372  if (fEvtGenFile.Contains(":")) // is the initial resonance provide as <decfile>.dec:iniRes ?
373  {
374  IniRes = fEvtGenFile(fEvtGenFile.Index(":")+1,1000);
375  fEvtGenFile = fEvtGenFile(0,fEvtGenFile.Index(":"));
376  }
377 
378  if (IniRes=="") // we need to search the decay file
379  {
380  std::ifstream fs(fEvtGenFile.Data());
381  char line[250];
382 
383  while (fs)
384  {
385  fs.getline(line,249);
386  TString s(line);
387  s.ReplaceAll("\r","");
388  if (IniRes=="" && s.Contains("Decay "))
389  {
390  if (s.Contains("#")) s=s(0,s.Index("#"));
391  s.ReplaceAll("Decay ","");
392  s.ReplaceAll(" ","");
393  IniRes = s;
394  }
395  }
396  fs.close();
397  }
398 
399  return IniRes;
400 }
PndMultiField * fField
Definition: sim_emc_apd.C:97
TString getInitialResonance(TString &fEvtGenFile)
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
TLorentzVector s
Definition: Pnd2DStar.C:50
Double_t par[3]
void getRange(TString par, double &min, double &max)
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
static const double mp
Definition: mzparameters.h:11
void SetUseFlatCov(bool v=true)
Definition: PndFastSim.h:64
void SetEventPrintFrequency(int freq)
Sets the frequency (accepted events) for printout (verbose>0) of accepted and generated events...
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)
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
int prod_fsim(TString prefix="", Int_t nEvents=100, TString inputGen="", Float_t pbeam=0.)
Primary generator with added event filtering capabilities.
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.
gDebug
Definition: sim_emc_apd.C:6
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
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 AddFilter(TString filterStr)
Registers a filter as a string to be parsed Each filter set consist of one filter definition or a num...
TLorentzVector fIni
Definition: full_core_ntp.C:29
void SetVerbose(Int_t verbose=12)
Set the level of commenting output.
Double_t rtime
Definition: hit_dirc.C:113
Double_t mult