Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_sPHENIX.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_sPHENIX.C
1 
3  const int nEvents = 10000,
4  const char * inputFile = "/sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/fieldmap/G4Hits_sPHENIX_e-_eta0_8GeV-0002.root",
5  const char * outputFile = "G4sPHENIXCells.root",
6  const char * embed_input_file = "/sphenix/sim/sim01/production/2016-07-12/sHijing/spacal2d/G4Hits_sPHENIX_sHijing-0-4.4fm.list"
7  )
8 {
9 
10  //===============
11  // Input options
12  //===============
13 
14  // Either:
15  // read previously generated g4-hits files, in this case it opens a DST and skips
16  // the simulations step completely. The G4Setup macro is only loaded to get information
17  // about the number of layers used for the cell reco code
18  //
19  // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
20  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
21  const bool readhits = false;
22  // Or:
23  // read files in HepMC format (typically output from event generators like hijing or pythia)
24  const bool readhepmc = false; // read HepMC files
25  // Or:
26  // Use particle generator
27  const bool runpythia8 = true;
28  const bool runpythia6 = false;
29  // And
30  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
31  // 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
32  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
33  const bool do_embedding = false;
34 
35  //======================
36  // What to run
37  //======================
38 
39  bool do_bbc = false;
40 
41  bool do_pipe = false;
42 
43  bool do_svtx = false;
44  bool do_svtx_cell = false;
45  bool do_svtx_track = false;
46  bool do_svtx_eval = false;
47 
48  bool do_preshower = false;
49 
50  bool do_cemc = false;
51  bool do_cemc_cell = false;
52  bool do_cemc_twr = false;
53  bool do_cemc_cluster = false;
54  bool do_cemc_eval = false;
55 
56  bool do_hcalin = false;
57  bool do_hcalin_cell = false;
58  bool do_hcalin_twr = false;
59  bool do_hcalin_cluster = false;
60  bool do_hcalin_eval = false;
61 
62  bool do_magnet = false;
63 
64  bool do_hcalout = false;
65  bool do_hcalout_cell = false;
66  bool do_hcalout_twr = false;
67  bool do_hcalout_cluster = false;
68  bool do_hcalout_eval = false;
69 
70  bool do_global = false;
71  bool do_global_fastsim = false;
72 
73  bool do_jet_reco = true;
74  bool do_jet_eval = false;
75 
76  bool do_dst_compress = false;
77 
78  //Option to convert DST to human command readable TTree for quick poke around the outputs
79  bool do_DSTReader = true;
80  //---------------
81  // Load libraries
82  //---------------
83 
84  gSystem->Load("libfun4all.so");
85  gSystem->Load("libg4detectors.so");
86  gSystem->Load("libphhepmc.so");
87  gSystem->Load("libg4testbench.so");
88  gSystem->Load("libg4hough.so");
89  gSystem->Load("libg4calo.so");
90  gSystem->Load("libg4eval.so");
91 
92  // establish the geometry and reconstruction setup
93  gROOT->LoadMacro("G4Setup_sPHENIX.C");
94  G4Init(do_svtx,do_preshower,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe);
95 
96  int absorberactive = 1; // set to 1 to make all absorbers active volumes
97  const string magfield = "0"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
98 // 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)
99  const float magfield_rescale = 1.4/1.5; // scale the map to a 1.4 T field
100 
101  //---------------
102  // Fun4All server
103  //---------------
104 
106  se->Verbosity(1);
107  // just if we set some flags somewhere in this macro
109  // By default every random number generator uses
110  // PHRandomSeed() which reads /dev/urandom to get its seed
111  // if the RANDOMSEED flag is set its value is taken as seed
112  // You ca neither set this to a random value using PHRandomSeed()
113  // which will make all seeds identical (not sure what the point of
114  // this would be:
115  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
116  // or set it to a fixed value so you can debug your code
117  // rc->set_IntFlag("RANDOMSEED", 12345);
118 
119  //-----------------
120  // Event generation
121  //-----------------
122 
123  if (readhits)
124  {
125  // Get the hits from a file
126  // The input manager is declared later
127 
128  if (do_embedding)
129  {
130  cout <<"Do not support read hits and embed background at the same time."<<endl;
131  exit(1);
132  }
133 
134  }
135  else if (readhepmc)
136  {
137  // this module is needed to read the HepMC records into our G4 sims
138  // but only if you read HepMC input files
139  HepMCNodeReader *hr = new HepMCNodeReader();
140  se->registerSubsystem(hr);
141  }
142  else if (runpythia8)
143  {
144  gSystem->Load("libPHPythia8.so");
145 
146 
147  PHPy8JetTrigger *theTrigger = new PHPy8JetTrigger();
148 // theTrigger->Verbosity(10);
149  theTrigger->SetEtaHighLow(-1, 1);
150  theTrigger->SetJetR(.4);
151  theTrigger->SetMinJetPt(25);
152 
153  PHPythia8* pythia8 = new PHPythia8();
154  // see coresoftware/generators/PHPythia8 for example config
155  pythia8->set_config_file("phpythia8.cfg");
156 
157  pythia8->beam_vertex_parameters(0,0,0,0,0,5);
158  pythia8->register_trigger(theTrigger);
159 // pythia8->set_trigger_AND();
160 
161  se->registerSubsystem(pythia8);
162  pythia8->print_config();
163 // pythia8->Verbosity(10);
164 
165  HepMCNodeReader *hr = new HepMCNodeReader();
166  se->registerSubsystem(hr);
167  }
168  else if (runpythia6)
169  {
170  gSystem->Load("libPHPythia6.so");
171 
172  PHPythia6 *pythia6 = new PHPythia6();
173  pythia6->set_config_file("phpythia6.cfg");
174  se->registerSubsystem(pythia6);
175 
176  HepMCNodeReader *hr = new HepMCNodeReader();
177  se->registerSubsystem(hr);
178  }
179  else
180  {
181  // toss low multiplicity dummy events
183  gen->add_particles("e-",1); // mu+,e+,proton,pi+,Upsilon
184  // gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
185  if (readhepmc || do_embedding) {
186  gen->set_reuse_existing_vertex(true);
187  gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
188  } else {
192  gen->set_vertex_distribution_mean(0.0,0.0,0.0);
193  gen->set_vertex_distribution_width(0.0,0.0,5.0);
194  }
196  gen->set_vertex_size_parameters(10.0,0.0);
197  gen->set_eta_range(-0.1, 0.1);
198  gen->set_phi_range(-1.0*TMath::Pi(), 1.0*TMath::Pi());
199  gen->set_pt_range(24, 24);
200  gen->Embed(1);
201  gen->Verbosity(0);
202  se->registerSubsystem(gen);
203  }
204 
205  if (!readhits)
206  {
207  //---------------------
208  // Detector description
209  //---------------------
210 
211  G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
212  do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe, magfield_rescale);
213  }
214 
215  //---------
216  // BBC Reco
217  //---------
218 
219  if (do_bbc)
220  {
221  gROOT->LoadMacro("G4_Bbc.C");
222  BbcInit();
223  Bbc_Reco();
224  }
225  //------------------
226  // Detector Division
227  //------------------
228 
229  if (do_svtx_cell) Svtx_Cells();
230 
231  if (do_cemc_cell) CEMC_Cells();
232 
233  if (do_hcalin_cell) HCALInner_Cells();
234 
235  if (do_hcalout_cell) HCALOuter_Cells();
236 
237  //-----------------------------
238  // CEMC towering and clustering
239  //-----------------------------
240 
241  if (do_cemc_twr) CEMC_Towers();
242  if (do_cemc_cluster) CEMC_Clusters();
243 
244  //-----------------------------
245  // HCAL towering and clustering
246  //-----------------------------
247 
248  if (do_hcalin_twr) HCALInner_Towers();
249  if (do_hcalin_cluster) HCALInner_Clusters();
250 
251  if (do_hcalout_twr) HCALOuter_Towers();
252  if (do_hcalout_cluster) HCALOuter_Clusters();
253 
254  if (do_dst_compress) ShowerCompress();
255 
256  //--------------
257  // SVTX tracking
258  //--------------
259 
260  if (do_svtx_track) Svtx_Reco();
261 
262  //-----------------
263  // Global Vertexing
264  //-----------------
265 
266  if (do_global)
267  {
268  gROOT->LoadMacro("G4_Global.C");
269  Global_Reco();
270  }
271 
272  else if (do_global_fastsim)
273  {
274  gROOT->LoadMacro("G4_Global.C");
275  Global_FastSim();
276  }
277 
278  //---------
279  // Jet reco
280  //---------
281 
282  if (do_jet_reco)
283  {
284  gROOT->LoadMacro("G4_Jets.C");
285  Jet_Reco();
286  }
287  //----------------------
288  // Simulation evaluation
289  //----------------------
290 
291  if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
292 
293  if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
294 
295  if (do_hcalin_eval) HCALInner_Eval("g4hcalin_eval.root");
296 
297  if (do_hcalout_eval) HCALOuter_Eval("g4hcalout_eval.root");
298 
299  if (do_jet_eval) Jet_Eval("g4jet_eval.root");
300 
301 
302 
303  //--------------
304  // IO management
305  //--------------
306 
307  if (readhits)
308  {
309  // Hits file
310  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
311  hitsin->fileopen(inputFile);
312  se->registerInputManager(hitsin);
313  }
314  if (do_embedding)
315  {
316  if (embed_input_file == NULL)
317  {
318  cout << "Missing embed_input_file! Exit";
319  exit(3);
320  }
321 
322  Fun4AllDstInputManager *in1 = new Fun4AllNoSyncDstInputManager("DSTinEmbed");
323  // in1->AddFile(embed_input_file); // if one use a single input file
324  in1->AddListFile(embed_input_file); // RecommendedL: if one use a text list of many input files
325  se->registerInputManager(in1);
326  }
327  if (readhepmc)
328  {
330  se->registerInputManager( in );
331  se->fileopen( in->Name().c_str(), inputFile );
332  }
333  else
334  {
335  // for single particle generators we just need something which drives
336  // the event loop, the Dummy Input Mgr does just that
338  se->registerInputManager( in );
339  }
340 
341  //temp lines for QA modules
342  {
343 
344  gSystem->Load("libqa_modules");
345 
346  if (do_jet_reco)
347  {
348  QAG4SimulationJet * calo_jet7 = new QAG4SimulationJet(
349  "AntiKt_Truth_r07",QAG4SimulationJet:: kProcessTruthSpectrum);
350  se->registerSubsystem(calo_jet7);
351 
352  QAG4SimulationJet * calo_jet7 = new QAG4SimulationJet(
353  "AntiKt_Truth_r04", QAG4SimulationJet::kProcessTruthSpectrum);
354  se->registerSubsystem(calo_jet7);
355 
356  QAG4SimulationJet * calo_jet7 = new QAG4SimulationJet(
357  "AntiKt_Truth_r02",QAG4SimulationJet:: kProcessTruthSpectrum);
358  se->registerSubsystem(calo_jet7);
359  }
360  }
361 
362  // HF jet trigger moudle
363  assert (gSystem->Load("libHFJetTruthGeneration") == 0);
364  {
365  if (do_jet_reco)
366  {
368  "HFJetTruthTrigger.root",5 , "AntiKt_Truth_r07");
369 // jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
370  se->registerSubsystem(jt);
371 
373  "HFJetTruthTrigger.root",5 , "AntiKt_Truth_r04");
374 // jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
375  se->registerSubsystem(jt);
376 
378  "HFJetTruthTrigger.root",5 , "AntiKt_Truth_r02");
379 // jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
380  se->registerSubsystem(jt);
381  }
382  }
383 
384  // HF jet analysis moudle
385  if (gSystem->Load("libslt") == 0)
386  {
387 
388 
389  if (do_jet_reco)
390  {
392  "AntiKt_Truth_r07");
393  se->registerSubsystem(slt);
394 
396  "AntiKt_Truth_r04");
397 // slt->Verbosity(1);
398  se->registerSubsystem(slt);
399 
401  "AntiKt_Truth_r02");
402  se->registerSubsystem(slt);
403  }
404  }
405 
406 
407  if (do_DSTReader)
408  {
409  //Convert DST to human command readable TTree for quick poke around the outputs
410  gROOT->LoadMacro("G4_DSTReader.C");
411 
412  G4DSTreader( outputFile, //
413  /*int*/ absorberactive ,
414  /*bool*/ do_svtx ,
415  /*bool*/ do_preshower ,
416  /*bool*/ do_cemc ,
417  /*bool*/ do_hcalin ,
418  /*bool*/ do_magnet ,
419  /*bool*/ do_hcalout ,
420  /*bool*/ do_cemc_twr ,
421  /*bool*/ do_hcalin_twr ,
422  /*bool*/ do_magnet ,
423  /*bool*/ do_hcalout_twr
424  );
425  }
426 
427 
428  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
429  if (do_dst_compress) DstCompress(out);
430  se->registerOutputManager(out);
431 
432  //-----------------
433  // Event processing
434  //-----------------
435  if (nEvents < 0)
436  {
437  return;
438  }
439  // if we run the particle generator and use 0 it'll run forever
440  if (nEvents == 0 && !readhits && !readhepmc)
441  {
442  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
443  cout << "it will run forever, so I just return without running anything" << endl;
444  return;
445  }
446 
447  se->run(nEvents);
448 
449  {
450 
451  gSystem->Load("libqa_modules");
452  QAHistManagerDef::saveQARootFile(string(outputFile) + "_qa.root");
453 
454  if (gSystem->Load("libslt") == 0)
455  {
457  string(outputFile) + "_slt.root");
458  }
459  }
460 
461  //-----
462  // Exit
463  //-----
464  gSystem->Exec("ps -o sid,ppid,pid,user,comm,vsize,rssize,time");
465  se->End();
466  std::cout << "All done" << std::endl;
467  delete se;
468  gSystem->Exit(0);
469 }