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  PHG4Reco* g4Reco = new PHG4Reco();
143  g4Reco->set_rapidity_coverage(1.1); // according to drawings
144 // uncomment to set QGSP_BERT_HP physics list for productions
145 // (default is QGSP_BERT for speed)
146  // g4Reco->SetPhysicsList("QGSP_BERT_HP");
147 
148  if (decayType != TPythia6Decayer::kAll) {
149  g4Reco->set_force_decay(decayType);
150  }
151 
152  double fieldstrength;
153  istringstream stringline(field);
154  stringline >> fieldstrength;
155  if (stringline.fail()) { // conversion to double fails -> we have a string
156 
157  if (field.find("sPHENIX.root") != string::npos) {
158  g4Reco->set_field_map(field, 1);
159  } else {
160  g4Reco->set_field_map(field, 2);
161  }
162  } else {
163  g4Reco->set_field(fieldstrength); // use const soleniodal field
164  }
165  g4Reco->set_field_rescale(magfield_rescale);
166 
167  double radius = 0.;
168 
169  //----------------------------------------
170  // PIPE
171  if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
172 
173  //----------------------------------------
174  // SVTX
175  if (do_svtx) radius = Svtx(g4Reco, radius, absorberactive);
176 
177  //----------------------------------------
178  // CEMC
179  //
180  if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
181  // if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
182 
183  //----------------------------------------
184  // HCALIN
185 
186  if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);
187 
188  //----------------------------------------
189  // MAGNET
190 
191  if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);
192 
193  //----------------------------------------
194  // HCALOUT
195 
196  if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);
197 
198  //----------------------------------------
199  // Forward tracking
200 
201  if ( do_FGEM )
202  FGEMSetup(g4Reco);
203 
204  if ( do_EGEM )
205  EGEMSetup(g4Reco);
206 
207  //----------------------------------------
208  // FEMC
209 
210  if ( do_FEMC )
211  FEMCSetup(g4Reco, absorberactive);
212 
213  //----------------------------------------
214  // FHCAL
215 
216  if ( do_FHCAL )
217  FHCALSetup(g4Reco, absorberactive);
218 
219  //----------------------------------------
220  // EEMC
221 
222  if ( do_EEMC )
223  EEMCSetup(g4Reco, absorberactive);
224 
225  //----------------------------------------
226  // PID
227 
228  if ( do_DIRC )
229  DIRCSetup(g4Reco);
230 
231  if ( do_RICH )
232  RICHSetup(g4Reco);
233 
234  if ( do_Aerogel )
235  AerogelSetup(g4Reco);
236 
237  // sPHENIX forward flux return(s)
238  PHG4CylinderSubsystem *flux_return_plus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
239  flux_return_plus->set_int_param("lengthviarapidity",0);
240  flux_return_plus->set_double_param("length",10.2);
241  flux_return_plus->set_double_param("radius",2.1);
242  flux_return_plus->set_double_param("thickness",263.5-5.0);
243  flux_return_plus->set_double_param("place_z",335.9);
244  flux_return_plus->set_string_param("material","G4_Fe");
245  flux_return_plus->SetActive(false);
246  flux_return_plus->SuperDetector("FLUXRET_ETA_PLUS");
247  flux_return_plus->OverlapCheck(overlapcheck);
248  g4Reco->registerSubsystem(flux_return_plus);
249 
250  PHG4CylinderSubsystem *flux_return_minus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
251  flux_return_minus->set_int_param("lengthviarapidity",0);
252  flux_return_minus->set_double_param("length",10.2);
253  flux_return_minus->set_double_param("radius",90.0);
254  flux_return_minus->set_double_param("place_z",-335.9);
255  flux_return_minus->set_double_param("thickness",263.5-5.0 - (90-2.1));
256  flux_return_minus->set_string_param("material","G4_Fe");
257  flux_return_minus->SetActive(false);
258  flux_return_minus->SuperDetector("FLUXRET_ETA_MINUS");
259  flux_return_minus->OverlapCheck(overlapcheck);
260  g4Reco->registerSubsystem(flux_return_minus);
261 
262  //----------------------------------------
263  // BLACKHOLE
264 
265  // swallow all particles coming out of the backend of sPHENIX
266  PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
267  blackhole->set_double_param("radius",radius + 100); // add 100 cm
268 
269  blackhole->set_int_param("lengthviarapidity",0);
270  blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
271  blackhole->BlackHole();
272  blackhole->set_double_param("thickness",0.1); // it needs some thickness
273  blackhole->SetActive(); // always see what leaks out
274  blackhole->OverlapCheck(overlapcheck);
275  g4Reco->registerSubsystem(blackhole);
276 
277  //----------------------------------------
278  // FORWARD BLACKHOLEs
279  // +Z
280  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
281  blackhole->SuperDetector("BH_FORWARD_PLUS");
282  blackhole->set_double_param("radius",0); // add 10 cm
283  blackhole->set_int_param("lengthviarapidity",0);
284  blackhole->set_double_param("length",0.1); // make it cover the world in length
285  blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1 - no_overlapp);
286  blackhole->BlackHole();
287  blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
288  blackhole->SetActive(); // always see what leaks out
289  blackhole->OverlapCheck(overlapcheck);
290  g4Reco->registerSubsystem(blackhole);
291 
292  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
293  blackhole->SuperDetector("BH_FORWARD_NEG");
294  blackhole->set_double_param("radius",0); // add 10 cm
295  blackhole->set_int_param("lengthviarapidity",0);
296  blackhole->set_double_param("length",0.1); // make it cover the world in length
297  blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1 + no_overlapp);
298  blackhole->BlackHole();
299  blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
300  blackhole->SetActive(); // always see what leaks out
301  blackhole->OverlapCheck(overlapcheck);
302  g4Reco->registerSubsystem(blackhole);
303 
304  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
305  g4Reco->registerSubsystem(truth);
306  se->registerSubsystem( g4Reco );
307 }
308 
309 
310 void ShowerCompress(int verbosity = 0) {
311 
312  gSystem->Load("libfun4all.so");
313  gSystem->Load("libg4eval.so");
314 
315  Fun4AllServer *se = Fun4AllServer::instance();
316 
317  PHG4DstCompressReco* compress = new PHG4DstCompressReco("PHG4DstCompressReco");
318  compress->AddHitContainer("G4HIT_PIPE");
319  compress->AddHitContainer("G4HIT_SVTXSUPPORT");
320  compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS");
321  compress->AddHitContainer("G4HIT_CEMC");
322  compress->AddHitContainer("G4HIT_ABSORBER_CEMC");
323  compress->AddHitContainer("G4HIT_CEMC_SPT");
324  compress->AddHitContainer("G4HIT_ABSORBER_HCALIN");
325  compress->AddHitContainer("G4HIT_HCALIN");
326  compress->AddHitContainer("G4HIT_HCALIN_SPT");
327  compress->AddHitContainer("G4HIT_MAGNET");
328  compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT");
329  compress->AddHitContainer("G4HIT_HCALOUT");
330  compress->AddHitContainer("G4HIT_BH_1");
331  compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS");
332  compress->AddHitContainer("G4HIT_BH_FORWARD_NEG");
333  compress->AddCellContainer("G4CELL_CEMC");
334  compress->AddCellContainer("G4CELL_HCALIN");
335  compress->AddCellContainer("G4CELL_HCALOUT");
336  compress->AddTowerContainer("TOWER_SIM_CEMC");
337  compress->AddTowerContainer("TOWER_RAW_CEMC");
338  compress->AddTowerContainer("TOWER_CALIB_CEMC");
339  compress->AddTowerContainer("TOWER_SIM_HCALIN");
340  compress->AddTowerContainer("TOWER_RAW_HCALIN");
341  compress->AddTowerContainer("TOWER_CALIB_HCALIN");
342  compress->AddTowerContainer("TOWER_SIM_HCALOUT");
343  compress->AddTowerContainer("TOWER_RAW_HCALOUT");
344  compress->AddTowerContainer("TOWER_CALIB_HCALOUT");
345 
346  compress->AddHitContainer("G4HIT_FEMC");
347  compress->AddHitContainer("G4HIT_ABSORBER_FEMC");
348  compress->AddHitContainer("G4HIT_FHCAL");
349  compress->AddHitContainer("G4HIT_ABSORBER_FHCAL");
350  compress->AddCellContainer("G4CELL_FEMC");
351  compress->AddCellContainer("G4CELL_FHCAL");
352  compress->AddTowerContainer("TOWER_SIM_FEMC");
353  compress->AddTowerContainer("TOWER_RAW_FEMC");
354  compress->AddTowerContainer("TOWER_CALIB_FEMC");
355  compress->AddTowerContainer("TOWER_SIM_FHCAL");
356  compress->AddTowerContainer("TOWER_RAW_FHCAL");
357  compress->AddTowerContainer("TOWER_CALIB_FHCAL");
358 
359  compress->AddHitContainer("G4HIT_EEMC");
360  compress->AddHitContainer("G4HIT_ABSORBER_EEMC");
361  compress->AddCellContainer("G4CELL_EEMC");
362  compress->AddTowerContainer("TOWER_SIM_EEMC");
363  compress->AddTowerContainer("TOWER_RAW_EEMC");
364  compress->AddTowerContainer("TOWER_CALIB_EEMC");
365 
366  se->registerSubsystem(compress);
367 
368  return;
369 }
370 
371 void DstCompress(Fun4AllDstOutputManager* out) {
372  if (out) {
373  out->StripNode("G4HIT_PIPE");
374  out->StripNode("G4HIT_SVTXSUPPORT");
375  out->StripNode("G4HIT_CEMC_ELECTRONICS");
376  out->StripNode("G4HIT_CEMC");
377  out->StripNode("G4HIT_ABSORBER_CEMC");
378  out->StripNode("G4HIT_CEMC_SPT");
379  out->StripNode("G4HIT_ABSORBER_HCALIN");
380  out->StripNode("G4HIT_HCALIN");
381  out->StripNode("G4HIT_HCALIN_SPT");
382  out->StripNode("G4HIT_MAGNET");
383  out->StripNode("G4HIT_ABSORBER_HCALOUT");
384  out->StripNode("G4HIT_HCALOUT");
385  out->StripNode("G4HIT_BH_1");
386  out->StripNode("G4HIT_BH_FORWARD_PLUS");
387  out->StripNode("G4HIT_BH_FORWARD_NEG");
388  out->StripNode("G4CELL_CEMC");
389  out->StripNode("G4CELL_HCALIN");
390  out->StripNode("G4CELL_HCALOUT");
391 
392  out->StripNode("G4HIT_FEMC");
393  out->StripNode("G4HIT_ABSORBER_FEMC");
394  out->StripNode("G4HIT_FHCAL");
395  out->StripNode("G4HIT_ABSORBER_FHCAL");
396  out->StripNode("G4CELL_FEMC");
397  out->StripNode("G4CELL_FHCAL");
398 
399  out->StripNode("G4HIT_EEMC");
400  out->StripNode("G4HIT_ABSORBER_EEMC");
401  out->StripNode("G4CELL_EEMC");
402  }
403 }