FairRoot/PandaRoot
Functions
simfast_opt.C File Reference

Go to the source code of this file.

Functions

int simfast_opt (TString Prefix, TString Decfile, Float_t Mom, Int_t nEvents=1000, TString Resonance="pbarpSystem0", int Pdgcode=11, TString DetOpt="MvdGem EmcBar Drc Dsc FwdSpec")
 

Function Documentation

int simfast_opt ( TString  Prefix,
TString  Decfile,
Float_t  Mom,
Int_t  nEvents = 1000,
TString  Resonance = "pbarpSystem0",
int  Pdgcode = 11,
TString  DetOpt = "MvdGem EmcBar Drc Dsc FwdSpec" 
)

Definition at line 25 of file simfast_opt.C.

References PndFastSim::AddDetector(), Bool_t, boxGen, ctime, Double_t, PndFastSim::EnableElectronBremsstrahlung(), PndFastSim::EnableSplitoffs(), fField, fIni, fRun, gDebug, PndFastSim::MergeNeutralClusters(), nEvents, primGen, printf(), rtime, PndFastSim::SetInvMassFilter(), PndFastSim::SetMultFilter(), PndEvtGenDirect::SetStoreTree(), PndDpmDirect::SetUnstable(), PndFastSim::SetUseFlatCov(), PndFastSim::SetVerbosity(), sqrt(), timer, and TString.

