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