EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_EICDetector.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_EICDetector.C
1 #ifndef MACRO_FUN4ALLG4EICDETECTOR_C
2 #define MACRO_FUN4ALLG4EICDETECTOR_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <DisplayOn.C>
7 #include <G4Setup_EICDetector.C>
8 #include <G4_Bbc.C>
9 #include <G4_CaloTrigger.C>
10 #include <G4_DSTReader_EICDetector.C>
11 #include <G4_FwdJets.C>
12 #include <G4_Global.C>
13 #include <G4_HIJetReco.C>
14 #include <G4_Input.C>
15 #include <G4_Jets.C>
16 #include <G4_Production.C>
17 #include <G4_QA_EIC.C>
18 
19 #include <qa_modules/QAHistManagerDef.h>
20 
21 #include <fun4all/Fun4AllDstOutputManager.h>
22 #include <fun4all/Fun4AllOutputManager.h>
23 #include <fun4all/Fun4AllServer.h>
24 
25 #include <phool/recoConsts.h>
26 
27 R__LOAD_LIBRARY(libfun4all.so)
28 
29 int Fun4All_G4_EICDetector(
30  const int nEvents = 1,
31  const string &inputFile = "/sphenix/data/data02/review_2017-08-02/single_particle/spacal2d/fieldmap/G4Hits_sPHENIX_e-_eta0_8GeV-0002.root",
32  const string &outputFile = "G4EICDetector.root",
33  const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
34  const int skip = 0,
35  const string &outdir = ".")
36 {
37  //---------------
38  // Fun4All server
39  //---------------
40  Fun4AllServer *se = Fun4AllServer::instance();
41  se->Verbosity(0);
42  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
43  //PHRandomSeed::Verbosity(1);
44 
45  // just if we set some flags somewhere in this macro
46  recoConsts *rc = recoConsts::instance();
47  // By default every random number generator uses
48  // PHRandomSeed() which reads /dev/urandom to get its seed
49  // if the RANDOMSEED flag is set its value is taken as initial seed
50  // which will produce identical results so you can debug your code
51  // rc->set_IntFlag("RANDOMSEED", 12345);
52 
53  //===============
54  // Input options
55  //===============
56 
57  // Either:
58  // read previously generated g4-hits files, in this case it opens a DST and skips
59  // the simulations step completely. The G4Setup macro is only loaded to get information
60  // about the number of layers used for the cell reco code
61  //
62  //Input::READHITS = true;
63  INPUTREADHITS::filename[0] = inputFile;
64 
65  // Or:
66  // Use one or more particle generators
67  // It is run if Input::<generator> is set to true
68  // all other options only play a role if it is active
69  // In case embedding into a production output, please double check your G4Setup_EICDetector.C and G4_*.C consistent with those in the production macro folder
70  // Input::EMBED = true;
71  INPUTEMBED::filename[0] = embed_input_file;
72  // Use Pythia 8
73  // Input::PYTHIA8 = true;
74 
75  // Use Pythia 6
76  // Input::PYTHIA6 = true;
77 
78  // Use Sartre
79  // Input::SARTRE = true;
80 
81  // Simple multi particle generator in eta/phi/pt ranges
82  Input::SIMPLE = true;
83  // Input::SIMPLE_VERBOSITY = 1;
84 
85  // Particle gun (same particles in always the same direction)
86  // Input::GUN = true;
87  Input::GUN_VERBOSITY = 0;
88 
89  // Upsilon generator
90  //Input::UPSILON = true;
91  Input::UPSILON_VERBOSITY = 0;
92 
93  // And/Or read generated particles from file
94 
95  // eic-smear output
96  // Input::READEIC = true;
97  INPUTREADEIC::filename = inputFile;
98 
99  // HepMC2 files
100  // Input::HEPMC = true;
101  Input::VERBOSITY = 0;
102  INPUTHEPMC::filename = inputFile;
103 
104  //-----------------
105  // Initialize the selected Input/Event generation
106  //-----------------
107  InputInit();
108  //--------------
109  // Set generator specific options
110  //--------------
111  // can only be set after InputInit() is called
112 
113  // Simple Input generator:
114  if (Input::SIMPLE)
115  {
116  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 1);
117  if (Input::HEPMC || Input::EMBED)
118  {
119  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
120  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
121  }
122  else
123  {
124  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
125  PHG4SimpleEventGenerator::Uniform,
126  PHG4SimpleEventGenerator::Uniform);
127  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
128  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 0.);
129  }
130  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-4, 4);
131  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
132  INPUTGENERATOR::SimpleEventGenerator[0]->set_p_range(0.1, 30.);
133  }
134  // Upsilons
135  if (Input::UPSILON)
136  {
137  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("mu", 0);
138  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
139  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
140  // Y species - select only one, last one wins
141  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
142  }
143  // particle gun
144  if (Input::GUN)
145  {
146  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
147  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
148  }
149  // pythia6
150  if (Input::PYTHIA6)
151  {
152  INPUTGENERATOR::Pythia6->set_config_file(string(getenv("CALIBRATIONROOT")) + "/Generators/phpythia6_ep.cfg");
153  }
154 
155  //--------------
156  // Set Input Manager specific options
157  //--------------
158  // can only be set after InputInit() is called
159 
160  if (Input::HEPMC)
161  {
162  INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 30, 0); //optional collision smear in space, time
163  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
164  // //optional choice of vertex distribution function in space, time
165  INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
170  //INPUTMANAGER::HepMCInputManager->set_embedding_id(2);
171  }
172 
173  // register all input generators with Fun4All
174  InputRegister();
175 
176  // set up production relatedstuff
177  // Enable::PRODUCTION = true;
178 
179  //======================
180  // Write the DST
181  //======================
182 
183  Enable::DSTOUT = true;
184  DstOut::OutputDir = outdir;
185  DstOut::OutputFile = outputFile;
186  Enable::DSTOUT_COMPRESS = false; // Compress DST files
187 
188  //Option to convert DST to human command readable TTree for quick poke around the outputs
189  // Enable::DSTREADER = true;
190 
191  //======================
192  // What to run
193  //======================
194  // Global options (enabled for all subsystems - if implemented)
195  Enable::ABSORBER = true;
196  // Enable::OVERLAPCHECK = true;
197  // Enable::VERBOSITY = 1;
198 
199  Enable::BBC = true;
200 
201  // whether to simulate the Be section of the beam pipe
202  Enable::PIPE = true;
203  // EIC beam pipe extension beyond the Be-section:
204  //G4PIPE::use_forward_pipes = true;
205 
206  Enable::EGEM = true;
207  Enable::FGEM = true;
208  // barrel tracker
209  Enable::BARREL = false;
210  Enable::FST = false;
211  // mvtx/tpc tracker
212  Enable::MVTX = true;
213  Enable::TPC = true;
214 
215 // Enable::TRACKING = true;
216 // Enable::TRACKING_EVAL = Enable::TRACKING && true;
217  G4TRACKING::DISPLACED_VERTEX = false; // this option exclude vertex in the track fitting and use RAVE to reconstruct primary and 2ndary vertexes
218  // projections to calorimeters
219  G4TRACKING::PROJECTION_CEMC = false;
220  G4TRACKING::PROJECTION_FEMC = false;
221  G4TRACKING::PROJECTION_FHCAL = false;
222 
223  Enable::CEMC = true;
224  // Enable::CEMC_ABSORBER = true;
225  Enable::CEMC_CELL = Enable::CEMC && true;
226  Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
227  Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
228 // Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
229  G4CEMC::TowerDigi = RawTowerDigitizer::kNo_digitization;
230 
231  Enable::HCALIN = true;
232  // Enable::HCALIN_ABSORBER = true;
233  Enable::HCALIN_CELL = Enable::HCALIN && true;
234  Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
235  Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
236 // Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
237  G4HCALIN::TowerDigi = RawTowerDigitizer::kNo_digitization;
238 
239  Enable::MAGNET = true;
240 
241  Enable::HCALOUT = true;
242  // Enable::HCALOUT_ABSORBER = true;
243  Enable::HCALOUT_CELL = Enable::HCALOUT && true;
244  Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
245  Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
246 // Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
247  G4HCALOUT::TowerDigi = RawTowerDigitizer::kNo_digitization;
248 
249  // EICDetector geometry - barrel
250 // DIRC occasionally produces lots of photons which allocates tons of memory
251 // which we cannot handle right now, needs fixing
252 // Enable::DIRC = true;
253 
254  // EICDetector geometry - 'hadron' direction
255  Enable::RICH = true;
256  Enable::AEROGEL = true;
257 
258  Enable::FEMC = true;
259  // Enable::FEMC_ABSORBER = true;
260  Enable::FEMC_TOWER = Enable::FEMC && true;
261  Enable::FEMC_CLUSTER = Enable::FEMC_TOWER && true;
262 // Enable::FEMC_EVAL = Enable::FEMC_CLUSTER && true;
263  G4FEMC::sampling_fraction = 0.4; // for pi
264 // G4FEMC::sampling_fraction = 0.249; // for e-
265 
266  Enable::FHCAL = true;
267  // Enable::FHCAL_ABSORBER = true;
268  Enable::FHCAL_TOWER = Enable::FHCAL && true;
269  Enable::FHCAL_CLUSTER = Enable::FHCAL_TOWER && true;
270 // Enable::FHCAL_EVAL = Enable::FHCAL_CLUSTER && true;
271 // set our sampling fraction based on electrons
272  G4FHCAL::sampling_fraction = 0.018;
273 
274  // EICDetector geometry - 'electron' direction
275  Enable::EEMC = true;
276  Enable::EEMC_TOWER = Enable::EEMC && true;
277  Enable::EEMC_CLUSTER = Enable::EEMC_TOWER && true;
278 // Enable::EEMC_EVAL = Enable::EEMC_CLUSTER && true;
279  G4EEMC::TowerDigi = RawTowerDigitizer::kNo_digitization;
280 
281  // Enable::PLUGDOOR = true;
282 
283  // Other options
284  // Enable::GLOBAL_RECO = true;
285 // Enable::GLOBAL_FASTSIM = true;
286 
287 // Enable::CALOTRIGGER = true && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
288 
289  // Select only one jet reconstruction- they currently use the same
290  // output collections on the node tree!
291 // Enable::JETS = true;
292  Enable::JETS_EVAL = Enable::JETS && true;
293 
294 // Enable::FWDJETS = true;
295  Enable::FWDJETS_EVAL = Enable::FWDJETS && true;
296 
297  // HI Jet Reco for jet simulations in Au+Au (default is false for
298  // single particle / p+p simulations, or for Au+Au simulations which
299  // don't care about jets)
300  Enable::HIJETS = false && Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
301 
302 // Enable::QA = true;
303 
304  // new settings using Enable namespace in GlobalVariables.C
305  Enable::BLACKHOLE = true;
306  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
307  //BlackHoleGeometry::visible = true;
308 
309  //---------------
310  // World Settings
311  //---------------
312  // G4WORLD::PhysicsList = "QGSP_BERT"; //FTFP_BERT_HP best for calo
313  // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
314 
315  //---------------
316  // Magnet Settings
317  //---------------
318 
319  // G4MAGNET::magfield = "1.4"; // alternatively to specify a constant magnetic field, give a float number, which will be translated to solenoidal field in T, if string use as fieldmap name (including path)
320  // G4MAGNET::magfield = "0"; // field off
321 // This is the 3d fieldmap setting (default)
322  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sphenix3dbigmapxyz.root"); // default map from the calibration database
323  // G4MAGNET::magfield_rescale = 1.; // in case you want to play with field
324 
325 // for old 2d map use these settings
326  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
327  // G4MAGNET::magfield_rescale = -1.4 / 1.5; // make consistent with expected Babar field strength of 1.4T
328 
329  //---------------
330  // Pythia Decayer
331  //---------------
332  // list of decay types in
333  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
334  // default is All:
335  // G4P6DECAYER::decayType = EDecayType::kAll;
336 
337  // Initialize the selected subsystems
338  G4Init();
339 
340  //---------------------
341  // GEANT4 Detector description
342  //---------------------
343 
344  // 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
345  if (!Input::READHITS)
346  {
347  G4Setup();
348  }
349 
350  //---------
351  // BBC Reco
352  //---------
353 
354  if (Enable::BBC)
355  {
356  BbcInit();
357  Bbc_Reco();
358  }
359 
360  //------------------
361  // Detector Division (only for barrel calorimeters)
362  //------------------
363 
364  if (Enable::CEMC_CELL) CEMC_Cells();
365 
366  if (Enable::HCALIN_CELL) HCALInner_Cells();
367 
368  if (Enable::HCALOUT_CELL) HCALOuter_Cells();
369 
370  //-----------------------------
371  // CEMC towering and clustering
372  //-----------------------------
373 
374  if (Enable::CEMC_TOWER) CEMC_Towers();
375  if (Enable::CEMC_CLUSTER) CEMC_Clusters();
376 
377  //-----------------------------
378  // HCAL towering and clustering
379  //-----------------------------
380 
381  if (Enable::HCALIN_TOWER) HCALInner_Towers();
382  if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
383 
384  if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
385  if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
386 
387  //-----------------------------
388  // e, h direction Calorimeter towering and clustering
389  //-----------------------------
390 
391  if (Enable::FEMC_TOWER) FEMC_Towers();
392  if (Enable::FEMC_CLUSTER) FEMC_Clusters();
393 
394  if (Enable::FHCAL_TOWER) FHCAL_Towers();
395  if (Enable::FHCAL_CLUSTER) FHCAL_Clusters();
396 
397  if (Enable::EEMC_TOWER) EEMC_Towers();
398  if (Enable::EEMC_CLUSTER) EEMC_Clusters();
399 
400  if (Enable::DSTOUT_COMPRESS) ShowerCompress();
401 
402  //--------------
403  // SVTX tracking
404  //--------------
405 
406  if (Enable::TRACKING) Tracking_Reco();
407 
408  //-----------------
409  // Global Vertexing
410  //-----------------
411 
412  if (Enable::GLOBAL_RECO)
413  {
414  Global_Reco();
415  }
416  else if (Enable::GLOBAL_FASTSIM)
417  {
418  Global_FastSim();
419  }
420 
421  //-----------------
422  // Calo Trigger Simulation
423  //-----------------
424 
425  if (Enable::CALOTRIGGER) CaloTrigger_Sim();
426 
427  //---------
428  // Jet reco
429  //---------
430 
431  if (Enable::JETS) Jet_Reco();
432 
433  if (Enable::HIJETS) HIJetReco();
434 
435  if (Enable::FWDJETS) Jet_FwdReco();
436 
437  if (Enable::QA) QAInit();
438 
439  string outputroot = outputFile;
440  string remove_this = ".root";
441  size_t pos = outputroot.find(remove_this);
442  if (pos != string::npos)
443  {
444  outputroot.erase(pos, remove_this.length());
445  }
446 
447  if (Enable::DSTREADER) G4DSTreader_EICDetector(outputroot + "_DSTReader.root");
448 
449  //----------------------
450  // Simulation evaluation
451  //----------------------
452  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4tracking_eval.root");
453 
454  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
455 
456  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
457 
458  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
459 
460  if (Enable::FEMC_EVAL) FEMC_Eval(outputroot + "_g4femc_eval.root");
461 
462  if (Enable::FHCAL_EVAL) FHCAL_Eval(outputroot + "_g4fhcal_eval.root");
463 
464  if (Enable::EEMC_EVAL) EEMC_Eval(outputroot + "_g4eemc_eval.root");
465 
466  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
467 
468  if (Enable::FWDJETS_EVAL) Jet_FwdEval(outputroot + "_g4fwdjet_eval.root");
469 
470  //--------------
471  // Set up Input Managers
472  //--------------
473 
474  InputManagers();
475 
476  //--------------
477  // Set up Output Manager
478  //--------------
479  if (Enable::PRODUCTION)
480  {
481  Production_CreateOutputDir();
482  }
483 
484  if (Enable::DSTOUT)
485  {
486  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
487  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
488  if (Enable::DSTOUT_COMPRESS) DstCompress(out);
489  se->registerOutputManager(out);
490  }
491 
492  //-----------------
493  // Event processing
494  //-----------------
495  if (nEvents < 0)
496  {
497  return 0;
498  }
499  // if we run any of the particle generators and use 0 it'll run forever
500  if (nEvents == 0 && !Input::READHITS && !Input::HEPMC && !Input::READEIC)
501  {
502  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
503  cout << "it will run forever, so I just return without running anything" << endl;
504  return 0;
505  }
506 
507  se->skip(skip);
508  se->run(nEvents);
509 
510  if (Enable::QA)
511  {
512  QAHistManagerDef::saveQARootFile(outputroot + "_qa.root");
513  }
514 
515  //-----
516  // Exit
517  //-----
518 
519  se->End();
520  std::cout << "All done" << std::endl;
521  delete se;
522  if (Enable::PRODUCTION)
523  {
524  Production_MoveOutput();
525  }
526  gSystem->Exit(0);
527  return 0;
528 }
529 #endif