Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_EICDetector_AnaTutorial.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_EICDetector_AnaTutorial.C
1 #ifndef MACRO_FUN4ALLG4EICDETECTOR_C
2 #define MACRO_FUN4ALLG4EICDETECTOR_C
3 
4 #include <anatutorial/AnaTutorial.h>
5 
6 #include <GlobalVariables.C>
7 
8 #include <DisplayOn.C>
9 #include <G4Setup_EICDetector.C>
10 #include <G4_Bbc.C>
12 #include <G4_FwdJets.C>
13 #include <G4_Global.C>
14 #include <G4_Input.C>
15 #include <G4_Jets.C>
16 #include <G4_Production.C>
17 #include <G4_User.C>
18 
19 #include <TROOT.h>
22 #include <fun4all/Fun4AllServer.h>
23 
24 #include <phool/recoConsts.h>
25 
26 R__LOAD_LIBRARY(libfun4all.so)
27 R__LOAD_LIBRARY(libanatutorial.so)
28 
30  const int nEvents = 1,
31  const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.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  //---------------
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
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  // if you use a filelist
65  // INPUTREADHITS::listfile[0] = inputFile;
66 
67  // Or:
68  // Use one or more particle generators
69  // It is run if Input::<generator> is set to true
70  // all other options only play a role if it is active
71  // 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
72  // Input::EMBED = true;
73  INPUTEMBED::filename[0] = embed_input_file;
74  // if you use a filelist
75  //INPUTEMBED::listfile[0] = embed_input_file;
76 
77  // Use Pythia 8
78  // Input::PYTHIA8 = true;
79 
80  // Use Pythia 6
81  // Input::PYTHIA6 = true;
82 
83  // Use Sartre
84  // Input::SARTRE = true;
85 
86  // Simple multi particle generator in eta/phi/pt ranges
87  Input::SIMPLE = true;
88  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
89  // Input::SIMPLE_VERBOSITY = 1;
90 
91  // Particle gun (same particles in always the same direction)
92  // Input::GUN = true;
93  // Input::GUN_NUMBER = 3; // if you need 3 of them
94  // Input::GUN_VERBOSITY = 0;
95 
96  // Upsilon generator
97  // Input::UPSILON = true;
98  // Input::UPSILON_NUMBER = 3; // if you need 3 of them
99  // Input::UPSILON_VERBOSITY = 0;
100 
101  // And/Or read generated particles from file
102 
103  // eic-smear output
104  // Input::READEIC = true;
105  INPUTREADEIC::filename = inputFile;
106 
107  // HepMC2 files
108  // Input::HEPMC = true;
109  Input::VERBOSITY = 0;
110  INPUTHEPMC::filename = inputFile;
111 
112  //-----------------
113  // Initialize the selected Input/Event generation
114  //-----------------
115  InputInit();
116  //--------------
117  // Set generator specific options
118  //--------------
119  // can only be set after InputInit() is called
120 
121  // Simple Input generator:
122  // if you run more than one of these Input::SIMPLE_NUMBER > 1
123  // add the settings for other with [1], next with [2]...
124  if (Input::SIMPLE)
125  {
126  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 5);
127  if (Input::HEPMC || Input::EMBED)
128  {
129  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
130  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
131  }
132  else
133  {
134  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
137  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
138  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 5.);
139  }
140  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-3, 3);
141  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
142  INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
143  }
144  // Upsilons
145  // if you run more than one of these Input::UPSILON_NUMBER > 1
146  // add the settings for other with [1], next with [2]...
147  if (Input::UPSILON)
148  {
149  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("mu", 0);
150  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
151  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
152  // Y species - select only one, last one wins
153  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
154  if (Input::HEPMC || Input::EMBED)
155  {
156  INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
157  INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
158  }
159  }
160  // particle gun
161  // if you run more than one of these Input::GUN_NUMBER > 1
162  // add the settings for other with [1], next with [2]...
163  if (Input::GUN)
164  {
165  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
166  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
167  }
168  // pythia6
169  if (Input::PYTHIA6)
170  {
171  INPUTGENERATOR::Pythia6->set_config_file(string(getenv("CALIBRATIONROOT")) + "/Generators/phpythia6_ep.cfg");
174  }
175  // pythia8
176  if (Input::PYTHIA8)
177  {
180  }
181  // Sartre
182  if (Input::SARTRE)
183  {
185  Input::ApplyEICBeamParameter(INPUTGENERATOR::Sartre);
186  }
187 
188  //--------------
189  // Set Input Manager specific options
190  //--------------
191  // can only be set after InputInit() is called
192 
193  if (Input::HEPMC)
194  {
197  // optional overriding beam parameters
198  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 30, 0); //optional collision smear in space, time
199  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
200  // //optional choice of vertex distribution function in space, time
201  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
206  //INPUTMANAGER::HepMCInputManager->set_embedding_id(2);
207  }
208 
209  // register all input generators with Fun4All
210  InputRegister();
211 
212  // Reads event generators in EIC smear files, which is registered in InputRegister
213  if (Input::READEIC)
214  {
217  }
218 
219  // set up production relatedstuff
220  // Enable::PRODUCTION = true;
221 
222  //======================
223  // Write the DST
224  //======================
225 
226  Enable::DSTOUT = false;
227  DstOut::OutputDir = outdir;
228  DstOut::OutputFile = outputFile;
229  Enable::DSTOUT_COMPRESS = false; // Compress DST files
230 
231  //Option to convert DST to human command readable TTree for quick poke around the outputs
232 // Enable::DSTREADER = true;
233 
234  // turn the display on (default off)
235  Enable::DISPLAY = false;
236 
237  //======================
238  // What to run
239  //======================
240  // Global options (enabled for all subsystems - if implemented)
241  // Enable::ABSORBER = true;
242  // Enable::OVERLAPCHECK = true;
243  // Enable::VERBOSITY = 1;
244 
245  // Enable::BBC = true;
246  Enable::BBCFAKE = true; // Smeared vtx and t0, use if you don't want real BBC in simulation
247 
248  // whether to simulate the Be section of the beam pipe
249  Enable::PIPE = true;
250  // EIC beam pipe extension beyond the Be-section:
251  G4PIPE::use_forward_pipes = false;
252  //EIC hadron far forward magnets and detectors. IP6 and IP8 are incompatible (pick either or);
253  Enable::HFARFWD_MAGNETS_IP6=true;
254  Enable::HFARFWD_VIRTUAL_DETECTORS_IP6=true;
255  Enable::HFARFWD_MAGNETS_IP8=false;
256  Enable::HFARFWD_VIRTUAL_DETECTORS_IP8=false;
257 
258  // gems
259  Enable::EGEM = true;
260  Enable::FGEM = true;
261  Enable::FGEM_ORIG = false; //5 forward gems; cannot be used with FST
262  // barrel tracker
263  Enable::BARREL = false;
264  //G4BARREL::SETTING::BARRELV6=true;
265  // fst
266  Enable::FST = true;
267  G4FST::SETTING::FST_TPC = true;
268  // mvtx/tpc tracker
269  Enable::MVTX = true;
270  Enable::TPC = true;
271  // Enable::TPC_ENDCAP = true;
272 
273  Enable::TRACKING = true;
274  Enable::TRACKING_EVAL = Enable::TRACKING && true;
275  G4TRACKING::DISPLACED_VERTEX = false; // this option exclude vertex in the track fitting and use RAVE to reconstruct primary and 2ndary vertexes
276  // projections to calorimeters
277  G4TRACKING::PROJECTION_EEMC = false;
278  G4TRACKING::PROJECTION_CEMC = false;
279  G4TRACKING::PROJECTION_FEMC = false;
280  G4TRACKING::PROJECTION_FHCAL = false;
281 
282  Enable::CEMC = true;
283  // Enable::CEMC_ABSORBER = true;
288 
289  Enable::HCALIN = true;
290  // Enable::HCALIN_ABSORBER = true;
295 
296  Enable::MAGNET = true;
297 
298  Enable::HCALOUT = true;
299  // Enable::HCALOUT_ABSORBER = true;
304 
305  // EICDetector geometry - barrel
306  Enable::DIRC = true;
307 
308  // EICDetector geometry - 'hadron' direction
309  Enable::RICH = true;
310  Enable::AEROGEL = true;
311 
312  Enable::FEMC = true;
313  // Enable::FEMC_ABSORBER = true;
314  Enable::FEMC_TOWER = Enable::FEMC && true;
315  Enable::FEMC_CLUSTER = Enable::FEMC_TOWER && true;
316  Enable::FEMC_EVAL = Enable::FEMC_CLUSTER && true;
317 
318  Enable::FHCAL = true;
319  // Enable::FHCAL_ABSORBER = true;
320  Enable::FHCAL_TOWER = Enable::FHCAL && true;
321  Enable::FHCAL_CLUSTER = Enable::FHCAL_TOWER && true;
322  Enable::FHCAL_EVAL = Enable::FHCAL_CLUSTER && true;
323 
324  // EICDetector geometry - 'electron' direction
325  Enable::EEMC = true;
326  Enable::EEMC_TOWER = Enable::EEMC && true;
327  Enable::EEMC_CLUSTER = Enable::EEMC_TOWER && true;
328  Enable::EEMC_EVAL = Enable::EEMC_CLUSTER && true;
329 
330  Enable::PLUGDOOR = true;
331 
332  // Other options
333  Enable::GLOBAL_RECO = true;
334  Enable::GLOBAL_FASTSIM = true;
335 
336  // Select only one jet reconstruction- they currently use the same
337  // output collections on the node tree!
338  Enable::JETS = true;
340 
341  Enable::FWDJETS = true;
342  Enable::FWDJETS_EVAL = Enable::FWDJETS && true;
343 
344  // new settings using Enable namespace in GlobalVariables.C
345  Enable::BLACKHOLE = true;
346  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
347  //BlackHoleGeometry::visible = true;
348 
349  //Enable::USER = true;
350 
351  //---------------
352  // World Settings
353  //---------------
354  // G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
355  // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
356 
357  //---------------
358  // Magnet Settings
359  //---------------
360 
361  // 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)
362  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
363  G4MAGNET::magfield_rescale = -1.4 / 1.5; // make consistent with expected Babar field strength of 1.4T
364 
365  //---------------
366  // Pythia Decayer
367  //---------------
368  // list of decay types in
369  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
370  // default is All:
371  // G4P6DECAYER::decayType = EDecayType::kAll;
372 
373  // Initialize the selected subsystems
374  G4Init();
375 
376  //---------------------
377  // GEANT4 Detector description
378  //---------------------
379 
380  // 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
381  if (!Input::READHITS)
382  {
383  G4Setup();
384  }
385 
386  //------------------
387  // Detector Division
388  //------------------
389 
390  if (Enable::BBC || Enable::BBCFAKE) Bbc_Reco();
391 
393 
395 
397 
398  //-----------------------------
399  // CEMC towering and clustering
400  //-----------------------------
401 
404 
405  //-----------------------------
406  // HCAL towering and clustering
407  //-----------------------------
408 
411 
414 
415  //-----------------------------
416  // e, h direction Calorimeter towering and clustering
417  //-----------------------------
418 
419  if (Enable::FEMC_TOWER) FEMC_Towers();
420  if (Enable::FEMC_CLUSTER) FEMC_Clusters();
421 
422  if (Enable::FHCAL_TOWER) FHCAL_Towers();
423  if (Enable::FHCAL_CLUSTER) FHCAL_Clusters();
424 
425  if (Enable::EEMC_TOWER) EEMC_Towers();
426  if (Enable::EEMC_CLUSTER) EEMC_Clusters();
427 
429 
430  //--------------
431  // SVTX tracking
432  //--------------
433 
434  if (Enable::TRACKING) Tracking_Reco();
435 
436  //-----------------
437  // Global Vertexing
438  //-----------------
439 
441  {
442  Global_Reco();
443  }
444  else if (Enable::GLOBAL_FASTSIM)
445  {
446  Global_FastSim();
447  }
448 
449  //---------
450  // Jet reco
451  //---------
452 
453  if (Enable::JETS) Jet_Reco();
454 
455  if (Enable::FWDJETS) Jet_FwdReco();
456 
457  string outputroot = outputFile;
458  string remove_this = ".root";
459  size_t pos = outputroot.find(remove_this);
460  if (pos != string::npos)
461  {
462  outputroot.erase(pos, remove_this.length());
463  }
464 
465  if (Enable::DSTREADER) G4DSTreader_EICDetector(outputroot + "_DSTReader.root");
466 
467  //----------------------
468  // Simulation evaluation
469  //----------------------
470  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4tracking_eval.root");
471 
472  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
473 
474  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
475 
476  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
477 
478  if (Enable::FEMC_EVAL) FEMC_Eval(outputroot + "_g4femc_eval.root");
479 
480  if (Enable::FHCAL_EVAL) FHCAL_Eval(outputroot + "_g4fhcal_eval.root");
481 
482  if (Enable::EEMC_EVAL) EEMC_Eval(outputroot + "_g4eemc_eval.root");
483 
484  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
485 
486  if (Enable::FWDJETS_EVAL) Jet_FwdEval(outputroot + "_g4fwdjet_eval.root");
487 
489 
490  AnaTutorial *anaTutorial = new AnaTutorial("anaTutorial", outputroot + "_anaTutorial.root");
491  anaTutorial->setMinJetPt(3.);
492  anaTutorial->Verbosity(0);
493  anaTutorial->analyzeTracks(true);
494  anaTutorial->analyzeClusters(true);
495  anaTutorial->analyzeJets(true);
496  anaTutorial->analyzeTruth(false);
497  se->registerSubsystem(anaTutorial);
498 
499  //--------------
500  // Set up Input Managers
501  //--------------
502 
503  InputManagers();
504 
505  //--------------
506  // Set up Output Manager
507  //--------------
508  if (Enable::PRODUCTION)
509  {
511  }
512 
513  if (Enable::DSTOUT)
514  {
515  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
516  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
518  se->registerOutputManager(out);
519  }
520 
521  //-----------------
522  // Event processing
523  //-----------------
524  if (Enable::DISPLAY)
525  {
526  DisplayOn();
527 
528  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
529  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
530 
531  cout << "-------------------------------------------------" << endl;
532  cout << "You are in event display mode. Run one event with" << endl;
533  cout << "se->run(1)" << endl;
534  cout << "Run Geant4 command with following examples" << endl;
535  gROOT->ProcessLine("displaycmd()");
536 
537  return 0;
538  }
539  // if we use a negative number of events we go back to the command line here
540  if (nEvents < 0)
541  {
542  return 0;
543  }
544  // if we run any of the particle generators and use 0 it'll run forever
545  if (nEvents == 0 && !Input::READHITS && !Input::HEPMC && !Input::READEIC)
546  {
547  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
548  cout << "it will run forever, so I just return without running anything" << endl;
549  return 0;
550  }
551 
552  se->skip(skip);
553  se->run(nEvents);
554 
555  //-----
556  // Exit
557  //-----
558 
559  se->End();
560  std::cout << "All done" << std::endl;
561  delete se;
562  if (Enable::PRODUCTION)
563  {
565  }
566  gSystem->Exit(0);
567  return 0;
568 }
569 #endif