Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_EICDetector_LQ.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_EICDetector_LQ.C
2  string n = "1000",
3  string ebeam="20",
4  string pbeam="250",
5  string seed="1",
6  string type="3pion"
7  )
8 {
9  // Set the number of TPC layer
10  const int n_TPC_layers = 40; // use 60 for backward compatibility only
11 
12  string inputFile = "/direct/phenix+u/spjeffas/temp/LQGENEP/LQeventgen/LQ_"+type+"_"+seed+"seed.1000event.root";
13 
14  //string inputFile = "/direct/phenix+u/spjeffas/analysis/EICAnalysis/eventgen/DST_p250_e20_"+seed+"seed_"+type+".root";
15 
16  string outputFile = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/g4sim/G4_Leptoquark_DST_p"+pbeam+"_e"+ebeam+"_"+n+"events_"+seed+"seed_SM.root";
17 
18  //Get parameter variables from parameter file
19 
20  int nEvents;
21  stringstream geek(n);
22  geek>>nEvents;
23 
24 
25  //===============
26  // Input options
27  //===============
28 
29  // Either:
30  // read previously generated g4-hits files, in this case it opens a DST and skips
31  // the simulations step completely. The G4Setup macro is only loaded to get information
32  // about the number of layers used for the cell reco code
33  //
34  // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
35  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
36  const bool readhits = false;
37  // Or:
38  // read files in HepMC format (typically output from event generators like hijing or pythia)
39  const bool readhepmc = false; // read HepMC files
40  // Or:
41  // read files in EICTree format generated by eicsmear package
42  const bool readeictree = false;
43  // Or:
44  // Use Pythia 8
45  const bool runpythia8 = false;
46  // Or:
47  // Use Pythia 6
48  const bool runpythia6 = true;
49  // Or:
50  // Use HEPGen
51  const bool runhepgen = false;
52  // Or:
53  // Use Sartre
54  const bool runsartre = false;
55 
56 
57 
58  // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
59  // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
60  const bool particles = false && !readhits;
61  // or gun/ very simple single particle gun generator
62  const bool usegun = false && !readhits;
63  // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also (note, careful to set Z vertex equal to Hijing events)
64  const bool upsilons = false && !readhits;
65 
66  //======================
67  // What to run
68  //======================
69 
70  // sPHENIX barrel
71  bool do_bbc = true;
72 
73  bool do_pipe = true;
74 
75  bool do_svtx = true;
76  bool do_svtx_cell = do_svtx && true;
77  bool do_svtx_track = do_svtx_cell && true;
78  bool do_svtx_eval = do_svtx_track && true;
79 
80  bool do_pstof = false;
81 
82  bool do_cemc = true;
83 
84  bool do_cemc_cell = true;
85  bool do_cemc_twr = true;
86  bool do_cemc_cluster = true;
87  bool do_cemc_eval = true;
88 
89  bool do_hcalin = true;
90  bool do_hcalin_cell = true;
91  bool do_hcalin_twr = true;
92  bool do_hcalin_cluster = true;
93  bool do_hcalin_eval = true;
94 
95  bool do_cemc_cell = do_cemc && true;
96  bool do_cemc_twr = do_cemc_cell && true;
97  bool do_cemc_cluster = do_cemc_twr && true;
98  bool do_cemc_eval = do_cemc_cluster && true;
99 
100  bool do_hcalin = true;
101  bool do_hcalin_cell = do_hcalin && true;
102  bool do_hcalin_twr = do_hcalin_cell && true;
103  bool do_hcalin_cluster = do_hcalin_twr && true;
104  bool do_hcalin_eval = do_hcalin_cluster && true;
105 
106 
107  bool do_magnet = true;
108 
109  bool do_hcalout = true;
110 
111  bool do_hcalout_cell = true;
112  bool do_hcalout_twr = true;
113  bool do_hcalout_cluster = true;
114  bool do_hcalout_eval = true;
115 
116  bool do_global = true;
117  bool do_global_fastsim = true;
118 
119  bool do_jet_reco = true;
120  bool do_jet_eval = true;
121 
122  bool do_fwd_jet_reco = true;
123  bool do_fwd_jet_eval = false;
124 
125  bool do_hcalout_cell = do_hcalout && true;
126  bool do_hcalout_twr = do_hcalout_cell && true;
127  bool do_hcalout_cluster = do_hcalout_twr && true;
128  bool do_hcalout_eval = do_hcalout_cluster && true;
129 
130 
131  // EICDetector geometry - barrel
132  bool do_DIRC = true;
133 
134  // EICDetector geometry - 'hadron' direction
135  bool do_FGEM = true;
136  bool do_FGEM_track = do_FGEM && true;
137 
138  bool do_RICH = true;
139  bool do_Aerogel = true;
140 
141  bool do_FEMC = true;
142  bool do_FEMC_cell = do_FEMC && true;
143  bool do_FEMC_twr = do_FEMC_cell && true;
144  bool do_FEMC_cluster = do_FEMC_twr && true;
145  bool do_FEMC_eval = do_FEMC_cluster && true;
146 
147  bool do_FHCAL = true;
148  bool do_FHCAL_cell = do_FHCAL && true;
149  bool do_FHCAL_twr = do_FHCAL_cell && true;
150  bool do_FHCAL_cluster = do_FHCAL_twr && true;
151  bool do_FHCAL_eval = do_FHCAL_cluster && true;
152 
153  // EICDetector geometry - 'electron' direction
154  bool do_EGEM = true;
155  bool do_EGEM_track = do_EGEM && true;
156 
157  bool do_EEMC = true;
158  bool do_EEMC_cell = do_EEMC && true;
159  bool do_EEMC_twr = do_EEMC_cell && true;
160  bool do_EEMC_cluster = do_EEMC_twr && true;
161  bool do_EEMC_eval = do_EEMC_cluster && true;
162 
163  //do leptoquark analysis modules
164  bool do_lepto_analysis = true;
165 
166  // Other options
167  bool do_global = true;
168  bool do_global_fastsim = false;
169 
170  bool do_calotrigger = false && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
171 
172  bool do_jet_reco = true;
173  bool do_jet_eval = do_jet_reco && true;
174 
175  bool do_fwd_jet_reco = true;
176  bool do_fwd_jet_eval = do_fwd_jet_reco && true;
177 
178  // HI Jet Reco for jet simulations in Au+Au (default is false for
179  // single particle / p+p simulations, or for Au+Au simulations which
180  // don't care about jets)
181  bool do_HIjetreco = false && do_jet_reco && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
182 
183  // Compress DST files
184  bool do_dst_compress = false;
185 
186  //Option to convert DST to human command readable TTree for quick poke around the outputs
187  bool do_DSTReader = false;
188 
189  //---------------
190  // Load libraries
191  //---------------
192 
193  gSystem->Load("libfun4all.so");
194  gSystem->Load("libg4detectors.so");
195  gSystem->Load("libphhepmc.so");
196  gSystem->Load("libg4testbench.so");
197  gSystem->Load("libg4hough.so");
198  gSystem->Load("libg4calo.so");
199  gSystem->Load("libg4eval.so");
200  gSystem->Load("libeicana.so");
201 
202 
203  // establish the geometry and reconstruction setup
204  gROOT->LoadMacro("G4Setup_EICDetector.C");
205  G4Init(do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,n_TPC_layers);
206 
207  int absorberactive = 0; // set to 1 to make all absorbers active volumes
208  // const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
209  const string magfield = "/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
210  const float magfield_rescale = 1.4/1.5; // scale the map to a 1.4 T field
211 
212  //---------------
213  // Fun4All server
214  //---------------
215 
217  se->Verbosity(0); // uncomment for batch production running with minimal output messages
218  // se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running
219 
220  // just if we set some flags somewhere in this macro
222  // By default every random number generator uses
223  // PHRandomSeed() which reads /dev/urandom to get its seed
224  // if the RANDOMSEED flag is set its value is taken as seed
225  // You can either set this to a random value using PHRandomSeed()
226  // which will make all seeds identical (not sure what the point of
227  // this would be:
228  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
229  // or set it to a fixed value so you can debug your code
230  // rc->set_IntFlag("RANDOMSEED", 12345);
231 
232  //-----------------
233  // Event generation
234  //-----------------
235 
236  if (readhits)
237  {
238  // Get the hits from a file
239  // The input manager is declared later
240  }
241  else if (readhepmc)
242  {
243  // this module is needed to read the HepMC records into our G4 sims
244  // but only if you read HepMC input files
245  HepMCNodeReader *hr = new HepMCNodeReader();
246  se->registerSubsystem(hr);
247  }
248  else if (readeictree)
249  {
250  // this module is needed to read the EICTree style records into our G4 sims
251  ReadEICFiles *eicr = new ReadEICFiles();
252  eicr->OpenInputFile(inputFile);
253 
254  se->registerSubsystem(eicr);
255  }
256  else if (runpythia8)
257  {
258  gSystem->Load("libPHPythia8.so");
259 
260  PHPythia8* pythia8 = new PHPythia8();
261  // see coresoftware/generators/PHPythia8 for example config
262  pythia8->set_config_file("/direct/phenix+u/spjeffas/coresoftware/generators/PHPythia8/phpythia8.cfg");
263  se->registerSubsystem(pythia8);
264 
265  HepMCNodeReader *hr = new HepMCNodeReader();
266  se->registerSubsystem(hr);
267  }
268  else if (runpythia6)
269  {
270  gSystem->Load("libPHPythia6.so");
271 
272  PHPythia6 *pythia6 = new PHPythia6();
273  // see coresoftware/generators/PHPythia6 for example config
274  pythia6->set_config_file("../../eventgen/config/phpythia6_ep.cfg");
275 
277  //trigger->SetParticleType(15);
278  trigger->SetQ2Min(1000);
279 
280  pythia6->register_trigger(trigger);
281 
282  se->registerSubsystem(pythia6);
283 
284  HepMCNodeReader *hr = new HepMCNodeReader();
285  se->registerSubsystem(hr);
286  }
287  else if (runhepgen)
288  {
289  gSystem->Load("libsHEPGen.so");
290 
291  sHEPGen *hepgen = new sHEPGen();
292  // see HEPGen source directory/share/vggdata for required .dat files
293  // see HEPGen source directory/share/datacards for required datacard files
294  hepgen->set_datacard_file("hepgen_dvcs.data");
295  hepgen->set_momentum_electron(-20);
296  hepgen->set_momentum_hadron(250);
297  se->registerSubsystem(hepgen);
298 
299  HepMCNodeReader *hr = new HepMCNodeReader();
300  se->registerSubsystem(hr);
301  }
302  else if (runsartre)
303  {
304  // see coresoftware/generators/PHSartre/README for setup instructions
305  // before running:
306  // setenv SARTRE_DIR /opt/sphenix/core/sartre-1.20_root-5.34.36
307  gSystem->Load("libPHSartre.so");
308 
309  PHSartre* mysartre = new PHSartre();
310  // see coresoftware/generators/PHSartre for example config
311  mysartre->set_config_file("sartre.cfg");
312 
313  // particle trigger to enhance forward J/Psi -> ee
314  PHSartreParticleTrigger* pTrig = new PHSartreParticleTrigger("MySartreTrigger");
315  pTrig->AddParticles(-11);
316  //pTrig->SetEtaHighLow(4.0,1.4);
317  pTrig->SetEtaHighLow(1.0,-1.1); // central arm
318  pTrig->PrintConfig();
319  mysartre->register_trigger((PHSartreGenTrigger *)pTrig);
320  se->registerSubsystem(mysartre);
321 
322  HepMCNodeReader *hr = new HepMCNodeReader();
323  se->registerSubsystem(hr);
324  }
325 
326  // If "readhepMC" is also set, the particles will be embedded in Hijing events
327  if(particles)
328  {
329  // toss low multiplicity dummy events
331 
332  //gen->add_particles("e-",5); // mu+,e+,proton,pi+,Upsilon
333  //gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
334  gen->add_particles("tau-",1); // mu-,e-,anti_proton,pi-
335  if (readhepmc) {
336  gen->set_reuse_existing_vertex(true);
337  gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
338  } else {
342  gen->set_vertex_distribution_mean(0.0,0.0,0.0);
343  gen->set_vertex_distribution_width(0.0,0.0,5.0);
344  }
346  gen->set_vertex_size_parameters(0.0,0.0);
347  gen->set_eta_range(0.1, 0.1);
348  //gen->set_eta_range(3.0, 3.0); //EICDetector FWD
349  gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
350  //gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
351  gen->set_p_range(30.0, 30.0);
352 
353  //gen->add_particles("pi-",1); // mu+,e+,proton,pi+,Upsilon
354  //gen->add_particles("pi+",100); // 100 pion option
355  if (readhepmc)
356  {
357  gen->set_reuse_existing_vertex(true);
358  gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
359  }
360  else
361  {
365  gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
366  gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
367  }
369  gen->set_vertex_size_parameters(0.0, 0.0);
370  gen->set_eta_range(-1.0, 1.0);
371  gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
372  //gen->set_pt_range(0.1, 50.0);
373  gen->set_pt_range(0.1, 20.0);
374 
375  gen->Embed(1);
376  gen->Verbosity(0);
377 
378  se->registerSubsystem(gen);
379  }
380  if (usegun)
381  {
382  // PHG4ParticleGun *gun = new PHG4ParticleGun();
383  // gun->set_name("anti_proton");
384  // gun->set_name("geantino");
385  // gun->set_vtx(0, 0, 0);
386  // gun->set_mom(10, 0, 0.01);
387  // gun->AddParticle("geantino",1.7776,-0.4335,0.);
388  // gun->AddParticle("geantino",1.7709,-0.4598,0.);
389  // gun->AddParticle("geantino",2.5621,0.60964,0.);
390  // gun->AddParticle("geantino",1.8121,0.253,0.);
391  // se->registerSubsystem(gun);
393  pgen->set_name("e-");
394  pgen->set_z_range(0,0);
395  pgen->set_eta_range(0.01,0.01);
396  pgen->set_mom_range(10,10);
397  pgen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
398  se->registerSubsystem(pgen);
399  }
400 
401  // If "readhepMC" is also set, the Upsilons will be embedded in Hijing events, if 'particles" is set, the Upsilons will be embedded in whatever particles are thrown
402  if(upsilons)
403  {
404  // run upsilons for momentum, dca performance, alone or embedded in Hijing
405 
407  vgen->add_decay_particles("e+","e-",0); // i = decay id
408  // event vertex
409  if (readhepmc || particles)
410  {
411  vgen->set_reuse_existing_vertex(true);
412  }
413  else
414  {
415  vgen->set_vtx_zrange(-10.0, +10.0);
416  }
417 
418  // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
419  vgen->set_rapidity_range(-1.0, +1.0);
420  vgen->set_pt_range(0.0, 10.0);
421 
422  int istate = 1;
423 
424  if(istate == 1)
425  {
426  // Upsilon(1S)
427  vgen->set_mass(9.46);
428  vgen->set_width(54.02e-6);
429  }
430  else if (istate == 2)
431  {
432  // Upsilon(2S)
433  vgen->set_mass(10.0233);
434  vgen->set_width(31.98e-6);
435  }
436  else
437  {
438  // Upsilon(3S)
439  vgen->set_mass(10.3552);
440  vgen->set_width(20.32e-6);
441  }
442 
443  vgen->Verbosity(0);
444  vgen->Embed(2);
445  se->registerSubsystem(vgen);
446 
447  cout << "Upsilon generator for istate = " << istate << " created and registered " << endl;
448  }
449 
450  if (!readhits)
451  {
452  //---------------------
453  // Detector description
454  //---------------------
455 
456  G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
457  do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,
458  do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,
459  magfield_rescale);
460 
461  }
462 
463  //---------
464  // BBC Reco
465  //---------
466 
467  if (do_bbc)
468  {
469  gROOT->LoadMacro("G4_Bbc.C");
470  BbcInit();
471  Bbc_Reco();
472  }
473 
474  //------------------
475  // Detector Division
476  //------------------
477 
478  if (do_svtx_cell) Svtx_Cells();
479 
480  if (do_cemc_cell) CEMC_Cells();
481 
482  if (do_hcalin_cell) HCALInner_Cells();
483 
484  if (do_hcalout_cell) HCALOuter_Cells();
485 
486  if (do_FEMC_cell) FEMC_Cells();
487 
488  if (do_FHCAL_cell) FHCAL_Cells();
489 
490  if (do_EEMC_cell) EEMC_Cells();
491 
492  //-----------------------------
493  // CEMC towering and clustering
494  //-----------------------------
495 
496  if (do_cemc_twr) CEMC_Towers();
497  if (do_cemc_cluster) CEMC_Clusters();
498 
499  //-----------------------------
500  // HCAL towering and clustering
501  //-----------------------------
502 
503  if (do_hcalin_twr) HCALInner_Towers();
504  if (do_hcalin_cluster) HCALInner_Clusters();
505 
506  if (do_hcalout_twr) HCALOuter_Towers();
507  if (do_hcalout_cluster) HCALOuter_Clusters();
508 
509  //-----------------------------
510  // e, h direction Calorimeter towering and clustering
511  //-----------------------------
512 
513  if (do_FEMC_twr) FEMC_Towers();
514  if (do_FEMC_cluster) FEMC_Clusters();
515 
516  if (do_FHCAL_twr) FHCAL_Towers();
517  if (do_FHCAL_cluster) FHCAL_Clusters();
518 
519  if (do_EEMC_twr) EEMC_Towers();
520  if (do_EEMC_cluster) EEMC_Clusters();
521 
522  if (do_dst_compress) ShowerCompress();
523 
524  //--------------
525  // SVTX tracking
526  //--------------
527 
528  if (do_svtx_track) Svtx_Reco();
529 
530  //--------------
531  // FGEM tracking
532  //--------------
533 
534  if(do_FGEM_track) FGEM_FastSim_Reco();
535 
536  //--------------
537  // EGEM tracking
538  //--------------
539 
540  if(do_EGEM_track) EGEM_FastSim_Reco();
541 
542  //-----------------
543  // Global Vertexing
544  //-----------------
545 
546  if (do_global)
547  {
548  gROOT->LoadMacro("G4_Global.C");
549  Global_Reco();
550  }
551 
552  else if (do_global_fastsim)
553  {
554  gROOT->LoadMacro("G4_Global.C");
555  Global_FastSim();
556  }
557 
558  //-----------------
559  // Calo Trigger Simulation
560  //-----------------
561 
562  if (do_calotrigger)
563  {
564  gROOT->LoadMacro("G4_CaloTrigger.C");
565  CaloTrigger_Sim();
566  }
567 
568 
569 
570 
571  //--------------
572  // IO management
573  //--------------
574 
575  if (readhits)
576  {
577  // Hits file
578  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
579  hitsin->fileopen(inputFile);
580  se->registerInputManager(hitsin);
581  }
582  if (readhepmc)
583  {
585  se->registerInputManager( in );
586  se->fileopen( in->Name().c_str(), inputFile );
587  }
588  else
589  {
590  // for single particle generators we just need something which drives
591  // the event loop, the Dummy Input Mgr does just that
593  se->registerInputManager( in );
594  }
595 
596  if (do_DSTReader)
597  {
598  //Convert DST to human command readable TTree for quick poke around the outputs
599  gROOT->LoadMacro("G4_DSTReader_EICDetector.C");
600 
601  G4DSTreader_EICDetector( outputFile, //
602  /*int*/ absorberactive ,
603  /*bool*/ do_svtx ,
604  /*bool*/ do_cemc ,
605  /*bool*/ do_hcalin ,
606  /*bool*/ do_magnet ,
607  /*bool*/ do_hcalout ,
608  /*bool*/ do_cemc_twr ,
609  /*bool*/ do_hcalin_twr ,
610  /*bool*/ do_magnet ,
611  /*bool*/ do_hcalout_twr,
612  /*bool*/ do_FGEM,
613  /*bool*/ do_EGEM,
614  /*bool*/ do_FHCAL,
615  /*bool*/ do_FHCAL_twr,
616  /*bool*/ do_FEMC,
617  /*bool*/ do_FEMC_twr,
618  /*bool*/ do_EEMC,
619  /*bool*/ do_EEMC_twr
620  );
621  }
622 
623  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
624  if (do_dst_compress) DstCompress(out);
625  se->registerOutputManager(out);
626 
627  //-----------------
628  // Event processing
629  //-----------------
630  if (nEvents < 0)
631  {
632  return;
633  }
634  // if we run the particle generator and use 0 it'll run forever
635  if (nEvents == 0 && !readhits && !readhepmc)
636  {
637  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
638  cout << "it will run forever, so I just return without running anything" << endl;
639  return;
640  }
641 
642  se->run(nEvents);
643 
644  //-----
645  // Exit
646  //-----
647 
648  se->End();
649  std::cout << "All done" << std::endl;
650  delete se;
651  gSystem->Exit(0);
652 }