Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
rundst_spiNo.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file rundst_spiNo.C
1 #ifndef MACRO_FUN4ALLG4CALO_C
2 #define MACRO_FUN4ALLG4CALO_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <DisplayOn.C>
7 #include <G4Setup_sPHENIX.C>
8 #include <G4_Mbd.C>
9 #include <G4_CaloTrigger.C>
10 #include <G4_DSTReader.C>
11 #include <G4_Global.C>
12 #include <G4_HIJetReco.C>
13 #include <G4_Input.C>
14 #include <G4_Jets.C>
15 #include <G4_ParticleFlow.C>
16 #include <G4_Production.C>
17 #include <G4_TopoClusterReco.C>
18 #include <G4_Tracking.C>
19 #include <G4_User.C>
20 
23 #include <fun4all/Fun4AllServer.h>
24 
25 //#include <litecaloeval/LiteCaloEval.h>
26 
27 
28 #include <phool/PHRandomSeed.h>
29 #include <phool/recoConsts.h>
30 
31 #include <litecaloeval/LiteCaloEval.h>
32 #include <calib_emc_pi0/CaloCalibEmc_Pi0.h>
33 
34 R__LOAD_LIBRARY(libfun4all.so)
35 R__LOAD_LIBRARY(libcaloCalibDBFile.so)
36 //R__LOAD_LIBRARY(libcalibCaloEmc_pi0.so)
37 R__LOAD_LIBRARY(libLiteCaloEvalTowSlope.so)
38 
39 
40 // For HepMC Hijing
41 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
42 
44  const int nEvents = 5,
45  // const string &inputFile0 = "DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root",
46  //const string &inputFile1 = "DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root",
47  const int mdc2_4_file_num = 20,
48  const string &outputFile = "dstnewoutput5_calo5",
49  // const string &inputFile0 = "/gpfs/mnt/gpfs02/sphenix/sim/sim01/sphnxpro/MDC2/sHijing_HepMC/fm_0_20/pileup/DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000003-00001.root", //"DST_CALO_G4HIT_sHijing_0_12fm-0000000001-00000.root",
50  // const string &inputFile1 = "/gpfs/mnt/gpfs02/sphenix/sim/sim01/sphnxpro/MDC2/sHijing_HepMC/fm_0_20/pileup/DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000003-00001.root", //"DST_VERTEX_sHijing_0_12fm-0000000001-00000.root",
51  // const string &outputFile = "G4sPHENIX_calo1",
52  const string &outdir = ".")
53 {
55  se->Verbosity(0);
56 
57 
58  // string inputFile0 = "DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root";
59  // string inputFile1 = "DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root";
60 
61  // string inputFile0 = "DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
62  //string inputFile0 = "DST_CALO_CLUSTER_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-00000.root";
63  // string inputFile1 = "DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
64 
65  string inputFile0 = "DST_CALO_G4HIT_single_pi0_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
66 
67  string inputFile1 = "DST_VERTEX_single_pi0_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
68 
69 
70  int ynum_int = 100000+ mdc2_4_file_num;
71  TString yn_tstr = "";
72  yn_tstr += ynum_int;
73  yn_tstr.Remove(0,1);
74  inputFile0 += yn_tstr.Data();
75  inputFile1 += yn_tstr.Data();
76 
77  inputFile0 += ".root";
78  inputFile1 += ".root";
79 
80  cout << "running over these files" << endl;
81  cout << inputFile0 << endl;
82  cout << inputFile1 << endl;
83 
84  //const string &outputFile = "newoutput1_calo1",
85 
86 
87  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
89  /*
90  long mtime = gSystem->Now();
91  TString fileOut = outputFile.c_str();
92  fileOut += (mtime / 100000 - 8542922)/3;
93  string outputFile2 = fileOut.Data();
94  */
95  string outputFile2 = outputFile.c_str();
96  outputFile2 = outputFile2 + ".root";
97 
98  // just if we set some flags somewhere in this macro
100  // By default every random number generator uses
101  // PHRandomSeed() which reads /dev/urandom to get its seed
102  // if the RANDOMSEED flag is set its value is taken as seed
103  // You can either set this to a random value using PHRandomSeed()
104  // which will make all seeds identical (not sure what the point of
105  // this would be:
106  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
107  // or set it to a fixed value so you can debug your code
108  // rc->set_IntFlag("RANDOMSEED", 12345);
109 
110  //===============
111  // Input options
112  //===============
113  // verbosity setting (applies to all input managers)
114  Input::VERBOSITY = 1;
115  // First enable the input generators
116  // Either:
117  // read previously generated g4-hits files, in this case it opens a DST and skips
118  // the simulations step completely. The G4Setup macro is only loaded to get information
119  // about the number of layers used for the cell reco code
120  Input::READHITS = true;
121  INPUTREADHITS::filename[0] = inputFile0;
122  INPUTREADHITS::filename[1] = inputFile1;
123 
124  //-----------------
125  // Initialize the selected Input/Event generation
126  //-----------------
127  // This creates the input generator(s)
128  InputInit();
129 
130  // register all input generators with Fun4All
131  InputRegister();
132 
133  // set up production relatedstuff
134  Enable::PRODUCTION = true;
135 
136  //======================
137  // Write the DST
138  //======================
139 
140  Enable::DSTOUT = true;
141  Enable::DSTOUT_COMPRESS = false;
142  DstOut::OutputDir = outdir;
143  //DstOut::OutputFile = outputFile;
144  DstOut::OutputFile = outputFile2;
145 
146 
147  //Option to convert DST to human command readable TTree for quick poke around the outputs
148  // Enable::DSTREADER = true;
149 
150  // turn the display on (default off)
151  Enable::DISPLAY = false;
152 
153  //======================
154  // What to run
155  //======================
156  // Global options (enabled for all enables subsystems - if implemented)
157  // Enable::ABSORBER = true;
158  // Enable::OVERLAPCHECK = true;
159  // Enable::VERBOSITY = 1;
160 
161  Enable::CEMC = true;
164  // Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
165 // Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
166 
167  Enable::HCALIN = true;
170  // Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
171 // Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
172 
173  Enable::HCALOUT = true;
176  // Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
177 // Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
178 
179 
180 // Enable::EPD = false;
181 
182 // Enable::GLOBAL_RECO = true;
183  // Enable::GLOBAL_FASTSIM = true;
184 
185 // Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
186 
187 // Enable::JETS = true;
188 // Enable::JETS_EVAL = Enable::JETS && true;
189 
190  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
191  // single particle / p+p-only simulations, or for p+Au / Au+Au
192  // simulations which don't particularly care about jets)
193 // Enable::HIJETS = false && Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
194 
195  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
196 // Enable::TOPOCLUSTER = true && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
197  // particle flow jet reconstruction - needs topoClusters!
198 // Enable::PARTICLEFLOW = true && Enable::TOPOCLUSTER;
199 
200  //---------------
201  // Magnet Settings
202  //---------------
203 
204  // const string 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)
205  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
206 // G4MAGNET::magfield_rescale = 1; // make consistent with expected Babar field strength of 1.4T
207 
208  //---------------
209  // Pythia Decayer
210  //---------------
211  // list of decay types in
212  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
213  // default is All:
214  // G4P6DECAYER::decayType = EDecayType::kAll;
215 
216  // Initialize the selected subsystems
217  G4Init();
218 
219  //---------------------
220  // GEANT4 Detector description
221  //---------------------
222  if (!Input::READHITS)
223  {
224  G4Setup();
225  }
226 
227  //------------------
228  // Detector Division
229  //------------------
230 
232 
237 
239 
241 
243 
244  //-----------------------------
245  // CEMC towering and clustering
246  //-----------------------------
247 
249  // if (Enable::CEMC_CLUSTER) CEMC_Clusters();
250 
251  //-----------------------------
252  // HCAL towering and clustering
253  //-----------------------------
254 
257 
260 
261  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
263 
265 
266  //--------------
267  // SVTX tracking
268  //--------------
273 
275  {
276  TrackingInit();
277  Tracking_Reco();
278  }
279  //-----------------
280  // Global Vertexing
281  //-----------------
282 
284  {
285  cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
286  gSystem->Exit(1);
287  }
289  {
290  Global_Reco();
291  }
292  else if (Enable::GLOBAL_FASTSIM)
293  {
294  Global_FastSim();
295  }
296 
297  //-----------------
298  // Calo Trigger Simulation
299  //-----------------
300 
302  {
303  CaloTrigger_Sim();
304  }
305 
306  //---------
307  // Jet reco
308  //---------
309 
310  if (Enable::JETS) Jet_Reco();
311  if (Enable::HIJETS) HIJetReco();
312 
314 
315  //----------------------
316  // Simulation evaluation
317  //----------------------
318  string outputroot = outputFile;
319  string remove_this = ".root";
320  size_t pos = outputroot.find(remove_this);
321  if (pos != string::npos)
322  {
323  outputroot.erase(pos, remove_this.length());
324  }
325 
326  //--------------
327  // Set up Input Managers
328  //--------------
329 
330  InputManagers();
331 
332  if (Enable::PRODUCTION)
333  {
335  }
336 
337  if (Enable::DSTOUT)
338  {
339  string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
340  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
341  out->AddNode("Sync");
342  out->AddNode("EventHeader");
343  out->AddNode("TOWER_SIM_HCALIN");
344  out->AddNode("TOWER_RAW_HCALIN");
345  out->AddNode("TOWER_CALIB_HCALIN");
346  out->AddNode("CLUSTER_HCALIN");
347  out->AddNode("TOWER_SIM_HCALOUT");
348  out->AddNode("TOWER_RAW_HCALOUT");
349  out->AddNode("TOWER_CALIB_HCALOUT");
350  out->AddNode("CLUSTER_HCALOUT");
351  out->AddNode("TOWER_SIM_CEMC");
352  out->AddNode("TOWER_RAW_CEMC");
353  out->AddNode("TOWER_CALIB_CEMC");
354  out->AddNode("TOWERINFO_SIM_CEMC");
355  out->AddNode("TOWERINFO_RAW_CEMC");
356  out->AddNode("TOWERINFO_CALIB_CEMC");
357  // out->AddNode("CLUSTER_CEMC");
358  //out->AddNode("CLUSTER_POS_COR_CEMC");
359 // leave the topo cluster here in case we run this during pass3
360  //out->AddNode("TOPOCLUSTER_ALLCALO");
361  //out->AddNode("TOPOCLUSTER_EMCAL");
362  //out->AddNode("TOPOCLUSTER_HCAL");
363  out->AddNode("GlobalVertexMap");
365  se->registerOutputManager(out);
366  }
367  //-----------------
368  // Event processing
369  //-----------------
370  if (Enable::DISPLAY)
371  {
372  DisplayOn();
373 
374  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
375  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
376 
377  cout << "-------------------------------------------------" << endl;
378  cout << "You are in event display mode. Run one event with" << endl;
379  cout << "se->run(1)" << endl;
380  cout << "Run Geant4 command with following examples" << endl;
381  gROOT->ProcessLine("displaycmd()");
382 
383  return 0;
384  }
385 
386  string outputfile = outputFile2 + "_ho4ho_eval.root";
387  string outputfile2 = outputFile2 + "_piemc.root";
388  string outputfile3 = outputFile2 + "_towslopemc.root";
389  string outputfile4 = outputFile2 + "_hcalin.root";
390  string outputfile5 = outputFile2 + "_hcalin_mod.root";
391  string outputfile6 = outputFile2 + "_emcal.root";
392 
393 
394  LiteCaloEval *eval3a = new LiteCaloEval("CEMCEVALUATOR", "CEMC", outputfile3.c_str());
395  // LiteCaloEval *eval3a = new LiteCaloEval("HOCEMCEVALUATOR", "HCALOUT", outputfile3.c_str());
396  // eval->Verbosity(verbosity);
397  eval3a->CaloType(LiteCaloEval::CEMC);
398  //eval3a->CaloType(LiteCaloEval::HCALOUT);
399  //eval3a->set_mode(1);
400  se->registerSubsystem(eval3a);
401 
402  /*
403  // LiteCaloEval *eval = new LiteCaloEval("CEMCEVALUATOR2", "CEMC", outputfile.c_str());
404  LiteCaloEval *eval = new LiteCaloEval("HOCEMCEVALUATOR2", "HCALOUT", outputfile.c_str());
405  // eval->Verbosity(verbosity);
406  // eval->CaloType(LiteCaloEval::CEMC);
407  eval->CaloType(LiteCaloEval::HCALOUT);
408  se->registerSubsystem(eval);
409 
410  LiteCaloEval *eval7 = new LiteCaloEval("CEMCEVALUATOR", "CEMC", outputfile6.c_str());
411  // LiteCaloEval *eval = new LiteCaloEval("HOCEMCEVALUATOR2", "HCALOUT", outputfile.c_str());
412  // eval->Verbosity(verbosity);
413  eval7->CaloType(LiteCaloEval::CEMC);
414  //eval->CaloType(LiteCaloEval::HCALOUT);
415  se->registerSubsystem(eval7);
416 
417 
418  */
419  /*
420  CaloCalibEmc_Pi0 *eval_pi1 = new CaloCalibEmc_Pi0("CEMC_CALIB_PI0", outputfile2);
421  // eval_pi1->set_mode(1);
422  // eval->Verbosity(verbosity);
423  se->registerSubsystem(eval_pi1);
424  */
425 
426 
427  /*
428  CaloCalibEmc_Pi0 *eval_pi2 = new CaloCalibEmc_Pi0("CEMC_CALIB_PI02", outputfile2);
429  // eval_pi2->set_mode(1);
430  // eval->Verbosity(verbosity);
431  se->registerSubsystem(eval_pi2);
432  cout << "successful registration of pi0 " << endl;
433  */
434 
435  /*
436  LiteCaloEval *eval3 = new LiteCaloEval("HcalIn_towslope", "HCALIN", outputfile4.c_str());
437  // eval->Verbosity(verbosity);
438  eval3->CaloType(LiteCaloEval::HCALIN);
439  se->registerSubsystem(eval3);
440 
441  LiteCaloEval *eval4 = new LiteCaloEval("HIN_towslope2", "HCALIN", outputfile5.c_str());
442  // eval->Verbosity(verbosity);
443  eval4->CaloType(LiteCaloEval::HCALIN);
444  eval4->set_mode(1);
445  se->registerSubsystem(eval4);
446 
447  */
448 
449  // if we use a negative number of events we go back to the command line here
450  if (nEvents < 0)
451  {
452  return 0;
453  }
454  // if we run the particle generator and use 0 it'll run forever
455  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS)
456  {
457  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
458  cout << "it will run forever, so I just return without running anything" << endl;
459  return 0;
460  }
461 
462  se->run(nEvents);
463 
464  //-----
465  // Exit
466  //-----
467 
468  se->End();
469  // se->PrintTimer();
470  //se->PrintMemoryTracker();
471  std::cout << "All done" << std::endl;
472  gSystem->Exit(0);
473  delete se;
474  /*
475  if (Enable::PRODUCTION)
476  {
477  Production_MoveOutput();
478  }
479  */
480 
481  return 0;
482 }
483 #endif