Belle II Software  release-06-01-15
GeoFarBeamLineCreator.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <ir/geometry/GeoFarBeamLineCreator.h>
10 
11 #include <geometry/Materials.h>
12 #include <geometry/CreatorFactory.h>
13 #include <geometry/utilities.h>
14 #include <framework/gearbox/Unit.h>
15 #include <ir/simulation/SensitiveDetector.h>
16 #include <simulation/background/BkgSensitiveDetector.h>
17 
18 #include <cmath>
19 #include <boost/format.hpp>
20 #include <boost/foreach.hpp>
21 #include <boost/algorithm/string.hpp>
22 
23 #include <G4LogicalVolume.hh>
24 #include <G4PVPlacement.hh>
25 
26 //Shapes
27 #include <G4Box.hh>
28 #include <G4Tubs.hh>
29 #include <G4Torus.hh>
30 #include <G4Polycone.hh>
31 #include <G4Trd.hh>
32 #include <G4IntersectionSolid.hh>
33 #include <G4SubtractionSolid.hh>
34 #include <G4UserLimits.hh>
35 
36 using namespace std;
37 using namespace boost;
38 
39 namespace Belle2 {
46  using namespace geometry;
47 
48  namespace ir {
49 
50  //-----------------------------------------------------------------
51  // Register the Creator
52  //-----------------------------------------------------------------
53 
54  geometry::CreatorFactory<GeoFarBeamLineCreator> GeoFarBeamLineFactory("FarBeamLineCreator");
55 
56  //-----------------------------------------------------------------
57  // Implementation
58  //-----------------------------------------------------------------
59 
60  GeoFarBeamLineCreator::GeoFarBeamLineCreator()
61  {
62  m_sensitive = new SensitiveDetector();
63  }
64 
65  GeoFarBeamLineCreator::~GeoFarBeamLineCreator()
66  {
67  delete m_sensitive;
68  }
69 
70  void GeoFarBeamLineCreator::createGeometry(G4LogicalVolume& topVolume, GeometryTypes)
71  {
72 
73  const int N = 2;
74 
75  double stepMax = 5.0 * Unit::mm;
76  int flag_limitStep = int(m_config.getParameter("LimitStepLength"));
77 
78 
79  //double unitFactor = 10.0;
80  const double unitFactor = Unit::cm / Unit::mm;
81 
82  map<string, FarBeamLineElement> elements;
83 
84 
85  //--------------
86  //- limits
87 
88  //--------------
89  //- TubeR
90 
91  FarBeamLineElement tubeR;
92 
93  //get parameters from .xml file
94  std::string prep = "TubeR.";
95 
96  int TubeR_N = int(m_config.getParameter(prep + "N"));
97 
98  std::vector<double> TubeR_Z(TubeR_N);
99  std::vector<double> TubeR_R(TubeR_N);
100  std::vector<double> TubeR_r(TubeR_N);
101 
102  for (int i = 0; i < TubeR_N; ++i) {
103  ostringstream ossZID;
104  ossZID << "Z" << i;
105 
106  ostringstream ossRID;
107  ossRID << "R" << i;
108 
109  ostringstream ossrID;
110  ossrID << "r" << i;
111 
112  TubeR_Z[i] = m_config.getParameter(prep + ossZID.str()) * unitFactor;
113  TubeR_R[i] = m_config.getParameter(prep + ossRID.str()) * unitFactor;
114  TubeR_r[i] = m_config.getParameter(prep + ossrID.str(), 0.0) * unitFactor;
115  }
116 
117  tubeR.transform = G4Translate3D(0.0, 0.0, 0.0);
118 
119  //define geometry
120  tubeR.geo = new G4Polycone("geo_TubeR_name", 0.0, 2 * M_PI, TubeR_N, &(TubeR_Z[0]), &(TubeR_r[0]), &(TubeR_R[0]));
121 
122  tubeR.logi = NULL;
123 
124  elements["TubeR"] = tubeR;
125 
126  //--------------
127  //- TubeL
128 
129  FarBeamLineElement tubeL;
130 
131  //get parameters from .xml file
132  prep = "TubeL.";
133 
134  int TubeL_N = int(m_config.getParameter(prep + "N"));
135 
136  std::vector<double> TubeL_Z(TubeL_N);
137  std::vector<double> TubeL_R(TubeL_N);
138  std::vector<double> TubeL_r(TubeL_N);
139 
140  for (int i = 0; i < TubeL_N; ++i) {
141  ostringstream ossZID;
142  ossZID << "Z" << i;
143 
144  ostringstream ossRID;
145  ossRID << "R" << i;
146 
147  ostringstream ossrID;
148  ossrID << "r" << i;
149 
150  TubeL_Z[i] = m_config.getParameter(prep + ossZID.str()) * unitFactor;
151  TubeL_R[i] = m_config.getParameter(prep + ossRID.str()) * unitFactor;
152  TubeL_r[i] = m_config.getParameter(prep + ossrID.str(), 0.0) * unitFactor;
153  }
154 
155  tubeL.transform = G4Translate3D(0.0, 0.0, 0.0);
156 
157  //define geometry
158  tubeL.geo = new G4Polycone("geo_TubeL_name", 0.0, 2 * M_PI, TubeL_N, &(TubeL_Z[0]), &(TubeL_r[0]), &(TubeL_R[0]));
159 
160  tubeL.logi = NULL;
161 
162  elements["TubeL"] = tubeL;
163 
164  std::vector<double> zero_r(N, 0.);
165 
166  std::vector<std::string> straightSections;
167  boost::split(straightSections, m_config.getParameterStr("Straight"), boost::is_any_of(" "));
168  for (const auto& name : straightSections) {
169  //--------------
170  //- Create straight element
171 
172  prep = name + ".";
173  string type = m_config.getParameterStr(prep + "type");
174 
175 
176  FarBeamLineElement polycone;
177 
178  std::vector<double> Polycone_Z(N);
179  std::vector<double> Polycone_R(N);
180  std::vector<double> Polycone_r(N);
181  Polycone_Z[0] = 0;
182  Polycone_Z[1] = m_config.getParameter(prep + "L") * unitFactor;
183  Polycone_R[0] = m_config.getParameter(prep + "R") * unitFactor;
184  Polycone_R[1] = m_config.getParameter(prep + "R") * unitFactor;
185  Polycone_r[0] = m_config.getParameter(prep + "r") * unitFactor;
186  Polycone_r[1] = m_config.getParameter(prep + "r") * unitFactor;
187 
188  double Polycone_X0 = m_config.getParameter(prep + "X0") * unitFactor;
189  double Polycone_Z0 = m_config.getParameter(prep + "Z0") * unitFactor;
190  double Polycone_PHI = m_config.getParameter(prep + "PHI");
191 
192  polycone.transform = G4Translate3D(Polycone_X0, 0.0, Polycone_Z0);
193  polycone.transform = polycone.transform * G4RotateY3D(Polycone_PHI / Unit::rad);
194 
195  //define geometry
196  string subtract = m_config.getParameterStr(prep + "Subtract", "");
197  string intersect = m_config.getParameterStr(prep + "Intersect", "");
198 
199  string geo_polyconexx_name = "geo_" + name + "xx_name";
200  string geo_polyconex_name = "geo_" + name + "x_name";
201  string geo_polycone_name = "geo_" + name + "_name";
202 
203  G4VSolid* geo_polyconexx(NULL), *geo_polycone(NULL);
204 
205  if (subtract != "" || intersect != "")
206  if (type == "pipe") // for pipes inner space will be created as vacuum
207  geo_polyconexx = new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
208  else
209  geo_polyconexx = new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
210  else if (type == "pipe") // for pipes inner space will be created as vacuum
211  geo_polycone = new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
212  else
213  geo_polycone = new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
214 
215 
216  if (subtract != "" && intersect != "") {
217  G4VSolid* geo_polyconex = new G4SubtractionSolid(geo_polyconex_name, geo_polyconexx, elements[subtract].geo,
218  polycone.transform.inverse()*elements[subtract].transform);
219  geo_polycone = new G4IntersectionSolid(geo_polycone_name, geo_polyconex, elements[intersect].geo,
220  polycone.transform.inverse()*elements[intersect].transform);
221  } else if (subtract != "")
222  geo_polycone = new G4SubtractionSolid(geo_polycone_name, geo_polyconexx, elements[subtract].geo,
223  polycone.transform.inverse()*elements[subtract].transform);
224  else if (intersect != "")
225  geo_polycone = new G4IntersectionSolid(geo_polycone_name, geo_polyconexx, elements[intersect].geo,
226  polycone.transform.inverse()*elements[intersect].transform);
227 
228  polycone.geo = geo_polycone;
229 
230  // define logical volume
231  string strMat_polycone = m_config.getParameterStr(prep + "Material");
232  G4Material* mat_polycone = Materials::get(strMat_polycone);
233  string logi_polycone_name = "logi_" + name + "_name";
234  G4LogicalVolume* logi_polycone = new G4LogicalVolume(polycone.geo, mat_polycone, logi_polycone_name);
235  setColor(*logi_polycone, "#CC0000");
236  setVisibility(*logi_polycone, false);
237 
238  polycone.logi = logi_polycone;
239 
240  //put volume
241  string phys_polycone_name = "phys_" + name + "_name";
242  new G4PVPlacement(polycone.transform, logi_polycone, phys_polycone_name, &topVolume, false, 0);
243 
244  elements[name] = polycone;
245 
246  double sum = 0.0;
247  for (int i = 0; i < N; ++i)
248  sum += Polycone_r[i]; // check that there is a space inside a pipe
249  if (type == "pipe" && sum != 0) { // add vacuum inside a pipe
250 
251  FarBeamLineElement vacuum;
252 
253  string nameVac = name + "Vac";
254 
255  //define geometry
256  string geo_vacuumxx_name = "geo_" + nameVac + "xx_name";
257  string geo_vacuum_name = "geo_" + nameVac + "_name";
258 
259  G4VSolid* geo_vacuumxx, *geo_vacuum;
260 
261  geo_vacuumxx = new G4Polycone(geo_vacuumxx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_r[0]));
262  geo_vacuum = new G4IntersectionSolid(geo_vacuumxx_name, geo_vacuumxx, geo_polycone);
263 
264  vacuum.geo = geo_vacuum;
265  vacuum.transform = polycone.transform;
266 
267  // define logical volume
268  G4Material* mat_vacuum = Materials::get("Vacuum");
269  string logi_vacuum_name = "logi_" + nameVac + "_name";
270  G4LogicalVolume* logi_vacuum = new G4LogicalVolume(vacuum.geo, mat_vacuum, logi_vacuum_name);
271  if (flag_limitStep) logi_vacuum->SetUserLimits(new G4UserLimits(stepMax));
272  setColor(*logi_vacuum, "#000000");
273  setVisibility(*logi_vacuum, false);
274 
275  vacuum.logi = logi_vacuum;
276 
277  //put volume
278  string phys_vacuum_name = "phys_" + nameVac + "_name";
279  //new G4PVPlacement(vacuum.transform, logi_vacuum, phys_vacuum_name, &topVolume, false, 0);
280  new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_polycone, false, 0);
281 
282  elements[nameVac] = vacuum;
283  }
284  }
285 
286 
287  std::vector<std::string> bendingSections;
288  boost::split(bendingSections, m_config.getParameterStr("Bending"), boost::is_any_of(" "));
289  for (const auto& name : bendingSections) {
290  //--------------
291  //- Create torus element
292 
293  prep = name + ".";
294  string type = m_config.getParameterStr(prep + "type");
295 
296  FarBeamLineElement torus;
297 
298  double torus_r = m_config.getParameter(prep + "r") * unitFactor;
299  double torus_R = m_config.getParameter(prep + "R") * unitFactor;
300  double torus_RT = m_config.getParameter(prep + "RT") * unitFactor;
301  double torus_X0 = m_config.getParameter(prep + "X0") * unitFactor;
302  double torus_Z0 = m_config.getParameter(prep + "Z0") * unitFactor;
303  double torus_SPHI = m_config.getParameter(prep + "SPHI");
304  double torus_DPHI = m_config.getParameter(prep + "DPHI");
305 
306  torus.transform = G4Translate3D(torus_X0, 0.0, torus_Z0);
307  torus.transform = torus.transform * G4RotateX3D(M_PI / 2 / Unit::rad);
308 
309  //define geometry
310  string subtract = m_config.getParameterStr(prep + "Subtract", "");
311  string intersect = m_config.getParameterStr(prep + "Intersect", "");
312 
313  string geo_torusxx_name = "geo_" + name + "xx_name";
314  string geo_torusx_name = "geo_" + name + "x_name";
315  string geo_torus_name = "geo_" + name + "_name";
316 
317  G4VSolid* geo_torusxx(NULL), *geo_torus(NULL);
318 
319  if (subtract != "" || intersect != "")
320  if (type == "pipe") // for pipes inner space will be created as vacuum
321  geo_torusxx = new G4Torus(geo_torusxx_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
322  else
323  geo_torusxx = new G4Torus(geo_torusxx_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
324  else if (type == "pipe") // for pipes inner space will be created as vacuum
325  geo_torus = new G4Torus(geo_torus_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
326  else
327  geo_torus = new G4Torus(geo_torus_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
328 
329  if (subtract != "" && intersect != "") {
330  G4VSolid* geo_torusx = new G4SubtractionSolid(geo_torusx_name, geo_torusxx, elements[subtract].geo,
331  torus.transform.inverse()*elements[subtract].transform);
332  geo_torus = new G4IntersectionSolid(geo_torus_name, geo_torusx, elements[intersect].geo,
333  torus.transform.inverse()*elements[intersect].transform);
334  } else if (subtract != "")
335  geo_torus = new G4SubtractionSolid(geo_torus_name, geo_torusxx, elements[subtract].geo,
336  torus.transform.inverse()*elements[subtract].transform);
337  else if (intersect != "")
338  geo_torus = new G4IntersectionSolid(geo_torus_name, geo_torusxx, elements[intersect].geo,
339  torus.transform.inverse()*elements[intersect].transform);
340 
341  torus.geo = geo_torus;
342 
343  // define logical volume
344  string strMat_torus = m_config.getParameterStr(prep + "Material");
345  G4Material* mat_torus = Materials::get(strMat_torus);
346  string logi_torus_name = "logi_" + name + "_name";
347  G4LogicalVolume* logi_torus = new G4LogicalVolume(torus.geo, mat_torus, logi_torus_name);
348  setColor(*logi_torus, "#CC0000");
349  setVisibility(*logi_torus, false);
350 
351  torus.logi = logi_torus;
352 
353  //put volume
354  string phys_torus_name = "phys_" + name + "_name";
355  new G4PVPlacement(torus.transform, logi_torus, phys_torus_name, &topVolume, false, 0);
356 
357  elements[name] = torus;
358 
359  if (type == "pipe" && torus_r != 0) { // add vacuum inside a pipe
360 
361  FarBeamLineElement vacuum;
362 
363  string nameVac = name + "Vac";
364 
365  //define geometry
366  string geo_vacuumxx_name = "geo_" + nameVac + "xx_name";
367  string geo_vacuum_name = "geo_" + nameVac + "_name";
368 
369  G4VSolid* geo_vacuumxx, *geo_vacuum;
370 
371  geo_vacuumxx = new G4Torus(geo_vacuumxx_name, 0.0, torus_r, torus_RT, torus_SPHI, torus_DPHI);
372  geo_vacuum = new G4IntersectionSolid(geo_vacuum_name, geo_vacuumxx, geo_torus);
373 
374  vacuum.geo = geo_vacuum;
375  vacuum.transform = torus.transform;
376 
377  // define logical volume
378  G4Material* mat_vacuum = Materials::get("Vacuum");
379  string logi_vacuum_name = "logi_" + nameVac + "_name";
380  G4LogicalVolume* logi_vacuum = new G4LogicalVolume(vacuum.geo, mat_vacuum, logi_vacuum_name);
381  if (flag_limitStep) logi_vacuum->SetUserLimits(new G4UserLimits(stepMax));
382  setColor(*logi_vacuum, "#000000");
383  setVisibility(*logi_vacuum, false);
384 
385  vacuum.logi = logi_vacuum;
386 
387  //put volume
388  string phys_vacuum_name = "phys_" + nameVac + "_name";
389  //new G4PVPlacement(vacuum.transform, logi_vacuum, phys_vacuum_name, &topVolume, false, 0);
390  new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_torus, false, 0);
391 
392  elements[nameVac] = vacuum;
393  }
394  }
395 
396 
397  //--------------
398  //- concrete wall
399 
400  //--------------
401  //- GateShield (gate shield)
402 
403  //get parameters from .xml file
404  prep = "GateShield.";
405 
406  double GateShield_X = m_config.getParameter(prep + "X") * unitFactor;
407  double GateShield_Y = m_config.getParameter(prep + "Y") * unitFactor;
408  double GateShield_Z = m_config.getParameter(prep + "Z") * unitFactor;
409  double TUN_X = m_config.getParameter(prep + "TUNX") * unitFactor;
410  double TUN_Y = m_config.getParameter(prep + "TUNY") * unitFactor;
411  double DET_Z = m_config.getParameter(prep + "DETZ") * unitFactor;
412  double DET_DZ = m_config.getParameter(prep + "DETDZ") * unitFactor;
413  double ROT = m_config.getParameter(prep + "ROT");
414 
415  G4Transform3D transform_DET = G4Translate3D(0.0, 0.0, DET_DZ);
416  G4Transform3D transform_ROT = G4Translate3D(0.0, 0.0, 0.0);
417  transform_ROT = transform_ROT * G4RotateY3D(ROT / Unit::rad);
418 
419  //define geometry
420  // wall is made from the box by excluding spaces for detector and tunnel
421  G4Box* geo_GateShieldxx = new G4Box("geo_GateShieldxx_name", GateShield_X, GateShield_Y, GateShield_Z);
422  G4Box* geo_TUN = new G4Box("geo_TUN_name", TUN_X, TUN_Y, GateShield_Z);
423  G4SubtractionSolid* geo_GateShieldx = new G4SubtractionSolid("geo_GateShieldx_name", geo_GateShieldxx, geo_TUN);
424  G4Box* geo_DET = new G4Box("geo_DET_name", GateShield_X, GateShield_Y, DET_Z);
425  G4SubtractionSolid* geo_GateShield = new G4SubtractionSolid("geo_GateShield_name", geo_GateShieldx, geo_DET, transform_DET);
426 
427  string strMat_GateShield = m_config.getParameterStr(prep + "Material");
428  G4Material* mat_GateShield = Materials::get(strMat_GateShield);
429  G4LogicalVolume* logi_GateShield = new G4LogicalVolume(geo_GateShield, mat_GateShield, "logi_GateShield_name");
430 
431  //put volume
432  setColor(*logi_GateShield, "#CC0000");
433  //setVisibility(*logi_GateShield, false);
434  new G4PVPlacement(transform_ROT, logi_GateShield, "phys_GateShield_name", &topVolume, false, 0);
435 
436 
437 
438  //--------------
439  //- Tube (virtual tube for radiation level study)
440 
441  //define geometry
442  G4Tubs* geo_Tube = new G4Tubs("geo_Tube_name", 3995 * CLHEP::mm, 4000 * CLHEP::mm, 29 * CLHEP::m, 0. * CLHEP::deg, 360.*CLHEP::deg);
443  G4Material* mat_Tube = Materials::get("G4_Si");
444  G4LogicalVolume* logi_Tube = new G4LogicalVolume(geo_Tube, mat_Tube, "logi_Tube_name");
445 
446  //put volume
447  setColor(*logi_Tube, "#CC0000");
448  //setVisibility(*logi_Tube, false);
449  bool radiation_study = false;
450  // cppcheck-suppress knownConditionTrueFalse
451  if (radiation_study) {
452  new G4PVPlacement(transform_ROT, logi_Tube, "phys_Tube_name", &topVolume, false, 0);
453  }
454 
455 
456  //--------------
457  //- polyethylene shields
458 
459  //--------------
460  //- PolyShieldR
461 
462  //get parameters from .xml file
463  prep = "PolyShieldR.";
464 
465  double PolyShieldR_Xp = m_config.getParameter(prep + "Xp") * unitFactor;
466  double PolyShieldR_Xm = m_config.getParameter(prep + "Xm") * unitFactor;
467  double PolyShieldR_Y = m_config.getParameter(prep + "Y") * unitFactor;
468  double PolyShieldR_Z = m_config.getParameter(prep + "Z") * unitFactor;
469  double PolyShieldR_DZ = m_config.getParameter(prep + "DZ") * unitFactor;
470  double PolyShieldR_r = m_config.getParameter(prep + "r") * unitFactor;
471  double PolyShieldR_dx = m_config.getParameter(prep + "dx") * unitFactor;
472 
473  double PolyShieldR_X = (PolyShieldR_Xp + PolyShieldR_Xm) / 2;
474  G4Transform3D transform_polyShieldR = G4Translate3D((PolyShieldR_Xp - PolyShieldR_Xm) / 2, 0.0, PolyShieldR_DZ);
475  G4Transform3D transform_polyShieldR_Hole = G4Translate3D(PolyShieldR_dx, 0.0, 0.0);
476 
477  //define geometry
478  G4Box* geo_polyShieldRx = new G4Box("geo_polyShieldRx_name", PolyShieldR_X, PolyShieldR_Y, PolyShieldR_Z);
479  G4Tubs* geo_polyShieldR_Hole = new G4Tubs("geo_polyShieldRxx_name", 0 * CLHEP::mm, PolyShieldR_r, PolyShieldR_Z, 0. * CLHEP::deg,
480  360.*CLHEP::deg);
481  G4SubtractionSolid* geo_polyShieldR
482  = new G4SubtractionSolid("geo_polyShieldR_name", geo_polyShieldRx, geo_polyShieldR_Hole, transform_polyShieldR_Hole);
483 
484  string strMat_polyShieldR = m_config.getParameterStr(prep + "Material");
485  G4Material* mat_polyShieldR = Materials::get(strMat_polyShieldR);
486  G4LogicalVolume* logi_polyShieldR = new G4LogicalVolume(geo_polyShieldR, mat_polyShieldR, "logi_polyShieldR_name");
487 
488  //put volume
489  setColor(*logi_polyShieldR, "#0000CC");
490  //setVisibility(*logi_polyShieldL, false);
491  new G4PVPlacement(transform_polyShieldR, logi_polyShieldR, "phys_polyShieldR_name", &topVolume, false, 0);
492 
493  //--------------
494  //- PolyShieldL
495 
496  //get parameters from .xml file
497  prep = "PolyShieldL.";
498 
499  double PolyShieldL_Xp = m_config.getParameter(prep + "Xp") * unitFactor;
500  double PolyShieldL_Xm = m_config.getParameter(prep + "Xm") * unitFactor;
501  double PolyShieldL_Y = m_config.getParameter(prep + "Y") * unitFactor;
502  double PolyShieldL_Z = m_config.getParameter(prep + "Z") * unitFactor;
503  double PolyShieldL_DZ = m_config.getParameter(prep + "DZ") * unitFactor;
504  double PolyShieldL_r = m_config.getParameter(prep + "r") * unitFactor;
505  double PolyShieldL_dx = m_config.getParameter(prep + "dx") * unitFactor;
506 
507  double PolyShieldL_X = (PolyShieldL_Xp + PolyShieldL_Xm) / 2;
508  G4Transform3D transform_polyShieldL = G4Translate3D((PolyShieldL_Xp - PolyShieldL_Xm) / 2, 0.0, PolyShieldL_DZ);
509  G4Transform3D transform_polyShieldL_Hole = G4Translate3D(PolyShieldL_dx, 0.0, 0.0);
510 
511  //define geometry
512  G4Box* geo_polyShieldLxx = new G4Box("geo_polyShieldLxx_name", PolyShieldL_X, PolyShieldL_Y, PolyShieldL_Z);
513  G4Tubs* geo_polyShieldL_Hole = new G4Tubs("geo_polyShieldLxxx_name", 0 * CLHEP::mm, PolyShieldL_r, PolyShieldL_Z, 0. * CLHEP::deg,
514  360.*CLHEP::deg);
515  G4SubtractionSolid* geo_polyShieldLx
516  = new G4SubtractionSolid("geo_polyShieldLx_name", geo_polyShieldLxx, geo_polyShieldL_Hole, transform_polyShieldL_Hole);
517  G4SubtractionSolid* geo_polyShieldL
518  = new G4SubtractionSolid("geo_polyShieldL_name", geo_polyShieldLx, geo_GateShield, transform_polyShieldL.inverse()*transform_ROT);
519 
520  string strMat_polyShieldL = m_config.getParameterStr(prep + "Material");
521  G4Material* mat_polyShieldL = Materials::get(strMat_polyShieldL);
522  G4LogicalVolume* logi_polyShieldL = new G4LogicalVolume(geo_polyShieldL, mat_polyShieldL, "logi_polyShieldL_name");
523 
524  //put volume
525  setColor(*logi_polyShieldL, "#0000CC");
526  //setVisibility(*logi_polyShieldL, false);
527  new G4PVPlacement(transform_polyShieldL, logi_polyShieldL, "phys_polyShieldL_name", &topVolume, false, 0);
528 
529  //--------------
530  //- concrete tunnel-end shields
531 
532  //--------------
533  //- ConcreteShieldR
534 
535  //get parameters from .xml file
536  prep = "ConcreteShieldR.";
537 
538  double ConcreteShieldR_X = m_config.getParameter(prep + "X") * unitFactor;
539  double ConcreteShieldR_Y = m_config.getParameter(prep + "Y") * unitFactor;
540  double ConcreteShieldR_Z = m_config.getParameter(prep + "Z") * unitFactor;
541  double ConcreteShieldR_DZ = m_config.getParameter(prep + "DZ") * unitFactor;
542  double ConcreteShieldR_x = m_config.getParameter(prep + "x") * unitFactor;
543  double ConcreteShieldR_y = m_config.getParameter(prep + "y") * unitFactor;
544  double ConcreteShieldR_dx = m_config.getParameter(prep + "dx") * unitFactor;
545  double ConcreteShieldR_dy = m_config.getParameter(prep + "dy") * unitFactor;
546 
547  G4Transform3D transform_ConcreteShieldR = G4Translate3D(0.0, 0.0, ConcreteShieldR_DZ);
548  transform_ConcreteShieldR = transform_ROT * transform_ConcreteShieldR;
549  G4Transform3D transform_ConcreteShieldR_Hole = G4Translate3D(ConcreteShieldR_dx, ConcreteShieldR_dy, 0.0);
550 
551 
552  //define geometry
553  G4Box* geo_ConcreteShieldRx = new G4Box("geo_ConcreteShieldRx_name", ConcreteShieldR_X, ConcreteShieldR_Y, ConcreteShieldR_Z);
554  G4Box* geo_ConcreteShieldR_Hole = new G4Box("geo_ConcreteShieldRxx_name", ConcreteShieldR_x, ConcreteShieldR_y, ConcreteShieldR_Z);
555  G4SubtractionSolid* geo_ConcreteShieldR = new G4SubtractionSolid("geo_ConcreteShieldR_name", geo_ConcreteShieldRx,
556  geo_ConcreteShieldR_Hole, transform_ConcreteShieldR_Hole);
557 
558  string strMat_ConcreteShieldR = m_config.getParameterStr(prep + "Material");
559  G4Material* mat_ConcreteShieldR = Materials::get(strMat_ConcreteShieldR);
560  G4LogicalVolume* logi_ConcreteShieldR = new G4LogicalVolume(geo_ConcreteShieldR, mat_ConcreteShieldR, "logi_ConcreteShieldR_name");
561 
562  //put volume
563  setColor(*logi_ConcreteShieldR, "#0000CC");
564  //setVisibility(*logi_ConcreteShieldR, false);
565  new G4PVPlacement(transform_ConcreteShieldR, logi_ConcreteShieldR, "phys_ConcreteShieldR_name", &topVolume, false, 0);
566 
567  //--------------
568  //- ConcreteShieldL
569 
570  //get parameters from .xml file
571  prep = "ConcreteShieldL.";
572 
573  double ConcreteShieldL_X = m_config.getParameter(prep + "X") * unitFactor;
574  double ConcreteShieldL_Y = m_config.getParameter(prep + "Y") * unitFactor;
575  double ConcreteShieldL_Z = m_config.getParameter(prep + "Z") * unitFactor;
576  double ConcreteShieldL_DZ = m_config.getParameter(prep + "DZ") * unitFactor;
577  double ConcreteShieldL_x = m_config.getParameter(prep + "x") * unitFactor;
578  double ConcreteShieldL_y = m_config.getParameter(prep + "y") * unitFactor;
579  double ConcreteShieldL_dx = m_config.getParameter(prep + "dx") * unitFactor;
580  double ConcreteShieldL_dy = m_config.getParameter(prep + "dy") * unitFactor;
581 
582  G4Transform3D transform_ConcreteShieldL = G4Translate3D(0.0, 0.0, ConcreteShieldL_DZ);
583  transform_ConcreteShieldL = transform_ROT * transform_ConcreteShieldL;
584  G4Transform3D transform_ConcreteShieldL_Hole = G4Translate3D(ConcreteShieldL_dx, ConcreteShieldL_dy, 0.0);
585 
586  //define geometry
587  G4Box* geo_ConcreteShieldLx = new G4Box("geo_ConcreteShieldLx_name", ConcreteShieldL_X, ConcreteShieldL_Y, ConcreteShieldL_Z);
588  G4Box* geo_ConcreteShieldL_Hole = new G4Box("geo_ConcreteShieldLxx_name", ConcreteShieldL_x, ConcreteShieldL_y, ConcreteShieldL_Z);
589  G4SubtractionSolid* geo_ConcreteShieldL = new G4SubtractionSolid("geo_ConcreteShieldL_name", geo_ConcreteShieldLx,
590  geo_ConcreteShieldL_Hole, transform_ConcreteShieldL_Hole);
591 
592  string strMat_ConcreteShieldL = m_config.getParameterStr(prep + "Material");
593  G4Material* mat_ConcreteShieldL = Materials::get(strMat_ConcreteShieldL);
594  G4LogicalVolume* logi_ConcreteShieldL = new G4LogicalVolume(geo_ConcreteShieldL, mat_ConcreteShieldL, "logi_ConcreteShieldL_name");
595 
596  //put volume
597  setColor(*logi_ConcreteShieldL, "#0000CC");
598  //setVisibility(*logi_ConcreteShieldL, false);
599  new G4PVPlacement(transform_ConcreteShieldL, logi_ConcreteShieldL, "phys_ConcreteShieldL_name", &topVolume, false, 0);
600 
601 
602  //---------------------------
603  // for dose simulation
604  //---------------------------
605 
606  //neutron shield (poly)
607  //logi_polyShieldL->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1001));
608  //logi_polyShieldR->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1002));
609 
610  //additional neutron shield (concrete)
611  //logi_ConcreteShieldL->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1003));
612  //logi_ConcreteShieldR->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1004));
613 
614  //gate shield (concrete)
615  //logi_GateShield->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1005));
616 
617  //virtual material outsire gate-shield
618 
619  // cppcheck-suppress knownConditionTrueFalse
620  if (radiation_study) {
621  logi_Tube->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1006));
622  }
623 
624 
625  //------------------
626  //- Collimators
627 
628  std::vector<std::string> collimators;
629  boost::split(collimators, m_config.getParameterStr("Collimator"), boost::is_any_of(" "));
630  for (const auto& name : collimators) {
631  //- Collimators consist of two independent jaws (trapezoids), identical in shape, positioned opposite to each other
632  //- Each jaw consists of copper body and high Z head
633 
634  prep = name + ".";
635 
636  string type = m_config.getParameterStr(prep + "type");
637  string motherVolume = m_config.getParameterStr(prep + "MotherVolume");
638  string motherVolumeVacuum = motherVolume + "Vac";
639 
640  // If zz < 0 (positioned at negative z) vertical collimator is flipped when rotated into Mother Volume system
641  G4Scale3D scale;
642  G4Rotate3D rotation;
643  G4Translate3D translation;
644  elements[motherVolumeVacuum].transform.getDecomposition(scale, rotation, translation);
645  double zz = rotation.zz();
646 
647  // d1, d2 are collimator jaws displacements from beam center, d1<0, d2>0
648  // Z is collimator position inside its Mother Volume
649  double collimator_d1 = m_config.getParameter(prep + "d1") * unitFactor;
650  double collimator_d2 = m_config.getParameter(prep + "d2") * unitFactor;
651  double collimator_fullH = m_config.getParameter(prep + "fullH") * unitFactor;
652  double collimator_headH = m_config.getParameter(prep + "headH") * unitFactor;
653  double collimator_minW = m_config.getParameter(prep + "minW") * unitFactor;
654  double collimator_maxW = m_config.getParameter(prep + "maxW") * unitFactor;
655  double collimator_th = m_config.getParameter(prep + "th") * unitFactor;
656  double collimator_Z = m_config.getParameter(prep + "Z") * unitFactor;
657 
658  B2WARNING("Collimator " << name << " displacement d1 is set to " << collimator_d1 << "mm (must be less than zero)");
659  B2WARNING("Collimator " << name << " displacement d2 is set to " << collimator_d2 << "mm (must be greater than zero)");
660 
661 
662  // Collimator heads
663 
664  // dx1,2 dy1,2 dz are trapezoid dimensions
665  double head_dx1;
666  double head_dx2;
667  double head_dy1;
668  double head_dy2;
669  double head_dz = collimator_headH / 2.0;
670  if (type == "vertical") {
671  head_dx1 = collimator_th / 2.0;
672  head_dx2 = collimator_th / 2.0;
673  head_dy1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
674  head_dy2 = collimator_minW / 2.0;
675  } else {
676  head_dx1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
677  head_dx2 = collimator_minW / 2.0;
678  head_dy1 = collimator_th / 2.0;
679  head_dy2 = collimator_th / 2.0;
680  }
681 
682  // storable elements
683  FarBeamLineElement collimator_head1;
684  FarBeamLineElement collimator_head2;
685 
686  // move collimator to position on beam line
687  G4Transform3D transform_head1 = G4Translate3D(0.0, 0.0, collimator_Z);
688  G4Transform3D transform_head2 = G4Translate3D(0.0, 0.0, collimator_Z);
689 
690  // rotate and move collimator jaws to their relative positions
691  if (type == "vertical") {
692  transform_head1 = transform_head1 * G4Translate3D(0.0, -head_dz + collimator_d1, 0.0);
693  transform_head1 = transform_head1 * G4RotateX3D(-M_PI / 2 / Unit::rad);
694 
695  transform_head2 = transform_head2 * G4Translate3D(0.0, head_dz + collimator_d2, 0.0);
696  transform_head2 = transform_head2 * G4RotateX3D(M_PI / 2 / Unit::rad);
697  } else {
698  if (zz > 0) {
699  transform_head1 = transform_head1 * G4Translate3D(-head_dz + collimator_d1, 0.0, 0.0);
700  transform_head1 = transform_head1 * G4RotateY3D(M_PI / 2 / Unit::rad);
701 
702  transform_head2 = transform_head2 * G4Translate3D(head_dz + collimator_d2, 0.0, 0.0);
703  transform_head2 = transform_head2 * G4RotateY3D(-M_PI / 2 / Unit::rad);
704  } else {
705  transform_head1 = transform_head1 * G4Translate3D(head_dz - collimator_d1, 0.0, 0.0);
706  transform_head1 = transform_head1 * G4RotateY3D(-M_PI / 2 / Unit::rad);
707 
708  transform_head2 = transform_head2 * G4Translate3D(-head_dz - collimator_d2, 0.0, 0.0);
709  transform_head2 = transform_head2 * G4RotateY3D(M_PI / 2 / Unit::rad);
710  }
711  }
712 
713  collimator_head1.transform = transform_head1;
714  collimator_head2.transform = transform_head2;
715 
716  // define geometry
717  string geo_headx_name = "geo_" + name + "_headx_name";
718 
719  string geo_head1_name = "geo_" + name + "_head1_name";
720  string geo_head2_name = "geo_" + name + "_head2_name";
721 
722  G4VSolid* geo_headx = new G4Trd(geo_headx_name, head_dx1, head_dx2, head_dy1, head_dy2, head_dz);
723 
724  G4VSolid* geo_head1 = new G4IntersectionSolid(geo_head1_name, geo_headx, elements[motherVolumeVacuum].geo,
725  collimator_head1.transform.inverse());
726  G4VSolid* geo_head2 = new G4IntersectionSolid(geo_head2_name, geo_headx, elements[motherVolumeVacuum].geo,
727  collimator_head2.transform.inverse());
728 
729  collimator_head1.geo = geo_head1;
730  collimator_head2.geo = geo_head2;
731 
732  // define logical volume
733  string strMat_head = m_config.getParameterStr(prep + "HeadMaterial");
734  G4Material* mat_head = Materials::get(strMat_head);
735  string logi_head1_name = "logi_" + name + "_head1_name";
736  string logi_head2_name = "logi_" + name + "_head2_name";
737  G4LogicalVolume* logi_head1 = new G4LogicalVolume(geo_head1, mat_head, logi_head1_name);
738  G4LogicalVolume* logi_head2 = new G4LogicalVolume(geo_head2, mat_head, logi_head2_name);
739  setColor(*logi_head1, "#CC0000");
740  setColor(*logi_head2, "#CC0000");
741  setVisibility(*logi_head1, false);
742  setVisibility(*logi_head2, false);
743 
744  // check if collimator is inside beam pipe
745  double volume_head1 = logi_head1->GetSolid()->GetCubicVolume();
746  double volume_head2 = logi_head2->GetSolid()->GetCubicVolume();
747 
748  collimator_head1.logi = logi_head1;
749  collimator_head2.logi = logi_head2;
750 
751  // put volume
752  string phys_head1_name = "phys_" + name + "_head1" + "_name";
753  string phys_head2_name = "phys_" + name + "_head2" + "_name";
754  if (volume_head1 != 0)
755  new G4PVPlacement(collimator_head1.transform, logi_head1, phys_head1_name, elements[motherVolumeVacuum].logi, false, 0);
756  if (volume_head2 != 0)
757  new G4PVPlacement(collimator_head2.transform, logi_head2, phys_head2_name, elements[motherVolumeVacuum].logi, false, 0);
758 
759  // to use it later in "intersect" and "subtract"
760  collimator_head1.transform = collimator_head1.transform * elements[motherVolumeVacuum].transform;
761  collimator_head2.transform = collimator_head2.transform * elements[motherVolumeVacuum].transform;
762 
763  string name_head1 = name + "_head1";
764  string name_head2 = name + "_head2";
765  elements[name_head1] = collimator_head1;
766  elements[name_head2] = collimator_head2;
767 
768 
769  // Collimator bodies
770 
771  // dx1,2 dy1,2 dz are trapezoid dimensions
772  double body_dx1;
773  double body_dx2;
774  double body_dy1;
775  double body_dy2;
776  double body_dz = (collimator_fullH - collimator_headH) / 2.0;
777  if (type == "vertical") {
778  body_dx1 = collimator_th / 2.0;
779  body_dx2 = collimator_th / 2.0;
780  body_dy1 = collimator_maxW / 2.0;
781  body_dy2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
782  } else {
783  body_dx1 = collimator_maxW / 2.0;
784  body_dx2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
785  body_dy1 = collimator_th / 2.0;
786  body_dy2 = collimator_th / 2.0;
787  }
788 
789  // storable elements
790  FarBeamLineElement collimator_body1;
791  FarBeamLineElement collimator_body2;
792 
793  // reuse head transfomation with additional shift
794  if (type == "vertical") {
795  collimator_body1.transform = G4Translate3D(0.0, -head_dz - body_dz, 0.0) * transform_head1;
796  collimator_body2.transform = G4Translate3D(0.0, head_dz + body_dz, 0.0) * transform_head2;
797  } else {
798  if (zz > 0) {
799  collimator_body1.transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head1;
800  collimator_body2.transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head2;
801  } else {
802  collimator_body1.transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head1;
803  collimator_body2.transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head2;
804  }
805  }
806 
807  // define geometry
808  string geo_bodyx_name = "geo_" + name + "_bodyx_name";
809 
810  string geo_body1_name = "geo_" + name + "_body1_name";
811  string geo_body2_name = "geo_" + name + "_body2_name";
812 
813  G4VSolid* geo_bodyx = new G4Trd(geo_bodyx_name, body_dx1, body_dx2, body_dy1, body_dy2, body_dz);
814 
815  G4VSolid* geo_body1 = new G4IntersectionSolid(geo_body1_name, geo_bodyx, elements[motherVolumeVacuum].geo,
816  collimator_body1.transform.inverse());
817  G4VSolid* geo_body2 = new G4IntersectionSolid(geo_body2_name, geo_bodyx, elements[motherVolumeVacuum].geo,
818  collimator_body2.transform.inverse());
819 
820  collimator_body1.geo = geo_body1;
821  collimator_body2.geo = geo_body2;
822 
823  // define logical volume
824  string strMat_body = m_config.getParameterStr(prep + "Material");
825  G4Material* mat_body = Materials::get(strMat_body);
826  string logi_body1_name = "logi_" + name + "_body1_name";
827  string logi_body2_name = "logi_" + name + "_body2_name";
828  G4LogicalVolume* logi_body1 = new G4LogicalVolume(geo_body1, mat_body, logi_body1_name);
829  G4LogicalVolume* logi_body2 = new G4LogicalVolume(geo_body2, mat_body, logi_body2_name);
830  setColor(*logi_body1, "#CC0000");
831  setColor(*logi_body2, "#CC0000");
832  setVisibility(*logi_body1, false);
833  setVisibility(*logi_body2, false);
834 
835  // check if collimator is inside beam pipe
836  double volume_body1 = logi_body1->GetSolid()->GetCubicVolume();
837  double volume_body2 = logi_body2->GetSolid()->GetCubicVolume();
838 
839  collimator_body1.logi = logi_body1;
840  collimator_body2.logi = logi_body2;
841 
842  // put volume
843  string phys_body1_name = "phys_" + name + "_body1" + "_name";
844  string phys_body2_name = "phys_" + name + "_body2" + "_name";
845  if (volume_body1 != 0)
846  new G4PVPlacement(collimator_body1.transform, logi_body1, phys_body1_name, elements[motherVolumeVacuum].logi, false, 0);
847  if (volume_body2 != 0)
848  new G4PVPlacement(collimator_body2.transform, logi_body2, phys_body2_name, elements[motherVolumeVacuum].logi, false, 0);
849 
850  // to use it later in "intersect" and "subtract"
851  collimator_body1.transform = collimator_body1.transform * elements[motherVolumeVacuum].transform;
852  collimator_body2.transform = collimator_body2.transform * elements[motherVolumeVacuum].transform;
853 
854  string name_body1 = name + "_body1";
855  string name_body2 = name + "_body2";
856  elements[name_body1] = collimator_body1;
857  elements[name_body2] = collimator_body2;
858  }
859 
860 
861  //-----------------------
862  //- Collimator shields
863 
864  std::vector<std::string> collimatorShields;
865  boost::split(collimatorShields, m_config.getParameterStr("CollimatorShield"), boost::is_any_of(" "));
866  for (const auto& name : collimatorShields) {
867  //- Collimator shield made as box with subtracted inner space
868 
869  prep = name + ".";
870 
871  double collShield_X0 = m_config.getParameter(prep + "X0") * unitFactor;
872  double collShield_Z0 = m_config.getParameter(prep + "Z0") * unitFactor;
873  double collShield_PHI = m_config.getParameter(prep + "PHI");
874  double collShield_W = m_config.getParameter(prep + "W") * unitFactor;
875  double collShield_H = m_config.getParameter(prep + "H") * unitFactor;
876  double collShield_L = m_config.getParameter(prep + "L") * unitFactor;
877  double collShield_d = m_config.getParameter(prep + "d") * unitFactor;
878 
879  // storable element
880  FarBeamLineElement collShield;
881 
882  // move collimator to its position on beam line
883  collShield.transform = G4Translate3D(collShield_X0, 0.0, collShield_Z0);
884  collShield.transform = collShield.transform * G4RotateY3D(collShield_PHI / Unit::rad);
885 
886  G4Transform3D transform_collShieldInner = G4Translate3D(0.0, -collShield_d / 2.0, 0.0);
887 
888  // define geometry
889  string geo_collShieldx_name = "geo_" + name + "x_name";
890  string geo_collShieldInner_name = "geo_" + name + "Inner" + "_name";
891  string geo_collShield_name = "geo_" + name + "_name";
892 
893  G4Box* geo_collShieldx = new G4Box(geo_collShieldx_name, collShield_W / 2.0, collShield_H / 2.0, collShield_L / 2.0);
894  G4Box* geo_collShieldInner = new G4Box(geo_collShieldInner_name, collShield_W / 2.0 - collShield_d,
895  collShield_H / 2.0 - collShield_d / 2.0, collShield_L / 2.0);
896  G4SubtractionSolid* geo_collShield = new G4SubtractionSolid(geo_collShield_name, geo_collShieldx, geo_collShieldInner,
897  transform_collShieldInner);
898 
899  collShield.geo = geo_collShield;
900 
901  // define logical volume
902  string strMat_collShield = m_config.getParameterStr(prep + "Material");
903  G4Material* mat_collShield = Materials::get(strMat_collShield);
904  string logi_collShield_name = "logi_" + name + "_name";
905  G4LogicalVolume* logi_collShield = new G4LogicalVolume(geo_collShield, mat_collShield, logi_collShield_name);
906  setColor(*logi_collShield, "#CC0000");
907  setVisibility(*logi_collShield, false);
908 
909  collShield.logi = logi_collShield;
910 
911  // put volume
912  string phys_collShield_name = "phys_" + name + "_name";
913  new G4PVPlacement(collShield.transform, logi_collShield, phys_collShield_name, &topVolume, false, 0);
914 
915  elements[name] = collShield;
916  }
917  }
918  }
920 }
The Class for BeamBackground Sensitive Detector.
The IR Sensitive Detector class.
int intersect(const TRGCDCLpar &lp1, const TRGCDCLpar &lp2, CLHEP::HepVector &v1, CLHEP::HepVector &v2)
intersection
Definition: Lpar.cc:249
void setVisibility(G4LogicalVolume &volume, bool visible)
Helper function to quickly set the visibility of a given volume.
Definition: utilities.cc:105
void setColor(G4LogicalVolume &volume, const std::string &color)
Set the color of a logical volume.
Definition: utilities.cc:97
GeometryTypes
Flag indiciating the type of geometry to be used.
Abstract base class for different kinds of events.
The struct for FarBeamLineElement.
G4LogicalVolume * logi
Logical volume.
G4Transform3D transform
Transformation.