EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_Input.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_Input.C
1 #ifndef MACRO_G4INPUT_C
2 #define MACRO_G4INPUT_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <phpythia6/PHPythia6.h>
7 
8 #include <phpythia8/PHPythia8.h>
9 
10 #include <g4main/HepMCNodeReader.h>
11 #include <g4main/PHG4IonGun.h>
15 #include <g4main/PHG4ParticleGun.h>
17 #include <g4main/ReadEICFiles.h>
18 
19 #include <fermimotionafterburner/FermimotionAfterburner.h>
20 
25 
26 #include <phsartre/PHSartre.h>
27 #include <phsartre/PHSartreParticleTrigger.h>
28 
33 #include <fun4all/Fun4AllServer.h>
34 
35 #include <set>
36 
37 R__LOAD_LIBRARY(libfun4all.so)
38 R__LOAD_LIBRARY(libg4testbench.so)
39 R__LOAD_LIBRARY(libPHPythia6.so)
40 R__LOAD_LIBRARY(libPHPythia8.so)
41 R__LOAD_LIBRARY(libPHSartre.so)
42 R__LOAD_LIBRARY(libFermimotionAfterburner.so)
43 
44 namespace Input
45 {
46  // Real Event generators
47  bool PYTHIA6 = false;
48  int PYTHIA6_EmbedId = 0;
49 
50  bool PYTHIA8 = false;
51  int PYTHIA8_EmbedId = 0;
52 
53  bool SARTRE = false;
54  int SARTRE_EmbedId = 0;
55 
56  // Single/multiple particle generators
57  bool DZERO = false;
58  int DZERO_NUMBER = 1;
59  int DZERO_VERBOSITY = 0;
60  std::set<int> DZERO_EmbedIds;
61 
62  bool GUN = false;
63  int GUN_NUMBER = 1;
64  int GUN_VERBOSITY = 0;
65  std::set<int> GUN_EmbedIds;
66 
67  bool IONGUN = false;
68  int IONGUN_NUMBER = 1;
70  std::set<int> IONGUN_EmbedIds;
71 
72  bool PGEN = false;
73  int PGEN_NUMBER = 1;
74  int PGEN_VERBOSITY = 0;
75  std::set<int> PGEN_EmbedIds;
76 
77  bool SIMPLE = false;
78  int SIMPLE_NUMBER = 1;
80 
81  int UPSILON_NUMBER = 1;
83  // other UPSILON settings which are also used elsewhere are in GlobalVariables.C
84 
85  double PILEUPRATE = 0.;
86  bool READHITS = false;
87  int VERBOSITY = 0;
88  int EmbedId = 1;
89 
93  {
94  if (HepMCGen == nullptr)
95  {
96  std::cout << "ApplysPHENIXBeamParameter(): Fatal Error - null input pointer HepMCGen" << std::endl;
97  exit(1);
98  }
99  HepMCGen->set_beam_direction_theta_phi(1e-3, 0, M_PI - 1e-3, 0); //2mrad x-ing of sPHENIX per sPH-TRG-2020-001
100 
102  100e-4, // approximation from past RICH data
103  100e-4, // approximation from past RICH data
104  7, // sPH-TRG-2020-001. Fig 3.2
105  20 / 29.9792); // 20cm collision length / speed of light in cm/ns
111  }
112 
117  {
118  if (HepMCGen == nullptr)
119  {
120  std::cout << "ApplyEICBeamParameter(): Fatal Error - null input pointer HepMCGen" << std::endl;
121  exit(1);
122  }
123 
124  //25mrad x-ing as in EIC CDR
125  const double EIC_hadron_crossing_angle = 25e-3;
126 
128  EIC_hadron_crossing_angle, // beamA_theta
129  0, // beamA_phi
130  M_PI, // beamB_theta
131  0 // beamB_phi
132  );
134  119e-6, 119e-6, // proton beam divergence horizontal & vertical, as in EIC CDR Table 1.1
135  211e-6, 152e-6 // electron beam divergence horizontal & vertical, as in EIC CDR Table 1.1
136  );
137 
138  // angular kick within a bunch as result of crab cavity
139  // using an naive assumption of transfer matrix from the cavity to IP,
140  // which is NOT yet validated with accelerator optics simulations!
141  const double z_hadron_cavity = 52e2; // CDR Fig 3.3
142  const double z_e_cavity = 38e2; // CDR Fig 3.2
144  -EIC_hadron_crossing_angle / 2. / z_hadron_cavity, 0,
145  -EIC_hadron_crossing_angle / 2. / z_e_cavity, 0);
146 
147  // calculate beam sigma width at IP as in EIC CDR table 1.1
148  const double sigma_p_h = sqrt(80 * 11.3e-7);
149  const double sigma_p_v = sqrt(7.2 * 1.0e-7);
150  const double sigma_p_l = 6;
151  const double sigma_e_h = sqrt(45 * 20.0e-7);
152  const double sigma_e_v = sqrt(5.6 * 1.3e-7);
153  const double sigma_e_l = 2;
154 
155  // combine two beam gives the collision sigma in z
156  const double collision_sigma_z = sqrt(sigma_p_l * sigma_p_l + sigma_e_l * sigma_e_l) / 2;
157  const double collision_sigma_t = collision_sigma_z / 29.9792; // speed of light in cm/ns
158 
160  sigma_p_h * sigma_e_h / sqrt(sigma_p_h * sigma_p_h + sigma_e_h * sigma_e_h), //x
161  sigma_p_v * sigma_e_v / sqrt(sigma_p_v * sigma_p_v + sigma_e_v * sigma_e_v), //y
162  collision_sigma_z, //z
163  collision_sigma_t); //t
169  }
170 } // namespace Input
171 
172 namespace INPUTHEPMC
173 {
174  string filename;
175  string listfile;
176  bool FLOW = false;
177  int FLOW_VERBOSITY = 0;
178  bool FERMIMOTION = false;
179 
180 } // namespace INPUTHEPMC
181 
182 namespace INPUTREADEIC
183 {
184  string filename;
185 } // namespace INPUTREADEIC
186 
187 namespace INPUTREADHITS
188 {
189  map<unsigned int, std::string> filename;
190  map<unsigned int, std::string> listfile;
191 } // namespace INPUTREADHITS
192 
193 namespace INPUTEMBED
194 {
195  map<unsigned int, std::string> filename;
196  map<unsigned int, std::string> listfile;
197 } // namespace INPUTEMBED
198 
199 namespace PYTHIA6
200 {
201  string config_file = string(getenv("CALIBRATIONROOT")) + "/Generators/phpythia6.cfg";
202 }
203 
204 namespace PYTHIA8
205 {
206  string config_file = string(getenv("CALIBRATIONROOT")) + "/Generators/phpythia8.cfg";
207 }
208 
209 namespace SARTRE
210 {
211  string config_file = string(getenv("CALIBRATIONROOT")) + "/Generators/sartre.cfg";
212 }
213 
214 namespace PILEUP
215 {
216  string pileupfile = "/sphenix/sim/sim01/sphnxpro/MDC1/sHijing_HepMC/data/sHijing_0_20fm-0000000001-00000.dat";
217  double TpcDriftVelocity = 8.0 / 1000.0;
218 } // namespace PILEUP
219 
220 // collection of pointers to particle generators we can grab in the Fun4All macro
221 namespace INPUTGENERATOR
222 {
223  std::vector<PHG4IonGun *> IonGun;
224  std::vector<PHG4ParticleGenerator *> ParticleGenerator;
225  std::vector<PHG4ParticleGeneratorD0 *> DZeroMesonGenerator;
226  std::vector<PHG4ParticleGeneratorVectorMeson *> VectorMesonGenerator;
227  std::vector<PHG4SimpleEventGenerator *> SimpleEventGenerator;
228  std::vector<PHG4ParticleGun *> Gun;
229  PHPythia6 *Pythia6 = nullptr;
230  PHPythia8 *Pythia8 = nullptr;
231  PHSartre *Sartre = nullptr;
234 } // namespace INPUTGENERATOR
235 
236 namespace INPUTMANAGER
237 {
240 } // namespace INPUTMANAGER
241 
242 void InputInit()
243 {
244  // first consistency checks - not all input generators play nice
245  // with each other
247  {
248  cout << "Reading Hits and Embedding into background at the same time is not supported" << endl;
249  gSystem->Exit(1);
250  }
252  {
253  cout << "Reading Hits and running G4 simultanously is not supported" << endl;
254  gSystem->Exit(1);
255  }
257  {
258  cout << "Pythia6 and Pythia8 cannot be run together - might be possible but needs R&D" << endl;
259  gSystem->Exit(1);
260  }
261 
263  {
264  cout << "Flow Afterburner and Pileup cannot be run simultanously" << endl;
265  gSystem->Exit(1);
266  }
267  // done with consistency checks, create generators in no specific order
268 
270  if (Input::PYTHIA6)
271  {
274 
277  Input::EmbedId++;
278  }
279  if (Input::PYTHIA8)
280  {
282  // see coresoftware/generators/PHPythia8 for example config
284 
287  Input::EmbedId++;
288  }
289  if (Input::SARTRE)
290  {
291  gSystem->Load("libPHSartre.so");
294  // particle trigger to enhance forward J/Psi -> ee
297  //INPUTGENERATOR::SartreTrigger->SetEtaHighLow(4.0,1.4);
298  INPUTGENERATOR::SartreTrigger->SetEtaHighLow(1.0, -1.1); // central arm
300 
303  Input::EmbedId++;
304  }
305  // single particle generators
306  if (Input::DZERO)
307  {
308  for (int i = 0; i < Input::DZERO_NUMBER; ++i)
309  {
310  std::string name = "DZERO_" + std::to_string(i);
312  dzero->Embed(Input::EmbedId);
314  Input::EmbedId++;
315  INPUTGENERATOR::DZeroMesonGenerator.push_back(dzero);
316  }
317  }
318  if (Input::GUN)
319  {
320  for (int i = 0; i < Input::GUN_NUMBER; ++i)
321  {
322  std::string name = "GUN_" + std::to_string(i);
323  PHG4ParticleGun *gun = new PHG4ParticleGun(name);
324  gun->Embed(Input::EmbedId);
326  Input::EmbedId++;
327  INPUTGENERATOR::Gun.push_back(gun);
328  }
329  }
330  if (Input::IONGUN)
331  {
332  for (int i = 0; i < Input::IONGUN_NUMBER; ++i)
333  {
334  std::string name = "IONGUN_" + std::to_string(i);
335  PHG4IonGun *iongun = new PHG4IonGun(name);
336  iongun->Embed(Input::EmbedId);
338  Input::EmbedId++;
339  INPUTGENERATOR::IonGun.push_back(iongun);
340  }
341  }
342  if (Input::PGEN)
343  {
344  for (int i = 0; i < Input::PGEN_NUMBER; ++i)
345  {
346  std::string name = "PGEN_" + std::to_string(i);
348  pgen->Embed(Input::EmbedId);
350  Input::EmbedId++;
351  INPUTGENERATOR::ParticleGenerator.push_back(pgen);
352  }
353  }
354  if (Input::SIMPLE)
355  {
356  for (int i = 0; i < Input::SIMPLE_NUMBER; ++i)
357  {
358  std::string name = "EVTGENERATOR_" + std::to_string(i);
360  simple->Embed(Input::EmbedId);
362  Input::EmbedId++;
363  INPUTGENERATOR::SimpleEventGenerator.push_back(simple);
364  }
365  }
366  if (Input::UPSILON)
367  {
368  for (int i = 0; i < Input::UPSILON_NUMBER; ++i)
369  {
370  std::string name = "UPSILON_" + std::to_string(i);
372  upsilon->Embed(Input::EmbedId);
374  Input::EmbedId++;
375  INPUTGENERATOR::VectorMesonGenerator.push_back(upsilon);
376  }
377  }
378 
379  // input managers for which we might need to set options
380  if (Input::HEPMC)
381  {
383  }
384  if (Input::PILEUPRATE > 0)
385  {
387  }
388 }
389 
391 {
393  if (Input::PYTHIA6)
394  {
396  }
397  if (Input::PYTHIA8)
398  {
400  }
401  if (Input::SARTRE)
402  {
405  }
406  if (Input::DZERO)
407  {
408  int verbosity = max(Input::DZERO_VERBOSITY, Input::VERBOSITY);
409  for (size_t icnt = 0; icnt < INPUTGENERATOR::DZeroMesonGenerator.size(); ++icnt)
410  {
411  INPUTGENERATOR::DZeroMesonGenerator[icnt]->Verbosity(verbosity);
413  }
414  }
415  if (Input::GUN)
416  {
417  int verbosity = max(Input::GUN_VERBOSITY, Input::VERBOSITY);
418  for (size_t icnt = 0; icnt < INPUTGENERATOR::Gun.size(); ++icnt)
419  {
420  INPUTGENERATOR::Gun[icnt]->Verbosity(verbosity);
422  }
423  }
424  if (Input::IONGUN)
425  {
426  int verbosity = max(Input::IONGUN_VERBOSITY, Input::VERBOSITY);
427  for (size_t icnt = 0; icnt < INPUTGENERATOR::IonGun.size(); ++icnt)
428  {
429  INPUTGENERATOR::IonGun[icnt]->Verbosity(verbosity);
431  }
432  }
433  if (Input::PGEN)
434  {
435  int verbosity = max(Input::PGEN_VERBOSITY, Input::VERBOSITY);
436  for (size_t icnt = 0; icnt < INPUTGENERATOR::ParticleGenerator.size(); ++icnt)
437  {
438  INPUTGENERATOR::ParticleGenerator[icnt]->Verbosity(verbosity);
440  }
441  }
442  if (Input::SIMPLE)
443  {
444  int verbosity = max(Input::SIMPLE_VERBOSITY, Input::VERBOSITY);
445  for (size_t icnt = 0; icnt < INPUTGENERATOR::SimpleEventGenerator.size(); ++icnt)
446  {
447  INPUTGENERATOR::SimpleEventGenerator[icnt]->Verbosity(verbosity);
449  }
450  }
451  if (Input::UPSILON)
452  {
453  for (size_t icnt = 0; icnt < INPUTGENERATOR::VectorMesonGenerator.size(); ++icnt)
454  {
457  {
458  INPUTGENERATOR::VectorMesonGenerator[icnt]->set_reuse_existing_vertex(true);
459  }
460  INPUTGENERATOR::VectorMesonGenerator[icnt]->Verbosity(verbosity);
462  }
463  }
464  if (Input::READEIC)
465  {
470  }
471  // here are the various utility modules which read particles and
472  // put them onto the G4 particle stack
474  {
475  if (Input::HEPMC)
476  {
477  // these need to be applied before the HepMCNodeReader since they
478  // work on the hepmc records
479  if (INPUTHEPMC::FLOW)
480  {
483  se->registerSubsystem(burn);
484  }
486  {
488  se->registerSubsystem(fermi);
489  }
490  }
491  // copy HepMC records into G4
492  HepMCNodeReader *hr = new HepMCNodeReader();
493  se->registerSubsystem(hr);
494  }
495 }
496 
498 {
500  if (Input::EMBED)
501  {
502  gSystem->Load("libg4dst.so");
503  if (!INPUTEMBED::filename.empty() && !INPUTEMBED::listfile.empty())
504  {
505  cout << "only filenames or filelists are supported, not mixtures" << endl;
506  gSystem->Exit(1);
507  }
508  if (INPUTEMBED::filename.empty() && INPUTEMBED::listfile.empty())
509  {
510  cout << "you need to give an input filenames or filelist" << endl;
511  gSystem->Exit(1);
512  }
513  for (auto iter = INPUTEMBED::filename.begin(); iter != INPUTEMBED::filename.end(); ++iter)
514  {
515  string mgrname = "DSTin" + to_string(iter->first);
516  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager(mgrname);
517  hitsin->fileopen(iter->second);
518  hitsin->Verbosity(Input::VERBOSITY);
519  hitsin->Repeat();
520  se->registerInputManager(hitsin);
521  }
522  for (auto iter = INPUTEMBED::listfile.begin(); iter != INPUTEMBED::listfile.end(); ++iter)
523  {
524  string mgrname = "DSTin" + to_string(iter->first);
525  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager(mgrname);
526  hitsin->AddListFile(iter->second);
527  hitsin->Verbosity(Input::VERBOSITY);
528  hitsin->Repeat();
529  se->registerInputManager(hitsin);
530  }
531  }
532  if (Input::HEPMC)
533  {
536  if (!INPUTHEPMC::filename.empty() && INPUTHEPMC::listfile.empty())
537  {
539  }
540  else if (!INPUTHEPMC::listfile.empty())
541  {
543  }
544  else
545  {
546  cout << "no filename INPUTHEPMC::filename or listfile INPUTHEPMC::listfile given" << endl;
547  gSystem->Exit(1);
548  }
549  }
550  else if (Input::READHITS)
551  {
552  gSystem->Load("libg4dst.so");
553  if (!INPUTREADHITS::filename.empty() && !INPUTREADHITS::listfile.empty())
554  {
555  cout << "only filenames or filelists are supported, not mixtures" << endl;
556  gSystem->Exit(1);
557  }
558  if (INPUTREADHITS::filename.empty() && INPUTREADHITS::listfile.empty())
559  {
560  cout << "you need to give an input filenames or filelist" << endl;
561  gSystem->Exit(1);
562  }
563  for (auto iter = INPUTREADHITS::filename.begin(); iter != INPUTREADHITS::filename.end(); ++iter)
564  {
565  string mgrname = "DSTin" + to_string(iter->first);
566  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager(mgrname);
567  hitsin->fileopen(iter->second);
568  hitsin->Verbosity(Input::VERBOSITY);
569  se->registerInputManager(hitsin);
570  }
571  for (auto iter = INPUTREADHITS::listfile.begin(); iter != INPUTREADHITS::listfile.end(); ++iter)
572  {
573  string mgrname = "DSTin" + to_string(iter->first);
574  Fun4AllInputManager *hitsin = new Fun4AllDstInputManager(mgrname);
575  hitsin->AddListFile(iter->second);
576  hitsin->Verbosity(Input::VERBOSITY);
577  se->registerInputManager(hitsin);
578  }
579  }
580  else
581  {
584  se->registerInputManager(in);
585  }
586  if (Input::PILEUPRATE > 0)
587  {
592  double time_window = 105.5 / PILEUP::TpcDriftVelocity;
593  INPUTMANAGER::HepMCPileupInputManager->set_time_window(-time_window, time_window);
595  }
596 }
597 #endif