ECCE @ 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>
9 #include <G4_EventEvaluator.C>
10 #include <G4_FwdJets.C>
11 #include <G4_Global.C>
12 #include <G4_Input.C>
13 #include <G4_Production.C>
14 #include <G4_User.C>
15 
16 #include <TROOT.h>
19 #include <fun4all/Fun4AllServer.h>
20 
21 #include <phool/recoConsts.h>
22 
23 #include <RooUnblindPrecision.h>
24 
25 R__LOAD_LIBRARY(libfun4all.so)
26 
28  const int nEvents = 1,
29  const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
30  const string &outputFile = "G4EICDetector.root",
31  const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
32  const int skip = 0,
33  const string &outdir = ".")
34 {
35  //---------------
36  // Fun4All server
37  //---------------
39  se->Verbosity(0);
40  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
41  //PHRandomSeed::Verbosity(1);
42 
43  // just if we set some flags somewhere in this macro
45  // By default every random number generator uses
46  // PHRandomSeed() which reads /dev/urandom to get its seed
47  // if the RANDOMSEED flag is set its value is taken as initial seed
48  // which will produce identical results so you can debug your code
49  // rc->set_IntFlag("RANDOMSEED", 12345);
50 
51  bool generate_seed = false;
52 
53  if (generate_seed)
54  {
55  size_t findSlash = inputFile.find_last_of("/");
56  string inputFileName = inputFile.substr(findSlash + 1, inputFile.size());
57 
58  RooRealVar dummyVal("dummy", "", 0);
59  RooUnblindPrecision blindVal("blindVal", "blindVal", inputFileName.c_str(), nEvents, skip + 1, dummyVal, kFALSE);
60  rc->set_IntFlag("RANDOMSEED", abs(ceil(blindVal.getVal() * 1e2)));
61  }
62 
63  //===============
64  // Input options
65  //===============
66 
67  // switching IPs by comment/uncommenting the following lines
68  // used for both beamline setting and for the event generator crossing boost
69  Enable::IP6 = true;
70  //Enable::IP8 = true;
71 
72 
73  //===============
74  // The following Ion energy and electron energy setting needs to be speficied
75  // The setting options for e-p high divergence setting (p energy x e energy):
76  // Option: 275x18, 275x10, 100x10, 100x5, 41x5
77  //
78  // The setting options for e-p high divergence setting (p energy x e energy):
79  // Option: 275x18, 275x10, 100x10, 100x5, 41x5
80  //
81  // The setting options for e-p high divergence setting (A energy x e energy):
82  // Option: 110x18, 110x10, 110x5, 41x5
83 
84  // Setting proton beam pipe energy. If you don't know what to set here, leave it at 275
86 
87  // Setting electron beam pipe energy. If you don't know what to set here, leave it at 18
89 
90  // Beam Scattering configuration setting specified by CDR
91  //
92  // Option 1: ep-high-acceptance
93  // Option 2: ep-high-divergence
94  // Option 3: eA
95  //
96  // Enable::BEAM_COLLISION_SETTING = "ep-high-divergence";
97  // If you don't know what to put here, set it to ep-high-divergence
98  //
99  // Enable::BEAM_COLLISION_SETTING = "eA";
100  Enable::BEAM_COLLISION_SETTING = "ep-high-divergence";
101 
102  // Either:
103  // read previously generated g4-hits files, in this case it opens a DST and skips
104  // the simulations step completely. The G4Setup macro is only loaded to get information
105  // about the number of layers used for the cell reco code
106  //
107  //Input::READHITS = true;
108  INPUTREADHITS::filename[0] = inputFile;
109  // if you use a filelist
110  // INPUTREADHITS::listfile[0] = inputFile;
111 
112  // Or:
113  // Use one or more particle generators
114  // It is run if Input::<generator> is set to true
115  // all other options only play a role if it is active
116  // 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
117  // Input::EMBED = true;
118  INPUTEMBED::filename[0] = embed_input_file;
119  // if you use a filelist
120  //INPUTEMBED::listfile[0] = embed_input_file;
121 
122  // Use Pythia 8
123  // Input::PYTHIA8 = true;
124 
125  // Use Pythia 6
126  // Input::PYTHIA6 = true;
127 
128  // Use Sartre
129  // Input::SARTRE = true;
130 
131  // Simple multi particle generator in eta/phi/pt ranges
132  Input::SIMPLE = true;
133  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
134  // Input::SIMPLE_VERBOSITY = 1;
135 
136  // Particle gun (same particles in always the same direction)
137  // Input::GUN = true;
138  // Input::GUN_NUMBER = 3; // if you need 3 of them
139  // Input::GUN_VERBOSITY = 0;
140 
141  // Particle ion gun
142  // Input::IONGUN = true;
143 
144  // Upsilon generator
145  // Input::UPSILON = true;
146  // Input::UPSILON_NUMBER = 3; // if you need 3 of them
147  // Input::UPSILON_VERBOSITY = 0;
148 
149  // And/Or read generated particles from file
150 
151  // eic-smear output
152  // Input::READEIC = true;
153  INPUTREADEIC::filename = inputFile;
154 
155  // HepMC2 files
156  // Input::HEPMC = true;
157  Input::VERBOSITY = 0;
158  INPUTHEPMC::filename = inputFile;
159 
160  //-----------------
161  // Initialize the selected Input/Event generation
162  //-----------------
163  InputInit();
164  //--------------
165  // Set generator specific options
166  //--------------
167  // can only be set after InputInit() is called
168 
169  // Simple Input generator:
170  // if you run more than one of these Input::SIMPLE_NUMBER > 1
171  // add the settings for other with [1], next with [2]...
172  if (Input::SIMPLE)
173  {
174  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 5);
175  if (Input::HEPMC || Input::EMBED)
176  {
177  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
178  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
179  }
180  else
181  {
182  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
185  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
186  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 5.);
187  }
188  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-3, 3);
189  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
190  INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
191  }
192  // Upsilons
193  // if you run more than one of these Input::UPSILON_NUMBER > 1
194  // add the settings for other with [1], next with [2]...
195  if (Input::UPSILON)
196  {
197  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("mu", 0);
198  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
199  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
200  // Y species - select only one, last one wins
201  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
202  if (Input::HEPMC || Input::EMBED)
203  {
204  INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
205  INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
206  }
207  }
208  // particle gun
209  // if you run more than one of these Input::GUN_NUMBER > 1
210  // add the settings for other with [1], next with [2]...
211  if (Input::GUN)
212  {
213  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
214  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
215  }
216 
217  if (Input::IONGUN)
218  {
219  float theta = -25e-3;
220 
221  INPUTGENERATOR::IonGun[0]->SetA(197);
222  INPUTGENERATOR::IonGun[0]->SetZ(79);
223  INPUTGENERATOR::IonGun[0]->SetCharge(79);
224  INPUTGENERATOR::IonGun[0]->SetMom(sin(theta)*110*197, 0,cos(theta)*110*197); // -25 mrad
225 
226  INPUTGENERATOR::IonGun[0]->Print();
227 
228  }
229 
230  // pythia6
231  if (Input::PYTHIA6)
232  {
233  INPUTGENERATOR::Pythia6->set_config_file(string(getenv("CALIBRATIONROOT")) + "/Generators/phpythia6_ep.cfg");
236  }
237  // pythia8
238  if (Input::PYTHIA8)
239  {
242  }
243  // Sartre
244  if (Input::SARTRE)
245  {
248  }
249 
250  //--------------
251  // Set Input Manager specific options
252  //--------------
253  // can only be set after InputInit() is called
254 
255  if (Input::HEPMC)
256  {
259  // optional overriding beam parameters
260  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 30, 0); //optional collision smear in space, time
261  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
262  // //optional choice of vertex distribution function in space, time
263  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
268  //INPUTMANAGER::HepMCInputManager->set_embedding_id(2);
269  }
270 
271  // register all input generators with Fun4All
272  InputRegister();
273 
274  // Reads event generators in EIC smear files, which is registered in InputRegister
275  if (Input::READEIC)
276  {
280  }
281 
282  // set up production relatedstuff
283  // Enable::PRODUCTION = true;
284 
285  //======================
286  // Write the DST
287  //======================
288 
289  // Enable::DSTOUT = true;
290  DstOut::OutputDir = outdir;
291  DstOut::OutputFile = outputFile;
292  Enable::DSTOUT_COMPRESS = true; // Compress DST files
293 
294  //Option to convert DST to human command readable TTree for quick poke around the outputs
295  // Enable::DSTREADER = true;
296 
297  // turn the display on (default off)
298  // Enable::DISPLAY = true;
299 
300  //======================
301  // What to run
302  //======================
303  // Global options (enabled for all subsystems - if implemented)
304  // Enable::ABSORBER = true;
305  // Enable::OVERLAPCHECK = true;
306  // Enable::VERBOSITY = 1;
307 
308  // whether to simulate the Be section of the beam pipe
309  Enable::PIPE = true;
310  // If need to disable EIC beam pipe extension beyond the Be-section:
312  //EIC hadron far forward magnets and detectors. IP6 and IP8 are incompatible (pick either or);
315 
318 
319  // gems
320  Enable::EGEM = false;
321  Enable::FGEM = false; // deactivated as it's replaced by a FTTL layer
322  // Enable::BGEM = true; // not yet defined in this model
323  Enable::RWELL = true;
324  // barrel tracker
326  // Enable::TrackingService_VERBOSITY = INT_MAX - 10;
327  Enable::BARREL = true;
328  // fst
329  Enable::FST = true;
330 
331  //AC-LGAD TOFs
332  Enable::FTTL = true;
333  Enable::ETTL = true;
334  Enable::CTTL = true;
335 
336  //mRPC TOFs
337  Enable::BTOF = false;
338  Enable::ETOF = false;
339  Enable::HTOF = false;
340  Enable::ETOF_GAS = Enable::ETOF && true;
341  Enable::HTOF_GAS = Enable::HTOF && true;
342 
343  Enable::TRACKING = true;
345  G4TRACKING::DISPLACED_VERTEX = true; // this option exclude vertex in the track fitting and use RAVE to reconstruct primary and 2ndary vertexes
346  // projections to calorimeters
356 
357  // enable barrel calos & magnet
358  Enable::BECAL = true;
359  Enable::HCALIN = true;
360  Enable::MAGNET = true;
361  Enable::HCALOUT = true;
362 
363  // EICDetector geometry - barrel
364  Enable::DIRC = true;
366 
367  Enable::BMMG = false;
368  // Enable::DIRC_VERBOSITY = 2;
369 
370  // EICDetector geometry - 'hadron' direction
371  Enable::RICH = true;
373 
374  Enable::TRD = false;
375  Enable::TRD_GAS = false;
376  // Enable::RICH_VERBOSITY = 2;
377 
378 
379  // enable forward calos
380  Enable::FEMC = true;
381  Enable::DRCALO = false;
382  Enable::LFHCAL = true;
383 
384  // EICDetector geometry - 'electron' direction
385  Enable::mRICH = true;
387  // Enable::mRICH_VERBOSITY = 2;
388 
389  // EICDetector geometry - 'electron' direction
390  Enable::EEMCH = true;
391  G4EEMCH::SETTING::USECUSTOMMAPUPDATED = true; // enable proper carbon structure
393  Enable::EHCAL = false;
394 
396 
397  Enable::PLUGDOOR = false;
398 
399  // Other options
400  Enable::GLOBAL_RECO = G4TRACKING::DISPLACED_VERTEX; // use reco vertex for global event vertex
401  Enable::GLOBAL_FASTSIM = true;
402 
403  // jet reconstruction
404  Enable::FWDJETS = true;
406 
407  // new settings using Enable namespace in GlobalVariables.C
408  Enable::BLACKHOLE = true;
409  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
410  //BlackHoleGeometry::visible = true;
411 
412  // ZDC
413  // Enable::ZDC = true;
414  // Enable::ZDC_DISABLE_BLACKHOLE = true;
415 
416  // B0
417  // B0 shape
418  // Enable::B0_DISABLE_HITPLANE = true;
419  // Enable::B0_FULLHITPLANE = true;
420  // Enable::B0_VAR_PIPE_HOLE = true;
421  // Enable::B0_CIRCLE_PIPE_HOLE = false;
422 
423  // B0 Tracking
424  // Enable::B0TRACKING = false; // For B0 Tracking
425  // Enable::B0TRACKING_EVAL = Enable::B0TRACKING && true; //For B0 Tracking
426 
427  // B0 ECAL
428  // Enable::B0ECALTOWERS = true; //To Construct Towers of B0ECal instead of one single volume
429  // Enable::B0ECAL = Enable::B0_DISABLE_HITPLANE && true;
430  // Enable::B0ECAL_CELL = Enable::B0ECAL && true;
431  // Enable::B0ECAL_TOWER = Enable::B0ECAL_CELL && true;
432  // Enable::B0ECAL_CLUSTER = Enable::B0ECAL_TOWER && true;
433  // Enable::B0ECAL_EVAL = Enable::B0ECAL_CLUSTER && true;
434 
435  // RP
436  // Enable::RP_DISABLE_HITPLANE = true;
437 
438  //Far Backward detectors
439 // Enable::BWD = true;
440 // Enable::BWDN[0]=true; // 1st Q2 tagger
441 // Enable::BWDN[1]=true; // 2nd Q2 tagger
442 // Enable::BWDN[2]=true; // 1st Lumi
443 // Enable::BWDN[3]=true; // 2nd Lumi (+)
444 // Enable::BWDN[4]=true; // 3rd Lumi (-)
445 // Enable::BWD_CELL = Enable::BWD && true;
446 // Enable::BWD_TOWER = Enable::BWD_CELL && true;
447 // Enable::BWD_CLUSTER = Enable::BWD_TOWER && true;
448 // Enable::BWD_EVAL = Enable::BWD_CLUSTER && true;
449 
450  //************************************************************
451  // details for calos: cells, towers, clusters
452  //************************************************************
457 
458  // Enable::HCALIN_ABSORBER = true;
463 
464  // Enable::HCALOUT_ABSORBER = true;
469 
470  // Enable::FEMC_ABSORBER = true;
474 
479 
480  Enable::LFHCAL_ABSORBER = false;
485 
489 
494 
495  // Enabling the event evaluator?
496  Enable::EVENT_EVAL = true;
500 
501  //Enable::USER = true;
502 
503  //---------------
504  // World Settings
505  //---------------
506  // G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
507  // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
508  // G4WORLD::WorldMaterial = "G4_Galactic"; // set to G4_GALACTIC for material scans
509 
510  //---------------
511  // Magnet Settings
512  //---------------
513 
514  // 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)
515  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
516  G4MAGNET::magfield_rescale = -1.4 / 1.5; // make consistent with expected Babar field strength of 1.4T
517 
518  //---------------
519  // Pythia Decayer
520  //---------------
521  // list of decay types in
522  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
523  // default is All:
524  // G4P6DECAYER::decayType = EDecayType::kAll;
525 
526  // Initialize the selected subsystems
527  G4Init();
528 
529  //---------------------
530  // GEANT4 Detector description
531  //---------------------
532 
533  // 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
534  if (!Input::READHITS)
535  {
536  G4Setup();
537  }
538 
539  //------------------
540  // Detector Division
541  //------------------
543 
545 
547 
548  //-----------------------------
549  // CEMC towering and clustering
550  //-----------------------------
551 
554 
555  //-----------------------------
556  // HCAL towering and clustering
557  //-----------------------------
558 
561 
564 
565  //-----------------------------
566  // e, h direction Calorimeter towering and clustering
567  //-----------------------------
568 
571 
574 
577 
580 
583 
586 
589 
592 
593  if (Enable::B0ECAL_TOWER) B0ECAL_Towers(); // For B0Ecal
594  if (Enable::B0ECAL_CLUSTER) B0ECAL_Clusters(); //For B0Ecal
595 
596  if (Enable::BWD_TOWER) BWD_Towers(); // For Bwd
597  if (Enable::BWD_CLUSTER) BWD_Clusters(); //For Bwd
598 
600 
601  //--------------
602  // Tracking and PID
603  //--------------
604 
606 
608 
610 
612 
613  //-----------------
614  // Global Vertexing
615  //-----------------
616 
618  {
619  Global_Reco();
620  }
621  else if (Enable::GLOBAL_FASTSIM)
622  {
623  Global_FastSim();
624  }
625 
626  //---------
627  // Jet reco
628  //---------
629 
631 
632  string outputroot = outputFile;
633  string remove_this = ".root";
634  size_t pos = outputroot.find(remove_this);
635  if (pos != string::npos)
636  {
637  outputroot.erase(pos, remove_this.length());
638  }
639 
640  if (Enable::DSTREADER) G4DSTreader_EICDetector(outputroot + "_DSTReader.root");
641 
642  //----------------------
643  // Simulation evaluation
644  //----------------------
645 
646  if (Enable::EVENT_EVAL) Event_Eval(outputroot + "_eventtree.root");
647 
648  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4tracking_eval.root");
649 
650  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
651 
652  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
653 
654  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
655 
656  if (Enable::FEMC_EVAL) FEMC_Eval(outputroot + "_g4femc_eval.root");
657 
658  if (Enable::FHCAL_EVAL) FHCAL_Eval(outputroot + "_g4fhcal_eval.root");
659 
660  if (Enable::EEMC_EVAL) EEMC_Eval(outputroot + "_g4eemc_eval.root");
661 
662  if (Enable::FFR_EVAL) FFR_Eval(outputroot + "_g4ffr_eval.root");
663 
664  if (Enable::B0ECAL_EVAL) B0ECAL_Eval(outputroot + "_g4b0ecal_eval_test.root"); // For B0Ecal
665 
666  if (Enable::BWD_EVAL) BWD_Eval(outputroot + "_g4bwd_eval_e0100_debug"); // For Bwd
667 
669 
671 
672  //--------------
673  // Set up Input Managers
674  //--------------
675 
676  InputManagers();
677 
678  //--------------
679  // Set up Output Manager
680  //--------------
681  if (Enable::PRODUCTION)
682  {
684  }
685 
686  if (Enable::DSTOUT)
687  {
688  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
689  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
691  se->registerOutputManager(out);
692  }
693 
694  //-----------------
695  // Event processing
696  //-----------------
697  if (Enable::DISPLAY)
698  {
699  DisplayOn();
700 
701  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
702  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
703 
704  cout << "-------------------------------------------------" << endl;
705  cout << "You are in event display mode. Run one event with" << endl;
706  cout << "se->run(1)" << endl;
707  cout << "Run Geant4 command with following examples" << endl;
708  gROOT->ProcessLine("displaycmd()");
709 
710  return 0;
711  }
712  // if we use a negative number of events we go back to the command line here
713  if (nEvents < 0)
714  {
715  return 0;
716  }
717  // if we run any of the particle generators and use 0 it'll run forever
718  if (nEvents == 0 && !Input::READHITS && !Input::HEPMC && !Input::READEIC)
719  {
720  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
721  cout << "it will run forever, so I just return without running anything" << endl;
722  return 0;
723  }
724 
725  se->skip(skip);
726  se->run(nEvents);
727 
728  //-----
729  // Exit
730  //-----
731 
732  se->End();
733  std::cout << "All done" << std::endl;
734  delete se;
735  if (Enable::PRODUCTION)
736  {
738  }
739  gSystem->Exit(0);
740  return 0;
741 }
742 #endif