Belle II Software  release-08-01-10
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_Y0 = m_config.getParameter(prep + "Y0", 0) * unitFactor;
190  double Polycone_Z0 = m_config.getParameter(prep + "Z0") * unitFactor;
191  double Polycone_PHI = m_config.getParameter(prep + "PHI");
192  double Polycone_PHIYZ = m_config.getParameter(prep + "PHIYZ", 0);
193 
194  polycone.transform = G4Translate3D(Polycone_X0, Polycone_Y0, Polycone_Z0);
195  polycone.transform = polycone.transform * G4RotateY3D(Polycone_PHI / Unit::rad);
196  if (Polycone_PHIYZ != 0)
197  polycone.transform = polycone.transform * G4RotateX3D(Polycone_PHIYZ / Unit::rad);
198 
199  //define geometry
200  string subtract = m_config.getParameterStr(prep + "Subtract", "");
201  string intersect = m_config.getParameterStr(prep + "Intersect", "");
202 
203  string geo_polyconexx_name = "geo_" + name + "xx_name";
204  string geo_polyconex_name = "geo_" + name + "x_name";
205  string geo_polycone_name = "geo_" + name + "_name";
206 
207  G4VSolid* geo_polyconexx(NULL), *geo_polycone(NULL);
208 
209  if (subtract != "" || intersect != "")
210  if (type == "pipe") // for pipes inner space will be created as vacuum
211  geo_polyconexx = new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
212  else
213  geo_polyconexx = new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
214  else if (type == "pipe") // for pipes inner space will be created as vacuum
215  geo_polycone = new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
216  else
217  geo_polycone = new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
218 
219 
220  if (subtract != "" && intersect != "") {
221  G4VSolid* geo_polyconex = new G4SubtractionSolid(geo_polyconex_name, geo_polyconexx, elements[subtract].geo,
222  polycone.transform.inverse()*elements[subtract].transform);
223  geo_polycone = new G4IntersectionSolid(geo_polycone_name, geo_polyconex, elements[intersect].geo,
224  polycone.transform.inverse()*elements[intersect].transform);
225  } else if (subtract != "")
226  geo_polycone = new G4SubtractionSolid(geo_polycone_name, geo_polyconexx, elements[subtract].geo,
227  polycone.transform.inverse()*elements[subtract].transform);
228  else if (intersect != "")
229  geo_polycone = new G4IntersectionSolid(geo_polycone_name, geo_polyconexx, elements[intersect].geo,
230  polycone.transform.inverse()*elements[intersect].transform);
231 
232  polycone.geo = geo_polycone;
233 
234  // define logical volume
235  string strMat_polycone = m_config.getParameterStr(prep + "Material");
236  G4Material* mat_polycone = Materials::get(strMat_polycone);
237  string logi_polycone_name = "logi_" + name + "_name";
238  G4LogicalVolume* logi_polycone = new G4LogicalVolume(polycone.geo, mat_polycone, logi_polycone_name);
239  setColor(*logi_polycone, "#CC0000");
240  setVisibility(*logi_polycone, false);
241 
242  polycone.logi = logi_polycone;
243 
244  //put volume
245  string phys_polycone_name = "phys_" + name + "_name";
246  new G4PVPlacement(polycone.transform, logi_polycone, phys_polycone_name, &topVolume, false, 0);
247 
248  elements[name] = polycone;
249 
250  double sum = 0.0;
251  for (int i = 0; i < N; ++i)
252  sum += Polycone_r[i]; // check that there is a space inside a pipe
253  if (type == "pipe" && sum != 0) { // add vacuum inside a pipe
254 
255  FarBeamLineElement vacuum;
256 
257  string nameVac = name + "Vac";
258 
259  //define geometry
260  string geo_vacuumxx_name = "geo_" + nameVac + "xx_name";
261  string geo_vacuum_name = "geo_" + nameVac + "_name";
262 
263  G4VSolid* geo_vacuumxx, *geo_vacuum;
264 
265  geo_vacuumxx = new G4Polycone(geo_vacuumxx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_r[0]));
266  geo_vacuum = new G4IntersectionSolid(geo_vacuumxx_name, geo_vacuumxx, geo_polycone);
267 
268  vacuum.geo = geo_vacuum;
269  vacuum.transform = polycone.transform;
270 
271  // define logical volume
272  G4Material* mat_vacuum = Materials::get("Vacuum");
273  string logi_vacuum_name = "logi_" + nameVac + "_name";
274  G4LogicalVolume* logi_vacuum = new G4LogicalVolume(vacuum.geo, mat_vacuum, logi_vacuum_name);
275  if (flag_limitStep) logi_vacuum->SetUserLimits(new G4UserLimits(stepMax));
276  setColor(*logi_vacuum, "#000000");
277  setVisibility(*logi_vacuum, false);
278 
279  vacuum.logi = logi_vacuum;
280 
281  //put volume
282  string phys_vacuum_name = "phys_" + nameVac + "_name";
283  //new G4PVPlacement(vacuum.transform, logi_vacuum, phys_vacuum_name, &topVolume, false, 0);
284  new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_polycone, false, 0);
285 
286  elements[nameVac] = vacuum;
287  }
288  }
289 
290 
291  std::vector<std::string> bendingSections;
292  boost::split(bendingSections, m_config.getParameterStr("Bending"), boost::is_any_of(" "));
293  for (const auto& name : bendingSections) {
294  //--------------
295  //- Create torus element
296 
297  prep = name + ".";
298  string type = m_config.getParameterStr(prep + "type");
299 
300  FarBeamLineElement torus;
301 
302  double torus_r = m_config.getParameter(prep + "r") * unitFactor;
303  double torus_R = m_config.getParameter(prep + "R") * unitFactor;
304  double torus_RT = m_config.getParameter(prep + "RT") * unitFactor;
305  double torus_X0 = m_config.getParameter(prep + "X0") * unitFactor;
306  double torus_Z0 = m_config.getParameter(prep + "Z0") * unitFactor;
307  double torus_SPHI = m_config.getParameter(prep + "SPHI");
308  double torus_DPHI = m_config.getParameter(prep + "DPHI");
309 
310  torus.transform = G4Translate3D(torus_X0, 0.0, torus_Z0);
311  torus.transform = torus.transform * G4RotateX3D(M_PI / 2 / Unit::rad);
312 
313  //define geometry
314  string subtract = m_config.getParameterStr(prep + "Subtract", "");
315  string intersect = m_config.getParameterStr(prep + "Intersect", "");
316 
317  string geo_torusxx_name = "geo_" + name + "xx_name";
318  string geo_torusx_name = "geo_" + name + "x_name";
319  string geo_torus_name = "geo_" + name + "_name";
320 
321  G4VSolid* geo_torusxx(NULL), *geo_torus(NULL);
322 
323  if (subtract != "" || intersect != "")
324  if (type == "pipe") // for pipes inner space will be created as vacuum
325  geo_torusxx = new G4Torus(geo_torusxx_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
326  else
327  geo_torusxx = new G4Torus(geo_torusxx_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
328  else if (type == "pipe") // for pipes inner space will be created as vacuum
329  geo_torus = new G4Torus(geo_torus_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
330  else
331  geo_torus = new G4Torus(geo_torus_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
332 
333  if (subtract != "" && intersect != "") {
334  G4VSolid* geo_torusx = new G4SubtractionSolid(geo_torusx_name, geo_torusxx, elements[subtract].geo,
335  torus.transform.inverse()*elements[subtract].transform);
336  geo_torus = new G4IntersectionSolid(geo_torus_name, geo_torusx, elements[intersect].geo,
337  torus.transform.inverse()*elements[intersect].transform);
338  } else if (subtract != "")
339  geo_torus = new G4SubtractionSolid(geo_torus_name, geo_torusxx, elements[subtract].geo,
340  torus.transform.inverse()*elements[subtract].transform);
341  else if (intersect != "")
342  geo_torus = new G4IntersectionSolid(geo_torus_name, geo_torusxx, elements[intersect].geo,
343  torus.transform.inverse()*elements[intersect].transform);
344 
345  torus.geo = geo_torus;
346 
347  // define logical volume
348  string strMat_torus = m_config.getParameterStr(prep + "Material");
349  G4Material* mat_torus = Materials::get(strMat_torus);
350  string logi_torus_name = "logi_" + name + "_name";
351  G4LogicalVolume* logi_torus = new G4LogicalVolume(torus.geo, mat_torus, logi_torus_name);
352  setColor(*logi_torus, "#CC0000");
353  setVisibility(*logi_torus, false);
354 
355  torus.logi = logi_torus;
356 
357  //put volume
358  string phys_torus_name = "phys_" + name + "_name";
359  new G4PVPlacement(torus.transform, logi_torus, phys_torus_name, &topVolume, false, 0);
360 
361  elements[name] = torus;
362 
363  if (type == "pipe" && torus_r != 0) { // add vacuum inside a pipe
364 
365  FarBeamLineElement vacuum;
366 
367  string nameVac = name + "Vac";
368 
369  //define geometry
370  string geo_vacuumxx_name = "geo_" + nameVac + "xx_name";
371  string geo_vacuum_name = "geo_" + nameVac + "_name";
372 
373  G4VSolid* geo_vacuumxx, *geo_vacuum;
374 
375  geo_vacuumxx = new G4Torus(geo_vacuumxx_name, 0.0, torus_r, torus_RT, torus_SPHI, torus_DPHI);
376  geo_vacuum = new G4IntersectionSolid(geo_vacuum_name, geo_vacuumxx, geo_torus);
377 
378  vacuum.geo = geo_vacuum;
379  vacuum.transform = torus.transform;
380 
381  // define logical volume
382  G4Material* mat_vacuum = Materials::get("Vacuum");
383  string logi_vacuum_name = "logi_" + nameVac + "_name";
384  G4LogicalVolume* logi_vacuum = new G4LogicalVolume(vacuum.geo, mat_vacuum, logi_vacuum_name);
385  if (flag_limitStep) logi_vacuum->SetUserLimits(new G4UserLimits(stepMax));
386  setColor(*logi_vacuum, "#000000");
387  setVisibility(*logi_vacuum, false);
388 
389  vacuum.logi = logi_vacuum;
390 
391  //put volume
392  string phys_vacuum_name = "phys_" + nameVac + "_name";
393  //new G4PVPlacement(vacuum.transform, logi_vacuum, phys_vacuum_name, &topVolume, false, 0);
394  new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_torus, false, 0);
395 
396  elements[nameVac] = vacuum;
397  }
398  }
399 
400 
401  //--------------------------------------------------------------------------------------------
402  //- Gate shields, end-of-tunnel concrete shields, polyethylene shields, collimator shields
403 
404  std::vector<std::string> shields;
405  boost::split(shields, m_config.getParameterStr("Shield"), boost::is_any_of(" "));
406  for (const auto& name : shields) {
407  prep = name + ".";
408 
409  //- Shield made as box with optional subtracted box-shaped inner space (hole)
410 
411  double shield_W = m_config.getParameter(prep + "W") * unitFactor;
412  double shield_H = m_config.getParameter(prep + "H") * unitFactor;
413  double shield_L = m_config.getParameter(prep + "L") * unitFactor;
414  double shield_X0 = m_config.getParameter(prep + "X0") * unitFactor;
415  double shield_Y0 = m_config.getParameter(prep + "Y0") * unitFactor;
416  double shield_Z0 = m_config.getParameter(prep + "Z0") * unitFactor;
417 
418  double shield_hole_W = m_config.getParameter(prep + "holeW", 0) * unitFactor;
419  double shield_hole_H = m_config.getParameter(prep + "holeH", 0) * unitFactor;
420  double shield_hole_L = m_config.getParameter(prep + "holeL", 0) * unitFactor;
421  double shield_hole_dX = m_config.getParameter(prep + "holeDX", 0) * unitFactor;
422  double shield_hole_dY = m_config.getParameter(prep + "holeDY", 0) * unitFactor;
423  double shield_hole_dZ = m_config.getParameter(prep + "holeDZ", 0) * unitFactor;
424 
425  double shield_PHI = m_config.getParameter(prep + "PHI");
426 
427  // storable element
428  FarBeamLineElement shield;
429 
430  shield.transform = G4Translate3D(shield_X0, shield_Y0, shield_Z0);
431  shield.transform = shield.transform * G4RotateY3D(shield_PHI / Unit::rad);
432 
433  G4Transform3D transform_shield_hole = G4Translate3D(shield_hole_dX, shield_hole_dY, shield_hole_dZ);
434 
435  //define geometry
436  string geo_shieldx_name = "geo_" + name + "x_name";
437  string geo_shield_hole_name = "geo_" + name + "_hole_name";
438  string geo_shield_name = "geo_" + name + "_name";
439 
440  if (shield_hole_W == 0 || shield_hole_H == 0 || shield_hole_L == 0) {
441  G4Box* geo_shield = new G4Box(geo_shield_name, shield_W / 2.0, shield_H / 2.0, shield_L / 2.0);
442 
443  shield.geo = geo_shield;
444  } else {
445  G4Box* geo_shieldx = new G4Box(geo_shieldx_name, shield_W / 2.0, shield_H / 2.0, shield_L / 2.0);
446  G4Box* geo_shield_hole = new G4Box(geo_shield_hole_name, shield_hole_W / 2.0, shield_hole_H / 2.0, shield_hole_L / 2.0);
447  G4SubtractionSolid* geo_shield = new G4SubtractionSolid(geo_shield_name, geo_shieldx, geo_shield_hole,
448  transform_shield_hole);
449 
450  shield.geo = geo_shield;
451  }
452 
453  string strMat_shield = m_config.getParameterStr(prep + "Material");
454  G4Material* mat_shield = Materials::get(strMat_shield);
455 
456  string logi_shield_name = "logi_" + name + "_name";
457  G4LogicalVolume* logi_shield = new G4LogicalVolume(shield.geo, mat_shield, logi_shield_name);
458 
459  shield.logi = logi_shield;
460 
461  //put volume
462  setColor(*logi_shield, "#0000CC");
463  //setVisibility(*logi_shield, false);
464  string phys_shield_name = "phys_" + name + "_name";
465  new G4PVPlacement(shield.transform, shield.logi, phys_shield_name, &topVolume, false, 0);
466 
467  elements[name] = shield;
468  }
469 
470 
471  //--------------
472  //- Tube (virtual tube for radiation level study)
473 
474  //define geometry
475  G4Tubs* geo_Tube = new G4Tubs("geo_Tube_name", 3995 * CLHEP::mm, 4000 * CLHEP::mm, 29 * CLHEP::m, 0. * CLHEP::deg, 360.*CLHEP::deg);
476  G4Material* mat_Tube = Materials::get("G4_Si");
477  G4LogicalVolume* logi_Tube = new G4LogicalVolume(geo_Tube, mat_Tube, "logi_Tube_name");
478 
479  //put volume
480  setColor(*logi_Tube, "#CC0000");
481  setVisibility(*logi_Tube, false);
482  bool radiation_study = false;
483  // cppcheck-suppress knownConditionTrueFalse
484  if (radiation_study && elements.count("GateShieldL")) {
485  new G4PVPlacement(elements["GateShieldL"].transform, logi_Tube, "phys_Tube_name", &topVolume, false, 0);
486  }
487 
488 
489  //---------------------------
490  // for dose simulation
491  //---------------------------
492 
493  // cppcheck-suppress knownConditionTrueFalse
494  if (radiation_study) {
495  //neutron shield (poly)
496  if (elements.count("PolyShieldL"))
497  elements["PolyShieldL"].logi->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1001));
498  if (elements.count("PolyShieldR"))
499  elements["PolyShieldR"].logi->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1002));
500 
501  //additional neutron shield (concrete)
502  if (elements.count("ConcreteShieldL"))
503  elements["ConcreteShieldL"].logi->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1003));
504  if (elements.count("ConcreteShieldR"))
505  elements["ConcreteShieldR"].logi->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1004));
506 
507  //gate shield (concrete)
508  if (elements.count("GateShieldL"))
509  elements["GateShieldL"].logi->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1005));
510  if (elements.count("GateShieldR"))
511  elements["GateShieldR"].logi->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1006));
512 
513  //virtual material outside gate-shield
514  logi_Tube->SetSensitiveDetector(new BkgSensitiveDetector("IR", 1007));
515  }
516 
517 
518  //------------------
519  //- Collimators
520 
521  std::vector<std::string> collimators;
522  boost::split(collimators, m_config.getParameterStr("Collimator"), boost::is_any_of(" "));
523  for (const auto& name : collimators) {
524  //- Collimators consist of two independent jaws (trapezoids), identical in shape, positioned opposite to each other
525  //- Each jaw consists of copper body and high Z head
526 
527  prep = name + ".";
528 
529  string type = m_config.getParameterStr(prep + "type");
530  string motherVolume = m_config.getParameterStr(prep + "MotherVolume");
531  string motherVolumeVacuum = motherVolume + "Vac";
532 
533  // If zz < 0 (positioned at negative z) vertical collimator is flipped when rotated into Mother Volume system
534  G4Scale3D scale;
535  G4Rotate3D rotation;
536  G4Translate3D translation;
537  elements[motherVolumeVacuum].transform.getDecomposition(scale, rotation, translation);
538  double zz = rotation.zz();
539 
540  // d1, d2 are collimator jaws displacements from beam center, d1<0, d2>0
541  // Z is collimator position inside its Mother Volume
542  double collimator_d1 = m_config.getParameter(prep + "d1") * unitFactor;
543  double collimator_d2 = m_config.getParameter(prep + "d2") * unitFactor;
544  double collimator_fullH = m_config.getParameter(prep + "fullH") * unitFactor;
545  double collimator_headH = m_config.getParameter(prep + "headH") * unitFactor;
546  double collimator_minW = m_config.getParameter(prep + "minW") * unitFactor;
547  double collimator_maxW = m_config.getParameter(prep + "maxW") * unitFactor;
548  double collimator_th = m_config.getParameter(prep + "th") * unitFactor;
549  double collimator_Z = m_config.getParameter(prep + "Z") * unitFactor;
550 
551  B2WARNING("Collimator " << name << " displacement d1 is set to " << collimator_d1 << "mm (must be negative)");
552  B2WARNING("Collimator " << name << " displacement d2 is set to " << collimator_d2 << "mm (must be positive)");
553 
554 
555  // Collimator heads
556 
557  // dx1,2 dy1,2 dz are trapezoid dimensions
558  double head_dx1;
559  double head_dx2;
560  double head_dy1;
561  double head_dy2;
562  double head_dz = collimator_headH / 2.0;
563  if (type == "vertical") {
564  head_dx1 = collimator_th / 2.0;
565  head_dx2 = collimator_th / 2.0;
566  head_dy1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
567  head_dy2 = collimator_minW / 2.0;
568  } else {
569  head_dx1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
570  head_dx2 = collimator_minW / 2.0;
571  head_dy1 = collimator_th / 2.0;
572  head_dy2 = collimator_th / 2.0;
573  }
574 
575  // storable elements
576  FarBeamLineElement collimator_head1;
577  FarBeamLineElement collimator_head2;
578 
579  // move collimator to position on beam line
580  G4Transform3D transform_head1 = G4Translate3D(0.0, 0.0, collimator_Z);
581  G4Transform3D transform_head2 = G4Translate3D(0.0, 0.0, collimator_Z);
582 
583  // rotate and move collimator jaws to their relative positions
584  if (type == "vertical") {
585  transform_head1 = transform_head1 * G4Translate3D(0.0, -head_dz + collimator_d1, 0.0);
586  transform_head1 = transform_head1 * G4RotateX3D(-M_PI / 2 / Unit::rad);
587 
588  transform_head2 = transform_head2 * G4Translate3D(0.0, head_dz + collimator_d2, 0.0);
589  transform_head2 = transform_head2 * G4RotateX3D(M_PI / 2 / Unit::rad);
590  } else {
591  if (zz > 0) {
592  transform_head1 = transform_head1 * G4Translate3D(-head_dz + collimator_d1, 0.0, 0.0);
593  transform_head1 = transform_head1 * G4RotateY3D(M_PI / 2 / Unit::rad);
594 
595  transform_head2 = transform_head2 * G4Translate3D(head_dz + collimator_d2, 0.0, 0.0);
596  transform_head2 = transform_head2 * G4RotateY3D(-M_PI / 2 / Unit::rad);
597  } else {
598  transform_head1 = transform_head1 * G4Translate3D(head_dz - collimator_d1, 0.0, 0.0);
599  transform_head1 = transform_head1 * G4RotateY3D(-M_PI / 2 / Unit::rad);
600 
601  transform_head2 = transform_head2 * G4Translate3D(-head_dz - collimator_d2, 0.0, 0.0);
602  transform_head2 = transform_head2 * G4RotateY3D(M_PI / 2 / Unit::rad);
603  }
604  }
605 
606  collimator_head1.transform = transform_head1;
607  collimator_head2.transform = transform_head2;
608 
609  // define geometry
610  string geo_headx_name = "geo_" + name + "_headx_name";
611 
612  string geo_head1_name = "geo_" + name + "_head1_name";
613  string geo_head2_name = "geo_" + name + "_head2_name";
614 
615  G4VSolid* geo_headx = new G4Trd(geo_headx_name, head_dx1, head_dx2, head_dy1, head_dy2, head_dz);
616 
617  G4VSolid* geo_head1 = new G4IntersectionSolid(geo_head1_name, geo_headx, elements[motherVolumeVacuum].geo,
618  collimator_head1.transform.inverse());
619  G4VSolid* geo_head2 = new G4IntersectionSolid(geo_head2_name, geo_headx, elements[motherVolumeVacuum].geo,
620  collimator_head2.transform.inverse());
621 
622  collimator_head1.geo = geo_head1;
623  collimator_head2.geo = geo_head2;
624 
625  // define logical volume
626  string strMat_head = m_config.getParameterStr(prep + "HeadMaterial");
627  G4Material* mat_head = Materials::get(strMat_head);
628  string logi_head1_name = "logi_" + name + "_head1_name";
629  string logi_head2_name = "logi_" + name + "_head2_name";
630  G4LogicalVolume* logi_head1 = new G4LogicalVolume(geo_head1, mat_head, logi_head1_name);
631  G4LogicalVolume* logi_head2 = new G4LogicalVolume(geo_head2, mat_head, logi_head2_name);
632  setColor(*logi_head1, "#CC0000");
633  setColor(*logi_head2, "#CC0000");
634  setVisibility(*logi_head1, false);
635  setVisibility(*logi_head2, false);
636 
637  // check if collimator is inside beam pipe
638  double volume_head1 = logi_head1->GetSolid()->GetCubicVolume();
639  double volume_head2 = logi_head2->GetSolid()->GetCubicVolume();
640 
641  collimator_head1.logi = logi_head1;
642  collimator_head2.logi = logi_head2;
643 
644  // put volume
645  string phys_head1_name = "phys_" + name + "_head1" + "_name";
646  string phys_head2_name = "phys_" + name + "_head2" + "_name";
647  if (volume_head1 != 0)
648  new G4PVPlacement(collimator_head1.transform, logi_head1, phys_head1_name, elements[motherVolumeVacuum].logi, false, 0);
649  if (volume_head2 != 0)
650  new G4PVPlacement(collimator_head2.transform, logi_head2, phys_head2_name, elements[motherVolumeVacuum].logi, false, 0);
651 
652  // to use it later in "intersect" and "subtract"
653  collimator_head1.transform = collimator_head1.transform * elements[motherVolumeVacuum].transform;
654  collimator_head2.transform = collimator_head2.transform * elements[motherVolumeVacuum].transform;
655 
656  string name_head1 = name + "_head1";
657  string name_head2 = name + "_head2";
658  elements[name_head1] = collimator_head1;
659  elements[name_head2] = collimator_head2;
660 
661 
662  // Collimator bodies
663 
664  // dx1,2 dy1,2 dz are trapezoid dimensions
665  double body_dx1;
666  double body_dx2;
667  double body_dy1;
668  double body_dy2;
669  double body_dz = (collimator_fullH - collimator_headH) / 2.0;
670  if (type == "vertical") {
671  body_dx1 = collimator_th / 2.0;
672  body_dx2 = collimator_th / 2.0;
673  body_dy1 = collimator_maxW / 2.0;
674  body_dy2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
675  } else {
676  body_dx1 = collimator_maxW / 2.0;
677  body_dx2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
678  body_dy1 = collimator_th / 2.0;
679  body_dy2 = collimator_th / 2.0;
680  }
681 
682  // storable elements
683  FarBeamLineElement collimator_body1;
684  FarBeamLineElement collimator_body2;
685 
686  // reuse head transfomation with additional shift
687  if (type == "vertical") {
688  collimator_body1.transform = G4Translate3D(0.0, -head_dz - body_dz, 0.0) * transform_head1;
689  collimator_body2.transform = G4Translate3D(0.0, head_dz + body_dz, 0.0) * transform_head2;
690  } else {
691  if (zz > 0) {
692  collimator_body1.transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head1;
693  collimator_body2.transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head2;
694  } else {
695  collimator_body1.transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head1;
696  collimator_body2.transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head2;
697  }
698  }
699 
700  // define geometry
701  string geo_bodyx_name = "geo_" + name + "_bodyx_name";
702 
703  string geo_body1_name = "geo_" + name + "_body1_name";
704  string geo_body2_name = "geo_" + name + "_body2_name";
705 
706  G4VSolid* geo_bodyx = new G4Trd(geo_bodyx_name, body_dx1, body_dx2, body_dy1, body_dy2, body_dz);
707 
708  G4VSolid* geo_body1 = new G4IntersectionSolid(geo_body1_name, geo_bodyx, elements[motherVolumeVacuum].geo,
709  collimator_body1.transform.inverse());
710  G4VSolid* geo_body2 = new G4IntersectionSolid(geo_body2_name, geo_bodyx, elements[motherVolumeVacuum].geo,
711  collimator_body2.transform.inverse());
712 
713  collimator_body1.geo = geo_body1;
714  collimator_body2.geo = geo_body2;
715 
716  // define logical volume
717  string strMat_body = m_config.getParameterStr(prep + "Material");
718  G4Material* mat_body = Materials::get(strMat_body);
719  string logi_body1_name = "logi_" + name + "_body1_name";
720  string logi_body2_name = "logi_" + name + "_body2_name";
721  G4LogicalVolume* logi_body1 = new G4LogicalVolume(geo_body1, mat_body, logi_body1_name);
722  G4LogicalVolume* logi_body2 = new G4LogicalVolume(geo_body2, mat_body, logi_body2_name);
723  setColor(*logi_body1, "#CC0000");
724  setColor(*logi_body2, "#CC0000");
725  setVisibility(*logi_body1, false);
726  setVisibility(*logi_body2, false);
727 
728  // check if collimator is inside beam pipe
729  double volume_body1 = logi_body1->GetSolid()->GetCubicVolume();
730  double volume_body2 = logi_body2->GetSolid()->GetCubicVolume();
731 
732  collimator_body1.logi = logi_body1;
733  collimator_body2.logi = logi_body2;
734 
735  // put volume
736  string phys_body1_name = "phys_" + name + "_body1" + "_name";
737  string phys_body2_name = "phys_" + name + "_body2" + "_name";
738  if (volume_body1 != 0)
739  new G4PVPlacement(collimator_body1.transform, logi_body1, phys_body1_name, elements[motherVolumeVacuum].logi, false, 0);
740  if (volume_body2 != 0)
741  new G4PVPlacement(collimator_body2.transform, logi_body2, phys_body2_name, elements[motherVolumeVacuum].logi, false, 0);
742 
743  // to use it later in "intersect" and "subtract"
744  collimator_body1.transform = collimator_body1.transform * elements[motherVolumeVacuum].transform;
745  collimator_body2.transform = collimator_body2.transform * elements[motherVolumeVacuum].transform;
746 
747  string name_body1 = name + "_body1";
748  string name_body2 = name + "_body2";
749  elements[name_body1] = collimator_body1;
750  elements[name_body2] = collimator_body2;
751  }
752  }
753  }
755 }
Belle2::BkgSensitiveDetector
The Class for BeamBackground Sensitive Detector.
Definition: BkgSensitiveDetector.h:25
Belle2::ir::SensitiveDetector
The IR Sensitive Detector class.
Definition: SensitiveDetector.h:41
Belle2::intersect
int intersect(const TRGCDCLpar &lp1, const TRGCDCLpar &lp2, CLHEP::HepVector &v1, CLHEP::HepVector &v2)
intersection
Definition: Lpar.cc:249
Belle2::geometry::setVisibility
void setVisibility(G4LogicalVolume &volume, bool visible)
Helper function to quickly set the visibility of a given volume.
Definition: utilities.cc:108
Belle2::geometry::setColor
void setColor(G4LogicalVolume &volume, const std::string &color)
Set the color of a logical volume.
Definition: utilities.cc:100
Belle2::geometry::GeometryTypes
GeometryTypes
Flag indiciating the type of geometry to be used.
Definition: GeometryManager.h:36
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
Belle2::ir::FarBeamLineElement
The struct for FarBeamLineElement.
Definition: GeoFarBeamLineCreator.h:30
Belle2::ir::FarBeamLineElement::logi
G4LogicalVolume * logi
Logical volume.
Definition: GeoFarBeamLineCreator.h:36
Belle2::ir::FarBeamLineElement::transform
G4Transform3D transform
Transformation.
Definition: GeoFarBeamLineCreator.h:32