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