EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_hFarFwdBeamLine_EIC.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_hFarFwdBeamLine_EIC.C
1 #ifndef MACRO_G4HFARFWDBEAMLINE_EIC_C
2 #define MACRO_G4HFARFWDBEAMLINE_EIC_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <g4detectors/BeamLineMagnetSubsystem.h>
7 #include <g4detectors/PHG4BlockSubsystem.h>
8 #include <g4detectors/PHG4ConeSubsystem.h>
9 #include <g4detectors/PHG4CylinderSubsystem.h>
10 
11 #include <eicg4zdc/EICG4ZDCHitTree.h>
12 #include <eicg4zdc/EICG4ZDCNtuple.h>
13 #include <eicg4zdc/EICG4ZDCSubsystem.h>
14 
15 #include <g4main/PHG4Reco.h>
16 
17 #include <TSystem.h>
18 
19 R__LOAD_LIBRARY(libEICG4ZDC.so)
20 
21 // This creates the Enable Flag to be used in the main steering macro
22 namespace Enable
23 {
24  bool HFARFWD_MAGNETS = false;
25  bool HFARFWD_VIRTUAL_DETECTORS = false;
26 
27  bool HFARFWD_PIPE = false;
28  bool HFARFWD_OVERLAPCHECK = false;
29  int HFARFWD_VERBOSITY = 0;
30  bool ZDC_DISABLE_BLACKHOLE = false;
31 
32  //enabled automatically in hFarFwdBeamLineInit(), unless overridden by user
33  bool HFARFWD_MAGNETS_IP6 = false;
34  bool HFARFWD_MAGNETS_IP8 = false;
35 
36  //enabled automatically in hFarFwdBeamLineInit(), unless overridden by user
37  bool HFARFWD_VIRTUAL_DETECTORS_IP6 = false;
38  bool HFARFWD_VIRTUAL_DETECTORS_IP8 = false;
39  bool IP6 = false;
40  bool IP8 = false;
41 
42 
43 } // namespace Enable
44 
45 namespace hFarFwdBeamLine
46 {
47  double starting_z = 450; //cm as center-forward interface
48  double enclosure_z_max = NAN;
49  double enclosure_r_max = NAN;
50  double enclosure_center = NAN;
51 
52  PHG4CylinderSubsystem *hFarFwdBeamLineEnclosure(nullptr);
53 
54  BeamLineMagnetSubsystem *B0Magnet = (nullptr);
55 } // namespace hFarFwdBeamLine
56 
57 void hFarFwdBeamLineInit()
58 {
59  Enable::HFARFWD_MAGNETS_IP6 |= Enable::HFARFWD_MAGNETS and Enable::IP6;
60  Enable::HFARFWD_MAGNETS_IP8 |= Enable::HFARFWD_MAGNETS and Enable::IP8;
61 
62  Enable::HFARFWD_VIRTUAL_DETECTORS_IP6 |= Enable::HFARFWD_VIRTUAL_DETECTORS and Enable::IP6;
63  Enable::HFARFWD_VIRTUAL_DETECTORS_IP8 |= Enable::HFARFWD_VIRTUAL_DETECTORS and Enable::IP8;
64 
65  if (Enable::HFARFWD_MAGNETS_IP6 && Enable::HFARFWD_MAGNETS_IP8)
66  {
67  cout << "You cannot have magnets for both IP6 and IP8 ON at the same time" << endl;
68  gSystem->Exit(1);
69  }
70 
71  if (Enable::HFARFWD_MAGNETS_IP6)
72  {
73  hFarFwdBeamLine::enclosure_z_max = 4500.;
74  BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, hFarFwdBeamLine::starting_z);
75  hFarFwdBeamLine::enclosure_r_max = 200.;
76  }
77 
78  if (Enable::HFARFWD_MAGNETS_IP8)
79  {
80  hFarFwdBeamLine::enclosure_z_max = 4800.;
81  BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -2050.);
82  hFarFwdBeamLine::enclosure_r_max = 150.;
83  }
84 
85  hFarFwdBeamLine::enclosure_center = 0.5 * (hFarFwdBeamLine::starting_z + hFarFwdBeamLine::enclosure_z_max);
86 
87  BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, hFarFwdBeamLine::enclosure_z_max);
88  BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, hFarFwdBeamLine::enclosure_r_max);
89 }
90 
91 void hFarFwdDefineMagnets(PHG4Reco *g4Reco)
92 {
93  bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARFWD_OVERLAPCHECK;
94  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
95 
96  hFarFwdBeamLine::hFarFwdBeamLineEnclosure = new PHG4CylinderSubsystem("hFarFwdBeamLineEnclosure");
97  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("place_z", hFarFwdBeamLine::enclosure_center);
98  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("radius", 0);
99  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("thickness", hFarFwdBeamLine::enclosure_r_max); // This is intentionally made large 25cm radius
100  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("length", hFarFwdBeamLine::enclosure_z_max - hFarFwdBeamLine::starting_z);
101  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_string_param("material", "G4_Galactic");
102  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_color(.5, .5, .5, 0.2);
103  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->OverlapCheck(overlapCheck);
104  if (verbosity)
105  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->Verbosity(verbosity);
106  g4Reco->registerSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
107 
108  string magFile;
109  if (Enable::HFARFWD_MAGNETS_IP6)
110  magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_h_farFwdBeamLineMagnets.dat";
111  else if (Enable::HFARFWD_MAGNETS_IP8)
112  magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip8_35mrad_h_farFwdBeamLineMagnets.dat";
113  else
114  {
115  cout << " You have to enable either the IP6 or IP8 Magnet configuration to define magnets! " << endl;
116  gSystem->Exit(1);
117  }
118 
119  // make magnet active volume if you want to study the hits
120  bool magnet_active = false;
121  int absorberactive = 0;
122 
123  // if you insert numbers it only displays those magnets, do not comment out the set declaration
124  set<int> magnetlist;
125  //magnetlist.insert(7);
126 
127  BeamLineMagnetSubsystem *bl = nullptr;
128  std::ifstream infile(magFile);
129  if (infile.is_open())
130  {
131  double biggest_z = 0.;
132  int imagnet = 0;
133  std::string line;
134  while (std::getline(infile, line))
135  {
136  if (!line.compare(0, 1, "B") ||
137  !line.compare(0, 1, "Q") ||
138  !line.compare(0, 1, "S"))
139  {
140  std::istringstream iss(line);
141  string magname;
142  double x;
143  double y;
144  double z;
145  double inner_radius_zin;
146  double inner_radius_zout;
147  double outer_magnet_diameter;
148  double length;
149  double angle;
150  double dipole_field_x;
151  double fieldgradient;
152  if (!(iss >> magname >> x >> y >> z >> inner_radius_zin >> inner_radius_zout >> outer_magnet_diameter >> length >> angle >> dipole_field_x >> fieldgradient))
153  {
154  cout << "coud not decode " << line << endl;
155  gSystem->Exit(1);
156  }
157  else
158  {
159  string magtype;
160  if (inner_radius_zin != inner_radius_zout)
161  {
162  cout << "inner radius at front of magnet " << inner_radius_zin
163  << " not equal radius at back of magnet " << inner_radius_zout
164  << " needs change in code (replace tube by cone for beamline)" << endl;
165  gSystem->Exit(1);
166  }
167  if (verbosity > 0)
168  {
169  cout << endl
170  << endl
171  << "\tID number " << imagnet << endl;
172  cout << "magname: " << magname << endl;
173  cout << "x: " << x << endl;
174  cout << "y: " << y << endl;
175  cout << "z: " << z << endl;
176  cout << "inner_radius_zin: " << inner_radius_zin << endl;
177  cout << "inner_radius_zout: " << inner_radius_zout << endl;
178  cout << "outer_magnet_diameter: " << outer_magnet_diameter << endl;
179  cout << "length: " << length << endl;
180  cout << "angle: " << angle << endl;
181  cout << "dipole_field_x: " << dipole_field_x << endl;
182  cout << "fieldgradient: " << fieldgradient << endl;
183  }
184  if (!magname.compare(0, 1, "B"))
185  {
186  magtype = "DIPOLE";
187  }
188  else if (!magname.compare(0, 1, "Q"))
189  {
190  magtype = "QUADRUPOLE";
191  }
192  else if (!magname.compare(0, 1, "S"))
193  {
194  magtype = "SEXTUPOLE";
195  }
196  else
197  {
198  cout << "cannot decode magnet name " << magname << endl;
199  gSystem->Exit(1);
200  }
201  // convert to our units (cm, deg)
202  x *= 100.;
203  y *= 100.;
204  z *= 100.;
205  length *= 100.;
206  inner_radius_zin *= 100.;
207  outer_magnet_diameter *= 100.;
208  angle = (angle / TMath::Pi() * 180.) / 1000.; // given in mrad
209 
210  if (magnetlist.empty() || magnetlist.find(imagnet) != magnetlist.end())
211  {
212  bl = new BeamLineMagnetSubsystem("BEAMLINEMAGNET", imagnet);
213  bl->set_double_param("field_y", dipole_field_x);
214  bl->set_double_param("fieldgradient", fieldgradient);
215  bl->set_string_param("magtype", magtype);
216  bl->set_double_param("length", length);
217  bl->set_double_param("place_x", x);// relative position to mother vol.
218  bl->set_double_param("place_y", y);// relative position to mother vol.
219  bl->set_double_param("place_z", z - hFarFwdBeamLine::enclosure_center);// relative position to mother vol.
220  bl->set_double_param("field_global_position_x", x);// abs. position to world for field manager
221  bl->set_double_param("field_global_position_y", y);// abs. position to world for field manager
222  bl->set_double_param("field_global_position_z", z);// abs. position to world for field manager
223  bl->set_double_param("rot_y", angle);
224  bl->set_double_param("field_global_rot_y", angle);// abs. rotation to world for field manager
225  bl->set_double_param("inner_radius", inner_radius_zin);
226  bl->set_double_param("outer_radius", outer_magnet_diameter / 2.);
227  bl->SetActive(magnet_active);
228  bl->BlackHole();
229  bl->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
230  if (absorberactive)
231  {
232  bl->SetAbsorberActive();
233  }
234  bl->OverlapCheck(overlapCheck);
235  bl->SuperDetector("BEAMLINEMAGNET");
236  if (verbosity)
237  bl->Verbosity(verbosity);
238  g4Reco->registerSubsystem(bl);
239 
240  // rag the B0 magnet
241  if (imagnet == 0)
242  hFarFwdBeamLine::B0Magnet = bl;
243  }
244  imagnet++;
245  if (fabs(z) + length > biggest_z)
246  {
247  biggest_z = fabs(z) + length;
248  }
249  }
250  }
251  }
252  infile.close();
253  }
254 }
255 
256 void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco)
257 {
258  bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARFWD_OVERLAPCHECK;
259  if (Enable::HFARFWD_VIRTUAL_DETECTORS_IP6 && Enable::HFARFWD_VIRTUAL_DETECTORS_IP8)
260  {
261  cout << "You cannot have detectors enabled for both IP6 and IP8 ON at the same time" << endl;
262  gSystem->Exit(1);
263  }
264 
265  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
266 
267  auto *detZDCsurrogate = new PHG4BlockSubsystem("zdcTruth");
268  const double detZDCsurrogate_size_z = 0.1;
269  detZDCsurrogate->SuperDetector("ZDCsurrogate");
270  detZDCsurrogate->set_double_param("place_x", 96.24);
271  detZDCsurrogate->set_double_param("place_y", 0);
272  detZDCsurrogate->set_double_param("place_z", 3750 - hFarFwdBeamLine::enclosure_center);
273  detZDCsurrogate->set_double_param("rot_y", -0.025 * TMath::RadToDeg());
274  detZDCsurrogate->set_double_param("size_x", 60);
275  detZDCsurrogate->set_double_param("size_y", 60);
276  detZDCsurrogate->set_double_param("size_z", detZDCsurrogate_size_z);
277  detZDCsurrogate->set_string_param("material", "G4_Si");
278  detZDCsurrogate->SetActive();
279  detZDCsurrogate->set_color(1, 0, 0, 0.5);
280  detZDCsurrogate->OverlapCheck(overlapCheck);
281  if (!Enable::ZDC_DISABLE_BLACKHOLE) detZDCsurrogate->BlackHole();
282  if (verbosity)
283  detZDCsurrogate->Verbosity(verbosity);
284  detZDCsurrogate->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
285  g4Reco->registerSubsystem(detZDCsurrogate);
286 
287  if (Enable::ZDC_DISABLE_BLACKHOLE)
288  {
289  EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC");
290  detZDC->SetActive();
291  detZDC->set_double_param("place_z", 3750. + detZDCsurrogate_size_z - hFarFwdBeamLine::enclosure_center);
292  detZDC->set_double_param("place_x", 96.24);
293  detZDC->set_double_param("rot_y", -0.025);
294  detZDC->OverlapCheck(overlapCheck);
295  detZDC->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
296  g4Reco->registerSubsystem(detZDC);
297  }
298 
299  const int offMomDetNr = 2;
300  const double om_zCent[offMomDetNr] = {3450, 3650};
301  const double om_xCent[offMomDetNr] = {162, 171};
302  for (int i = 0; i < offMomDetNr; i++)
303  {
304  auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i), i);
305  detOM->SuperDetector("offMomTruth");
306  detOM->set_double_param("place_x", om_xCent[i]);
307  detOM->set_double_param("place_y", 0);
308  detOM->set_double_param("place_z", om_zCent[i] - hFarFwdBeamLine::enclosure_center);
309  detOM->set_double_param("rot_y", -0.045 * TMath::RadToDeg());
310  detOM->set_double_param("size_x", 50);
311  detOM->set_double_param("size_y", 35);
312  detOM->set_double_param("size_z", 0.03);
313  detOM->set_string_param("material", "G4_Si");
314  detOM->SetActive();
315  if (verbosity)
316  detOM->Verbosity(verbosity);
317  detOM->OverlapCheck(overlapCheck);
318  detOM->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
319  g4Reco->registerSubsystem(detOM);
320  }
321 
322  const int rpDetNr = 2;
323  const double rp_zCent[rpDetNr] = {2600, 2800};
324  const double rp_xCent[rpDetNr] = {84.49, 93.59};
325  for (int i = 0; i < rpDetNr; i++)
326  {
330  //
331  // auto *detRP = new PHG4BlockSubsystem(Form("rpTruth_%d",i));
333  // detRP->SuperDetector(Form("RomanPots_%d",i));
334  // detRP->set_double_param("place_x",rp_xCent[i]);
335  // detRP->set_double_param("place_y",0);
336  // detRP->set_double_param("place_z",rp_zCent[i]);
337  // detRP->set_double_param("rot_y",-0.025*TMath::RadToDeg());
338  // detRP->set_double_param("size_x",25);
339  // detRP->set_double_param("size_y",10);
340  // detRP->set_double_param("size_z",0.03);
341  // detRP->set_string_param("material","G4_Si");
342 
346 
347  auto *detRP = new PHG4CylinderSubsystem(Form("rpTruth_%d", i), i);
348  detRP->SuperDetector("rpTruth");
349  detRP->set_double_param("place_x", rp_xCent[i]);
350  detRP->set_double_param("place_y", 0);
351  detRP->set_double_param("place_z", rp_zCent[i] - hFarFwdBeamLine::enclosure_center);
352  detRP->set_double_param("rot_y", -0.025 * TMath::RadToDeg());
353  detRP->set_double_param("radius", 0);
354  detRP->set_double_param("thickness", 25); // This is intentionally made large 25cm radius
355  detRP->set_double_param("length", 0.03);
356  detRP->set_string_param("material", "G4_Si");
357  detRP->OverlapCheck(overlapCheck);
358  detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
359 
360  detRP->SetActive();
361  if (verbosity)
362  detRP->Verbosity(verbosity);
363  g4Reco->registerSubsystem(detRP);
364  }
365 
366  const int b0DetNr = 4;
367  const double b0Mag_zCent = 590;
368  const double b0Mag_zLen = 120;
369  for (int i = 0; i < b0DetNr; i++)
370  {
371  auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), i);
372  detB0->SuperDetector("b0Truth");
373  detB0->set_double_param("radius", 0);
374  detB0->set_double_param("thickness", 20);
375  detB0->set_double_param("length", 0.1);
376  detB0->set_string_param("material", "G4_Si");
377  detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2)); // relative to B0 magnet
378  detB0->SetActive(true);
379  if (verbosity)
380  detB0->Verbosity(verbosity);
381  detB0->OverlapCheck(overlapCheck);
382  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
383  g4Reco->registerSubsystem(detB0);
384  }
385 }
386 
387 void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco)
388 {
389  cout << __PRETTY_FUNCTION__ << " : IP8 setup is not yet validated!" << endl;
390  gSystem->Exit(1);
391 
392  bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARFWD_OVERLAPCHECK;
393  if (Enable::HFARFWD_VIRTUAL_DETECTORS_IP6 && Enable::HFARFWD_VIRTUAL_DETECTORS_IP8)
394  {
395  cout << "You cannot have detectors enabled for both IP6 and IP8 ON at the same time" << endl;
396  gSystem->Exit(1);
397  }
398 
399  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
400 
401  const int offMomDetNr = 3;
402  const double om_zCent[offMomDetNr] = {4250, 4400, 4550};
403  const double om_xCent[offMomDetNr] = {100, 100, 100};
404  for (int i = 0; i < offMomDetNr; i++)
405  {
406  auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i));
407  detOM->set_double_param("place_x", om_xCent[i]);
408  detOM->set_double_param("place_y", 0);
409  detOM->set_double_param("place_z", om_zCent[i]);
410  detOM->set_double_param("rot_y", -0.029 * TMath::RadToDeg());
411  detOM->set_double_param("size_x", 100);
412  detOM->set_double_param("size_y", 100);
413  detOM->set_double_param("size_z", 0.03);
414  detOM->set_string_param("material", "G4_Si");
415  detOM->SetActive();
416  if (verbosity)
417  detOM->Verbosity(verbosity);
418  g4Reco->registerSubsystem(detOM);
419  }
420 
421  auto *detZDCsurrogate = new PHG4BlockSubsystem("zdcTruth");
422  const double detZDCsurrogate_size_z = 0.1;
423  detZDCsurrogate->SuperDetector("ZDCsurrogate");
424  detZDCsurrogate->set_double_param("place_x", 127.8);
425  detZDCsurrogate->set_double_param("place_y", 0);
426  detZDCsurrogate->set_double_param("place_z", 3650);
427  detZDCsurrogate->set_double_param("rot_y", -0.035 * TMath::RadToDeg());
428  detZDCsurrogate->set_double_param("size_x", 60);
429  detZDCsurrogate->set_double_param("size_y", 60);
430  detZDCsurrogate->set_double_param("size_z", detZDCsurrogate_size_z);
431  detZDCsurrogate->set_string_param("material", "G4_Si");
432  detZDCsurrogate->SetActive();
433  detZDCsurrogate->OverlapCheck(overlapCheck);
434  detZDCsurrogate->set_color(1, 0, 0, 0.5);
435  if (!Enable::ZDC_DISABLE_BLACKHOLE) detZDCsurrogate->BlackHole();
436  if (verbosity)
437  detZDCsurrogate->Verbosity(verbosity);
438  g4Reco->registerSubsystem(detZDCsurrogate);
439 
440  EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC");
441  detZDC->SetActive();
442  detZDC->set_double_param("place_z", 3650. + detZDCsurrogate_size_z);
443  detZDC->set_double_param("place_x", 127.8);
444  detZDC->set_double_param("rot_y", -0.035);
445  detZDC->OverlapCheck(overlapCheck);
446  g4Reco->registerSubsystem(detZDC);
447 
448  const int rpDetNr = 4;
449  const double rp_zCent[rpDetNr] = {2200, 2500, 2800, 3100};
450  const double rp_xCent[rpDetNr] = {75, 75, 75, 75};
451  for (int i = 0; i < rpDetNr; i++)
452  {
453  auto *detRP = new PHG4BlockSubsystem(Form("rpTruth_%d", i));
454  // detRP->SuperDetector("RomanPots");
455  detRP->SuperDetector(Form("RomanPots_%d", i));
456  detRP->set_double_param("place_x", rp_xCent[i]);
457  detRP->set_double_param("place_y", 0);
458  detRP->set_double_param("place_z", rp_zCent[i]);
459  detRP->set_double_param("rot_y", -0.0215 * TMath::RadToDeg());
460  detRP->set_double_param("size_x", 100);
461  detRP->set_double_param("size_y", 100);
462  detRP->set_double_param("size_z", 0.03);
463  detRP->set_string_param("material", "G4_Si");
464  detRP->SetActive();
465  if (verbosity)
466  detRP->Verbosity(verbosity);
467  g4Reco->registerSubsystem(detRP);
468  }
469 
470  const int b0DetNr = 4;
471  const double b0Mag_zCent = 610;
472  const double b0Mag_zLen = 120;
473  for (int i = 0; i < b0DetNr; i++)
474  {
475  auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), 0);
476  //detB0->SuperDetector("B0detectors");
477  detB0->set_double_param("radius", 0);
478  detB0->set_double_param("thickness", 20);
479  detB0->set_double_param("length", 0.1);
480  detB0->set_string_param("material", "G4_Si");
481  detB0->set_double_param("place_x", 21.2);
482  detB0->set_double_param("place_y", 0);
483  detB0->set_double_param("place_z", (b0Mag_zCent - b0Mag_zLen / 2) + b0Mag_zLen / (b0DetNr - 1) * i);
484  detB0->SetActive(true);
485  if (verbosity)
486  detB0->Verbosity(verbosity);
487  g4Reco->registerSubsystem(detB0);
488  }
489 }
490 
491 void hFarFwdDefineBeamPipe(PHG4Reco *g4Reco)
492 {
493  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
494  //exit window
495  PHG4CylinderSubsystem *exitWin = new PHG4CylinderSubsystem("exitWin", 0);
496  exitWin->set_double_param("radius", 3.2);
497  exitWin->set_double_param("thickness", 11.8);
498  exitWin->set_double_param("length", 0.15);
499  exitWin->set_double_param("rot_y", -0.025 * TMath::RadToDeg());
500  exitWin->set_string_param("material", "G4_STAINLESS-STEEL");
501  exitWin->set_double_param("place_x", 12.5);
502  exitWin->set_double_param("place_y", 0);
503  exitWin->set_double_param("place_z", 500);
504  exitWin->SetActive(false);
505  g4Reco->registerSubsystem(exitWin);
506 
507  //B0 magnet pipe
508  PHG4CylinderSubsystem *pipeB0 = new PHG4CylinderSubsystem("beamPipeB0", 0);
509  pipeB0->set_double_param("radius", 2.8);
510  pipeB0->set_double_param("thickness", 0.25);
511  pipeB0->set_double_param("length", 195);
512  pipeB0->set_double_param("rot_y", -0.025 * TMath::RadToDeg());
513  pipeB0->set_string_param("material", "G4_Al");
514  pipeB0->set_double_param("place_x", 14.748);
515  pipeB0->set_double_param("place_y", 0);
516  pipeB0->set_double_param("place_z", 590);
517  pipeB0->SetActive(false);
518  g4Reco->registerSubsystem(pipeB0);
519 
520  //Quad pipes
521  const int nSecQ = 5; //B0apf, Q1apf, Q1bpf, Q2pf, B1pf
522  const string nm[nSecQ] = {"B0apf", "Q1apf", "Q1bpf", "Q2pf", "B1pf"};
523  const double qlen[nSecQ] = {160, 150, 220, 440, 330};
524  const double qir[nSecQ] = {4, 5.1, 7, 12, 12.2};
525  const double qor[nSecQ] = {4.1, 5.2, 7.2, 12.2, 12.4};
526  const double qrot[nSecQ] = {25, 19.5, 15, 15, 34}; //mrad
527  const double qxC[nSecQ] = {19.8, 24.47, 30.05, 39.5, 48};
528  const double qyC[nSecQ] = {0, 0, 0, 0, 0};
529  const double qzC[nSecQ] = {770, 922.8, 1106.3, 1416.7, 1806.7};
530  for (int i = 0; i < nSecQ; i++)
531  {
532  PHG4CylinderSubsystem *pipe = new PHG4CylinderSubsystem(Form("beamPipe%s", nm[i].c_str()), 0);
533  pipe->set_double_param("radius", qir[i]);
534  pipe->set_double_param("thickness", qor[i] - qir[i]);
535  pipe->set_double_param("length", qlen[i]);
536  pipe->set_double_param("rot_y", -qrot[i] / 1000 * TMath::RadToDeg());
537  pipe->set_string_param("material", "G4_Al");
538  pipe->set_double_param("place_x", qxC[i]);
539  pipe->set_double_param("place_y", qyC[i]);
540  pipe->set_double_param("place_z", qzC[i]);
541  pipe->SetActive(false);
542  // pipe->SetActive(true);
543  g4Reco->registerSubsystem(pipe);
544  }
545 
546  //Electron pipe
547  PHG4CylinderSubsystem *pipeElectron = new PHG4CylinderSubsystem("beamPipeElectron", 0);
548  pipeElectron->set_double_param("radius", 1);
549  pipeElectron->set_double_param("thickness", 1);
550  pipeElectron->set_double_param("length", 3000);
551  pipeElectron->set_double_param("rot_y", -0.025 * TMath::RadToDeg());
552  pipeElectron->set_string_param("material", "G4_Al");
553  pipeElectron->set_double_param("place_x", 0);
554  pipeElectron->set_double_param("place_y", 0);
555  pipeElectron->set_double_param("place_z", 2000);
556  pipeElectron->SetActive(false);
557  //g4Reco->registerSubsystem(pipeElectron);
558 
559  //ZDC pipe
560  PHG4CylinderSubsystem *pipeZDC = new PHG4CylinderSubsystem("beamPipeZDC", 0);
561  pipeZDC->set_double_param("radius", 16.5);
562  pipeZDC->set_double_param("thickness", 0.1);
563  pipeZDC->set_double_param("length", 170);
564  pipeZDC->set_double_param("rot_y", -0.025 * TMath::RadToDeg());
565  pipeZDC->set_string_param("material", "G4_Al");
566  pipeZDC->set_double_param("place_x", 59);
567  pipeZDC->set_double_param("place_y", 0);
568  pipeZDC->set_double_param("place_z", 2041.59);
569  pipeZDC->SetActive(false);
570  // pipeZDC->SetActive(true);
571  g4Reco->registerSubsystem(pipeZDC);
572 
573  //Roman Pot pipe
574  const int nSec = 2;
575  const double len[nSec] = {850, 1150};
576  const double ir1[nSec] = {17, 17};
577  const double or1[nSec] = {17.1, 17.1};
578  const double ir2[nSec] = {17, 7};
579  const double or2[nSec] = {17.1, 7.1};
580  const double xC[nSec] = {83, 130};
581  const double yC[nSec] = {0, 0};
582  const double zC[nSec] = {2550, 3550};
583  for (int i = 0; i < nSec; i++)
584  {
585  PHG4ConeSubsystem *pipe = new PHG4ConeSubsystem(Form("beamPipeRP%d", i), 0);
586  pipe->set_string_param("material", "G4_STAINLESS-STEEL");
587  pipe->set_double_param("place_x", xC[i]);
588  pipe->set_double_param("place_y", yC[i]);
589  pipe->set_double_param("place_z", zC[i]);
590  pipe->set_double_param("length", len[i] / 2);
591  pipe->set_double_param("rmin1", ir1[i]);
592  pipe->set_double_param("rmin2", ir2[i]);
593  pipe->set_double_param("rmax1", or1[i]);
594  pipe->set_double_param("rmax2", or2[i]);
595  pipe->set_double_param("rot_y", -0.047 * TMath::RadToDeg());
596  g4Reco->registerSubsystem(pipe);
597  }
598 }
599 
600 #endif