ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_AnaTutorial_sPHENIX_Jets.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_AnaTutorial_sPHENIX_Jets.C
1 #ifndef MACRO_FUN4ALLG4SPHENIX_C
2 #define MACRO_FUN4ALLG4SPHENIX_C
3 
4 #include <anatutorial/AnaTutorial.h>
5 
6 #include <GlobalVariables.C>
7 
8 #include <DisplayOn.C>
9 #include <G4Setup_sPHENIX.C>
10 #include <G4_Bbc.C>
11 #include <G4_CaloTrigger.C>
12 #include <G4_DSTReader.C>
13 #include <G4_Global.C>
14 #include <G4_HIJetReco.C>
15 #include <G4_Input.C>
16 #include <G4_Jets.C>
17 #include <G4_KFParticle.C>
18 #include <G4_ParticleFlow.C>
19 #include <G4_Production.C>
20 #include <G4_TopoClusterReco.C>
21 #include <G4_Tracking.C>
22 #include <G4_User.C>
23 #include <QA.C>
24 
27 #include <fun4all/Fun4AllServer.h>
28 
29 #include <phool/PHRandomSeed.h>
30 #include <phool/recoConsts.h>
31 
32 R__LOAD_LIBRARY(libfun4all.so)
33 R__LOAD_LIBRARY(libanatutorial.so)
34 
35 // For HepMC Hijing
36 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
37 
39  const int nEvents = 1,
40  const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
41  const string &outputFile = "G4sPHENIX.root",
42  const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
43  const int skip = 0,
44  const string &outdir = ".")
45 {
47  se->Verbosity(0);
48 
49  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
51 
52  // just if we set some flags somewhere in this macro
54  // By default every random number generator uses
55  // PHRandomSeed() which reads /dev/urandom to get its seed
56  // if the RANDOMSEED flag is set its value is taken as seed
57  // You can either set this to a random value using PHRandomSeed()
58  // which will make all seeds identical (not sure what the point of
59  // this would be:
60  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
61  // or set it to a fixed value so you can debug your code
62  // rc->set_IntFlag("RANDOMSEED", 12345);
63 
64  //===============
65  // Input options
66  //===============
67  // verbosity setting (applies to all input managers)
68  Input::VERBOSITY = 0;
69  // First enable the input generators
70  // Either:
71  // read previously generated g4-hits files, in this case it opens a DST and skips
72  // the simulations step completely. The G4Setup macro is only loaded to get information
73  // about the number of layers used for the cell reco code
74  // Input::READHITS = true;
75  INPUTREADHITS::filename[0] = inputFile;
76  // if you use a filelist
77  // INPUTREADHITS::listfile[0] = inputFile;
78  // Or:
79  // Use particle generator
80  // And
81  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
82  // 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
83  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
84  // Input::EMBED = true;
85  INPUTEMBED::filename[0] = embed_input_file;
86  // if you use a filelist
87  //INPUTEMBED::listfile[0] = embed_input_file;
88 
89  Input::SIMPLE = false;
90  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
91  // Input::SIMPLE_VERBOSITY = 1;
92 
93  // Input::PYTHIA6 = true;
94 
95  Input::PYTHIA8 = true;
96 
97  // Input::GUN = true;
98  // Input::GUN_NUMBER = 3; // if you need 3 of them
99  // Input::GUN_VERBOSITY = 1;
100 
101  //D0 generator
102  //Input::DZERO = false;
103  //Input::DZERO_VERBOSITY = 0;
104  //Lambda_c generator //Not ready yet
105  //Input::LAMBDAC = false;
106  //Input::LAMBDAC_VERBOSITY = 0;
107  // Upsilon generator
108  //Input::UPSILON = true;
109  //Input::UPSILON_NUMBER = 3; // if you need 3 of them
110  //Input::UPSILON_VERBOSITY = 0;
111 
112  // Input::HEPMC = true;
113  INPUTHEPMC::filename = inputFile;
114 
115  // Event pile up simulation with collision rate in Hz MB collisions.
116  //Input::PILEUPRATE = 100e3;
117 
118  //-----------------
119  // Initialize the selected Input/Event generation
120  //-----------------
121  // This creates the input generator(s)
122  InputInit();
123 
124  //--------------
125  // Set generator specific options
126  //--------------
127  // can only be set after InputInit() is called
128 
129  // Simple Input generator:
130  // if you run more than one of these Input::SIMPLE_NUMBER > 1
131  // add the settings for other with [1], next with [2]...
132  if (Input::SIMPLE)
133  {
134  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 5);
135  if (Input::HEPMC || Input::EMBED)
136  {
137  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
138  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
139  }
140  else
141  {
142  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
145  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
146  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 5.);
147  }
148  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
149  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
150  INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
151  }
152  // Upsilons
153  // if you run more than one of these Input::UPSILON_NUMBER > 1
154  // add the settings for other with [1], next with [2]...
155  if (Input::UPSILON)
156  {
157  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
158  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
159  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
160  // Y species - select only one, last one wins
161  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
162  if (Input::HEPMC || Input::EMBED)
163  {
164  INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
165  INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
166  }
167  }
168  // particle gun
169  // if you run more than one of these Input::GUN_NUMBER > 1
170  // add the settings for other with [1], next with [2]...
171  if (Input::GUN)
172  {
173  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
174  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
175  }
176 
177  // pythia6
178  if (Input::PYTHIA6)
179  {
182  }
183  // pythia8
184  if (Input::PYTHIA8)
185  {
188  }
189 
190  //--------------
191  // Set Input Manager specific options
192  //--------------
193  // can only be set after InputInit() is called
194 
195  if (Input::HEPMC)
196  {
199 
200  // optional overriding beam parameters
201  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0); //optional collision smear in space, time
202  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
203  // //optional choice of vertex distribution function in space, time
204  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
209  //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
210  if (Input::PILEUPRATE > 0)
211  {
212  // Copy vertex settings from foreground hepmc input
214  // and then modify the ones you want to be different
215  // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
216  }
217  }
218  if (Input::PILEUPRATE > 0)
219  {
222  }
223  // register all input generators with Fun4All
224  InputRegister();
225 
226  // set up production relatedstuff
227  // Enable::PRODUCTION = true;
228 
229  //======================
230  // Write the DST
231  //======================
232 
233  //Enable::DSTOUT = true;
234  Enable::DSTOUT_COMPRESS = false;
235  DstOut::OutputDir = outdir;
236  DstOut::OutputFile = outputFile;
237 
238  //Option to convert DST to human command readable TTree for quick poke around the outputs
239  // Enable::DSTREADER = true;
240 
241  // turn the display on (default off)
242  Enable::DISPLAY = false;
243 
244  //======================
245  // What to run
246  //======================
247 
248  // QA, main switch
249  Enable::QA = true;
250 
251  // Global options (enabled for all enables subsystems - if implemented)
252  // Enable::ABSORBER = true;
253  // Enable::OVERLAPCHECK = true;
254  // Enable::VERBOSITY = 1;
255 
256  // Enable::BBC = true;
257  Enable::BBCFAKE = true; // Smeared vtx and t0, use if you don't want real BBC in simulation
258 
259  Enable::PIPE = true;
260  Enable::PIPE_ABSORBER = true;
261 
262  // central tracking
263  Enable::MVTX = true;
264  Enable::MVTX_CELL = Enable::MVTX && true;
265  Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
266  Enable::MVTX_QA = Enable::MVTX_CLUSTER and Enable::QA && true;
267 
268  Enable::INTT = true;
269  Enable::INTT_CELL = Enable::INTT && true;
270  Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
271  Enable::INTT_QA = Enable::INTT_CLUSTER and Enable::QA && true;
272 
273  Enable::TPC = true;
274  Enable::TPC_ABSORBER = true;
275  Enable::TPC_CELL = Enable::TPC && true;
276  Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
277  Enable::TPC_QA = Enable::TPC_CLUSTER and Enable::QA && true;
278 
279  //Enable::MICROMEGAS = true;
282 
283  Enable::TRACKING_TRACK = true;
284  Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && true;
285  Enable::TRACKING_QA = Enable::TRACKING_TRACK and Enable::QA && true;
286 
287  // cemc electronics + thin layer of W-epoxy to get albedo from cemc
288  // into the tracking, cannot run together with CEMC
289  // Enable::CEMCALBEDO = true;
290 
291  Enable::CEMC = true;
292  Enable::CEMC_ABSORBER = true;
297  Enable::CEMC_QA = Enable::CEMC_CLUSTER and Enable::QA && true;
298 
299  Enable::HCALIN = true;
305  Enable::HCALIN_QA = Enable::HCALIN_CLUSTER and Enable::QA && true;
306 
307  Enable::MAGNET = true;
309 
310  Enable::HCALOUT = true;
316  Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER and Enable::QA && true;
317 
318  // forward EMC
319  //Enable::FEMC = true;
320  Enable::FEMC_ABSORBER = true;
323  Enable::FEMC_EVAL = Enable::FEMC_CLUSTER and Enable::QA && true;
324 
325  Enable::EPD = false;
326 
328  //Enable::PLUGDOOR = true;
330 
331  Enable::GLOBAL_RECO = true;
332  //Enable::GLOBAL_FASTSIM = true;
333  //Enable::KFPARTICLE = true;
334  //Enable::KFPARTICLE_VERBOSITY = 1;
335  //Enable::KFPARTICLE_TRUTH_MATCH = true;
336  //Enable::KFPARTICLE_SAVE_NTUPLE = true;
337 
338  Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
339 
340  Enable::JETS = true;
342  Enable::JETS_QA = Enable::JETS and Enable::QA && true;
343 
344  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
345  // single particle / p+p-only simulations, or for p+Au / Au+Au
346  // simulations which don't particularly care about jets)
348 
349  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
351  // particle flow jet reconstruction - needs topoClusters!
353 
354  // new settings using Enable namespace in GlobalVariables.C
355  Enable::BLACKHOLE = true;
356  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
357  //BlackHoleGeometry::visible = true;
358 
359  // run user provided code (from local G4_User.C)
360  //Enable::USER = true;
361 
362  //---------------
363  // World Settings
364  //---------------
365  // G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
366  // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
367 
368  //---------------
369  // Magnet Settings
370  //---------------
371 
372  // 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)
373  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
374  G4MAGNET::magfield_rescale = -1.4 / 1.5; // make consistent with expected Babar field strength of 1.4T
375 
376  //---------------
377  // Pythia Decayer
378  //---------------
379  // list of decay types in
380  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
381  // default is All:
382  // G4P6DECAYER::decayType = EDecayType::kAll;
383 
384  // Initialize the selected subsystems
385  G4Init();
386 
387  //---------------------
388  // GEANT4 Detector description
389  //---------------------
390  if (!Input::READHITS)
391  {
392  G4Setup();
393  }
394 
395  //------------------
396  // Detector Division
397  //------------------
398 
400 
401  if (Enable::MVTX_CELL) Mvtx_Cells();
402  if (Enable::INTT_CELL) Intt_Cells();
403  if (Enable::TPC_CELL) TPC_Cells();
405 
407 
409 
411 
412  //-----------------------------
413  // CEMC towering and clustering
414  //-----------------------------
415 
418 
419  //-----------------------------
420  // HCAL towering and clustering
421  //-----------------------------
422 
425 
426  if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
428 
429  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
431 
434 
435  //--------------
436  // SVTX tracking
437  //--------------
438  if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
439  if (Enable::INTT_CLUSTER) Intt_Clustering();
440  if (Enable::TPC_CLUSTER) TPC_Clustering();
442 
443  if (Enable::TRACKING_TRACK)
444  {
445  TrackingInit();
446  Tracking_Reco();
447  }
448  //-----------------
449  // Global Vertexing
450  //-----------------
451 
453  {
454  cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
455  gSystem->Exit(1);
456  }
458  {
459  Global_Reco();
460  }
461  else if (Enable::GLOBAL_FASTSIM)
462  {
463  Global_FastSim();
464  }
465 
466  //-----------------
467  // Calo Trigger Simulation
468  //-----------------
469 
470  if (Enable::CALOTRIGGER)
471  {
472  CaloTrigger_Sim();
473  }
474 
475  //---------
476  // Jet reco
477  //---------
478 
479  if (Enable::JETS) Jet_Reco();
480  if (Enable::HIJETS) HIJetReco();
481 
483 
484  //----------------------
485  // Simulation evaluation
486  //----------------------
487  string outputroot = outputFile;
488  string remove_this = ".root";
489  size_t pos = outputroot.find(remove_this);
490  if (pos != string::npos)
491  {
492  outputroot.erase(pos, remove_this.length());
493  }
494 
495  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
496 
497  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
498 
499  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
500 
501  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
502 
503  if (Enable::FEMC_EVAL) FEMC_Eval(outputroot + "_g4femc_eval.root");
504 
505  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
506 
507  if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
508 
510 
511  AnaTutorial *anaTutorial = new AnaTutorial("anaTutorial", outputroot + "_anaTutorial.root");
512  anaTutorial->setMinJetPt(25.);
513  anaTutorial->Verbosity(0);
514  anaTutorial->analyzeTracks(false);
515  anaTutorial->analyzeClusters(false);
516  anaTutorial->analyzeJets(true);
517  anaTutorial->analyzeTruth(false);
518  se->registerSubsystem(anaTutorial);
519 
520  //======================
521  // Run KFParticle on evt
522  //======================
525  //if (Enable::KFPARTICLE && Input::LAMBDAC) KFParticle_Lambdac_Reco();
526 
527  //----------------------
528  // Standard QAs
529  //----------------------
530 
531  if (Enable::CEMC_QA) CEMC_QA();
534 
535  if (Enable::JETS_QA) Jet_QA();
536 
537  if (Enable::MVTX_QA) Mvtx_QA();
538  if (Enable::INTT_QA) Intt_QA();
539  if (Enable::TPC_QA) TPC_QA();
540  if (Enable::TRACKING_QA) Tracking_QA();
541 
542  if (Enable::TRACKING_QA and Enable::CEMC_QA and Enable::HCALIN_QA and Enable::HCALOUT_QA) QA_G4CaloTracking();
543 
544  //--------------
545  // Set up Input Managers
546  //--------------
547 
548  InputManagers();
549 
550  if (Enable::PRODUCTION)
551  {
553  }
554 
555  if (Enable::DSTOUT)
556  {
557  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
558  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
560  {
561  ShowerCompress();
562  DstCompress(out);
563  }
564  se->registerOutputManager(out);
565  }
566  //-----------------
567  // Event processing
568  //-----------------
569  if (Enable::DISPLAY)
570  {
571  DisplayOn();
572 
573  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
574  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
575 
576  cout << "-------------------------------------------------" << endl;
577  cout << "You are in event display mode. Run one event with" << endl;
578  cout << "se->run(1)" << endl;
579  cout << "Run Geant4 command with following examples" << endl;
580  gROOT->ProcessLine("displaycmd()");
581 
582  return 0;
583  }
584 
585  // if we use a negative number of events we go back to the command line here
586  if (nEvents < 0)
587  {
588  return 0;
589  }
590  // if we run the particle generator and use 0 it'll run forever
591  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS)
592  {
593  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
594  cout << "it will run forever, so I just return without running anything" << endl;
595  return 0;
596  }
597 
598  se->skip(skip);
599  se->run(nEvents);
600 
601  //-----
602  // QA output
603  //-----
604 
605  if (Enable::QA) QA_Output(outputroot + "_qa.root");
606 
607  //-----
608  // Exit
609  //-----
610 
611  se->End();
612  std::cout << "All done" << std::endl;
613  delete se;
614  if (Enable::PRODUCTION)
615  {
617  }
618 
619  gSystem->Exit(0);
620  return 0;
621 }
622 #endif