Analysis Software
Documentation for sPHENIX simulation 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::max_z + 2*no_overlapp;
41  double enclosure_z_max = 2050. + (800-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  if (Enable::BEAMLINE_BLACKHOLE) pipe->BlackHole();
255  pipe->set_color(1,0,0,1.);
256  pipe->set_double_param("radius", qir[i]);
257  pipe->set_double_param("thickness", qor[i] - qir[i]);
258  pipe->set_double_param("length", qlen[i]);
259  pipe->set_double_param("rot_y", qrot[i]);
260  pipe->set_string_param("material", "G4_STAINLESS-STEEL");
261  pipe->set_double_param("place_x", PosFlip(qxC[i]));
262  pipe->set_double_param("place_y", qyC[i]);
263  if (qzC[i] > 0)
264  {
265  pipe->set_double_param("place_z", qzC[i] - G4BEAMLINE::enclosure_center); // relative position to mother vol.
266  }
267  else
268  {
269  pipe->set_double_param("place_z", qzC[i] + G4BEAMLINE::enclosure_center); // relative position to mother vol.
270  }
271  if (AbsorberActive) pipe->SetActive();
272  pipe->SuperDetector("PIPE");
273  if (qzC[i] > 0)
274  {
275  pipe->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
276  }
277  else
278  {
279  pipe->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
280  }
281 
282  pipe->OverlapCheck(OverlapCheck);
283  g4Reco->registerSubsystem(pipe);
284  }
285 
286  //Roman Pot pipe
287  const int nSec = 2;
288  const double len[nSec] = {20.87, 20.87};
289  const double ir1[nSec] = {7.14, 14.60};
290  const double or1[nSec] = {7.77, 15.24};
291  const double ir2[nSec] = {14.60, 7.14};
292  const double or2[nSec] = {15.24, 7.77};
293  const double xC[nSec] = {0, 0};
294  const double yC[nSec] = {0, 0};
295  const double zC[nSec] = {1394.25, -1394.25};
296  for (int i = 0; i < nSec; i++)
297  {
298  string name = "beamPipeRP" + to_string(i);
299  PHG4ConeSubsystem *pipe = new PHG4ConeSubsystem(name, G4BEAMLINE::roman_pot_pipe_id_offset + i);
300  pipe->set_string_param("material", "G4_STAINLESS-STEEL");
301  pipe->set_double_param("place_x", PosFlip(xC[i]));
302  pipe->set_double_param("place_y", yC[i]);
303  if (zC[i] > 0)
304  {
305  pipe->set_double_param("place_z", zC[i] - G4BEAMLINE::enclosure_center);
306  }
307  else
308  {
309  pipe->set_double_param("place_z", zC[i] + G4BEAMLINE::enclosure_center);
310  }
311  pipe->set_double_param("length", len[i]);
312  pipe->set_double_param("rmin1", ir1[i]);
313  pipe->set_double_param("rmin2", ir2[i]);
314  pipe->set_double_param("rmax1", or1[i]);
315  pipe->set_double_param("rmax2", or2[i]);
316  if (zC[i] > 0)
317  {
318  pipe->SetMotherSubsystem(G4BEAMLINE::ForwardBeamLineEnclosure);
319  }
320  else
321  {
322  pipe->SetMotherSubsystem(G4BEAMLINE::BackwardBeamLineEnclosure);
323  }
324  if (AbsorberActive) pipe->SetActive();
325  if (Enable::BEAMLINE_BLACKHOLE) pipe->BlackHole();
326  pipe->set_color(1,0,0,1.);
327  pipe->SuperDetector("PIPE");
328  pipe->OverlapCheck(OverlapCheck);
329  g4Reco->registerSubsystem(pipe);
330  }
331 }
332 
333 float PosFlip(float pos)
334 {
335  return pos;
336 };
337 float AngleFlip(float angle)
338 {
339  return angle;
340 };
341 float MagFieldFlip(float Bfield)
342 {
343  return Bfield;
344 };
345 
346 #endif