Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_SinglePhotonGen_sPHENIX.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_SinglePhotonGen_sPHENIX.C
1 #ifndef MACRO_FUN4ALLG4SPHENIX_C
2 #define MACRO_FUN4ALLG4SPHENIX_C
3 
4 #include <../src/singlePhotonClusterAna.h>
5 #include <GlobalVariables.C>
6 
7 #include <DisplayOn.C>
8 #include <G4Setup_sPHENIX.C>
9 #include <G4_Bbc.C>
10 #include <G4_CaloTrigger.C>
11 #include <G4_Centrality.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(libsinglePhotonClusterAna.so)
34 // For HepMC Hijing
35 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
36 
38  const int nEvents = 1,
39  const int pid = 22,
40  const string &outputFile = "G4sPHENIX.root",
41  const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.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  // conditions DB flags
66  //===============
67  // tag
68  rc->set_StringFlag("XPLOAD_TAG","sPHENIX_ExampleGT_1");
69  // database config
70  rc->set_StringFlag("XPLOAD_CONFIG","sPHENIX_cdb");
71  // 64 bit timestamp
72  rc->set_uint64Flag("TIMESTAMP",12345678912345);
73 
74 
75  //===============
76  // Input options
77  //===============
78  // verbosity setting (applies to all input managers)
79  Input::VERBOSITY = 0;
80  // First enable the input generators
81  // Either:
82  // read previously generated g4-hits files, in this case it opens a DST and skips
83  // the simulations step completely. The G4Setup macro is only loaded to get information
84  // about the number of layers used for the cell reco code
85  // Input::READHITS = true;
86  INPUTREADHITS::filename[0] = inputFile;
87  // if you use a filelist
88  // INPUTREADHITS::listfile[0] = inputFile;
89  // Or:
90  // Use particle generator
91  // And
92  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
93  // 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
94  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
95  // Input::EMBED = true;
96  INPUTEMBED::filename[0] = embed_input_file;
97  // if you use a filelist
98  //INPUTEMBED::listfile[0] = embed_input_file;
99 
100  Input::SIMPLE = true;
101  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
103 
104  // Input::PYTHIA6 = true;
105 
106  // Input::PYTHIA8 = true;
107 
108  // Input::GUN = true;
109  // Input::GUN_NUMBER = 3; // if you need 3 of them
110  // Input::GUN_VERBOSITY = 1;
111 
112  //D0 generator
113  //Input::DZERO = false;
114  //Input::DZERO_VERBOSITY = 0;
115  //Lambda_c generator //Not ready yet
116  //Input::LAMBDAC = false;
117  //Input::LAMBDAC_VERBOSITY = 0;
118  // Upsilon generator
119  //Input::UPSILON = true;
120  //Input::UPSILON_NUMBER = 3; // if you need 3 of them
121  //Input::UPSILON_VERBOSITY = 0;
122 
123  // Input::HEPMC = true;
124  INPUTHEPMC::filename = inputFile;
125 
126  // Event pile up simulation with collision rate in Hz MB collisions.
127  //Input::PILEUPRATE = 100e3;
128 
129  //-----------------
130  // Initialize the selected Input/Event generation
131  //-----------------
132  // This creates the input generator(s)
133  InputInit();
134 
135  //--------------
136  // Set generator specific options
137  //--------------
138  // can only be set after InputInit() is called
139 
140  // Simple Input generator:
141  // if you run more than one of these Input::SIMPLE_NUMBER > 1
142  // add the settings for other with [1], next with [2]...
143  if (Input::SIMPLE)
144  {
145  INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(pid, 1);
146 
147  if (Input::HEPMC || Input::EMBED)
148  {
149  INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
150  INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
151  }
152  else
153  {
154  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
157  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
158  /* INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 0.); */
159  INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0.01, 0.01, 5.);
160  }
161  INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1.2, 1.2);
162  INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
163  INPUTGENERATOR::SimpleEventGenerator[0]->set_p_range(2.0, 20.0);
164  }
165  // Upsilons
166  // if you run more than one of these Input::UPSILON_NUMBER > 1
167  // add the settings for other with [1], next with [2]...
168  if (Input::UPSILON)
169  {
170  INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
171  INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
172  INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
173  // Y species - select only one, last one wins
174  INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
175  if (Input::HEPMC || Input::EMBED)
176  {
177  INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
178  INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
179  }
180  }
181  // particle gun
182  // if you run more than one of these Input::GUN_NUMBER > 1
183  // add the settings for other with [1], next with [2]...
184  if (Input::GUN)
185  {
186  INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
187  INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
188  }
189 
190  // pythia6
191  if (Input::PYTHIA6)
192  {
195  }
196  // pythia8
197  if (Input::PYTHIA8)
198  {
201  }
202 
203  //--------------
204  // Set Input Manager specific options
205  //--------------
206  // can only be set after InputInit() is called
207 
208  if (Input::HEPMC)
209  {
212 
213  // optional overriding beam parameters
214  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0); //optional collision smear in space, time
215  // INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
216  // //optional choice of vertex distribution function in space, time
217  //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
222  //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
223  if (Input::PILEUPRATE > 0)
224  {
225  // Copy vertex settings from foreground hepmc input
227  // and then modify the ones you want to be different
228  // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
229  }
230  }
231  if (Input::PILEUPRATE > 0)
232  {
235  }
236  // register all input generators with Fun4All
237  InputRegister();
238 
239  // set up production relatedstuff
240  // Enable::PRODUCTION = true;
241 
242  //======================
243  // Write the DST
244  //======================
245 
246  //Enable::DSTOUT = true;
247  Enable::DSTOUT_COMPRESS = false;
248  DstOut::OutputDir = outdir;
249  DstOut::OutputFile = outputFile;
250 
251  //Option to convert DST to human command readable TTree for quick poke around the outputs
252  // Enable::DSTREADER = true;
253 
254  // turn the display on (default off)
255  //Enable::DISPLAY = true;
256 
257  //======================
258  // What to run
259  //======================
260 
261  // QA, main switch
262  Enable::QA = false;
263 
264  // Global options (enabled for all enables subsystems - if implemented)
265  // Enable::ABSORBER = true;
266  // Enable::OVERLAPCHECK = true;
267  // Enable::VERBOSITY = 1;
268 
269  /* Enable::BBC = true; */
270  // Enable::BBC_SUPPORT = true; // save hist in bbc support structure
271  Enable::BBCFAKE = true; // Smeared vtx and t0, use if you don't want real BBC in simulation
272 
273  Enable::PIPE = true;
274  //Enable::PIPE_ABSORBER = true;
275 
276  // central tracking
277  Enable::MVTX = true;
278  //Enable::MVTX_CELL = Enable::MVTX && true;
281  /* Enable::TrackingService = false; */
282 
283  Enable::INTT = true;
284 // Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
285 // Enable::INTT_SUPPORT = true; // enable global support structure readout
286  // Enable::INTT_CELL = Enable::INTT && true;
288  Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;
289 
290  Enable::TPC = true;
291  // Enable::TPC_ABSORBER = true;
292  // Enable::TPC_CELL = Enable::TPC && true;
294  Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;
295 
296  Enable::MICROMEGAS = true;
297  //Enable::MICROMEGAS_CELL = fasle;// = Enable::MICROMEGAS && true;
299  Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;
300 
301  //Enable::TRACKING_TRACK = true;
302  Enable::TRACKING_EVAL = false;//= Enable::TRACKING_TRACK && true;
303  Enable::TRACKING_QA = false;//xs= Enable::TRACKING_TRACK && Enable::QA && true;
304 
305  // cemc electronics + thin layer of W-epoxy to get albedo from cemc
306  // into the tracking, cannot run together with CEMC
307  // Enable::CEMCALBEDO = true;
308 
309  Enable::CEMC = true;
310  Enable::CEMC_ABSORBER = true;
314  Enable::CEMC_EVAL = false;//Enable::CEMC_CLUSTER && true;
315  Enable::CEMC_QA = false;//Enable::CEMC_CLUSTER && Enable::QA && true;
316 
317  Enable::HCALIN =false;
323  Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && true;
324 
325  Enable::MAGNET = false;
326  Enable::MAGNET_ABSORBER = false;
327 
328  Enable::HCALOUT = false;
334  Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && true;
335 
336  Enable::EPD = false;
337 
338  //Enable::BEAMLINE = true;
339 // Enable::BEAMLINE_ABSORBER = true; // makes the beam line magnets sensitive volumes
340 // Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
342 
343  Enable::ZDC = false;
344 // Enable::ZDC_ABSORBER = true;
345 // Enable::ZDC_SUPPORT = true;
346  Enable::ZDC_TOWER = Enable::ZDC && true;
348 
350  //Enable::PLUGDOOR = true;
351  //Enable::PLUGDOOR_ABSORBER = true;
353 
354  Enable::GLOBAL_RECO = true;
355  //Enable::GLOBAL_FASTSIM = true;
356 
357  //Enable::KFPARTICLE = true;
358  //Enable::KFPARTICLE_VERBOSITY = 1;
359  //Enable::KFPARTICLE_TRUTH_MATCH = true;
360  //Enable::KFPARTICLE_SAVE_NTUPLE = true;
361 
363 
364  Enable::JETS = false;
366  Enable::JETS_QA = Enable::JETS && Enable::QA && true;
367 
368  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
369  // single particle / p+p-only simulations, or for p+Au / Au+Au
370  // simulations which don't particularly care about jets)
372 
373  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
375  // particle flow jet reconstruction - needs topoClusters!
377  // centrality reconstruction
378  Enable::CENTRALITY = true;
379 
380  // new settings using Enable namespace in GlobalVariables.C
381  Enable::BLACKHOLE = true;
382  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
383  //BlackHoleGeometry::visible = true;
384 
385  // run user provided code (from local G4_User.C)
386  //Enable::USER = true;
387 
388  //---------------
389  // World Settings
390  //---------------
391  // G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
392  // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
393 
394  //---------------
395  // Magnet Settings
396  //---------------
397 
398  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root"); // default map from the calibration database
399  // 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)
400 // G4MAGNET::magfield_rescale = 1.; // make consistent with expected Babar field strength of 1.4T
401 
402  //---------------
403  // Pythia Decayer
404  //---------------
405  // list of decay types in
406  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
407  // default is All:
408  // G4P6DECAYER::decayType = EDecayType::kAll;
409 
410  // Initialize the selected subsystems
411  G4Init();
412 
413  //---------------------
414  // GEANT4 Detector description
415  //---------------------
416  if (!Input::READHITS)
417  {
418  G4Setup();
419  }
420 
421  //------------------
422  // Detector Division
423  //------------------
424 
425  if (Enable::BBC || Enable::BBCFAKE) Bbc_Reco();
426 
431 
433 
435 
437 
438  //-----------------------------
439  // CEMC towering and clustering
440  //-----------------------------
441 
444 
445  //-----------------------------
446  // HCAL towering and clustering
447  //-----------------------------
448 
451 
452  if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
454 
455  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
457 
458  //--------------
459  // SVTX tracking
460  //--------------
462  {
463  TrackingInit();
464  }
469 
471  {
472  Tracking_Reco();
473  }
474  //-----------------
475  // Global Vertexing
476  //-----------------
477 
479  {
480  cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
481  gSystem->Exit(1);
482  }
484  {
485  Global_Reco();
486  }
487  else if (Enable::GLOBAL_FASTSIM)
488  {
489  Global_FastSim();
490  }
491 
492  //-----------------
493  // Centrality Determination
494  //-----------------
495 
496  if (Enable::CENTRALITY)
497  {
498  Centrality();
499  }
500 
501  //-----------------
502  // Calo Trigger Simulation
503  //-----------------
504 
506  {
507  CaloTrigger_Sim();
508  }
509 
510  //---------
511  // Jet reco
512  //---------
513 
514  if (Enable::JETS) Jet_Reco();
515  if (Enable::HIJETS) HIJetReco();
516 
518 
519  //----------------------
520  // Simulation evaluation
521  //----------------------
522  string outputroot = outputFile;
523  string remove_this = ".root";
524  size_t pos = outputroot.find(remove_this);
525  if (pos != string::npos)
526  {
527  outputroot.erase(pos, remove_this.length());
528  }
529 
530  if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
531 
532  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
533 
534  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
535 
536  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
537 
538  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
539 
540  if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
541 
543 
544 
545 
546 
547 
548  singlePhotonClusterAna *eval = new singlePhotonClusterAna("dummy", outputroot + "_singlePhotonClusterAna.root");
549  se->registerSubsystem(eval);
550 
551 
552  //======================
553  // Run KFParticle on evt
554  //======================
557 
558  //----------------------
559  // Standard QAs
560  //----------------------
561 
562  if (Enable::CEMC_QA) CEMC_QA();
565 
566  if (Enable::JETS_QA) Jet_QA();
567 
568  if (Enable::MVTX_QA) Mvtx_QA();
569  if (Enable::INTT_QA) Intt_QA();
570  if (Enable::TPC_QA) TPC_QA();
573 
575 
576  //--------------
577  // Set up Input Managers
578  //--------------
579 
580  InputManagers();
581 
582  if (Enable::PRODUCTION)
583  {
585  }
586 
587  if (Enable::DSTOUT)
588  {
589  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
590  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
592  {
593  ShowerCompress();
594  DstCompress(out);
595  }
596  se->registerOutputManager(out);
597  }
598  //-----------------
599  // Event processing
600  //-----------------
601  if (Enable::DISPLAY)
602  {
603  DisplayOn();
604 
605  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
606  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
607 
608  cout << "-------------------------------------------------" << endl;
609  cout << "You are in event display mode. Run one event with" << endl;
610  cout << "se->run(1)" << endl;
611  cout << "Run Geant4 command with following examples" << endl;
612  gROOT->ProcessLine("displaycmd()");
613 
614  return 0;
615  }
616 
617  // if we use a negative number of events we go back to the command line here
618  if (nEvents < 0)
619  {
620  return 0;
621  }
622  // if we run the particle generator and use 0 it'll run forever
623  // for embedding it runs forever if the repeat flag is set
624  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
625  {
626  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
627  cout << "it will run forever, so I just return without running anything" << endl;
628  return 0;
629  }
630 
631  se->skip(skip);
632  se->run(nEvents);
633 
634  //-----
635  // QA output
636  //-----
637 
638  if (Enable::QA) QA_Output(outputroot + "_qa.root");
639 
640  //-----
641  // Exit
642  //-----
643 
644  se->End();
645  std::cout << "All done" << std::endl;
646  delete se;
647  if (Enable::PRODUCTION)
648  {
650  }
651 
652  gSystem->Exit(0);
653  return 0;
654 }
655 #endif