Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Setup_EICDetector.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Setup_EICDetector.C
1 
2 double no_overlapp = 0.0001; // added to radii to avoid overlapping volumes
3 bool overlapcheck = false; // set to true if you want to check for overlaps
4 
5 void G4Init(bool do_svtx = true,
6  bool do_cemc = true,
7  bool do_hcalin = true,
8  bool do_magnet = true,
9  bool do_hcalout = true,
10  bool do_pipe = true,
11  bool do_FGEM = true,
12  bool do_EGEM = true,
13  bool do_FEMC = true,
14  bool do_FHCAL = true,
15  bool do_EEMC = true,
16  bool do_DIRC = true,
17  bool do_RICH = true,
18  bool do_Aerogel = true,
19  int n_TPC_layers = 40)
20 {
21 
22  // load detector/material macros and execute Init() function
23 
24  if (do_pipe)
25  {
26  gROOT->LoadMacro("G4_Pipe.C");
27  PipeInit();
28  }
29  if (do_svtx)
30  {
31  gROOT->LoadMacro("G4_Svtx_maps_ladders+intt_ladders+tpc_KalmanPatRec.C");
32  SvtxInit(n_TPC_layers);
33  }
34 
35  if (do_cemc)
36  {
37  gROOT->LoadMacro("G4_CEmc_Spacal.C");
38  CEmcInit(72); // make it 2*2*2*3*3 so we can try other combinations
39  }
40 
41  if (do_hcalin)
42  {
43  gROOT->LoadMacro("G4_HcalIn_ref.C");
44  HCalInnerInit();
45  }
46 
47  if (do_magnet)
48  {
49  gROOT->LoadMacro("G4_Magnet.C");
50  MagnetInit();
51  }
52  if (do_hcalout)
53  {
54  gROOT->LoadMacro("G4_HcalOut_ref.C");
55  HCalOuterInit();
56  }
57 
58  if (do_FGEM)
59  {
60  gROOT->LoadMacro("G4_FGEM_EIC.C");
61  FGEM_Init();
62  }
63 
64  if (do_EGEM)
65  {
66  gROOT->LoadMacro("G4_EGEM_EIC.C");
67  EGEM_Init();
68  }
69 
70  if (do_FEMC)
71  {
72  gROOT->LoadMacro("G4_FEMC.C");
73  FEMCInit();
74  }
75 
76  if (do_FHCAL)
77  {
78  gROOT->LoadMacro("G4_FHCAL.C");
79  FHCALInit();
80  }
81 
82  if (do_EEMC)
83  {
84  gROOT->LoadMacro("G4_EEMC.C");
85  EEMCInit();
86  }
87 
88  if (do_DIRC)
89  {
90  gROOT->LoadMacro("G4_DIRC.C");
91  DIRCInit();
92  }
93 
94  if (do_RICH)
95  {
96  gROOT->LoadMacro("G4_RICH.C");
97  RICHInit();
98  }
99 
100  if (do_Aerogel)
101  {
102  gROOT->LoadMacro("G4_Aerogel.C");
103  AerogelInit();
104  }
105 
106 
107 }
108 
109 
110 int G4Setup(const int absorberactive = 0,
111  const string &field ="1.5",
112  const EDecayType decayType = TPythia6Decayer::kAll,
113  const bool do_svtx = true,
114  const bool do_cemc = true,
115  const bool do_hcalin = true,
116  const bool do_magnet = true,
117  const bool do_hcalout = true,
118  const bool do_pipe = true,
119  const bool do_FGEM = true,
120  const bool do_EGEM = true,
121  const bool do_FEMC = false,
122  const bool do_FHCAL = false,
123  const bool do_EEMC = true,
124  const bool do_DIRC = true,
125  const bool do_RICH = true,
126  const bool do_Aerogel = true,
127  const float magfield_rescale = 1.0) {
128 
129  //---------------
130  // Load libraries
131  //---------------
132 
133  gSystem->Load("libg4detectors.so");
134  gSystem->Load("libg4testbench.so");
135 
136  //---------------
137  // Fun4All server
138  //---------------
139 
140  Fun4AllServer *se = Fun4AllServer::instance();
141 
142  // read-in HepMC events to Geant4 if there is any
143  HepMCNodeReader *hr = new HepMCNodeReader();
144  se->registerSubsystem(hr);
145 
146  PHG4Reco* g4Reco = new PHG4Reco();
147  g4Reco->set_rapidity_coverage(1.1); // according to drawings
148 // uncomment to set QGSP_BERT_HP physics list for productions
149 // (default is QGSP_BERT for speed)
150  // g4Reco->SetPhysicsList("QGSP_BERT_HP");
151 
152  if (decayType != TPythia6Decayer::kAll) {
153  g4Reco->set_force_decay(decayType);
154  }
155 
156  double fieldstrength;
157  istringstream stringline(field);
158  stringline >> fieldstrength;
159  if (stringline.fail()) { // conversion to double fails -> we have a string
160 
161  if (field.find("sPHENIX.root") != string::npos) {
162  g4Reco->set_field_map(field, 1);
163  } else {
164  g4Reco->set_field_map(field, 2);
165  }
166  } else {
167  g4Reco->set_field(fieldstrength); // use const soleniodal field
168  }
169  g4Reco->set_field_rescale(magfield_rescale);
170 
171  double radius = 0.;
172 
173  //----------------------------------------
174  // PIPE
175  if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
176 
177  //----------------------------------------
178  // SVTX
179  if (do_svtx) radius = Svtx(g4Reco, radius, absorberactive);
180 
181  //----------------------------------------
182  // CEMC
183  //
184  if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
185  // if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
186 
187  //----------------------------------------
188  // HCALIN
189 
190  if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);
191 
192  //----------------------------------------
193  // MAGNET
194 
195  if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);
196 
197  //----------------------------------------
198  // HCALOUT
199 
200  if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);
201 
202  //----------------------------------------
203  // Forward tracking
204 
205  if ( do_FGEM )
206  FGEMSetup(g4Reco);
207 
208  if ( do_EGEM )
209  EGEMSetup(g4Reco);
210 
211  //----------------------------------------
212  // FEMC
213 
214  if ( do_FEMC )
215  FEMCSetup(g4Reco, absorberactive);
216 
217  //----------------------------------------
218  // FHCAL
219 
220  if ( do_FHCAL )
221  FHCALSetup(g4Reco, absorberactive);
222 
223  //----------------------------------------
224  // EEMC
225 
226  if ( do_EEMC )
227  EEMCSetup(g4Reco, absorberactive);
228 
229  //----------------------------------------
230  // PID
231 
232  if ( do_DIRC )
233  DIRCSetup(g4Reco);
234 
235  if ( do_RICH )
236  RICHSetup(g4Reco);
237 
238  if ( do_Aerogel )
239  AerogelSetup(g4Reco);
240 
241  // sPHENIX forward flux return(s)
242  PHG4CylinderSubsystem *flux_return_plus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
243  flux_return_plus->set_int_param("lengthviarapidity",0);
244  flux_return_plus->set_double_param("length",10.2);
245  flux_return_plus->set_double_param("radius",2.1);
246  flux_return_plus->set_double_param("thickness",263.5-5.0);
247  flux_return_plus->set_double_param("place_z",335.9);
248  flux_return_plus->set_string_param("material","G4_Fe");
249  flux_return_plus->SetActive(false);
250  flux_return_plus->SuperDetector("FLUXRET_ETA_PLUS");
251  flux_return_plus->OverlapCheck(overlapcheck);
252  g4Reco->registerSubsystem(flux_return_plus);
253 
254  PHG4CylinderSubsystem *flux_return_minus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
255  flux_return_minus->set_int_param("lengthviarapidity",0);
256  flux_return_minus->set_double_param("length",10.2);
257  flux_return_minus->set_double_param("radius",90.0);
258  flux_return_minus->set_double_param("place_z",-335.9);
259  flux_return_minus->set_double_param("thickness",263.5-5.0 - (90-2.1));
260  flux_return_minus->set_string_param("material","G4_Fe");
261  flux_return_minus->SetActive(false);
262  flux_return_minus->SuperDetector("FLUXRET_ETA_MINUS");
263  flux_return_minus->OverlapCheck(overlapcheck);
264  g4Reco->registerSubsystem(flux_return_minus);
265 
266  //----------------------------------------
267  // BLACKHOLE
268 
269  // swallow all particles coming out of the backend of sPHENIX
270  PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
271  blackhole->set_double_param("radius",radius + 100); // add 100 cm
272 
273  blackhole->set_int_param("lengthviarapidity",0);
274  blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
275  blackhole->BlackHole();
276  blackhole->set_double_param("thickness",0.1); // it needs some thickness
277  blackhole->SetActive(); // always see what leaks out
278  blackhole->OverlapCheck(overlapcheck);
279  g4Reco->registerSubsystem(blackhole);
280 
281  //----------------------------------------
282  // FORWARD BLACKHOLEs
283  // +Z
284  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
285  blackhole->SuperDetector("BH_FORWARD_PLUS");
286  blackhole->set_double_param("radius",0); // add 10 cm
287  blackhole->set_int_param("lengthviarapidity",0);
288  blackhole->set_double_param("length",0.1); // make it cover the world in length
289  blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1 - no_overlapp);
290  blackhole->BlackHole();
291  blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
292  blackhole->SetActive(); // always see what leaks out
293  blackhole->OverlapCheck(overlapcheck);
294  g4Reco->registerSubsystem(blackhole);
295 
296  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
297  blackhole->SuperDetector("BH_FORWARD_NEG");
298  blackhole->set_double_param("radius",0); // add 10 cm
299  blackhole->set_int_param("lengthviarapidity",0);
300  blackhole->set_double_param("length",0.1); // make it cover the world in length
301  blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1 + no_overlapp);
302  blackhole->BlackHole();
303  blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
304  blackhole->SetActive(); // always see what leaks out
305  blackhole->OverlapCheck(overlapcheck);
306  g4Reco->registerSubsystem(blackhole);
307 
308  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
309  g4Reco->registerSubsystem(truth);
310  se->registerSubsystem( g4Reco );
311 }
312 
313 
314 void ShowerCompress(int verbosity = 0) {
315 
316  gSystem->Load("libfun4all.so");
317  gSystem->Load("libg4eval.so");
318 
319  Fun4AllServer *se = Fun4AllServer::instance();
320 
321  PHG4DstCompressReco* compress = new PHG4DstCompressReco("PHG4DstCompressReco");
322  compress->AddHitContainer("G4HIT_PIPE");
323  compress->AddHitContainer("G4HIT_SVTXSUPPORT");
324  compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS");
325  compress->AddHitContainer("G4HIT_CEMC");
326  compress->AddHitContainer("G4HIT_ABSORBER_CEMC");
327  compress->AddHitContainer("G4HIT_CEMC_SPT");
328  compress->AddHitContainer("G4HIT_ABSORBER_HCALIN");
329  compress->AddHitContainer("G4HIT_HCALIN");
330  compress->AddHitContainer("G4HIT_HCALIN_SPT");
331  compress->AddHitContainer("G4HIT_MAGNET");
332  compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT");
333  compress->AddHitContainer("G4HIT_HCALOUT");
334  compress->AddHitContainer("G4HIT_BH_1");
335  compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS");
336  compress->AddHitContainer("G4HIT_BH_FORWARD_NEG");
337  compress->AddCellContainer("G4CELL_CEMC");
338  compress->AddCellContainer("G4CELL_HCALIN");
339  compress->AddCellContainer("G4CELL_HCALOUT");
340  compress->AddTowerContainer("TOWER_SIM_CEMC");
341  compress->AddTowerContainer("TOWER_RAW_CEMC");
342  compress->AddTowerContainer("TOWER_CALIB_CEMC");
343  compress->AddTowerContainer("TOWER_SIM_HCALIN");
344  compress->AddTowerContainer("TOWER_RAW_HCALIN");
345  compress->AddTowerContainer("TOWER_CALIB_HCALIN");
346  compress->AddTowerContainer("TOWER_SIM_HCALOUT");
347  compress->AddTowerContainer("TOWER_RAW_HCALOUT");
348  compress->AddTowerContainer("TOWER_CALIB_HCALOUT");
349 
350  compress->AddHitContainer("G4HIT_FEMC");
351  compress->AddHitContainer("G4HIT_ABSORBER_FEMC");
352  compress->AddHitContainer("G4HIT_FHCAL");
353  compress->AddHitContainer("G4HIT_ABSORBER_FHCAL");
354  compress->AddCellContainer("G4CELL_FEMC");
355  compress->AddCellContainer("G4CELL_FHCAL");
356  compress->AddTowerContainer("TOWER_SIM_FEMC");
357  compress->AddTowerContainer("TOWER_RAW_FEMC");
358  compress->AddTowerContainer("TOWER_CALIB_FEMC");
359  compress->AddTowerContainer("TOWER_SIM_FHCAL");
360  compress->AddTowerContainer("TOWER_RAW_FHCAL");
361  compress->AddTowerContainer("TOWER_CALIB_FHCAL");
362 
363  compress->AddHitContainer("G4HIT_EEMC");
364  compress->AddHitContainer("G4HIT_ABSORBER_EEMC");
365  compress->AddCellContainer("G4CELL_EEMC");
366  compress->AddTowerContainer("TOWER_SIM_EEMC");
367  compress->AddTowerContainer("TOWER_RAW_EEMC");
368  compress->AddTowerContainer("TOWER_CALIB_EEMC");
369 
370  se->registerSubsystem(compress);
371 
372  return;
373 }
374 
375 void DstCompress(Fun4AllDstOutputManager* out) {
376  if (out) {
377  out->StripNode("G4HIT_PIPE");
378  out->StripNode("G4HIT_SVTXSUPPORT");
379  out->StripNode("G4HIT_CEMC_ELECTRONICS");
380  out->StripNode("G4HIT_CEMC");
381  out->StripNode("G4HIT_ABSORBER_CEMC");
382  out->StripNode("G4HIT_CEMC_SPT");
383  out->StripNode("G4HIT_ABSORBER_HCALIN");
384  out->StripNode("G4HIT_HCALIN");
385  out->StripNode("G4HIT_HCALIN_SPT");
386  out->StripNode("G4HIT_MAGNET");
387  out->StripNode("G4HIT_ABSORBER_HCALOUT");
388  out->StripNode("G4HIT_HCALOUT");
389  out->StripNode("G4HIT_BH_1");
390  out->StripNode("G4HIT_BH_FORWARD_PLUS");
391  out->StripNode("G4HIT_BH_FORWARD_NEG");
392  out->StripNode("G4CELL_CEMC");
393  out->StripNode("G4CELL_HCALIN");
394  out->StripNode("G4CELL_HCALOUT");
395 
396  out->StripNode("G4HIT_FEMC");
397  out->StripNode("G4HIT_ABSORBER_FEMC");
398  out->StripNode("G4HIT_FHCAL");
399  out->StripNode("G4HIT_ABSORBER_FHCAL");
400  out->StripNode("G4CELL_FEMC");
401  out->StripNode("G4CELL_FHCAL");
402 
403  out->StripNode("G4HIT_EEMC");
404  out->StripNode("G4HIT_ABSORBER_EEMC");
405  out->StripNode("G4CELL_EEMC");
406  }
407 }