Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_sPHENIX_ForTrackCutStudy.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_sPHENIX_ForTrackCutStudy.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_Centrality.C>
11 #include <G4_DSTReader.C>
12 #include <G4_Global.C>
13 #include <G4_HIJetReco.C>
14 #include <G4_Input.C>
15 #include <G4_Jets.C>
16 #include <G4_KFParticle.C>
17 #include <G4_ParticleFlow.C>
18 #include <G4_Production.C>
19 #include <G4_TopoClusterReco.C>
20 #include <G4_Tracking.C>
21 #include <G4_User.C>
22 #include <QA.C>
23 
24 #include <ffamodules/FlagHandler.h>
25 #include <ffamodules/HeadReco.h>
26 #include <ffamodules/SyncReco.h>
27 #include <ffamodules/XploadInterface.h>
28 
31 #include <fun4all/Fun4AllServer.h>
32 
33 #include <phool/PHRandomSeed.h>
34 #include <phool/recoConsts.h>
35 
36 R__LOAD_LIBRARY(libfun4all.so)
37 R__LOAD_LIBRARY(libffamodules.so)
38 
39 // For HepMC Hijing
40 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
41 
43  const int nEvents = 1,
44  const string &inputFile = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/fm_0_20/g4hits/G4Hits_sHijing_0_20fm-0000000040-00848.root",
45  const string &outputFile = "sPhenixG4_test00019.root",
46  const string &embed_input_file = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/fm_0_20/g4hits/G4Hits_sHijing_0_20fm-0000000040-00848.root",
47  const int skip = 0,
48  const string &outdir = ".")
49 {
51  se->Verbosity(0);
52 
53  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
55 
56  // just if we set some flags somewhere in this macro
58  // By default every random number generator uses
59  // PHRandomSeed() which reads /dev/urandom to get its seed
60  // if the RANDOMSEED flag is set its value is taken as seed
61  // You can either set this to a random value using PHRandomSeed()
62  // which will make all seeds identical (not sure what the point of
63  // this would be:
64  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
65  // or set it to a fixed value so you can debug your code
66  // rc->set_IntFlag("RANDOMSEED", 12345);
67 
68 
69  //===============
70  // Input options
71  //===============
72  // verbosity setting (applies to all input managers)
73  Input::VERBOSITY = 0;
74  // First enable the input generators
75  // Either:
76  // read previously generated g4-hits files, in this case it opens a DST and skips
77  // the simulations step completely. The G4Setup macro is only loaded to get information
78  // about the number of layers used for the cell reco code
79  // Input::READHITS = true;
80  INPUTREADHITS::filename[0] = inputFile;
81  // if you use a filelist
82  // INPUTREADHITS::listfile[0] = inputFile;
83  // Or:
84  // Use particle generator
85  // And
86  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
87  // 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
88  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
89  // Input::EMBED = true;
90  INPUTEMBED::filename[0] = embed_input_file;
91  // if you use a filelist
92  //INPUTEMBED::listfile[0] = embed_input_file;
93 
94  Input::SIMPLE = true;
95  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
96  // Input::SIMPLE_VERBOSITY = 1;
97 
98  // Input::PYTHIA6 = true;
99 
100  // Input::PYTHIA8 = true;
101 
102  // Input::GUN = true;
103  // Input::GUN_NUMBER = 3; // if you need 3 of them
104  // Input::GUN_VERBOSITY = 1;
105 
106  //D0 generator
107  //Input::DZERO = false;
108  //Input::DZERO_VERBOSITY = 0;
109  //Lambda_c generator //Not ready yet
110  //Input::LAMBDAC = false;
111  //Input::LAMBDAC_VERBOSITY = 0;
112  // Upsilon generator
113  //Input::UPSILON = true;
114  //Input::UPSILON_NUMBER = 3; // if you need 3 of them
115  //Input::UPSILON_VERBOSITY = 0;
116 
117  // Input::HEPMC = true;
118  INPUTHEPMC::filename = inputFile;
119 
120  // Event pile up simulation with collision rate in Hz MB collisions.
121  //Input::PILEUPRATE = 100e3;
122 
123  //-----------------
124  // Initialize the selected Input/Event generation
125  //-----------------
126  // This creates the input generator(s)
127  InputInit();
128 
129  //--------------
130  // Set generator specific options
131  //--------------
132  // can only be set after InputInit() is called
133 
134  // Simple Input generator:
135  // if you run more than one of these Input::SIMPLE_NUMBER > 1
136  // add the settings for other with [1], next with [2]...
137  if (Input::SIMPLE)
138  {
139  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 5);
140  if (Input::HEPMC || Input::EMBED)
141  {
142  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
143  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
144  }
145  else
146  {
147  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
150  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
151  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0.01, 0.01, 5.);
152  }
153  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
154  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
155  INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
156  }
157  // Upsilons
158  // if you run more than one of these Input::UPSILON_NUMBER > 1
159  // add the settings for other with [1], next with [2]...
160  if (Input::UPSILON)
161  {
162  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
163  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
164  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
165  // Y species - select only one, last one wins
166  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
167  if (Input::HEPMC || Input::EMBED)
168  {
169  INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
170  INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
171  }
172  }
173  // particle gun
174  // if you run more than one of these Input::GUN_NUMBER > 1
175  // add the settings for other with [1], next with [2]...
176  if (Input::GUN)
177  {
178  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
179  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
180  }
181 
182  // pythia6
183  if (Input::PYTHIA6)
184  {
187  }
188  // pythia8
189  if (Input::PYTHIA8)
190  {
193  }
194 
195  //--------------
196  // Set Input Manager specific options
197  //--------------
198  // can only be set after InputInit() is called
199 
200  if (Input::HEPMC)
201  {
204 
205  // optional overriding beam parameters
206  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0); //optional collision smear in space, time
207  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
208  // //optional choice of vertex distribution function in space, time
209  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
214  //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
215  if (Input::PILEUPRATE > 0)
216  {
217  // Copy vertex settings from foreground hepmc input
219  // and then modify the ones you want to be different
220  // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
221  }
222  }
223  if (Input::PILEUPRATE > 0)
224  {
227  }
228  // register all input generators with Fun4All
229  InputRegister();
230 
231  if (! Input::READHITS)
232  {
233  rc->set_IntFlag("RUNNUMBER",1);
234 
235  SyncReco *sync = new SyncReco();
236  se->registerSubsystem(sync);
237 
238  HeadReco *head = new HeadReco();
239  se->registerSubsystem(head);
240  }
241 // Flag Handler is always needed to read flags from input (if used)
242 // and update our rc flags with them. At the end it saves all flags
243 // again on the DST in the Flags node under the RUN node
244  FlagHandler *flag = new FlagHandler();
245  se->registerSubsystem(flag);
246 
247  // set up production relatedstuff
248  // Enable::PRODUCTION = true;
249 
250  //======================
251  // Write the DST
252  //======================
253 
254  //Enable::DSTOUT = true;
255  Enable::DSTOUT_COMPRESS = false;
256  DstOut::OutputDir = outdir;
257  DstOut::OutputFile = outputFile;
258 
259  //Option to convert DST to human command readable TTree for quick poke around the outputs
260  // Enable::DSTREADER = true;
261 
262  // turn the display on (default off)
263  //Enable::DISPLAY = true;
264 
265  //======================
266  // What to run
267  //======================
268 
269  // QA, main switch
270  //Enable::QA = true;
271  Enable::QA = false; // no need [Derek, 12.03.2022]
272 
273  // Global options (enabled for all enables subsystems - if implemented)
274  // Enable::ABSORBER = true;
275  // Enable::OVERLAPCHECK = true;
276  // Enable::VERBOSITY = 1;
277 
278  // Enable::BBC = true;
279  // Enable::BBC_SUPPORT = true; // save hist in bbc support structure
280  Enable::BBCFAKE = true; // Smeared vtx and t0, use if you don't want real BBC in simulation
281 
282  Enable::PIPE = true;
283  Enable::PIPE_ABSORBER = true;
284 
285  // central tracking
286  Enable::MVTX = true;
290 
291  Enable::INTT = true;
292 // Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
293 // Enable::INTT_SUPPORT = true; // enable global support structure readout
296  Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;
297 
298  Enable::TPC = true;
299  Enable::TPC_ABSORBER = true;
300  Enable::TPC_CELL = Enable::TPC && true;
302  Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;
303 
304  Enable::MICROMEGAS = true;
307  Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;
308 
311  Enable::TRACKING_QA = Enable::TRACKING_TRACK && Enable::QA && true;
312 
313  // cemc electronics + thin layer of W-epoxy to get albedo from cemc
314  // into the tracking, cannot run together with CEMC
315  // Enable::CEMCALBEDO = true;
316 
317  Enable::CEMC = true;
318  Enable::CEMC_ABSORBER = true;
322  //Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
323  Enable::CEMC_EVAL = false; // no need [Derek, 12.03.2022]
324  Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && true;
325 
326  Enable::HCALIN = true;
331  //Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
332  Enable::HCALIN_EVAL = false; // no need [Derek, 12.03.2022]
333  Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && true;
334 
335  Enable::MAGNET = true;
337 
338  Enable::HCALOUT = true;
343  //Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
344  Enable::HCALOUT_EVAL = false; // no need [Derek, 12.03.2022]
345  Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && true;
346 
347  Enable::EPD = true;
348 
349  Enable::BEAMLINE = true;
350 // Enable::BEAMLINE_ABSORBER = true; // makes the beam line magnets sensitive volumes
351 // Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
352  Enable::ZDC = true;
353 // Enable::ZDC_ABSORBER = true;
354 // Enable::ZDC_SUPPORT = true;
355  Enable::ZDC_TOWER = Enable::ZDC && true;
356  //Enable::ZDC_EVAL = Enable::ZDC_TOWER && true;
357  Enable::ZDC_EVAL = false; // no need [Derek, 12.03.2022]
358 
360  //Enable::PLUGDOOR = true;
362 
363  Enable::GLOBAL_RECO = (Enable::BBCFAKE || Enable::TRACKING_TRACK) && true;
364  //Enable::GLOBAL_FASTSIM = true;
365 
366  //Enable::KFPARTICLE = true;
367  //Enable::KFPARTICLE_VERBOSITY = 1;
368  //Enable::KFPARTICLE_TRUTH_MATCH = true;
369  //Enable::KFPARTICLE_SAVE_NTUPLE = true;
370 
372 
374  //Enable::JETS_EVAL = Enable::JETS && true;
375  Enable::JETS_EVAL = false; // no need [Derek, 12.03.2022]
376  Enable::JETS_QA = Enable::JETS && Enable::QA && true;
377 
378  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
379  // single particle / p+p-only simulations, or for p+Au / Au+Au
380  // simulations which don't particularly care about jets)
381  Enable::HIJETS = Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
382 
383  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
384  Enable::TOPOCLUSTER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
385  // particle flow jet reconstruction - needs topoClusters!
387  // centrality reconstruction
388  Enable::CENTRALITY = true;
389 
390  // new settings using Enable namespace in GlobalVariables.C
391  Enable::BLACKHOLE = true;
392  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
393  //Enable::BLACKHOLE_FORWARD_SAVEHITS = false; // disable forward/backward hits
394  //BlackHoleGeometry::visible = true;
395 
396  // run user provided code (from local G4_User.C)
397  //Enable::USER = true;
398 
399  //===============
400  // conditions DB flags
401  //===============
402  //Enable::XPLOAD = true;
403  // tag
404  rc->set_StringFlag("XPLOAD_TAG",XPLOAD::tag);
405  // database config
406  rc->set_StringFlag("XPLOAD_CONFIG",XPLOAD::config);
407  // 64 bit timestamp
408  rc->set_uint64Flag("TIMESTAMP",XPLOAD::timestamp);
409 
410  //---------------
411  // World Settings
412  //---------------
413  // G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
414  // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
415 
416  //---------------
417  // Magnet Settings
418  //---------------
419 
420  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root"); // default map from the calibration database
421  // G4MAGNET::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)
422 // G4MAGNET::magfield_rescale = 1.; // make consistent with expected Babar field strength of 1.4T
423 
424  //---------------
425  // Pythia Decayer
426  //---------------
427  // list of decay types in
428  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
429  // default is All:
430  // G4P6DECAYER::decayType = EDecayType::kAll;
431 
432  // Initialize the selected subsystems
433  G4Init();
434 
435  //---------------------
436  // GEANT4 Detector description
437  //---------------------
438  if (!Input::READHITS)
439  {
440  G4Setup();
441  }
442 
443  //------------------
444  // Detector Division
445  //------------------
446 
447  if (Enable::BBC || Enable::BBCFAKE) Bbc_Reco();
448 
453 
455 
457 
459 
460  //-----------------------------
461  // CEMC towering and clustering
462  //-----------------------------
463 
466 
467  //-----------------------------
468  // HCAL towering and clustering
469  //-----------------------------
470 
473 
474  if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
476 
477  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
479 
480  //--------------
481  // SVTX tracking
482  //--------------
484  {
485  TrackingInit();
486  }
491 
493  {
494  Tracking_Reco();
495  }
496  //-----------------
497  // Global Vertexing
498  //-----------------
499 
501  {
502  cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
503  gSystem->Exit(1);
504  }
506  {
507  Global_Reco();
508  }
509  else if (Enable::GLOBAL_FASTSIM)
510  {
511  Global_FastSim();
512  }
513 
514  //-----------------
515  // Centrality Determination
516  //-----------------
517 
518  if (Enable::CENTRALITY)
519  {
520  Centrality();
521  }
522 
523  //-----------------
524  // Calo Trigger Simulation
525  //-----------------
526 
528  {
529  CaloTrigger_Sim();
530  }
531 
532  //---------
533  // Jet reco
534  //---------
535 
536  if (Enable::JETS) Jet_Reco();
537  if (Enable::HIJETS) HIJetReco();
538 
540 
541  //----------------------
542  // Simulation evaluation
543  //----------------------
544  string outputroot = outputFile;
545  string remove_this = ".root";
546  size_t pos = outputroot.find(remove_this);
547  if (pos != string::npos)
548  {
549  outputroot.erase(pos, remove_this.length());
550  }
551 
552  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
553 
554  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
555 
556  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
557 
558  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
559 
560  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
561 
562  if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
563 
565 
566  //======================
567  // Run KFParticle on evt
568  //======================
571 
572  //----------------------
573  // Standard QAs
574  //----------------------
575 
576  if (Enable::CEMC_QA) CEMC_QA();
579 
580  if (Enable::JETS_QA) Jet_QA();
581 
582  if (Enable::MVTX_QA) Mvtx_QA();
583  if (Enable::INTT_QA) Intt_QA();
584  if (Enable::TPC_QA) TPC_QA();
587 
589 
590  //--------------
591  // Set up Input Managers
592  //--------------
593 
594  InputManagers();
595 
596  if (Enable::PRODUCTION)
597  {
599  }
600 
601  if (Enable::DSTOUT)
602  {
603  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
604  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
606  {
607  ShowerCompress();
608  DstCompress(out);
609  }
610  se->registerOutputManager(out);
611  }
612  //-----------------
613  // Event processing
614  //-----------------
615  if (Enable::DISPLAY)
616  {
617  DisplayOn();
618 
619  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
620  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
621 
622  cout << "-------------------------------------------------" << endl;
623  cout << "You are in event display mode. Run one event with" << endl;
624  cout << "se->run(1)" << endl;
625  cout << "Run Geant4 command with following examples" << endl;
626  gROOT->ProcessLine("displaycmd()");
627 
628  return 0;
629  }
630 
631  // if we use a negative number of events we go back to the command line here
632  if (nEvents < 0)
633  {
634  return 0;
635  }
636  // if we run the particle generator and use 0 it'll run forever
637  // for embedding it runs forever if the repeat flag is set
638  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
639  {
640  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
641  cout << "it will run forever, so I just return without running anything" << endl;
642  return 0;
643  }
644 
645  se->skip(skip);
646  se->run(nEvents);
647 
648  //-----
649  // QA output
650  //-----
651 
652  if (Enable::QA) QA_Output(outputroot + "_qa.root");
653 
654  //-----
655  // Exit
656  //-----
657 
658  se->End();
659  std::cout << "All done" << std::endl;
660  delete se;
661  if (Enable::PRODUCTION)
662  {
664  }
665 
666  gSystem->Exit(0);
667  return 0;
668 }
669 #endif