Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Jin_BJet.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Jin_BJet.C
1 int Jin_BJet(
2  const int nEvents = 1000,
3  char* inputFile,
4  char* outputFile,
5  int which_tracking = 0
6  )
7 {
8  //===============
9  // Input options
10  //===============
11 
12  //char *inputFile;
13  //char *outputFile;
14  char *embed_input_file;
15 
16  bool do_svtx = true;
17  bool do_svtx_cell = true;
18  bool do_svtx_track = true;
19  bool do_svtx_eval = true;
20 
21  bool do_bbc = false;
22  bool do_global = false;
23 
24  bool do_jet_reco = true;
25 
26  embed_input_file = NULL;
27 
28  int which_reco = 1;
29 
30  // Either:
31  // read previously generated g4-hits files, in this case it opens a DST and skips
32  // the simulations step completely. The G4Setup macro is only loaded to get information
33  // about the number of layers used for the cell reco code
34  const bool readhits = false;
35  // Or:
36  // read files in HepMC format (typically output from event generators like hijing or pythia)
37  const bool readhepmc = true; // read HepMC files
38  // Or:
39  // Use particle generator
40  const bool runpythia = false;
41 
42  //======================
43  // What to run
44  //======================
45 
46 
47  bool do_pipe = true;
48 
49  bool do_preshower = false;
50 
51  bool do_cemc = false;
52  bool do_cemc_cell = false;
53  bool do_cemc_twr = false;
54  bool do_cemc_cluster = false;
55  bool do_cemc_eval = false;
56 
57  bool do_hcalin = false;
58  bool do_hcalin_cell = false;
59  bool do_hcalin_twr = false;
60  bool do_hcalin_cluster = false;
61  bool do_hcalin_eval = false;
62 
63  bool do_magnet = true;
64 
65  bool do_hcalout = false;
66  bool do_hcalout_cell = false;
67  bool do_hcalout_twr = false;
68  bool do_hcalout_cluster = false;
69  bool do_hcalout_eval = false;
70 
71  bool do_global_fastsim = false;
72 
73  bool do_jet_eval = false;
74 
75  bool do_dst_compress = false;
76 
77  //Option to convert DST to human command readable TTree for quick poke around the outputs
78  bool do_DSTReader = false;
79  //---------------
80  // Load libraries
81  //---------------
82 
83  gSystem->Load("libfun4all.so");
84  gSystem->Load("libg4detectors.so");
85  gSystem->Load("libphhepmc.so");
86  gSystem->Load("libg4testbench.so");
87  gSystem->Load("libg4hough.so");
88  gSystem->Load("libg4calo.so");
89  gSystem->Load("libg4eval.so");
90  gSystem->Load("libg4vertex.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,which_tracking);
95 
96  int absorberactive = 1; // set to 1 to make all absorbers active volumes
97  // const string magfield = "1.5"; // 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 
105  Fun4AllServer *se = Fun4AllServer::instance();
106  // se->Verbosity(10);
107  // just if we set some flags somewhere in this macro
108  recoConsts *rc = recoConsts::instance();
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 (embed_input_file)
129  {
130  cout <<"Do not support read hits and embed background"<<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 
141 
142  if (which_reco == 1)
143  hr->Embed( 17 );
144 
145  se->registerSubsystem(hr);
146  }
147  else if (runpythia)
148  {
149  gSystem->Load("libPHPythia8.so");
150 
151  PHPythia8* pythia8 = new PHPythia8();
152  // see coresoftware/generators/PHPythia8 for example config
153  pythia8->set_config_file("phpythia8.cfg");
154  se->registerSubsystem(pythia8);
155 
156  HepMCNodeReader *hr = new HepMCNodeReader();
157  se->registerSubsystem(hr);
158  }
159  else
160  {
161  // toss low multiplicity dummy events
162  PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
163  gen->add_particles(inputFile,1); // mu+,e+,proton,pi+,Upsilon
164  // gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
165  if (readhepmc || embed_input_file) {
166  gen->set_reuse_existing_vertex(true);
167  gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
168  } else {
169  gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
170  PHG4SimpleEventGenerator::Uniform,
171  PHG4SimpleEventGenerator::Uniform);
172  gen->set_vertex_distribution_mean(0.0,0.0,0.0);
173  gen->set_vertex_distribution_width(0.0,0.0,0.0);
174  }
175  gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
176  gen->set_vertex_size_parameters(0.0,0.0);
177  gen->set_eta_range(-0, .1);
178  // gen->set_eta_range(.9, 1);
179  gen->set_phi_range(-1.0*TMath::Pi(), 1.0*TMath::Pi());
180  gen->set_p_range(6, 6);
181  // gen->set_p_range(6, 6);
182  gen->Embed(1);
183  gen->Verbosity(0);
184  se->registerSubsystem(gen);
185  }
186 
187  if (!readhits)
188  {
189  //---------------------
190  // Detector description
191  //---------------------
192 
193  G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
194  do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe, magfield_rescale);
195  }
196 
197  //---------
198  // BBC Reco
199  //---------
200 
201  if (do_bbc)
202  {
203  gROOT->LoadMacro("G4_Bbc.C");
204  BbcInit();
205  Bbc_Reco();
206  }
207  //------------------
208  // Detector Division
209  //------------------
210 
211  if (do_svtx_cell) Svtx_Cells();
212 
213  if (do_cemc_cell) CEMC_Cells();
214 
215  if (do_hcalin_cell) HCALInner_Cells();
216 
217  if (do_hcalout_cell) HCALOuter_Cells();
218 
219  //-----------------------------
220  // CEMC towering and clustering
221  //-----------------------------
222 
223  if (do_cemc_twr) CEMC_Towers();
224  if (do_cemc_cluster) CEMC_Clusters();
225 
226  //-----------------------------
227  // HCAL towering and clustering
228  //-----------------------------
229 
230  if (do_hcalin_twr) HCALInner_Towers();
231  if (do_hcalin_cluster) HCALInner_Clusters();
232 
233  if (do_hcalout_twr) HCALOuter_Towers();
234  if (do_hcalout_cluster) HCALOuter_Clusters();
235 
236  if (do_dst_compress) ShowerCompress();
237 
238  //--------------
239  // SVTX tracking
240  //--------------
241 
242  if (do_svtx_track) Svtx_Reco();
243 
244  //-----------------
245  // Global Vertexing
246  //-----------------
247 
248  if (do_global)
249  {
250  gROOT->LoadMacro("G4_Global.C");
251  Global_Reco();
252  }
253 
254  else if (do_global_fastsim)
255  {
256  gROOT->LoadMacro("G4_Global.C");
257  Global_FastSim();
258  }
259 
260  //---------
261  // Jet reco
262  //---------
263 
264  if (do_jet_reco)
265  {
266 
267  gROOT->LoadMacro("G4_Jets.C");
268 
269  //if (embed_input_file) gROOT->LoadMacro("G4_Jets.C(0,true)");
270  //else gROOT->LoadMacro("G4_Jets.C(0,false)");
271  Jet_Reco();
272  }
273  //----------------------
274  // Simulation evaluation
275  //----------------------
276 
277  if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
278 
279  if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
280 
281  if (do_hcalin_eval) HCALInner_Eval("g4hcalin_eval.root");
282 
283  if (do_hcalout_eval) HCALOuter_Eval("g4hcalout_eval.root");
284 
285  if (do_jet_eval) Jet_Eval("g4jet_eval.root");
286 
287  //--------------
288  // IO management
289  //--------------
290 
291  if (readhits)
292  {
293  // Hits file
294  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
295  hitsin->fileopen(inputFile);
296  se->registerInputManager(hitsin);
297  }
298  if (embed_input_file)
299  {
300 
301  Fun4AllDstInputManager *in1 = new Fun4AllNoSyncDstInputManager("DSTin1");
302  in1->AddFile(embed_input_file);
303  //in1->AddListFile(embed_input_file);
304  se->registerInputManager(in1);
305 
306  }
307  if (readhepmc)
308  {
309  Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
310  se->registerInputManager( in );
311  se->fileopen( in->Name().c_str(), inputFile );
312  }
313  else
314  {
315  // for single particle generators we just need something which drives
316  // the event loop, the Dummy Input Mgr does just that
317  Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
318  se->registerInputManager( in );
319  }
320 
321  // HF jet trigger moudle
322  assert (gSystem->Load("libHFJetTruthGeneration") == 0);
323  {
324  if (do_jet_reco)
325  {
326  HFJetTruthTrigger * jt = new HFJetTruthTrigger(
327  "HFJetTruthTrigger.root",5 , "AntiKt_Truth_r07");
328  //jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
329  se->registerSubsystem(jt);
330 
331  HFJetTruthTrigger * jt = new HFJetTruthTrigger(
332  "HFJetTruthTrigger.root",5 , "AntiKt_Truth_r04");
333  //jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
334  jt->set_pt_min(20);
335  se->registerSubsystem(jt);
336 
337  HFJetTruthTrigger * jt = new HFJetTruthTrigger(
338  "HFJetTruthTrigger.root",5 , "AntiKt_Truth_r02");
339  //jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
340  se->registerSubsystem(jt);
341  }
342  }
343 
344  if (do_DSTReader)
345  {
346  //Convert DST to human command readable TTree for quick poke around the outputs
347  gROOT->LoadMacro("G4_DSTReader.C");
348 
349  G4DSTreader( outputFile, //
350  /*int*/ absorberactive ,
351  /*bool*/ do_svtx ,
352  /*bool*/ do_preshower ,
353  /*bool*/ do_cemc ,
354  /*bool*/ do_hcalin ,
355  /*bool*/ do_magnet ,
356  /*bool*/ do_hcalout ,
357  /*bool*/ do_cemc_twr ,
358  /*bool*/ do_hcalin_twr ,
359  /*bool*/ do_magnet ,
360  /*bool*/ do_hcalout_twr
361  );
362  }
363 
364 
365  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
366  if (do_dst_compress) DstCompress(out);
367  se->registerOutputManager(out);
368 
369  //-----------------
370  // Event processing
371  //-----------------
372  if (nEvents < 0)
373  {
374  return;
375  }
376  // if we run the particle generator and use 0 it'll run forever
377  if (nEvents == 0 && !readhits && !readhepmc)
378  {
379  cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
380  cout << "it will run forever, so I just return without running anything" << endl;
381  return;
382  }
383  se->run(nEvents);
384 
385  //-----
386  // Exit
387  //-----
388 
389  se->End();
390  std::cout << "All done" << std::endl;
391  delete se;
392  gSystem->Exit(0);
393 }