Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_W.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_W.C
1 int
2 Fun4All_G4_W(const int nEvents = 100, const double momentum = 5, const char *outfile = "test.root")
3 {
4 
5  gSystem->Load("libfun4all");
6  gSystem->Load("libg4detectors");
7  gSystem->Load("libg4testbench");
8  gSystem->Load("libg4eval");
9  gSystem->Load("libg4histos");
10 
12  // Make the Server
14  Fun4AllServer *se = Fun4AllServer::instance();
15  se->Verbosity(1);
16 
17 // gSystem->Load("libPHPythia8.so");
18 //
19 // PHPythia8* pythia8 = new PHPythia8();
20 // // see coresoftware/generators/PHPythia8 for example config
21 // pythia8->set_config_file("./phpythia8_510.cfg");
22 // se->registerSubsystem(pythia8);
23 //
24 // HepMCNodeReader *hr = new HepMCNodeReader();
25 // se->registerSubsystem(hr);
26 
27 // // particle gun
28 // PHG4ParticleGun *gun = new PHG4ParticleGun("PGUN");
29 // // gun->set_name("anti_proton");
30 // gun->set_name("mu-");
31 // // gun->set_name("proton");
32 // gun->set_vtx(-0,0,0);
33 // gun->set_mom(10,0,0.);
34 // se->registerSubsystem(gun);
35 
36  // Fun4All G4 module
37  PHG4Reco* g4Reco = new PHG4Reco();
38  // no magnetic field
39  g4Reco->set_field(0);
40  g4Reco->set_field_rescale(1);
41  // size of the world - every detector has to fit in here
42  g4Reco->SetWorldSizeX(500);
43  g4Reco->SetWorldSizeY(500);
44  g4Reco->SetWorldSizeZ(500);
45  // shape of our world - it is a box
46  g4Reco->SetWorldShape("G4BOX");
47  // this is what our world is filled with
48  g4Reco->SetWorldMaterial("G4_AIR");
49  // Geant4 Physics list to use
50  g4Reco->SetPhysicsList("QGSP_BERT");
51 
53  // Event generator
55  // toss low multiplicity dummy events
56  PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
57  gen->add_particles("e-", 1); // mu+,e+,proton,pi+,Upsilon
58  // gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
59 
60  gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
61  PHG4SimpleEventGenerator::Uniform, PHG4SimpleEventGenerator::Uniform);
62  gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
63  gen->set_vertex_distribution_width(0.0, 0.0, 5.0);
64  gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
65  gen->set_vertex_size_parameters(0.0, 0.0);
66  gen->set_eta_range(-0.1, 0.1);
67  gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
68  gen->set_p_range(momentum, momentum);
69  gen->Embed(1);
70  gen->Verbosity(0);
71  se->registerSubsystem(gen);
72 
74  // W-shashlyk
76  // radiation length/layer = 0.12/0.3504 + 0.03/1.436 + 0.15/41.31 = 0.37
77  // 20 Radiation length = 54 layers
78  // Total length = 16 cm
79  double scintiwidth = 0.15;
80  double tungstenwidth = 0.12;
81  double copperwidth = 0.015;
82  double no_overlapp = 0.0001;
83  double radius = 95;
84  int Max_cemc_layer = 54;
85  int absorberactive = 1;
86  int overlapcheck = 1;
87 
88  PHG4CylinderSubsystem *cemc;
89 
90  for (int ilayer = 0; ilayer <= Max_cemc_layer; ilayer++)
91  {
92 
93  cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", 3 * ilayer + 0);
94  cemc->set_double_param("radius", radius);
95  cemc->set_string_param("material", "G4_Cu");
96  cemc->set_double_param("thickness", copperwidth);
97  if (absorberactive)
98  cemc->SetActive();
99  cemc->OverlapCheck(overlapcheck);
100  g4Reco->registerSubsystem(cemc);
101  cemc->SuperDetector("ABSORBER_CEMC");
102 
103  radius += copperwidth;
104  radius += no_overlapp;
105 
106  cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", 3 * ilayer + 1);
107  cemc->set_double_param("radius", radius);
108  cemc->set_string_param("material", "G4_W");
109  cemc->set_double_param("thickness", tungstenwidth);
110  if (absorberactive)
111  cemc->SetActive();
112  cemc->OverlapCheck(overlapcheck);
113  g4Reco->registerSubsystem(cemc);
114  cemc->SuperDetector("ABSORBER_CEMC");
115 
116  radius += tungstenwidth;
117  radius += no_overlapp;
118 
119  cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", 3 * ilayer + 2);
120  cemc->set_double_param("radius", radius);
121  cemc->set_string_param("material", "G4_Cu");
122  cemc->set_double_param("thickness", copperwidth);
123  if (absorberactive)
124  cemc->SetActive();
125  cemc->OverlapCheck(overlapcheck);
126  g4Reco->registerSubsystem(cemc);
127  cemc->SuperDetector("ABSORBER_CEMC");
128 
129  radius += copperwidth;
130  radius += no_overlapp;
131 
132  cemc = new PHG4CylinderSubsystem("CEMC", 3 * ilayer + 0);
133  cemc->set_double_param("radius", radius);
134  cemc->set_string_param("material", "G4_POLYSTYRENE");
135  cemc->set_double_param("thickness", scintiwidth);
136  if (absorberactive)
137  cemc->SetActive();
138  cemc->OverlapCheck(overlapcheck);
139  g4Reco->registerSubsystem(cemc);
140  cemc->SuperDetector("CEMC");
141 
142  radius += scintiwidth;
143  radius += no_overlapp;
144  }
145 
146  //----------------------------------------
147  // BLACKHOLE
148 
149  // swallow all particles coming out of the backend of sPHENIX
150  PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
151  blackhole->set_double_param("radius", radius + 10); // add 10 cm
152 
153  blackhole->set_int_param("lengthviarapidity", 0);
154  blackhole->set_double_param("length", g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
155  blackhole->BlackHole();
156  blackhole->set_double_param("thickness", 0.1); // it needs some thickness
157  blackhole->SetActive(); // always see what leaks out
158  blackhole->OverlapCheck(overlapcheck);
159  g4Reco->registerSubsystem(blackhole);
160 
161  //----------------------------------------
162  // FORWARD BLACKHOLEs
163  // +Z
164  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
165  blackhole->SuperDetector("BH_FORWARD_PLUS");
166  blackhole->set_double_param("radius", 0); // add 10 cm
167  blackhole->set_int_param("lengthviarapidity", 0);
168  blackhole->set_double_param("length", 0.1); // make it cover the world in length
169  blackhole->set_double_param("place_z",
170  g4Reco->GetWorldSizeZ() / 2. - 0.1 - no_overlapp);
171  blackhole->BlackHole();
172  blackhole->set_double_param("thickness", radius - no_overlapp); // it needs some thickness
173  blackhole->SetActive(); // always see what leaks out
174  blackhole->OverlapCheck(overlapcheck);
175  g4Reco->registerSubsystem(blackhole);
176 
177  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
178  blackhole->SuperDetector("BH_FORWARD_NEG");
179  blackhole->set_double_param("radius", 0); // add 10 cm
180  blackhole->set_int_param("lengthviarapidity", 0);
181  blackhole->set_double_param("length", 0.1); // make it cover the world in length
182  blackhole->set_double_param("place_z",
183  -g4Reco->GetWorldSizeZ() / 2. + 0.1 + no_overlapp);
184  blackhole->BlackHole();
185  blackhole->set_double_param("thickness", radius - no_overlapp); // it needs some thickness
186  blackhole->SetActive(); // always see what leaks out
187  blackhole->OverlapCheck(overlapcheck);
188  g4Reco->registerSubsystem(blackhole);
189 
190  //----------------------------------------
191  // PHG4TruthSubsystem
192 
193  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
194  g4Reco->registerSubsystem(truth);
195 
196  se->registerSubsystem(g4Reco);
197 
199  // Output
201 
202 // if (outfile)
203 // {
204 // Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",outfile);
205 // se->registerOutputManager(out);
206 //
207 // }
208  if (outfile)
209  {
210  // save a comprehensive evaluation file
211  PHG4DSTReader* ana = new PHG4DSTReader(
212  string(outfile) + string("_DSTReader.root"));
213  ana->set_save_particle(true);
214  ana->set_load_all_particle(false);
215  ana->set_load_active_particle(false);
216  ana->set_save_vertex(true);
217  if (nEvents > 0 && nEvents < 2)
218  {
219  ana->Verbosity(2);
220  }
221  ana->AddNode("CEMC");
222  se->registerSubsystem(ana);
223 
224  gSystem->Load("libqa_modules");
225 
226  QAG4SimulationCalorimeter * qa = new QAG4SimulationCalorimeter("CEMC",
227  QAG4SimulationCalorimeter::kProcessG4Hit);
228  se->registerSubsystem(qa);
229 
230  }
231 
232  // input - we need a dummy to drive the event loop
233  Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
234  se->registerInputManager(in);
235 
236  //-----------------
237  // Event processing
238  //-----------------
239  if (nEvents < 0)
240  {
241  return;
242  }
243  // if we run the particle generator and use 0 it'll run forever
244  if (nEvents == 0)
245  {
246  cout
247  << "using 0 for number of events is a bad idea when using particle generators"
248  << endl;
249  cout << "it will run forever, so I just return without running anything"
250  << endl;
251  return;
252  }
253 
254  se->run(nEvents);
255 
256 // if (!readhits)
257 // {
258 // // Geometry export:
259 // PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
260 // g4->Dump_GDML("sPHENIX.gdml");
261 // }
262  //-----
263  // Exit
264  //-----
265  if (outfile)
266  {
267  gSystem->Load("libqa_modules");
268  QAHistManagerDef::saveQARootFile(string(outfile) + "_qa.root");
269 
270  }
271 
272  se->End();
273  std::cout << "All done" << std::endl;
274  delete se;
275  gSystem->Exit(0);
276 
277  return;
278 
279 }
280 
281 void
282 G4Cmd(const char * cmd)
283 {
284  Fun4AllServer *se = Fun4AllServer::instance();
285  PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
286  g4->ApplyCommand(cmd);
287 }
288 
289 PHG4ParticleGun *
290 get_gun(const char *name = "PGUN")
291 {
292  Fun4AllServer *se = Fun4AllServer::instance();
293  PHG4ParticleGun *pgun = (PHG4ParticleGun *) se->getSubsysReco(name);
294  return pgun;
295 }
296