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