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