Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_Calo.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_Calo.C
1 #ifndef MACRO_FUN4ALLG4CALO_C
2 #define MACRO_FUN4ALLG4CALO_C
3 
4 #include <GlobalVariables.C>
5 #include <../../pi0Efficiency.h>
6 
7 #include <DisplayOn.C>
8 #include <G4Setup_sPHENIX.C>
9 #include <G4_Bbc.C>
10 #include <G4_CaloTrigger.C>
11 #include <G4_DSTReader.C>
12 #include <G4_Global.C>
13 #include <G4_HIJetReco.C>
14 #include <G4_Input.C>
15 #include <G4_Jets.C>
16 #include <G4_ParticleFlow.C>
17 #include <G4_Production.C>
18 #include <G4_TopoClusterReco.C>
19 #include <G4_Tracking.C>
20 #include <G4_User.C>
21 
22 #include <ffamodules/FlagHandler.h>
23 #include <ffamodules/XploadInterface.h>
24 
27 #include <fun4all/Fun4AllServer.h>
28 
29 #include <phool/PHRandomSeed.h>
30 #include <phool/recoConsts.h>
31 
32 R__LOAD_LIBRARY(libffamodules.so)
33 R__LOAD_LIBRARY(libfun4all.so)
34 R__LOAD_LIBRARY(libpi0Efficiency.so)
35 
36 // For HepMC Hijing
37 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
38 
40  const int nEvents = 1,
41  const string &inputFile0 = "g4hits.list",
42  //const string &inputFile1 = "DST_VERTEX_sHijing_0_12fm-0000000001-00000.root",
43  const string &outputFile = "G4sPHENIX_calo.root",
44  const string &outdir = ".")
45 {
47  se->Verbosity(0);
48 
49  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
51 
52  // just if we set some flags somewhere in this macro
54  // By default every random number generator uses
55  // PHRandomSeed() which reads /dev/urandom to get its seed
56  // if the RANDOMSEED flag is set its value is taken as seed
57  // You can either set this to a random value using PHRandomSeed()
58  // which will make all seeds identical (not sure what the point of
59  // this would be:
60  // rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
61  // or set it to a fixed value so you can debug your code
62  // rc->set_IntFlag("RANDOMSEED", 12345);
63  Enable::XPLOAD = true;
64  rc->set_StringFlag("XPLOAD_TAG",XPLOAD::tag);
65  rc->set_StringFlag("XPLOAD_CONFIG",XPLOAD::config);
66  rc->set_uint64Flag("TIMESTAMP",XPLOAD::timestamp);
67 
68  //===============
69  // Input options
70  //===============
71  // verbosity setting (applies to all input managers)
72  Input::VERBOSITY = 1;
73  // First enable the input generators
74  // Either:
75  // read previously generated g4-hits files, in this case it opens a DST and skips
76  // the simulations step completely. The G4Setup macro is only loaded to get information
77  // about the number of layers used for the cell reco code
78  Input::READHITS = true;
79  INPUTREADHITS::listfile[0] = inputFile0;
80  //INPUTREADHITS::filename[0] = inputFile0;
81  //INPUTREADHITS::filename[1] = inputFile1;
82 
83  //-----------------
84  // Initialize the selected Input/Event generation
85  //-----------------
86  // This creates the input generator(s)
87  InputInit();
88 
89  // register all input generators with Fun4All
90  InputRegister();
91 
92 // register the flag handling
93  FlagHandler *flag = new FlagHandler();
94  se->registerSubsystem(flag);
95 
96  // set up production relatedstuff
97  Enable::PRODUCTION = true;
98 
99  //======================
100  // Write the DST
101  //======================
102 
103  Enable::DSTOUT = true;
104  Enable::DSTOUT_COMPRESS = false;
105  DstOut::OutputDir = outdir;
106  DstOut::OutputFile = outputFile;
107 
108  //Option to convert DST to human command readable TTree for quick poke around the outputs
109  // Enable::DSTREADER = true;
110 
111  // turn the display on (default off)
112  Enable::DISPLAY = false;
113 
114  //======================
115  // What to run
116  //======================
117  // Global options (enabled for all enables subsystems - if implemented)
118  // Enable::ABSORBER = true;
119  // Enable::OVERLAPCHECK = true;
120  // Enable::VERBOSITY = 1;
121 
122  Enable::CEMC = true;
127 
128  // Enable::HCALIN = true;
129 // Enable::HCALIN_CELL = Enable::HCALIN && true;
130 // Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
131 // Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
132 // // Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
133 
134 // Enable::HCALOUT = true;
135 // Enable::HCALOUT_CELL = Enable::HCALOUT && true;
136 // Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
137 // Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
138 // // Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
139 
140 
141 // Enable::EPD = false;
142 
143  Enable::GLOBAL_RECO = true;
144  // Enable::GLOBAL_FASTSIM = true;
145 
146 // Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
147 
148 // Enable::JETS = true;
149 // Enable::JETS_EVAL = Enable::JETS && true;
150 
151  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
152  // single particle / p+p-only simulations, or for p+Au / Au+Au
153  // simulations which don't particularly care about jets)
154 // Enable::HIJETS = false && Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
155 
156  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
157 // Enable::TOPOCLUSTER = true && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
158  // particle flow jet reconstruction - needs topoClusters!
159 // Enable::PARTICLEFLOW = true && Enable::TOPOCLUSTER;
160 
161  //---------------
162  // Magnet Settings
163  //---------------
164 
165  // 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)
166  // G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); // default map from the calibration database
167 // G4MAGNET::magfield_rescale = 1; // make consistent with expected Babar field strength of 1.4T
168 
169  //---------------
170  // Pythia Decayer
171  //---------------
172  // list of decay types in
173  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
174  // default is All:
175  // G4P6DECAYER::decayType = EDecayType::kAll;
176 
177  // Initialize the selected subsystems
178  G4Init();
179 
180  //---------------------
181  // GEANT4 Detector description
182  //---------------------
183  if (!Input::READHITS)
184  {
185  G4Setup();
186  }
187 
188  //------------------
189  // Detector Division
190  //------------------
191 
192  if (Enable::BBC || Enable::BBCFAKE) Bbc_Reco();
193 
198 
200 
202 
204 
205  //-----------------------------
206  // CEMC towering and clustering
207  //-----------------------------
208 
211 
212  //-----------------------------
213  // HCAL towering and clustering
214  //-----------------------------
215 
218 
221 
222  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
224 
226 
227  //--------------
228  // SVTX tracking
229  //--------------
234 
236  {
237  TrackingInit();
238  Tracking_Reco();
239  }
240  //-----------------
241  // Global Vertexing
242  //-----------------
243 
245  {
246  cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
247  gSystem->Exit(1);
248  }
250  {
251  Global_Reco();
252  }
253  else if (Enable::GLOBAL_FASTSIM)
254  {
255  Global_FastSim();
256  }
257 
258  //-----------------
259  // Calo Trigger Simulation
260  //-----------------
261 
263  {
264  CaloTrigger_Sim();
265  }
266 
267  //---------
268  // Jet reco
269  //---------
270 
271  if (Enable::JETS) Jet_Reco();
272  if (Enable::HIJETS) HIJetReco();
273 
275 
276  //----------------------
277  // Simulation evaluation
278  //----------------------
279  string outputroot = outputFile;
280  string remove_this = ".root";
281  size_t pos = outputroot.find(remove_this);
282  if (pos != string::npos)
283  {
284  outputroot.erase(pos, remove_this.length());
285  }
286 
287  //--------------
288  // Set up Input Managers
289  //--------------
290 
291  InputManagers();
292 
293  if (Enable::PRODUCTION)
294  {
296  }
297 
298  if (Enable::DSTOUT)
299  {
300  string FullOutFile = DstOut::OutputFile;
301  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
302  out->AddNode("Sync");
303  out->AddNode("EventHeader");
304  out->AddNode("TOWER_SIM_HCALIN");
305  out->AddNode("TOWER_RAW_HCALIN");
306  out->AddNode("TOWER_CALIB_HCALIN");
307  out->AddNode("CLUSTER_HCALIN");
308  out->AddNode("TOWER_SIM_HCALOUT");
309  out->AddNode("TOWER_RAW_HCALOUT");
310  out->AddNode("TOWER_CALIB_HCALOUT");
311  out->AddNode("CLUSTER_HCALOUT");
312  out->AddNode("TOWER_SIM_CEMC");
313  out->AddNode("TOWER_RAW_CEMC");
314  out->AddNode("TOWER_CALIB_CEMC");
315  out->AddNode("CLUSTER_CEMC");
316  out->AddNode("CLUSTER_POS_COR_CEMC");
317 // leave the topo cluster here in case we run this during pass3
318  out->AddNode("TOPOCLUSTER_ALLCALO");
319  out->AddNode("TOPOCLUSTER_EMCAL");
320  out->AddNode("TOPOCLUSTER_HCAL");
321  out->AddNode("GlobalVertexMap");
323  se->registerOutputManager(out);
324  }
325  //-----------------
326  // Event processing
327  //-----------------
328  if (Enable::DISPLAY)
329  {
330  DisplayOn();
331 
332  gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
333  gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
334 
335  cout << "-------------------------------------------------" << endl;
336  cout << "You are in event display mode. Run one event with" << endl;
337  cout << "se->run(1)" << endl;
338  cout << "Run Geant4 command with following examples" << endl;
339  gROOT->ProcessLine("displaycmd()");
340 
341  return 0;
342  }
343 
344  // if we use a negative number of events we go back to the command line here
345  if (nEvents < 0)
346  {
347  return 0;
348  }
349  // if we run the particle generator and use 0 it'll run forever
350  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS)
351  {
352  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
353  cout << "it will run forever, so I just return without running anything" << endl;
354  return 0;
355  }
356  pi0Efficiency *eval = new pi0Efficiency("dummy", outputroot + "_pi0Efficiency.root");
357  se->registerSubsystem(eval);
358 
359  se->run(nEvents);
360 
361  //-----
362  // Exit
363  //-----
364 
365  //XploadInterface::instance()->Print(); // print used DB files
366  se->End();
367  //se->PrintTimer();
368  std::cout << "All done" << std::endl;
369  delete se;
370  // if (Enable::PRODUCTION)
371  // {
372  // Production_MoveOutput();
373  // }
374 
375  gSystem->Exit(0);
376  return 0;
377 }
378 #endif