16 if (prefix==
"" || inputGen==
"" || pbeam==0.)
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";
41 double X = (pbeam*pbeam-2*mp*
mp)/(2*mp);
42 pbeam =
sqrt(X*X-mp*mp);
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);
53 Bool_t enableSplitoff =
true;
54 Bool_t mergeNeutrals =
true;
55 Bool_t electronBrems =
true;
58 Bool_t usePndEventFilter =
false;
62 TString BaseDir = gSystem->Getenv(
"VMCWORKDIR");
63 TString splitpars = BaseDir+
"/fastsim/splitpars.dat";
67 TString OutputFile = prefix+
"_fsim.root";
78 FairRunSim *
fRun =
new FairRunSim();
79 fRun->SetOutputFile(OutputFile.Data());
80 fRun->SetGenerateRunInfo(kFALSE);
83 FairLogger::GetLogger()->SetLogToFile(kFALSE);
91 if (!usePndEventFilter) primGen->
SetVerbose(0);
92 fRun->SetGenerator(primGen);
93 fRun->SetName(
"TGeant3");
99 if (usePndEventFilter)
102 cout <<
"Using PndEventFilter"<<endl;
104 primGen->
AddFilter(
"(t+-;4..) && M(e+ e-; m[3.1,0.6])");
114 if (inputGen.BeginsWith(
"DPM",TString::kIgnoreCase) && !inputGen.EndsWith(
".dec",TString::kIgnoreCase))
118 if (inputGen==
"dpm1") mode = 1;
119 if (inputGen==
"dpm2") mode = 2;
134 primGen->AddGenerator(Dpm);
140 else if (inputGen.BeginsWith(
"FTF",TString::kIgnoreCase) && !inputGen.EndsWith(
".dec",TString::kIgnoreCase))
144 primGen->AddGenerator(Ftf);
150 else if (inputGen.BeginsWith(
"BOX",TString::kIgnoreCase) && !inputGen.EndsWith(
".dec",TString::kIgnoreCase))
170 inputGen.ReplaceAll(
"box",
"");
171 inputGen.ReplaceAll(
" ",
"");
174 while (inputGen.Contains(
":"))
176 TString curpar = inputGen(0,inputGen.Index(
":"));
177 inputGen = inputGen(inputGen.Index(
":")+1,1000);
178 curpar.ReplaceAll(
"[",
"("); curpar.ReplaceAll(
"]",
")");
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);
188 cout <<
"BOX generator range: type["<<BoxType<<
","<<BoxMult<<
"] p["<<BoxMomMin<<
","<<BoxMomMax
189 <<
"] tht["<<BoxThtMin<<
","<<BoxThtMax<<
"]"<<(BoxCosTht?
"*":
"")<<
" phi["<<BoxPhiMin<<
","<<BoxPhiMax<<
"]"<<endl;
195 boxGen->SetPRange(BoxMomMin,BoxMomMax);
196 boxGen->SetPhiRange(BoxPhiMin, BoxPhiMax);
197 boxGen->SetThetaRange(BoxThtMin, BoxThtMax);
199 if (BoxCosTht) boxGen->SetCosTheta();
201 boxGen->SetXYZ(0., 0., 0.);
202 primGen->AddGenerator(boxGen);
211 Resonance.ReplaceAll(
"pbp",
"pbarpSystem");
215 primGen->AddGenerator(EvtGen);
219 primGen->DoTracking(kFALSE);
223 fRun->SetField(fField);
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");
270 fastSim->
AddDetector(
"ScVtxMvd",
"thtMin=5. thtMax=145. ptmin=0.1 vtxRes=0.005 efficiency=1.");
271 fastSim->
AddDetector(
"ScVtxNoMvd",
"thtMin=0. thtMax=5. ptmin=0.0 vtxRes=0.05 efficiency=1.");
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");
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");
295 fastSim->
AddDetector(
"SttPid",
"thtMin=7.8 thtMax=159.5 ptmin=0.1 dEdxRes=0.15 efficiency=1.");
297 fastSim->
AddDetector(
"MvdPid",
"thtMin=5. thtMax=133.6 ptmin=0.1 dEdxResMulti=1. efficiency=1.");
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");
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");
315 fRun->AddTask(fastSim);
328 FairSystemInfo sysInfo;
329 Float_t maxMemory=sysInfo.GetMaxMemory();
333 Float_t cpuUsage=ctime/
rtime;
336 cout <<
"[INFO ] Macro call : prod_fsim.C(\""<<prefix<<
"\", "<<
nEvents<<
", \""<<inputGen<<
"\", "<<pbeam<<
")" <<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;
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)
bool AddDetector(std::string name, std::string params="")
void getRange(TString par, double &min, double &max)
void SetUseFlatCov(bool v=true)
void SetEventPrintFrequency(int freq)
Sets the frequency (accepted events) for printout (verbose>0) of accepted and generated events...
FairPrimaryGenerator * primGen
void SetUnstable(int pdg)
bool MergeNeutralClusters(bool merge=true, double par=0.389)
void SetStoreTree(Bool_t store=true)
void EnableElectronBremsstrahlung(bool brems=true)
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.
bool EnableSplitoffs(std::string fname="splitpars.dat")
FairBoxGenerator * boxGen
void SetVerbosity(int vb)
void AddFilter(TString filterStr)
Registers a filter as a string to be parsed Each filter set consist of one filter definition or a num...
void SetVerbose(Int_t verbose=12)
Set the level of commenting output.