EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_BeamLine.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_BeamLine.C
1 #ifndef MACRO_G4BEAMLINE_C
2 #define MACRO_G4BEAMLINE_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <G4_Pipe.C>
7 
8 #include <g4detectors/BeamLineMagnetSubsystem.h>
9 #include <g4detectors/PHG4BlockSubsystem.h>
10 #include <g4detectors/PHG4ConeSubsystem.h>
11 #include <g4detectors/PHG4CylinderSubsystem.h>
12 #include <g4detectors/PHG4ZDCSubsystem.h>
13 #include <g4detectors/PHG4DetectorSubsystem.h>
14 
15 #include <g4main/PHG4Reco.h>
16 
17 #include <TSystem.h>
18 
19 R__LOAD_LIBRARY(libg4detectors.so)
20 
21 float PosFlip(float pos);
22 float AngleFlip(float angle);
23 float MagFieldFlip(float Bfield);
24 
25 // This creates the Enable Flag to be used in the main steering macro
26 namespace Enable
27 {
28  bool BEAMLINE = false;
29  bool BEAMLINE_ABSORBER = false;
30  bool BEAMLINE_BLACKHOLE = false;
31  bool BEAMLINE_OVERLAPCHECK = false;
32  int BEAMLINE_VERBOSITY = 0;
33 
34 } // namespace Enable
35 
36 namespace G4BEAMLINE
37 {
38 // the beampipes seem to add 2 no_overlaps - needs to be looked at
39 // but this z position takes care of our current overlap issues
40  double starting_z = G4PIPE::be_pipe_length / 2. + G4PIPE::al_pipe_length + G4PIPE::al_pipe_cone_length + G4PIPE::al_pipe_ext_length + 2*no_overlapp;
41  double enclosure_z_max = 2050. + (700-starting_z);
42  double enclosure_r_max = 30.; // 30cm radius to cover magnets
43  double enclosure_center = 0.5 * (starting_z + enclosure_z_max);
44  double skin_thickness = 0.; // if center of magnet iron is black hole - thickness of Fe surrounding it
45  int pipe_id_offset = 100;
46  int roman_pot_pipe_id_offset = 200;
47  PHG4CylinderSubsystem *ForwardBeamLineEnclosure(nullptr);
48  PHG4CylinderSubsystem *BackwardBeamLineEnclosure(nullptr);
49 
50 } // namespace BeamLine
51 
52 void BeamLineInit()
53 {
54  BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -G4BEAMLINE::enclosure_z_max);
55  BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4BEAMLINE::enclosure_z_max);
56  BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4BEAMLINE::enclosure_r_max);
57 }
58 
59 void BeamLineDefineMagnets(PHG4Reco *g4Reco)
60 {
61  bool overlapCheck = Enable::OVERLAPCHECK || Enable::BEAMLINE_OVERLAPCHECK;
62  bool AbsorberActive = Enable::ABSORBER || Enable::BEAMLINE_ABSORBER;
63 
64  int verbosity = std::max(Enable::VERBOSITY, Enable::BEAMLINE_VERBOSITY);
65 
66  G4BEAMLINE::ForwardBeamLineEnclosure = new PHG4CylinderSubsystem("ForwardBeamLineEnclosure");
67  G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("place_z", G4BEAMLINE::enclosure_center );
68  G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("radius", 0);
69  G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("thickness", G4BEAMLINE::enclosure_r_max);
70  G4BEAMLINE::ForwardBeamLineEnclosure->set_double_param("length", G4BEAMLINE::enclosure_z_max - G4BEAMLINE::starting_z);
71  G4BEAMLINE::ForwardBeamLineEnclosure->set_string_param("material", "G4_Galactic");
72  G4BEAMLINE::ForwardBeamLineEnclosure->set_color(.5, .5, .5, 0.2);
73  G4BEAMLINE::ForwardBeamLineEnclosure->OverlapCheck(overlapCheck);
74  if (verbosity) G4BEAMLINE::ForwardBeamLineEnclosure->Verbosity(verbosity);
75  g4Reco->registerSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
76 
77  G4BEAMLINE::BackwardBeamLineEnclosure = new PHG4CylinderSubsystem("BackwardBeamLineEnclosure");
78  G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("place_z", -G4BEAMLINE::enclosure_center);
79  G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("radius", 0);
80  G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("thickness", G4BEAMLINE::enclosure_r_max); // This is intentionally made large 25cm radius
81  G4BEAMLINE::BackwardBeamLineEnclosure->set_double_param("length", G4BEAMLINE::enclosure_z_max - G4BEAMLINE::starting_z);
82  G4BEAMLINE::BackwardBeamLineEnclosure->set_string_param("material", "G4_Galactic");
83  G4BEAMLINE::BackwardBeamLineEnclosure->set_color(.5, .5, .5, 0.2);
84  G4BEAMLINE::BackwardBeamLineEnclosure->OverlapCheck(overlapCheck);
85  if (verbosity) G4BEAMLINE::BackwardBeamLineEnclosure->Verbosity(verbosity);
86  g4Reco->registerSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
87 
88  string magFile;
89  magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/D0DXMagnets.dat";
90 
91  // if you insert numbers it only displays those magnets, do not comment out the set declaration
92  set<int> magnetlist;
93  //magnetlist.insert(7);
94 
95  BeamLineMagnetSubsystem *bl = nullptr;
96  std::ifstream infile(magFile);
97  if (infile.is_open())
98  {
99  double biggest_z = 0.;
100  int imagnet = 0;
101  std::string line;
102  while (std::getline(infile, line))
103  {
104  if (!line.compare(0, 1, "B") ||
105  !line.compare(0, 1, "Q") ||
106  !line.compare(0, 1, "S"))
107  {
108  std::istringstream iss(line);
109  string magname;
110  double x;
111  double y;
112  double z;
113  double inner_radius_zin;
114  double inner_radius_zout;
115  double outer_magnet_diameter;
116  double length;
117  double angle;
118  double dipole_field_x;
119  double fieldgradient;
120  if (!(iss >> magname >> x >> y >> z >> inner_radius_zin >> inner_radius_zout >> outer_magnet_diameter >> length >> angle >> dipole_field_x >> fieldgradient))
121  {
122  cout << "could not decode " << line << endl;
123  gSystem->Exit(1);
124  }
125  else
126  {
127  //------------------------
128 
129  string magtype;
130  if (inner_radius_zin != inner_radius_zout)
131  {
132  cout << "inner radius at front of magnet " << inner_radius_zin
133  << " not equal radius at back of magnet " << inner_radius_zout
134  << " needs change in code (replace tube by cone for beamline)" << endl;
135  gSystem->Exit(1);
136  }
137  if (verbosity > 0)
138  {
139  cout << endl
140  << endl
141  << "\tID number " << imagnet << endl;
142  cout << "magname: " << magname << endl;
143  cout << "x: " << x << endl;
144  cout << "y: " << y << endl;
145  cout << "z: " << z << endl;
146  cout << "inner_radius_zin: " << inner_radius_zin << endl;
147  cout << "inner_radius_zout: " << inner_radius_zout << endl;
148  cout << "outer_magnet_diameter: " << outer_magnet_diameter << endl;
149  cout << "length: " << length << endl;
150  cout << "angle: " << angle << endl;
151  cout << "dipole_field_x: " << dipole_field_x << endl;
152  cout << "fieldgradient: " << fieldgradient << endl;
153  }
154  if (!magname.compare(0, 1, "B"))
155  {
156  magtype = "DIPOLE";
157  }
158  else if (!magname.compare(0, 1, "Q"))
159  {
160  magtype = "QUADRUPOLE";
161  }
162  else if (!magname.compare(0, 1, "S"))
163  {
164  magtype = "SEXTUPOLE";
165  }
166  else
167  {
168  cout << "cannot decode magnet name " << magname << endl;
169  gSystem->Exit(1);
170  }
171  // convert to our units (cm, deg)
172  x *= 100.;
173  y *= 100.;
174  z *= 100.;
175  length *= 100.;
176  inner_radius_zin *= 100.;
177  outer_magnet_diameter *= 100.;
178  angle = (angle / M_PI * 180.) / 1000.; // given in mrad
179 
180  //------------------------
181 
182  if (magnetlist.empty() || magnetlist.find(imagnet) != magnetlist.end())
183  {
184  bl = new BeamLineMagnetSubsystem("BEAMLINEMAGNET", imagnet);
185  bl->set_double_param("field_y", MagFieldFlip(dipole_field_x));
186  bl->set_double_param("fieldgradient", MagFieldFlip(fieldgradient));
187  bl->set_string_param("magtype", magtype);
188  bl->set_double_param("length", length);
189  bl->set_double_param("place_x", PosFlip(x)); // relative position to mother vol.
190  bl->set_double_param("place_y", y); // relative position to mother vol.
191  if (z > 0)
192  {
193  bl->set_double_param("place_z", z - G4BEAMLINE::enclosure_center); // relative position to mother vol.
194  }
195  else
196  {
197  bl->set_double_param("place_z", z + G4BEAMLINE::enclosure_center); // relative position to mother vol.
198  }
199  bl->set_double_param("field_global_position_x", PosFlip(x)); // abs. position to world for field manager
200  bl->set_double_param("field_global_position_y", y); // abs. position to world for field manager
201  bl->set_double_param("field_global_position_z", z); // abs. position to world for field manager
202  bl->set_double_param("rot_y", AngleFlip(angle));
203  bl->set_double_param("field_global_rot_y", AngleFlip(angle)); // abs. rotation to world for field manager
204  bl->set_double_param("inner_radius", inner_radius_zin);
205  bl->set_double_param("outer_radius", outer_magnet_diameter / 2.);
206  bl->set_double_param("skin_thickness",G4BEAMLINE::skin_thickness);
207  bl->SetActive(AbsorberActive);
208  bl->SetAbsorberActive(AbsorberActive);
209  if (Enable::BEAMLINE_BLACKHOLE) bl->BlackHole(); // turn magnets into black holes
210  if (z > 0)
211  {
212  bl->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
213  }
214  else
215  {
216  bl->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
217  }
218  bl->OverlapCheck(overlapCheck);
219  bl->SuperDetector("BEAMLINEMAGNET");
220  if (verbosity) bl->Verbosity(verbosity);
221  g4Reco->registerSubsystem(bl);
222  }
223  imagnet++;
224  if (fabs(z) + length > biggest_z)
225  {
226  biggest_z = fabs(z) + length;
227  }
228  }
229  }
230  }
231  infile.close();
232  }
233 }
234 
235 void BeamLineDefineBeamPipe(PHG4Reco *g4Reco)
236 {
237  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::BEAMLINE_OVERLAPCHECK;
238  bool AbsorberActive = Enable::ABSORBER || Enable::BEAMLINE_ABSORBER;
239  int verbosity = std::max(Enable::VERBOSITY, Enable::BEAMLINE_VERBOSITY);
240 
241  const int ntube = 10;
242  const string nm[ntube] = {"B00", "B01", "B10", "B11", "B20", "B21", "B30", "B31", "B32", "B33"};
243  const double qlen[ntube] = {233.8, 233.8, 118.7, 118.7, 217.16, 217.16, 182.5, 182.5, 182.5, 182.5};
244  const double qir[ntube] = {6.08, 6.08, 14.60, 14.60, 20.0, 20.0, 6.07, 6.07, 6.07, 6.07};
245  const double qor[ntube] = {6.35, 6.35, 15.24, 15.24, 20.96, 20.96, 6.35, 6.35, 6.35, 6.35};
246  const double qrot[ntube] = {0, 0, 0, 0, 0, 0, 1.074, -1.074, -1.074, 1.074}; //degree
247  const double qxC[ntube] = {0, 0, 0, 0, 0, 0, 12.820, -12.820, 12.820, -12.820};
248  const double qyC[ntube] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
249  const double qzC[ntube] = {863.1, -863.1, 1474.470, -1474.470, 1642.4, -1642.4, 1843.2, 1843.2, -1843.2, -1843.2};
250  for (int i = 0; i < ntube; i++)
251  {
252  string name = "beamPipe" + nm[i];
253  PHG4CylinderSubsystem *pipe = new PHG4CylinderSubsystem(name, G4BEAMLINE::pipe_id_offset + i);
254  pipe->set_double_param("radius", qir[i]);
255  pipe->set_double_param("thickness", qor[i] - qir[i]);
256  pipe->set_double_param("length", qlen[i]);
257  pipe->set_double_param("rot_y", qrot[i]);
258  pipe->set_string_param("material", "G4_STAINLESS-STEEL");
259  pipe->set_double_param("place_x", PosFlip(qxC[i]));
260  pipe->set_double_param("place_y", qyC[i]);
261  if (qzC[i] > 0)
262  {
263  pipe->set_double_param("place_z", qzC[i] - G4BEAMLINE::enclosure_center); // relative position to mother vol.
264  }
265  else
266  {
267  pipe->set_double_param("place_z", qzC[i] + G4BEAMLINE::enclosure_center); // relative position to mother vol.
268  }
269  if (AbsorberActive) pipe->SetActive();
270  pipe->SuperDetector("PIPE");
271  if (qzC[i] > 0)
272  {
273  pipe->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
274  }
275  else
276  {
277  pipe->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
278  }
279 
280  pipe->OverlapCheck(OverlapCheck);
281  g4Reco->registerSubsystem(pipe);
282  }
283 
284  //Roman Pot pipe
285  const int nSec = 2;
286  const double len[nSec] = {20.87, 20.87};
287  const double ir1[nSec] = {7.14, 14.60};
288  const double or1[nSec] = {7.77, 15.24};
289  const double ir2[nSec] = {14.60, 7.14};
290  const double or2[nSec] = {15.24, 7.77};
291  const double xC[nSec] = {0, 0};
292  const double yC[nSec] = {0, 0};
293  const double zC[nSec] = {1394.25, -1394.25};
294  for (int i = 0; i < nSec; i++)
295  {
296  string name = "beamPipeRP" + to_string(i);
297  PHG4ConeSubsystem *pipe = new PHG4ConeSubsystem(name, G4BEAMLINE::roman_pot_pipe_id_offset + i);
298  pipe->set_string_param("material", "G4_STAINLESS-STEEL");
299  pipe->set_double_param("place_x", PosFlip(xC[i]));
300  pipe->set_double_param("place_y", yC[i]);
301  if (zC[i] > 0)
302  {
303  pipe->set_double_param("place_z", zC[i] - G4BEAMLINE::enclosure_center);
304  }
305  else
306  {
307  pipe->set_double_param("place_z", zC[i] + G4BEAMLINE::enclosure_center);
308  }
309  pipe->set_double_param("length", len[i]);
310  pipe->set_double_param("rmin1", ir1[i]);
311  pipe->set_double_param("rmin2", ir2[i]);
312  pipe->set_double_param("rmax1", or1[i]);
313  pipe->set_double_param("rmax2", or2[i]);
314  if (zC[i] > 0)
315  {
316  pipe->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
317  }
318  else
319  {
320  pipe->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
321  }
322  if (AbsorberActive) pipe->SetActive();
323  pipe->SuperDetector("PIPE");
324  pipe->OverlapCheck(OverlapCheck);
325  g4Reco->registerSubsystem(pipe);
326  }
327 }
328 
329 float PosFlip(float pos)
330 {
331  return pos;
332 };
333 float AngleFlip(float angle)
334 {
335  return angle;
336 };
337 float MagFieldFlip(float Bfield)
338 {
339  return Bfield;
340 };
341 
342 #endif