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