28 {
29  // Prevent generator from throwing a lot of warnings
30  TLorentzVector fIni(0,0,Mom,0.938272+sqrt(Mom*Mom+0.938272*0.938272));
31  TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",fIni.M(),kFALSE,0.1,0, "",88888);
32  TDatabasePDG::Instance()->AddParticle("pbarpSystem0","pbarpSystem0",fIni.M(),kFALSE,0.1,0, "",88880);
33 
34  //-----Evaluate Detector Setup ---------------------------------------
35  bool SwMvdGem = false;
36  bool SwEmcBar = false;
37  bool SwDrc = false;
38  bool SwDsc = false;
39  bool SwFwdSpec = false;
40 
41  if (DetOpt.Contains("MvdGem") || DetOpt.Contains("1") ) SwMvdGem = true;
42  if (DetOpt.Contains("EmcBar") || DetOpt.Contains("2") ) SwEmcBar = true;
43  if (DetOpt.Contains("Drc") || DetOpt.Contains("3") ) SwDrc = true;
44  if (DetOpt.Contains("Dsc") || DetOpt.Contains("4") ) SwDsc = true;
45  if (DetOpt.Contains("FwdSpec") || DetOpt.Contains("5") ) SwFwdSpec = true;
46 
47  //----- Switches for Simulation Options ------------------------------
48  Bool_t enableSplitoff = false; // create e.-m. and hadronic split offs
49  Bool_t mergeNeutrals = false; // merge neutrals (for merged pi0s)
50  Bool_t electronBrems = false; // bremsstrahlung loss for electrons
51  Bool_t useEventFilter = false; // enable event filter. *** Needs configuration (see below) ***
52 
53  //-----General settings-----------------------------------------------
54  TString BaseDir = gSystem->Getenv("VMCWORKDIR");
55  TString splitpars = BaseDir+"/fsim/splitpars.dat";
56  gRandom->SetSeed();
57 
58  //-----User Settings:-------------------------------------------------
59  TString OutputFile = Prefix+"_fast.root";
60  gDebug = 0;
61 
62  // choose your event generator
63  Bool_t UseEvtGenDirect = kTRUE;
64  Bool_t UseDpm = kFALSE;
65  Bool_t UseBoxGenerator = kFALSE;
66 
67  // use DPM generator; default: inelastic @ pbarmom = mom
68  if (Decfile=="DPM")
69  {
70  UseEvtGenDirect = kFALSE;
71  UseDpm = kTRUE;
72  }
73 
74  // use BOX generator; default: single mu-, 0<tht<180, 0<phi<360, 0.1<p<mom
75  if (Decfile=="BOX")
76  {
77  UseEvtGenDirect = kFALSE;
78  UseBoxGenerator = kTRUE;
79  }
80 
81  Double_t MomMin = 0.1; // minimum momentum for box generator
82 
83  // for negative values of Mom and use of Box generator -> Just generate tracks with p=-Mom
84  if (Mom<0)
85  {
86  Mom = -Mom;
87  MomMin = Mom;
88  }
89  Double_t MomMax = Mom; // maximum " "
90 
91  // Start a stop watch
92  TStopwatch timer;
93  timer.Start();
94 
95  // Create the Simulation run manager
96  // --------------------------------
97  FairRunSim *fRun = new FairRunSim();
98  fRun->SetOutputFile(OutputFile.Data());
99  fRun->SetGenerateRunInfo(kFALSE);
100 
101  FairLogger::GetLogger()->SetLogToFile(kFALSE);
102 
103 
104  // Create and Set Event Generator
105  // -------------------------------
106  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
107  fRun->SetGenerator(primGen);
108  fRun->SetName("TGeant3");
109 
110  if(UseBoxGenerator)
111  { // Box Generator
112  FairBoxGenerator* boxGen = new FairBoxGenerator(Pdgcode, 1); // 211 = pion; 1 = multipl.
113  boxGen->SetPRange(MomMin,MomMax); // GeV/c
114  boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree]
115  boxGen->SetThetaRange(0., 180.); // Polar angle in lab system range [degree]
116  boxGen->SetXYZ(0., 0., 0.); //cm
117  primGen->AddGenerator(boxGen);
118  }
119  if(UseDpm)
120  {
121  PndDpmDirect *Dpm= new PndDpmDirect(Mom,0); // 0 = inelastic, 1 = inelastic & elastic, 2 = elastic
122  Dpm->SetUnstable(111); // pi0
123  Dpm->SetUnstable(310); // K_S0
124  Dpm->SetUnstable(3122); // Lambda
125  Dpm->SetUnstable(-3122); // anti-Lambda
126  Dpm->SetUnstable(221); // eta
127  primGen->AddGenerator(Dpm);
128  }
129  if(UseEvtGenDirect)
130  {
131  PndEvtGenDirect *EvtGen = new PndEvtGenDirect(Resonance, Decfile.Data(), Mom);
132  EvtGen->SetStoreTree(kTRUE);
133  primGen->AddGenerator(EvtGen);
134  }
135 
136  // ------------- switch off the transport of particles
137  primGen->DoTracking(kFALSE);
138 
139  //---------------------Create and Set the Field(s)----------
140  PndMultiField *fField= new PndMultiField("FULL");
141  fRun->SetField(fField);
142 
143  //Setup the Fast Simulation Task
144  //-----------------------------
145  PndFastSim* fastSim = new PndFastSim();
146 
147  // increasing verbosity increases the amount of console output (mainly for debugging)
148  fastSim->SetVerbosity(0);
149 
150  //set event filters
151  if (useEventFilter)
152  {
153  // Filters are:
154  // -----------
155  // fastSim->SetMultFilter(type, min, max);
156  // requires min <= mult <= max
157 
158  // available types are:
159 
160  // "+" : positive charged particles
161  // "-" : negative charged particles
162  // "gam" : gammas
163  // "pi0" : pi0 candidates ( -> 2 gammas); mass window 0.135 +- 0.03 GeV
164  // "eta" : eta candidates ( -> 2 gammas); mass window 0.547 +- 0.04 GeV
165  // "ks" : K_S candidates ( -> pi+ pi-); mass window 0.497 +- 0.04 GeV
166 
167  fastSim->SetMultFilter("+", 2,1000); // at least 2 trk+
168  fastSim->SetMultFilter("-", 2,1000); // at least 2 trk-
169  fastSim->SetMultFilter("gam", 0, 4); // at most 4 gammas
170 
171  // fastSim->SetInvMassFilter(comb, m_min, m_max, mult);
172 
173  // requires at least mult combined candidates with m_min < m < m_max
174 
175  // comb is a TString describing the combinatoric
176  // - particle codes are: e+ e- mu+ mu- pi+ pi- k+ k- p+ p- gam pi0 ks eta
177  // - codes must be separated with a single blank
178  // - for charged final states only the mass is set; no pdg code selection is done!
179  // - optional a 'cc' added at the end of also takes into account charge conjugation
180 
181  // Examples:
182  // - ("k+ k-", 0.98, 1.1, 2) : forms K+ K- candidate and requires >=2 in the given window
183  // - ("ks k+ pi- cc", 2.8, 3.2,1 ) : forms ks k+ pi- / ks k- pi+ cands and req. at least one in window
184 
185  fastSim->SetInvMassFilter("e+ e-",2.8,3.3,1); // look for J/psi -> e+ e- candidate
186  }
187 
188  // enable the merging of neutrals if they have similar direction
189  fastSim->MergeNeutralClusters(mergeNeutrals);
190 
191  // enable bremsstahlung loss for electrons
192  fastSim->EnableElectronBremsstrahlung(electronBrems);
193 
194  //enable the producting of parametrized neutral (hadronic) split offs
195  // generate electro-magnetic / hadronic split offs in the EMC? switch off when running w/o EMC
196 
197  if (enableSplitoff)
198  fastSim->EnableSplitoffs(splitpars.Data());
199 
200  fastSim->SetUseFlatCov(true);
201  // -----------------------------------------------------------------------------------
202  //Tracking: Set up in parts of theta coverage. All modelled by PndFsmSimpleTracker.
203  // Mind: Numbers on resolution (pRes,thtRes,phiRes) and efficiency are guessed
204  // -----------------------------------------------------------------------------------
205  if (SwMvdGem) // MVD and GEM are enabled; combined tracking available
206  {
207  // - (Full Panda Tracking: STT MVD GEM FTS)
208 
209  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");
210  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");
211  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");
212  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");
213  }
214  else // MVD and GEM are disabled; only STT tracking in central region
215  {
216  // - STT alone:
217  fastSim->AddDetector("ScSttAlone", "thtMin=133.6 thtMax=159.5 ptmin=0.1 pmin=0.0 pRes=0.04 thtRes=0.006 phiRes=0.007 efficiency=0.25");
218  fastSim->AddDetector("ScSttAlone2", "thtMin=20.9 thtMax=133.6 ptmin=0.1 pmin=0.0 pRes=0.04 thtRes=0.006 phiRes=0.007 efficiency=0.80");
219  fastSim->AddDetector("ScSttAlone3", "thtMin=7.8 thtMax=20.9 ptmin=0.1 pmin=0.0 pRes=0.04 thtRes=0.006 phiRes=0.007 efficiency=0.25");
220  }
221 
222  if (SwFwdSpec) // Fwd spectrometer enabled -> use Fwd tracking system
223  {
224  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");
225  }
226 
227  // -----------------------------------------------------------------------------------
228  // Vertexing
229  // -----------------------------------------------------------------------------------
230  if (SwMvdGem) // MVD and GEM are enabled -> better vertexing in central region
231  {
232  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
233  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
234  }
235  else // MVD and GEM are disabled -> no good vertexing at all
236  {
237  fastSim->AddDetector("ScVtxNoMvd", "thtMin=0. thtMax=160. ptmin=0.1 vtxRes=0.1 efficiency=1."); // efficiency=1: all tracks found in trackers will get a vertex information
238  }
239  // -----------------------------------------------------------------------------------
240  // EM Calorimeters w/ default parameters
241  // (don't have to be set, just to list the available parameters
242  // -----------------------------------------------------------------------------------
243 
244  fastSim->AddDetector("EmcFwCap", "thtMin=10.0 thtMax=22.0 Emin=0.01 dist=2.5");
245  fastSim->AddDetector("EmcBwCap", "thtMin=142.0 thtMax=160.0 Emin=0.01 dist=0.7");
246 
247  if (SwEmcBar)
248  {
249  // EmcBarrel also allows to set phiMin and phiMax and can be added multiple times as EmcBarrel1, EmcBarrel2, etc.
250  // Should be made constistent with EmcPidBarrel below
251  fastSim->AddDetector("EmcBarrel","thtMin=22.0 thtMax=142.0 Emin=0.01 barrelRadius=0.5");
252  }
253 
254  if (SwFwdSpec) // Fwd spectrometer enabled -> use Fwd EMC
255  {
256  fastSim->AddDetector("EmcFS", "thtMin=0.05 thtMax=10.0 aPar=0.013 bPar=0.0283 Emin=0.01 dist=8.2");
257  }
258 
259  // -----------------------------------------------------------------------------------
260  // PID
261  // -----------------------------------------------------------------------------------
262 
263  // PID detectors being always in: STT, MUO Barrel, EMC FwdCap, EMC BwdCap
264  //Note: A dEdX parametrization from 2008
265  fastSim->AddDetector("SttPid","thtMin=7.8 thtMax=159.5 ptmin=0.1 dEdxRes=0.15 efficiency=1.");
266  fastSim->AddDetector("ScMdtPidBarrel", "thtMin=10.0 thtMax=130.0 pmin=0.5 efficiency=0.95 misId=0.01");
267  fastSim->AddDetector("ScEmcPidFwCap", "thtMin=10.0 thtMax=22.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
268  fastSim->AddDetector("ScEmcPidBwCap", "thtMin=142.0 thtMax=160.0 ptmin=0.0 pmin=0.0 efficiency=1.0");
269 
270  if (SwMvdGem) // MVD and GEM are enabled -> MVD PID available
271  {
272  //Note: A Bethe-Bloch-Landau-Gauss Prametrization from 2008
273  fastSim->AddDetector("MvdPid","thtMin=5. thtMax=133.6 ptmin=0.1 dEdxResMulti=1. efficiency=1.");
274  }
275 
276  if (SwEmcBar) // EMC Barrel enable -> EMC barrel PID available
277  {
278  fastSim->AddDetector("ScEmcPidBarrel", "thtMin=22.0 thtMax=142.0 ptmin=0.2 pmin=0.0 efficiency=1.0");
279  }
280 
281  if (SwDrc) // Barrel DIRC enabled
282  {
283  fastSim->AddDetector("DrcBarrel","thtMin=22.0 thtMax=140.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
284  }
285 
286  if (SwDsc) // Disc DIRC enabled
287  {
288  fastSim->AddDetector("DrcDisc","thtMin=5.0 thtMax=22.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
289  }
290 
291  if (SwFwdSpec) // Fwd spectrometer enabled -> use RICH, FwdMUO and EMC FS
292  {
293  fastSim->AddDetector("ScEmcPidFS", "thtMin=0.5 thtMax=10.0 ptmin=0.0 pmin=0.5 efficiency=1.0");
294  fastSim->AddDetector("Rich","angleXMax=5.0 angleYMax=10.0 dthtc=0.01 nPhotMin=5 effNPhotons=0.075");
295  fastSim->AddDetector("ScMdtPidForward","thtMin=0.0 thtMax=10.0 pmin=0.5 efficiency=0.95 misId=0.01");
296  }
297 
298 
299  fRun->AddTask(fastSim);
300  //------------------------- Initialize the RUN -----------------
301  fRun->Init();
302  //------------------------- Run the Simulation -----------------
303  fRun->Run(nEvents);
304  //------------------------Print some info and exit----------------
305  timer.Stop();
306  Double_t rtime = timer.RealTime();
307  Double_t ctime = timer.CpuTime();
308  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
309  return 0;
310 }
PndMultiField * fField
Definition: sim_emc_apd.C:97
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
bool AddDetector(std::string name, std::string params="")
Definition: PndFastSim.cxx:313
void SetUseFlatCov(bool v=true)
Definition: PndFastSim.h:64
FairPrimaryGenerator * primGen
Definition: sim_emc_apd.C:81
FairRunAna * fRun
Definition: hit_dirc.C:58
void SetUnstable(int pdg)
bool MergeNeutralClusters(bool merge=true, double par=0.389)
Definition: PndFastSim.h:61
void SetInvMassFilter(TString filter, double min, double max, int mult=1)
Definition: PndFastSim.cxx:373
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
void SetMultFilter(TString type, int min, int max=1000)
Definition: PndFastSim.cxx:347
void SetVerbosity(int vb)
Definition: PndFastSim.h:59
TLorentzVector fIni
Definition: full_core_ntp.C:29
Double_t rtime
Definition: hit_dirc.C:113