Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_sPHENIX_photons.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_sPHENIX_photons.C
1 #include <iostream>
2 using namespace std;
4  const int nEvents = 2,
5  const int jobnum = 0,
6  const float photenergy = 10,
7  const int doembedding = 0,
8  const char *inputFile = "/sphenix/data/data02/review_2017-08-02/single_particle/spacal2d/fieldmap/G4Hits_sPHENIX_e-_eta0_8GeV-0002.root",
9  const char *outputFile = "G4sPHENIX.root",
10  const char *embed_input_file = "/sphenix/data/data02/review_2017-08-02/sHijing/fm_0-4.list")
11 {
12  // Set the number of TPC layer
13  const int n_TPC_layers = 40; // use 60 for backward compatibility only
14 
15  //===============
16  // Input options
17  //===============
18 
19  // Either:
20  // read previously generated g4-hits files, in this case it opens a DST and skips
21  // the simulations step completely. The G4Setup macro is only loaded to get information
22  // about the number of layers used for the cell reco code
23  //
24  // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
25  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
26  const bool readhits = false;
27  // Or:
28  // read files in HepMC format (typically output from event generators like hijing or pythia)
29  const bool readhepmc = false && doembedding; // read HepMC files
30  // Or:
31  // Use pythia
32  const bool runpythia8 = false;
33  const bool runpythia6 = false;
34  //
35  // **** And ****
36  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
37  // 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
38  // E.g. /sphenix/data/data02/review_2017-08-02/
39  const bool do_embedding = doembedding && false;
40 
41  // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
42  // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
43  const bool particles = true && !readhits;
44  // or gun/ very simple single particle gun generator
45  const bool usegun = false && !readhits;
46  // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also (note, careful to set Z vertex equal to Hijing events)
47  const bool upsilons = false && !readhits;
48  // Event pile up simulation with collision rate in Hz MB collisions.
49  // Note please follow up the macro to verify the settings for beam parameters
50  const double pileup_collision_rate = 0; // 100e3 for 100kHz nominal AuAu collision rate.
51 
52  //======================
53  // What to run
54  //======================
55 
56  bool do_bbc = true;
57 
58  bool do_pipe = true;
59 
60  bool do_svtx = true;
61  bool do_svtx_cell = do_svtx && true;
62  bool do_svtx_track = do_svtx_cell && false;
63  bool do_svtx_eval = do_svtx_track && false;
64 
65  bool do_pstof = false;
66 
67  bool do_cemc = true;
68  bool do_cemc_cell = do_cemc && true;
69  bool do_cemc_twr = do_cemc_cell && true;
70  bool do_cemc_cluster = do_cemc_twr && true;
71  bool do_cemc_eval = do_cemc_cluster && false;
72 
73  bool do_hcalin = true;
74  bool do_hcalin_cell = do_hcalin && true;
75  bool do_hcalin_twr = do_hcalin_cell && true;
76  bool do_hcalin_cluster = do_hcalin_twr && true;
77  bool do_hcalin_eval = do_hcalin_cluster && false;
78 
79  bool do_magnet = true;
80 
81  bool do_hcalout = true;
82  bool do_hcalout_cell = do_hcalout && true;
83  bool do_hcalout_twr = do_hcalout_cell && true;
84  bool do_hcalout_cluster = do_hcalout_twr && true;
85  bool do_hcalout_eval = do_hcalout_cluster && false;
86 
88  bool do_plugdoor = false;
89 
90  bool do_global = true;
91  bool do_global_fastsim = true;
92 
93  bool do_calotrigger = true && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
94 
95  bool do_jet_reco = false;
96  bool do_jet_eval = do_jet_reco && false;
97 
98  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
99  // single particle / p+p-only simulations, or for p+Au / Au+Au
100  // simulations which don't particularly care about jets)
101  bool do_HIjetreco = false && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
102 
103  bool do_dst_compress = false;
104 
105  //Option to convert DST to human command readable TTree for quick poke around the outputs
106  bool do_DSTReader = false;
107  //---------------
108  // Load libraries
109  //---------------
110 
111  gSystem->Load("libfun4all.so");
112  gSystem->Load("libg4detectors.so");
113  gSystem->Load("libphhepmc.so");
114  gSystem->Load("libg4testbench.so");
115  gSystem->Load("libg4hough.so");
116  gSystem->Load("libg4eval.so");
117 
118  // establish the geometry and reconstruction setup
119  gROOT->LoadMacro("G4Setup_sPHENIX.C");
120  G4Init(do_svtx, do_pstof, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe, do_plugdoor, n_TPC_layers);
121 
122  int absorberactive = 1; // set to 1 to make all absorbers active volumes
123  // const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
124  const string magfield = "/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
125  const float magfield_rescale = -1.4 / 1.5; // scale the map to a 1.4 T field
126 
127  //---------------
128  // Fun4All server
129  //---------------
130 
132  se->Verbosity(0);
133  // just if we set some flags somewhere in this macro
135  // By default every random number generator uses
136  // PHRandomSeed() which reads /dev/urandom to get its seed
137  // if the RANDOMSEED flag is set its value is taken as seed
138  // You ca neither set this to a random value using PHRandomSeed()
139  // which will make all seeds identical (not sure what the point of
140  // this would be:
141  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
142  // or set it to a fixed value so you can debug your code
143  // rc->set_IntFlag("RANDOMSEED", 12345);
144 
145  //-----------------
146  // Event generation
147  //-----------------
148 
149  if (readhits)
150  {
151  // Get the hits from a file
152  // The input manager is declared later
153 
154  if (do_embedding)
155  {
156  cout << "Do not support read hits and embed background at the same time." << endl;
157  exit(1);
158  }
159  }
160  else
161  {
162  // running Geant4 stage. First load event generators.
163 
164  if (readhepmc)
165  {
166  // place holder. Additional action is performed in later stage at the input manager level
167  }
168 
169  if (runpythia8)
170  {
171  gSystem->Load("libPHPythia8.so");
172 
173  PHPythia8 *pythia8 = new PHPythia8();
174  // see coresoftware/generators/PHPythia8 for example config
175  pythia8->set_config_file("phpythia8.cfg");
176  if (readhepmc)
177  pythia8->set_reuse_vertex(0); // reuse vertex of subevent with embedding ID of 0
178  // pythia8->set_vertex_distribution_width(0,0,10,0); // additional vertex smearing if needed, more vertex options available
179  se->registerSubsystem(pythia8);
180  }
181 
182  if (runpythia6)
183  {
184  gSystem->Load("libPHPythia6.so");
185 
186  PHPythia6 *pythia6 = new PHPythia6();
187  pythia6->set_config_file("phpythia6.cfg");
188  if (readhepmc)
189  pythia6->set_reuse_vertex(0); // reuse vertex of subevent with embedding ID of 0
190  // pythia6->set_vertex_distribution_width(0,0,10,0); // additional vertex smearing if needed, more vertex options available
191  se->registerSubsystem(pythia6);
192  }
193 
194  // If "readhepMC" is also set, the particles will be embedded in Hijing events
195  if (particles)
196  {
197  // toss low multiplicity dummy events
199  gen->add_particles(22, 1); // mu+,e+,proton,pi+,Upsilon
200  //gen->add_particles("pi+",100); // 100 pion option
201  if (readhepmc || do_embedding || runpythia8 || runpythia6)
202  {
203  gen->set_reuse_existing_vertex(true);
204  gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
205  }
206  else
207  {
211  gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
212  gen->set_vertex_distribution_width(0.0, 0.0, 5.0);
213  }
215  gen->set_vertex_size_parameters(0.0, 0.0);
216  gen->set_eta_range(0.5, 0.7);
217  gen->set_phi_range(0.2,0.4);
218 
219  gen->set_p_range(photenergy, photenergy);
220  if(doembedding)
221  gen->Embed(1);
222  gen->Verbosity(0);
223 
224  se->registerSubsystem(gen);
225  }
226 
227  if (usegun)
228  {
229  PHG4ParticleGun *gun = new PHG4ParticleGun();
230  // gun->set_name("anti_proton");
231  gun->set_name("geantino");
232  gun->set_vtx(0, 0, 0);
233  gun->set_mom(10, 0, 0.01);
234  // gun->AddParticle("geantino",1.7776,-0.4335,0.);
235  // gun->AddParticle("geantino",1.7709,-0.4598,0.);
236  // gun->AddParticle("geantino",2.5621,0.60964,0.);
237  // gun->AddParticle("geantino",1.8121,0.253,0.);
238  // se->registerSubsystem(gun);
240  pgen->set_name("geantino");
241  pgen->set_z_range(0, 0);
242  pgen->set_eta_range(0.01, 0.01);
243  pgen->set_mom_range(10, 10);
244  pgen->set_phi_range(5.3 / 180. * TMath::Pi(), 5.7 / 180. * TMath::Pi());
245  se->registerSubsystem(pgen);
246  }
247 
248  // 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
249  if (upsilons)
250  {
251  // run upsilons for momentum, dca performance, alone or embedded in Hijing
252 
254  vgen->add_decay_particles("e+", "e-", 0); // i = decay id
255  // event vertex
256  if (readhepmc || do_embedding || particles || runpythia8 || runpythia6)
257  {
258  vgen->set_reuse_existing_vertex(true);
259  }
260  else
261  {
262  vgen->set_vtx_zrange(-10.0, +10.0);
263  }
264 
265  // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
266  vgen->set_rapidity_range(-1.0, +1.0);
267  vgen->set_pt_range(0.0, 10.0);
268 
269  int istate = 1;
270 
271  if (istate == 1)
272  {
273  // Upsilon(1S)
274  vgen->set_mass(9.46);
275  vgen->set_width(54.02e-6);
276  }
277  else if (istate == 2)
278  {
279  // Upsilon(2S)
280  vgen->set_mass(10.0233);
281  vgen->set_width(31.98e-6);
282  }
283  else
284  {
285  // Upsilon(3S)
286  vgen->set_mass(10.3552);
287  vgen->set_width(20.32e-6);
288  }
289 
290  vgen->Verbosity(0);
291  vgen->Embed(3);
292  se->registerSubsystem(vgen);
293 
294  cout << "Upsilon generator for istate = " << istate << " created and registered " << endl;
295  }
296  }
297 
298  if (!readhits)
299  {
300  //---------------------
301  // Detector description
302  //---------------------
303 
304  G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
305  do_svtx, do_pstof, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe,do_plugdoor, magfield_rescale);
306  }
307 
308  //---------
309  // BBC Reco
310  //---------
311 
312  if (do_bbc)
313  {
314  gROOT->LoadMacro("G4_Bbc.C");
315  BbcInit();
316  Bbc_Reco();
317  }
318  //------------------
319  // Detector Division
320  //------------------
321 
322  if (do_svtx_cell) Svtx_Cells();
323 
324  if (do_cemc_cell) CEMC_Cells();
325 
326  if (do_hcalin_cell) HCALInner_Cells();
327 
328  if (do_hcalout_cell) HCALOuter_Cells();
329 
330  //-----------------------------
331  // CEMC towering and clustering
332  //-----------------------------
333 
334  if (do_cemc_twr) CEMC_Towers();
335  if (do_cemc_cluster) CEMC_Clusters();
336 
337  //-----------------------------
338  // HCAL towering and clustering
339  //-----------------------------
340 
341  if (do_hcalin_twr) HCALInner_Towers();
342  if (do_hcalin_cluster) HCALInner_Clusters();
343 
344  if (do_hcalout_twr) HCALOuter_Towers();
345  if (do_hcalout_cluster) HCALOuter_Clusters();
346 
347  if (do_dst_compress) ShowerCompress();
348 
349  //--------------
350  // SVTX tracking
351  //--------------
352 
353  if (do_svtx_track) Svtx_Reco();
354 
355  //-----------------
356  // Global Vertexing
357  //-----------------
358 
359  if (do_global)
360  {
361  gROOT->LoadMacro("G4_Global.C");
362  Global_Reco();
363  }
364 
365  else if (do_global_fastsim)
366  {
367  gROOT->LoadMacro("G4_Global.C");
368  Global_FastSim();
369  }
370 
371  //-----------------
372  // Calo Trigger Simulation
373  //-----------------
374 
375  if (do_calotrigger)
376  {
377  gROOT->LoadMacro("G4_CaloTrigger.C");
378  CaloTrigger_Sim();
379  }
380 
381  //---------
382  // Jet reco
383  //---------
384 
385  if (do_jet_reco)
386  {
387  gROOT->LoadMacro("G4_Jets.C");
388  Jet_Reco();
389  }
390 
391  if (do_HIjetreco)
392  {
393  gROOT->LoadMacro("G4_HIJetReco.C");
394  HIJetReco();
395  }
396 
397  //----------------------
398  // Simulation evaluation
399  //----------------------
400 
401  if (do_svtx_eval) Svtx_Eval(string(outputFile) + "_g4svtx_eval.root");
402 
403  if (do_cemc_eval) CEMC_Eval(string(outputFile) + "_g4cemc_eval.root");
404 
405  if (do_hcalin_eval) HCALInner_Eval(string(outputFile) + "_g4hcalin_eval.root");
406 
407  if (do_hcalout_eval) HCALOuter_Eval(string(outputFile) + "_g4hcalout_eval.root");
408 
409  if (do_jet_eval) Jet_Eval(string(outputFile) + "_g4jet_eval.root");
410 
411  gSystem->Load("libPhotons.so");
412  Photons *photons = new Photons(outputFile);
413  photons->Verbosity(0);
414 
415  photons->set_cluspt_mincut(0.5);
416  photons->set_eta_lowhigh(-1.05,1.05);
417  if(doembedding)
418  photons->use_embedding(1);
419  if(!doembedding)
420  photons->use_embedding(0);
421  se->registerSubsystem(photons);
422 
423 
424  gSystem->ListLibraries();
425 
426 
427 
428  //--------------
429  // IO management
430  //--------------
431 
432  if (readhits)
433  {
434  //meta-lib for DST objects used in simulation outputs
435  gSystem->Load("libg4dst.so");
436 
437  // Hits file
438  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
439  hitsin->fileopen(inputFile);
440  se->registerInputManager(hitsin);
441  }
442 
443  if (do_embedding)
444  {
445  if (embed_input_file == NULL)
446  {
447  cout << "Missing embed_input_file! Exit";
448  exit(3);
449  }
450 
451  //meta-lib for DST objects used in simulation outputs
452  gSystem->Load("libg4dst.so");
453 
454  Fun4AllDstInputManager *in1 = new Fun4AllNoSyncDstInputManager("DSTinEmbed");
455  // in1->AddFile(embed_input_file); // if one use a single input file
456  in1->AddListFile(embed_input_file); // RecommendedL: if one use a text list of many input files
457  se->registerInputManager(in1);
458  }
459 
460  if (readhepmc)
461  {
462  //meta-lib for DST objects used in simulation outputs
463  gSystem->Load("libg4dst.so");
464 
466  se->registerInputManager(in);
467  se->fileopen(in->Name().c_str(), inputFile);
468  //in->set_vertex_distribution_width(100e-4,100e-4,30,0);//optional collision smear in space, time
469  //in->set_vertex_distribution_mean(0,0,1,0);//optional collision central position shift in space, time
470  // //optional choice of vertex distribution function in space, time
471  //in->set_vertex_distribution_function(PHHepMCGenHelper::Gaus,PHHepMCGenHelper::Gaus,PHHepMCGenHelper::Uniform,PHHepMCGenHelper::Gaus);
476  //in->set_embedding_id(2);
477  }
478  else
479  {
480  // for single particle generators we just need something which drives
481  // the event loop, the Dummy Input Mgr does just that
483  se->registerInputManager(in);
484  }
485 
486  if (pileup_collision_rate > 0)
487  {
488  // pile up simulation.
489  // add random beam collisions following a collision diamond and rate from a HepMC stream
490  Fun4AllHepMCPileupInputManager *pileup = new Fun4AllHepMCPileupInputManager("HepMCPileupInput");
491  se->registerInputManager(pileup);
492 
493  const string pileupfile("/sphenix/sim/sim01/sHijing/sHijing_0-12fm.dat");
494  pileup->AddFile(pileupfile); // HepMC events used in pile up collisions. You can add multiple files, and the file list will be reused.
495  //pileup->set_vertex_distribution_width(100e-4,100e-4,30,5);//override collision smear in space time
496  //pileup->set_vertex_distribution_mean(0,0,0,0);//override collision central position shift in space time
497  pileup->set_collision_rate(pileup_collision_rate);
498 
499  double time_window_minus = -35000;
500  double time_window_plus = 35000;
501 
502  if (do_svtx)
503  {
504  // double TPCDriftVelocity = 6.0 / 1000.0; // cm/ns, which is loaded from G4_SVTX*.C macros
505  time_window_minus = -105.5 / TPCDriftVelocity; // ns
506  time_window_plus = 105.5 / TPCDriftVelocity; // ns;
507  }
508  pileup->set_time_window(time_window_minus, time_window_plus); // override timing window in ns
509  cout << "Collision pileup enabled using file " << pileupfile << " with collision rate " << pileup_collision_rate
510  << " and time window " << time_window_minus << " to " << time_window_plus << endl;
511  }
512 
513  if (do_DSTReader)
514  {
515  //Convert DST to human command readable TTree for quick poke around the outputs
516  gROOT->LoadMacro("G4_DSTReader.C");
517 
518  G4DSTreader(outputFile, //
519  /*int*/ absorberactive,
520  /*bool*/ do_svtx,
521  /*bool*/ do_pstof,
522  /*bool*/ do_cemc,
523  /*bool*/ do_hcalin,
524  /*bool*/ do_magnet,
525  /*bool*/ do_hcalout,
526  /*bool*/ do_cemc_twr,
527  /*bool*/ do_hcalin_twr,
528  /*bool*/ do_magnet,
529  /*bool*/ do_hcalout_twr);
530  }
531 
532  // Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
533  // if (do_dst_compress) DstCompress(out);
534  // se->registerOutputManager(out);
535 
536  //-----------------
537  // Event processing
538  //-----------------
539  if (nEvents < 0)
540  {
541  return;
542  }
543  // if we run the particle generator and use 0 it'll run forever
544  if (nEvents == 0 && !readhits && !readhepmc)
545  {
546  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
547  cout << "it will run forever, so I just return without running anything" << endl;
548  return;
549  }
550 
551  se->run(nEvents);
552 
553  //-----
554  // Exit
555  //-----
556 
557  se->End();
558  std::cout << "All done" << std::endl;
559  delete se;
560  gSystem->Exit(0);
561 }