ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_sPHENIX.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_sPHENIX.C
1 #ifndef MACRO_FUN4ALLG4SPHENIX_C
2 #define MACRO_FUN4ALLG4SPHENIX_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <DisplayOn.C>
7 #include <G4Setup_sPHENIX.C>
8 #include <G4_Bbc.C>
9 #include <G4_CaloTrigger.C>
10 #include <G4_DSTReader.C>
11 #include <G4_Global.C>
12 #include <G4_HIJetReco.C>
13 #include <G4_Input.C>
14 #include <G4_Jets.C>
15 #include <G4_ParticleFlow.C>
16 #include <G4_Production.C>
17 #include <G4_TopoClusterReco.C>
18 #include <G4_Tracking.C>
19 #include <G4_User.C>
20 
23 #include <fun4all/Fun4AllServer.h>
24 
25 #include <phool/PHRandomSeed.h>
26 #include <phool/recoConsts.h>
27 
28 R__LOAD_LIBRARY(libfun4all.so)
29 
30 // For HepMC Hijing
31 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
32 
34  const int nEvents = 1,
35  const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
36  const string &outputFile = "G4sPHENIX.root",
37  const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
38  const int skip = 0,
39  const string &outdir = ".")
40 {
42  se->Verbosity(0);
43 
44  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
46 
47  // just if we set some flags somewhere in this macro
49  // By default every random number generator uses
50  // PHRandomSeed() which reads /dev/urandom to get its seed
51  // if the RANDOMSEED flag is set its value is taken as seed
52  // You can either set this to a random value using PHRandomSeed()
53  // which will make all seeds identical (not sure what the point of
54  // this would be:
55  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
56  // or set it to a fixed value so you can debug your code
57  // rc->set_IntFlag("RANDOMSEED", 12345);
58 
59  //===============
60  // Input options
61  //===============
62  // verbosity setting (applies to all input managers)
63  Input::VERBOSITY = 0;
64  // First enable the input generators
65  // Either:
66  // read previously generated g4-hits files, in this case it opens a DST and skips
67  // the simulations step completely. The G4Setup macro is only loaded to get information
68  // about the number of layers used for the cell reco code
69  // Input::READHITS = true;
70  INPUTREADHITS::filename = inputFile;
71 
72  // Or:
73  // Use particle generator
74  // And
75  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
76  // In case embedding into a production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
77  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
78  // Input::EMBED = true;
79  INPUTEMBED::filename = embed_input_file;
80 
81  Input::SIMPLE = true;
82  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
83  // Input::SIMPLE_VERBOSITY = 1;
84 
85  // Input::PYTHIA6 = true;
86 
87  // Input::PYTHIA8 = true;
88 
89  // Input::GUN = true;
90  // Input::GUN_NUMBER = 3; // if you need 3 of them
91  // Input::GUN_VERBOSITY = 1;
92 
93  // Upsilon generator
94  // Input::UPSILON = true;
95  // Input::UPSILON_NUMBER = 3; // if you need 3 of them
96  // Input::UPSILON_VERBOSITY = 0;
97 
98  // Input::HEPMC = true;
99  INPUTHEPMC::filename = inputFile;
100 
101  // Event pile up simulation with collision rate in Hz MB collisions.
102  //Input::PILEUPRATE = 100e3;
103 
104  //-----------------
105  // Initialize the selected Input/Event generation
106  //-----------------
107  // This creates the input generator(s)
108  InputInit();
109 
110  //--------------
111  // Set generator specific options
112  //--------------
113  // can only be set after InputInit() is called
114 
115  // Simple Input generator:
116  // if you run more than one of these Input::SIMPLE_NUMBER > 1
117  // add the settings for other with [1], next with [2]...
118  if (Input::SIMPLE)
119  {
120  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 5);
121  if (Input::HEPMC || Input::EMBED)
122  {
123  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
124  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
125  }
126  else
127  {
128  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
131  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
132  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 5.);
133  }
134  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
135  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
136  INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
137  }
138  // Upsilons
139  // if you run more than one of these Input::UPSILON_NUMBER > 1
140  // add the settings for other with [1], next with [2]...
141  if (Input::UPSILON)
142  {
143  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
144  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
145  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
146  // Y species - select only one, last one wins
147  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
148  }
149  // particle gun
150  // if you run more than one of these Input::GUN_NUMBER > 1
151  // add the settings for other with [1], next with [2]...
152  if (Input::GUN)
153  {
154  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
155  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
156  }
157 
158  //--------------
159  // Set Input Manager specific options
160  //--------------
161  // can only be set after InputInit() is called
162 
163  if (Input::HEPMC)
164  {
165  INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0); //optional collision smear in space, time
166  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
167  // //optional choice of vertex distribution function in space, time
173  //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
174  if (Input::PILEUPRATE > 0)
175  {
176  // Copy vertex settings from foreground hepmc input
178  // and then modify the ones you want to be different
179  // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
180  }
181  }
182  // register all input generators with Fun4All
183  InputRegister();
184 
185  // set up production relatedstuff
186  // Enable::PRODUCTION = true;
187 
188  //======================
189  // Write the DST
190  //======================
191 
192  // Enable::DSTOUT = true;
193  Enable::DSTOUT_COMPRESS = false;
194  DstOut::OutputDir = outdir;
195  DstOut::OutputFile = outputFile;
196 
197  //Option to convert DST to human command readable TTree for quick poke around the outputs
198  // Enable::DSTREADER = true;
199 
200  // turn the display on (default off)
201  Enable::DISPLAY = false;
202 
203  //======================
204  // What to run
205  //======================
206  // Global options (enabled for all enables subsystems - if implemented)
207  // Enable::ABSORBER = true;
208  // Enable::OVERLAPCHECK = true;
209  // Enable::VERBOSITY = 1;
210 
211  // Enable::BBC = true;
212  Enable::BBCFAKE = true; // Smeared vtx and t0, use if you don't want real BBC in simulation
213 
214  Enable::PIPE = true;
215  Enable::PIPE_ABSORBER = true;
216 
217  // central tracking
218  Enable::MVTX = true;
219  Enable::MVTX_CELL = Enable::MVTX && true;
220  Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
221 
222  Enable::INTT = true;
223  Enable::INTT_CELL = Enable::INTT && true;
224  Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
225 
226  Enable::TPC = true;
227  Enable::TPC_ABSORBER = true;
228  Enable::TPC_CELL = Enable::TPC && true;
229  Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
230 
231  //Enable::MICROMEGAS = true;
234 
235  Enable::TRACKING_TRACK = true;
236  Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && true;
237 
238  // cemc electronics + thin layer of W-epoxy to get albedo from cemc
239  // into the tracking, cannot run together with CEMC
240  // Enable::CEMCALBEDO = true;
241 
242  Enable::CEMC = true;
243  Enable::CEMC_ABSORBER = true;
248 
249  Enable::HCALIN = true;
255 
256  Enable::MAGNET = true;
258 
259  Enable::HCALOUT = true;
265 
266  // forward EMC
267  //Enable::FEMC = true;
268  Enable::FEMC_ABSORBER = true;
273 
274  Enable::EPD = false;
275 
277  //Enable::PLUGDOOR = true;
279 
280  Enable::GLOBAL_RECO = true;
281  // Enable::GLOBAL_FASTSIM = true;
282 
283  Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
284 
285  Enable::JETS = true;
287 
288  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
289  // single particle / p+p-only simulations, or for p+Au / Au+Au
290  // simulations which don't particularly care about jets)
292 
293  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
295  // particle flow jet reconstruction - needs topoClusters!
297 
298  // new settings using Enable namespace in GlobalVariables.C
299  Enable::BLACKHOLE = true;
300  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
301  //BlackHoleGeometry::visible = true;
302 
303  // run user provided code (from local G4_User.C)
304  //Enable::USER = true;
305 
306  //---------------
307  // World Settings
308  //---------------
309  // G4WORLD::PhysicsList = "QGSP_BERT"; //FTFP_BERT_HP best for calo
310  G4WORLD::WorldMaterial = "G4_Galactic"; // set to G4_Galactic for material scans
311 
312  //---------------
313  // Magnet Settings
314  //---------------
315 
316  // const string magfield = "1.5"; // 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)
317  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
318  G4MAGNET::magfield_rescale = -1.4 / 1.5; // make consistent with expected Babar field strength of 1.4T
319 
320  //---------------
321  // Pythia Decayer
322  //---------------
323  // list of decay types in
324  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
325  // default is All:
326  // G4P6DECAYER::decayType = EDecayType::kAll;
327 
328  // Initialize the selected subsystems
329  G4Init();
330 
331  //---------------------
332  // GEANT4 Detector description
333  //---------------------
334  if (!Input::READHITS)
335  {
336  G4Setup();
337  }
338 
339  //------------------
340  // Detector Division
341  //------------------
342 
344 
345  if (Enable::MVTX_CELL) Mvtx_Cells();
346  if (Enable::INTT_CELL) Intt_Cells();
347  if (Enable::TPC_CELL) TPC_Cells();
349 
351 
353 
355 
357 
358  //-----------------------------
359  // CEMC towering and clustering
360  //-----------------------------
361 
364 
365  //-----------------------------
366  // HCAL towering and clustering
367  //-----------------------------
368 
371 
372  if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
374 
375  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
377 
380 
382 
383  //--------------
384  // SVTX tracking
385  //--------------
386  if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
387  if (Enable::INTT_CLUSTER) Intt_Clustering();
388  if (Enable::TPC_CLUSTER) TPC_Clustering();
390 
391  if (Enable::TRACKING_TRACK)
392  {
393  TrackingInit();
394  Tracking_Reco();
395  }
396  //-----------------
397  // Global Vertexing
398  //-----------------
399 
401  {
402  cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
403  gSystem->Exit(1);
404  }
406  {
407  Global_Reco();
408  }
409  else if (Enable::GLOBAL_FASTSIM)
410  {
411  Global_FastSim();
412  }
413 
414  //-----------------
415  // Calo Trigger Simulation
416  //-----------------
417 
418  if (Enable::CALOTRIGGER)
419  {
420  CaloTrigger_Sim();
421  }
422 
423  //---------
424  // Jet reco
425  //---------
426 
427  if (Enable::JETS) Jet_Reco();
428  if (Enable::HIJETS) HIJetReco();
429 
431 
432  //----------------------
433  // Simulation evaluation
434  //----------------------
435  string outputroot = outputFile;
436  string remove_this = ".root";
437  size_t pos = outputroot.find(remove_this);
438  if (pos != string::npos)
439  {
440  outputroot.erase(pos, remove_this.length());
441  }
442 
443  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
444 
445  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
446 
447  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
448 
449  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
450 
451  if (Enable::FEMC_EVAL) FEMC_Eval(outputroot + "_g4femc_eval.root");
452 
453  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
454 
455  if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
456 
458 
459  //--------------
460  // Set up Input Managers
461  //--------------
462 
463  InputManagers();
464 
465  if (Enable::PRODUCTION)
466  {
468  }
469 
470  if (Enable::DSTOUT)
471  {
472  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
473  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
475  se->registerOutputManager(out);
476  }
477  //-----------------
478  // Event processing
479  //-----------------
480  if (Enable::DISPLAY)
481  {
482  DisplayOn();
483 
484  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
485  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
486 
487  cout << "-------------------------------------------------" << endl;
488  cout << "You are in event display mode. Run one event with" << endl;
489  cout << "se->run(1)" << endl;
490  cout << "Run Geant4 command with following examples" << endl;
491  gROOT->ProcessLine("displaycmd()");
492 
493  return 0;
494  }
495 
496  // if we use a negative number of events we go back to the command line here
497  if (nEvents < 0)
498  {
499  return 0;
500  }
501  // if we run the particle generator and use 0 it'll run forever
502  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS)
503  {
504  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
505  cout << "it will run forever, so I just return without running anything" << endl;
506  return 0;
507  }
508 
509  se->skip(skip);
510  se->run(nEvents);
511 
512  //-----
513  // Exit
514  //-----
515 
516  se->End();
517  std::cout << "All done" << std::endl;
518  delete se;
519  if (Enable::PRODUCTION)
520  {
522  }
523 
524  gSystem->Exit(0);
525  return 0;
526 }
527 #endif