Belle II Software  release-05-01-25
barrel.cc
1 #include <ecl/geometry/GeoECLCreator.h>
2 #include "G4LogicalVolume.hh"
3 #include "G4PVPlacement.hh"
4 #include <G4Tubs.hh>
5 #include <G4Box.hh>
6 #include <G4IntersectionSolid.hh>
7 #include <G4SubtractionSolid.hh>
8 #include <G4Region.hh>
9 #include <G4Trd.hh>
10 #include <G4TwoVector.hh>
11 #include <G4ExtrudedSolid.hh>
12 #include <G4PVReplica.hh>
13 
14 #include <iostream>
15 #include "CLHEP/Matrix/Vector.h"
16 #include "CLHEP/Matrix/Matrix.h"
17 #include "G4Vector3D.hh"
18 #include "G4Point3D.hh"
19 
20 #include "ecl/geometry/BelleLathe.h"
21 #include "ecl/geometry/BelleCrystal.h"
22 #include "ecl/geometry/shapes.h"
23 #include <geometry/Materials.h>
24 
25 using namespace std;
26 using namespace Belle2;
27 using namespace Belle2::geometry;
28 using namespace ECL;
29 
30 namespace {
31  Point_t p1[] = {{0.0, 0.0}, {263.4, 0.0}, {263.4, 197.2}, {248.6, 206.6}, {212.0, 152.6}, {145.2, 199.0}, {108.5, 146.6}, {46.4, 192.1}, {7.9, 140.5}, { -50.0, 185.2}, { -90.3, 134.2}, { -143.8, 178.3}, { -170.7, 145.0}};
32  Point_t p2[] = {{0.0, 0.0}, {384.8, 0.0}, {181.9, 172.5}, {136.8, 122.3}, {90.7, 163.9}, {46.5, 115.6}, {4.1, 156.3}, { -42.5, 109.9}, { -80.9, 149.5}, { -110.9, 119.8}};
33  Point_t p3[] = {{0.0, 0.0}, {343.2, 0.0}, {209.6, 144.2}, {158.0, 98.1}, {125.8, 135.0}, {74.3, 92.4}, {46.1, 128.4}, { -6.9, 87.3}, { -32.5, 121.3}, { -67.3, 95.8}};
34  Point_t p4[] = {{0.0, 0.0}, {311.5, 0.0}, {229.5, 116.7}, {170.3, 76.5}, {150.6, 107.6}, {92.8, 71.9}, {75.0, 101.2}, {16.1, 67.5}, {1.0, 95.2}, { -37.8, 75.6}};
35  Point_t p5[] = {{0.0, 0.0}, {289.2, 0.0}, {243.1, 92.2}, {177.4, 60.9}, {166.6, 85.0}, {103.3, 57.8}, {94.3, 79.8}, {30.0, 55.5}, {22.9, 75.3}, { -19.3, 61.8}};
36  Point_t p6[] = {{0.0, 0.0}, {276.1, 0.0}, {253.3, 73.3}, {185.2, 53.4}, {180.7, 69.1}, {110.5, 51.9}, {107.5, 65.3}, {38.7, 51.8}, {36.7, 62.8}, { -7.5, 56.8}};
37  Point_t p7[] = {{0.0, 0.0}, {253.1, 0.0}, {244.8, 62.7}, {171.0, 54.0}, {170.2, 60.3}, {100.8, 55.6}, {100.3, 59.6}, {38.2, 58.1}, {38.1, 59.5}, {0.0, 60.9}};
38  Point_t p11[] = {{0.0, 0.0}, {314.1, 0.0}, {351.9, 75.6}, {312.4, 96.9}, {297.5, 69.4}, {238.2, 102.9}, {221.4, 73.2}, {163.1, 108.8}, {143.5, 77.4}, {87.1, 114.8}, {64.3, 81.7}, {9.3, 121.5}, {0.0, 108.4}};
39  Point_t* pp[] = {p1, p2, p3, p4, p5, p6, p7, p11};
40  int npp[] = {13, 10, 10, 10, 10, 10, 10, 13};
41 
42  G4VSolid* get_crystal_support(int n)
43  {
44  std::vector<G4TwoVector> p; p.reserve(npp[n]);
45  for (int i = 0; i < npp[n]; i++)
46  p.push_back(G4TwoVector(pp[n][i].x, pp[n][i].y));
47 
48  string name("barrel_crystal_support");
49  name += to_string(n);
50 
51  G4TwoVector off(0, 0);
52  return new G4ExtrudedSolid(name, p, 1.5, off, 1, off, 1);
53  }
54 }
55 
56 void Belle2::ECL::GeoECLCreator::barrel(G4LogicalVolume& _top)
57 {
58  G4LogicalVolume* top = &_top;
59 
60  // vector<cplacement_t> bp = load_placements("/ecl/data/crystal_placement_barrel.dat");
61  vector<cplacement_t> bp = load_placements(m_sap, ECLParts::barrel);
62 
63  vector<cplacement_t>::iterator fp = find_if(bp.begin(), bp.end(), [](const cplacement_t& p) {
64  const int ECL_barrel_part = 1001;
65  return p.nshape == ECL_barrel_part;
66  });
67  // global transformation before placing the whole forward part in the top logical volume
68  G4Transform3D gT = (fp == bp.end()) ? G4Translate3D(0, 0, 0) : get_transform(*fp);
69  if (fp != bp.end()) bp.erase(fp); // now not needed
70 
71  bool b_crystals = true;
72  bool b_forward_support_legs = true;
73  bool b_backward_support_legs = true;
74  bool b_forward_support_ring = true;
75  bool b_backward_support_ring = true;
76  bool b_forward_support_wall = true;
77  bool b_backward_support_wall = true;
78  bool b_crystal_support = true;
79  bool b_preamplifier = true;
80  bool b_septum_walls = true;
81  bool b_inner_support_wall = true;
82  bool b_support_ribs = true;
83  bool b_outer_plates = true;
84  bool b_forward_part5 = true;
85  bool b_crystal_holder = true;
86  bool b_pipe = true;
87 
88  double phi0 = M_PI / 2 - M_PI / 72, dphi = 2 * M_PI; //2*2*M_PI/72;
89 
90  int nseg = 72;
91 
92  int overlap = m_overlap;
93 
94  G4LogicalVolume* sectorlogical; // tilted sector
95  {
96  double R = 1251.6;
97  double zp = 4.9, zm = 2.5;
98  // zr_t bint[] = {{-1000.8-zm, 0}, {1982.6 + zp,0}, {2288.,198.2 - zp*tand(32.982)}, {2288.,404 + 12}, {-1219.,404 + 12}, {-1219., 404 -112.6-zm*tand(52.897)-3.2}};
99  zr_t bint[] = {{ -1000.8 - zm, 0}, {1982.6 + zp, 0}, {2288., 198.2 - zp * tand(32.982)}, {2288., 1582 - R}, {2380, 1582 - R}, {2380, 404 + 12}, { -1330., 404 + 12}, { -1330., 1543 + 84 - R}, { -1219., 1543 + 84 - R}, { -1219., 404 - 112.6 - zm * tand(52.897) - 3.2}};
100  for (unsigned int i = 0; i < sizeof(bint) / sizeof(zr_t); i++) bint[i].r += R;
101  std::vector<zr_t> contourb(bint, bint + sizeof(bint) / sizeof(zr_t));
102  G4VSolid* sect0 = new BelleLathe("sect0", phi0, 2 * M_PI, contourb);
103 
104  G4Vector3D n0(cos(phi0), sin(phi0), 0);
105  // shift r0 by half of septum wall thickness so from one side of the sector the wall will be completely inside the sector volume
106  G4Point3D r0 = R * n0 + 0.5 / 2 * (G4RotateZ3D(M_PI / 2) * n0), rt0 = r0 + G4RotateZ3D(M_PI / 144 - M_PI / 72) * (G4Vector3D(0, 417,
107  0)), r1 = G4RotateZ3D(2 * M_PI / 72) * r0, rt1 = G4RotateZ3D(2 * M_PI / 72) * rt0;
108  G4Transform3D t0 = G4RotateZ3D(-M_PI / 144) * G4Translate3D(-r0);
109  G4Point3D tr0 = t0 * r0, trt0 = t0 * rt0, tr1 = t0 * r1, trt1 = t0 * rt1;
110  G4Vector3D tn1 = trt1 - tr1;
111  double v0 = -(tr1.y() - tr0.y()) / tn1.y();
112  double v1 = -(tr1.y() - trt0.y()) / tn1.y();
113 
114  G4Point3D u0 = tr1 + v0 * tn1, u1 = tr1 + v1 * tn1;
115  G4VSolid* sect1 = new G4Trd("sect1", -u0.x() / 2, (trt0 - u1).x() / 2, 2000, 2000, trt0.y() / 2);
116  G4IntersectionSolid* sector = new G4IntersectionSolid("sector", sect0, sect1,
117  G4Translate3D(r0)*G4RotateZ3D(M_PI / 144)*G4Translate3D(u0.x() / 2, trt0.y() / 2, (2288 - 1219.) / 2)*G4RotateX3D(-M_PI / 2));
118  sectorlogical = new G4LogicalVolume(sector, Materials::get("G4_AIR"), "ECLBarrelSectorLogical", 0, 0, 0);
119  sectorlogical->SetVisAttributes(att("air"));
120 
121  // Set up region for production cuts
122  G4Region* aRegion = new G4Region("ECLBarrelSector");
123  sectorlogical->SetRegion(aRegion);
124  aRegion->AddRootLogicalVolume(sectorlogical);
125 
126  for (int i = 0; i < nseg; i++) {
127  double phi = i * M_PI / 36 - M_PI / 2;
128  new G4PVPlacement(gT * G4RotateZ3D(phi), sectorlogical, suf("ECLBarrelSectorPhysical", i), top, false, i, overlap);
129  }
130  }
131 
132  // phi septum wall
133  // cppcheck-suppress knownConditionTrueFalse
134  if (b_septum_walls) {
135  double zh = (0.5 - 0.001) / 2;
136  G4ThreeVector psw[] = {G4ThreeVector(-1000.8, 0, -zh), G4ThreeVector(1982.6, 0, -zh), G4ThreeVector(2288, 198.2, -zh),
137  G4ThreeVector(2288, 404, -zh), G4ThreeVector(-1219., 404, -zh), G4ThreeVector(-1219., 404 - 112.6, -zh),
138  G4ThreeVector(-1000.8, 0, zh), G4ThreeVector(1982.6, 0, zh), G4ThreeVector(2288, 198.2, zh),
139  G4ThreeVector(2288, 404, zh), G4ThreeVector(-1219., 404, zh), G4ThreeVector(-1219., 404 - 112.6, zh)
140  };
141  G4VSolid* pswsolid = new BelleCrystal("pswsolid", 6, psw);
142  G4LogicalVolume* pswlogical = new G4LogicalVolume(pswsolid, Materials::get("A5052"), "pswlogical", 0, 0, 0);
143  pswlogical->SetVisAttributes(att("alum"));
144  new G4PVPlacement(G4RotateZ3D(2.5 / 180 * M_PI)*G4Translate3D(0, 1251.6 + 0.01,
145  0)*G4RotateZ3D(1.25 / 180 * M_PI)*G4RotateY3D(-M_PI / 2), pswlogical, "phiwallphysical", sectorlogical, false, 0, overlap);
146  }
147 
148  // cppcheck-suppress knownConditionTrueFalse
149  if (b_forward_support_ring) {
150  zr_t vc4[] = {{0, 1449.2}, {0, 1582}, {90, 1582}, {90, 1582 - 75.7}};
151  std::vector<zr_t> contour4(vc4, vc4 + 4);
152  G4VSolid* part4solid = new BelleLathe("part4solid", phi0, dphi, contour4);
153  G4LogicalVolume* part4logical = new G4LogicalVolume(part4solid, Materials::get("SUS304"), "part4logical", 0, 0, 0);
154  part4logical->SetVisAttributes(att("iron"));
155  new G4PVPlacement(gT * G4Translate3D(0, 0, 2290), part4logical, "part4physical", top, false, 0, overlap);
156  }
157 
158  // cppcheck-suppress knownConditionTrueFalse
159  if (b_backward_support_ring) {
160  zr_t vc6[] = {{0, 1543}, { -105, 1543}, { -105, 1543 + 84}, {0, 1543 + 84}};
161  std::vector<zr_t> contour6(vc6, vc6 + 4);
162  G4VSolid* part6solid = new BelleLathe("part6solid", phi0, dphi, contour6);
163  G4LogicalVolume* part6logical = new G4LogicalVolume(part6solid, Materials::get("SUS304"), "part6logical", 0, 0, 0);
164  part6logical->SetVisAttributes(att("iron"));
165  new G4PVPlacement(gT * G4Translate3D(0, 0, -1225), part6logical, "part6physical", top, false, 0, overlap);
166  }
167 
168  // cppcheck-suppress knownConditionTrueFalse
169  if (b_forward_support_wall) {
170  zr_t vc2[] = {{0, 1180.1}, {0, 1249.5}, {54.3, 1249.5}, {445.8, 1503.6}, {445.8 + 4 * sind(32.982), 1503.6 - 4 * cosd(32.982)}, {54.3 + 4 * sind(32.982 / 2), 1249.5 - 4}, {4, 1249.5 - 4}, {4, 1180.1}};
171  std::vector<zr_t> contour2(vc2, vc2 + 8);
172  G4VSolid* part2solid = new BelleLathe("part2solid", phi0, dphi, contour2);
173  G4LogicalVolume* part2logical = new G4LogicalVolume(part2solid, Materials::get("A5083"), "part2logical", 0, 0, 0);
174  part2logical->SetVisAttributes(att("alum"));
175  new G4PVPlacement(gT * G4Translate3D(0, 0, 1930), part2logical, "part2physical", top, false, 0, overlap);
176  }
177 
178  // cppcheck-suppress knownConditionTrueFalse
179  if (b_backward_support_wall) {
180  zr_t vc3[] = {{0, 1180.1}, { -4, 1180.1}, { -4, 1249.5 - 4}, { -61.9 - 4 * sind(52.897 / 2), 1249.5 - 4}, { -285, 1539.1 - 4 * sind(52.897 / 2)}, { -285, 1618.8}, { -285 + 4, 1618.8}, { -285 + 4, 1539.1}, { -61.9, 1249.5}, {0, 1249.5}};
181  std::vector<zr_t> contour3(vc3, vc3 + 10);
182  G4VSolid* part3solid = new BelleLathe("part3solid", phi0, dphi, contour3);
183  G4LogicalVolume* part3logical = new G4LogicalVolume(part3solid, Materials::get("A5083"), "part3logical", 0, 0, 0);
184  part3logical->SetVisAttributes(att("alum"));
185  new G4PVPlacement(gT * G4Translate3D(0, 0, -940), part3logical, "part3physical", top, false, 0, overlap);
186  }
187 
188  // cppcheck-suppress knownConditionTrueFalse
189  if (b_inner_support_wall) {
190  G4VSolid* part1solid = new G4Tubs("part1solid", 1250.1, 1250.1 + 1.5, (988.5 + 1972.5) / 2, phi0, dphi);
191  G4LogicalVolume* part1logical = new G4LogicalVolume(part1solid, Materials::get("A5083"), "part1logical", 0, 0, 0);
192  part1logical->SetVisAttributes(att("alum2"));
193  new G4PVPlacement(gT * G4Translate3D(0, 0, (-988.5 + 1972.5) / 2), part1logical, "part1physical", top, false, 0, overlap);
194  }
195 
196  if (b_support_ribs || b_outer_plates || b_forward_part5) {
197  // double zh = (1219.+2288.)/2, dz = (-1219.+2288.)/2;
198  double zh = (1330. + 2380.) / 2, dz = (-1330. + 2380.) / 2;
199  if (b_support_ribs) {
200  G4ThreeVector vc9[] = {G4ThreeVector(0, 0, -zh), G4ThreeVector(25, 0, -zh), G4ThreeVector(25, 24.5 - 25.0 * tand(2.5), -zh), G4ThreeVector(0, 24.5, -zh),
201  G4ThreeVector(0, 0, zh), G4ThreeVector(25, 0, zh), G4ThreeVector(25, 24.5 - 25.0 * tand(2.5), zh), G4ThreeVector(0, 24.5, zh)
202  };
203  G4VSolid* part9solid = new BelleCrystal("part9solid", 4, vc9);
204  G4LogicalVolume* part9logical = new G4LogicalVolume(part9solid, Materials::get("SUS304"), "part9logical", 0, 0, 0);
205  part9logical->SetVisAttributes(att("silv"));
206  G4Transform3D tpart9a = G4RotateZ3D(2.5 / 180 * M_PI) * G4Translate3D(0, 1251.6,
207  0) * G4RotateZ3D(1.25 / 180 * M_PI) * G4Translate3D(0.5 / 2, 404 - 25 - 3.7 + 0.2, dz);
208  G4Transform3D tpart9b = G4RotateZ3D(-2.5 / 180 * M_PI) * G4Translate3D(0, 1251.6,
209  0) * G4RotateZ3D(1.25 / 180 * M_PI) * G4Translate3D(-0.5 / 2, 404 - 25 - 3.7 + 0.5, dz) * G4RotateY3D(M_PI);
210  new G4PVPlacement(tpart9a, part9logical, "part9aphysical", sectorlogical, false, 0, overlap);
211  new G4PVPlacement(tpart9b, part9logical, "part9bphysical", sectorlogical, false, 0, overlap);
212  }
213 
214  if (b_outer_plates) {
215  G4VSolid* part10solid = new G4Box("part10solid", 143. / 2, 8. / 2, zh);
216  G4LogicalVolume* part10logical = new G4LogicalVolume(part10solid, Materials::get("SUS304"), "part10logical", 0, 0, 0);
217  part10logical->SetVisAttributes(att("iron2"));
218  G4Transform3D tpart10 = G4Translate3D(0, 1251.6, 0) * G4Translate3D(-404 * tand(1.25), 404 - 1.1,
219  dz) * G4RotateZ3D(0.3 / 180 * M_PI);
220  new G4PVPlacement(tpart10, part10logical, "part10physical", sectorlogical, false, 0, overlap);
221  }
222 
223  // G4VSolid *part5solid = new G4Box("part5solid", 60./2, 45./2, 90./2);
224  // G4LogicalVolume* part5logical = new G4LogicalVolume(part5solid, Materials::get("SUS304"), "part5logical", 0, 0, 0);
225  // part5logical->SetVisAttributes(att("iron2"));
226  // new G4PVPlacement(G4RotateZ3D(2.796/180*M_PI)*G4Translate3D(0,1582+45./2,2380.-90./2), part5logical, "part5physical", sectorlogical, false, 1, overlap);
227 
228 
229  if (b_forward_part5) { // slice the element to fit into the sector
230  double R = 1251.6, phi1 = phi0 + 1.25 * M_PI / 180, phi2 = (2.796 + 90 - 5) / 180 * M_PI;
231  G4Vector3D n0(cos(phi0), sin(phi0), 0), n1(cos(phi1), sin(phi1), 0), n2(cos(phi2), sin(phi2), 0);
232  // shift r0 by half of septum wall thickness so from one side of the sector the wall will be completely inside the sector volume
233  G4Point3D r0 = R * n0 + 0.5 / 2 * (G4RotateZ3D(M_PI / 2) * n0);
234 
235  // n2*(r0 + t*n1 - n2*1582) = 0;
236  double t0 = -n2 * (r0 - n2 * 1582) / (n2 * n1);
237  double t1 = -n2 * (r0 - n2 * (1582 + 45)) / (n2 * n1);
238 
239  G4Point3D x0 = r0 + t0 * n1;
240  G4Point3D x1 = (t1 - t0) * n1;
241 
242  Point_t c5a[] = {{0, 0}, {30 * cosd(5 - 2.796), -30 * sind(5 - 2.796)}, {30 * cosd(5 - 2.796) + 45 * cos(phi2), -30 * sind(5 - 2.796) + 45 * sin(phi2)}, {x1.x(), x1.y()}};
243  G4ThreeVector vc5a[8];
244  for (int i = 0; i < 4; i++) vc5a[i] = G4ThreeVector(c5a[i].x, c5a[i].y, -90. / 2);
245  for (int i = 0; i < 4; i++) vc5a[i + 4] = G4ThreeVector(c5a[i].x, c5a[i].y, 90. / 2);
246  G4VSolid* part5asolid = new BelleCrystal("part5asolid", 4, vc5a);
247  G4LogicalVolume* part5alogical = new G4LogicalVolume(part5asolid, Materials::get("SUS304"), "part5alogical", 0, 0, 0);
248  part5alogical->SetVisAttributes(att("iron2"));
249  new G4PVPlacement(G4RotateZ3D(2 * M_PI / 72)*G4Translate3D(x0.x(), x0.y(), 2380. - 90. / 2), part5alogical, "part5aphysical",
250  sectorlogical, false, 0, overlap);
251 
252  Point_t c5b[] = {{ -30 * cosd(5 - 2.796), 30 * sind(5 - 2.796)}, {0, 0}, {x1.x(), x1.y()}, { -30 * cosd(5 - 2.796) + 44.8 * cos(phi2), 30 * sind(5 - 2.796) + 44.8 * sin(phi2)}};
253  G4ThreeVector vc5b[8];
254  for (int i = 0; i < 4; i++) vc5b[i] = G4ThreeVector(c5b[i].x, c5b[i].y, -90. / 2);
255  for (int i = 0; i < 4; i++) vc5b[i + 4] = G4ThreeVector(c5b[i].x, c5b[i].y, 90. / 2);
256  G4VSolid* part5bsolid = new BelleCrystal("part5bsolid", 4, vc5b);
257  G4LogicalVolume* part5blogical = new G4LogicalVolume(part5bsolid, Materials::get("SUS304"), "part5blogical", 0, 0, 0);
258  part5blogical->SetVisAttributes(att("iron2"));
259  new G4PVPlacement(G4Translate3D(x0.x(), x0.y(), 2380. - 90. / 2), part5blogical, "part5bphysical", sectorlogical, false, 0,
260  overlap);
261  }
262 
263  }
264 
265  // cppcheck-suppress knownConditionTrueFalse
266  if (b_forward_support_legs) {
267  G4ThreeVector vc11a[] = {G4ThreeVector(0, 0, -63. / 2), G4ThreeVector(35, 24.5, -63. / 2), G4ThreeVector(35, 158 - 35, -63. / 2), G4ThreeVector(0, 158 - 35, -63. / 2),
268  G4ThreeVector(0, 0, 63. / 2), G4ThreeVector(35, 24.5, 63. / 2), G4ThreeVector(35, 158 - 35, 63. / 2), G4ThreeVector(0, 158 - 35, 63. / 2)
269  };
270  G4VSolid* part11asolid = new BelleCrystal("part11asolid", 4, vc11a);
271  G4VSolid* part11bsolid = new G4Box("part11bsolid", 63. / 2, 35. / 2, 100. / 2);
272  G4LogicalVolume* part11alogical = new G4LogicalVolume(part11asolid, Materials::get("SUS304"), "part11alogical", 0, 0, 0);
273  part11alogical->SetVisAttributes(att("iron2"));
274  G4LogicalVolume* part11blogical = new G4LogicalVolume(part11bsolid, Materials::get("SUS304"), "part11blogical", 0, 0, 0);
275  part11blogical->SetVisAttributes(att("iron"));
276 
277  double Ro = 158.4;
278  zr_t vsf[] = {{0, 0}, {35, 24.5}, {35, 158 - 35}, {100, 158 - 35}, {100, Ro}, {0, Ro}};
279  for (zr_t* it = vsf; it != vsf + 6; it++) {it->r += 1668 - 158;}
280  std::vector<zr_t> csf(vsf, vsf + 6);
281  G4VSolid* sf = new BelleLathe("sf", 0, 2 * M_PI, csf);
282  G4LogicalVolume* sfl = new G4LogicalVolume(sf, Materials::get("G4_AIR"), "supportfwd", 0, 0, 0);
283  sfl->SetVisAttributes(att("air"));
284  new G4PVPlacement(gT * G4Translate3D(0, 0, 2380), sfl, "supportfwdphysical", top, false, 0, overlap);
285 
286  G4VSolid* sfs = new BelleLathe("sfs", -M_PI / 72, 2 * M_PI / 72, csf);
287  G4LogicalVolume* sfsl = new G4LogicalVolume(sfs, Materials::get("G4_AIR"), "supportfwdsector", 0, 0, 0);
288  sfsl->SetVisAttributes(att("air"));
289 
290  new G4PVPlacement(G4Translate3D(1668 - 158, 0, 0)*G4RotateZ3D(-M_PI / 2)*G4RotateY3D(-M_PI / 2), part11alogical, "part11aphysical",
291  sfsl, false, 0, overlap);
292  new G4PVPlacement(G4Translate3D(1668 - 158, 0, 0)*G4RotateZ3D(-M_PI / 2)*G4Translate3D(0, (158. - 35.) + 35. / 2, 100. / 2),
293  part11blogical, "part11bphysical", sfsl, false, 0, overlap);
294  new G4PVReplica("supportfwdrep", sfsl, sfl, kPhi, nseg, 2 * M_PI / 72, (2.796 - 2.5)*M_PI / 180);
295  }
296 
297  // cppcheck-suppress knownConditionTrueFalse
298  if (b_backward_support_legs) {
299  G4VSolid* part12asolid = new G4Box("part12asolid", 90. / 2, (118. - 35.) / 2, 35. / 2);
300  G4VSolid* part12bsolid = new G4Box("part12bsolid", 90. / 2, 35. / 2, 120. / 2);
301  G4LogicalVolume* part12alogical = new G4LogicalVolume(part12asolid, Materials::get("SUS304"), "part12alogical", 0, 0, 0);
302  part12alogical->SetVisAttributes(att("iron2"));
303  G4LogicalVolume* part12blogical = new G4LogicalVolume(part12bsolid, Materials::get("SUS304"), "part12blogical", 0, 0, 0);
304  part12blogical->SetVisAttributes(att("iron"));
305 
306  double Rb = 118.7;
307  zr_t vsb[] = {{0, 0}, {0, Rb}, { -120, Rb}, { -120, 118 - 35}, { -35, 118 - 35}, { -35, 0}};
308  for (zr_t* it = vsb; it != vsb + 6; it++) {it->r += 1668 - 118;}
309  std::vector<zr_t> csb(vsb, vsb + 6);
310  G4VSolid* sb = new BelleLathe("sb", 0, 2 * M_PI, csb);
311  G4LogicalVolume* sbl = new G4LogicalVolume(sb, Materials::get("G4_AIR"), "supportbkw", 0, 0, 0);
312  sbl->SetVisAttributes(att("air"));
313  new G4PVPlacement(gT * G4Translate3D(0, 0, -1330), sbl, "supportbkwphysical", top, false, 0, overlap);
314 
315  G4VSolid* sbs = new BelleLathe("sbs", -M_PI / 72, 2 * M_PI / 72, csb);
316  G4LogicalVolume* sbsl = new G4LogicalVolume(sbs, Materials::get("G4_AIR"), "supportbkwsector", 0, 0, 0);
317  sbsl->SetVisAttributes(att("air"));
318 
319  new G4PVPlacement(G4Translate3D(1668 - 118, 0, 0)*G4RotateZ3D(-M_PI / 2)*G4Translate3D(0, (118. - 35.) / 2, -35. / 2),
320  part12alogical, "part12aphysical", sbsl, false, 0, overlap);
321  new G4PVPlacement(G4Translate3D(1668 - 118, 0, 0)*G4RotateZ3D(-M_PI / 2)*G4Translate3D(0, (118. - 35.) + 35. / 2, -120. / 2),
322  part12blogical, "part12bphysical", sbsl, false, 0, overlap);
323 
324  new G4PVReplica("supportbkwrep", sbsl, sbl, kPhi, nseg, 2 * M_PI / 72, (2.796 - 2.5)*M_PI / 180);
325  }
326 
327  // cppcheck-suppress knownConditionTrueFalse
328  if (b_septum_walls) {
329  struct zwall_t { double zlow, zup;};
330  zwall_t zs[] = {{ -682.1, -883.9}, { -445.9, -571.7}, { -220.8, -274.2}, {0, 0}, {220.8, 274.2}, {445.9, 571.7}, {682.1, 883.9}, {936.6, 1220.5}, {1217.1, 1591.2}, {1532., 2007.4}};
331 
332  const double dz = 0.5 / 2; // half of theta fin thickness
333 
334  G4Vector3D n0(sind(90 + 2.5 + 1.25), -cosd(90 + 2.5 + 1.25), 0); // normal vector of the left sector wall pointing inside the sector
335  G4Point3D r0 = G4Point3D(-1251.6 * sind(2.5), 0, 0) + dz * n0; // point in the left sector wall shifted by phi fin
336  G4Vector3D n1(-sind(90 - 2.5 + 1.25), cosd(90 - 2.5 + 1.25),
337  0); // normal vector of the right sector wall pointing inside the sector
338  G4Point3D r1 = G4Point3D(1251.6 * sind(2.5), 0, 0) + dz * n1; // point in the right sector wall shifted by phi fin
339 
340  G4Vector3D n2(0, 1, 0); G4Point3D r2(0, 0, 0), r3(0, 310, 0);
341 
342  for (unsigned int i = 0; i < sizeof(zs) / sizeof(zwall_t); i++) {
343  double th = atan2(404, -(zs[i].zlow - zs[i].zup));
344 
345  G4Vector3D nf(0, -cos(th), sin(th));
346  // take into account finite thickness of the theta fin
347  G4Point3D rf(0, 0, copysign(dz / sin(th), cos(th)));
348 
349  // calculate three plane intersection point
350  auto inter = [](const G4Vector3D & l_n0, const G4Point3D & l_r0,
351  const G4Vector3D & l_n1, const G4Point3D & l_r1,
352  const G4Vector3D & l_n2, const G4Point3D & l_r2) -> G4Point3D {
353  CLHEP::HepMatrix A(3, 3);
354  CLHEP::HepVector B(3);
355  A[0][0] = l_n0.x(), A[0][1] = l_n0.y(), A[0][2] = l_n0.z();
356  A[1][0] = l_n1.x(), A[1][1] = l_n1.y(), A[1][2] = l_n1.z();
357  A[2][0] = l_n2.x(), A[2][1] = l_n2.y(), A[2][2] = l_n2.z();
358 
359  B[0] = l_r0 * l_n0;
360  B[1] = l_r1 * l_n1;
361  B[2] = l_r2 * l_n2;
362 
363  CLHEP::HepVector r = A.inverse() * B;
364  G4Point3D res(r[0], r[1], r[2]);
365  return res;
366  };
367 
368  G4Transform3D tfin(G4RotateX3D(-th + M_PI / 2)), tfini(G4RotateX3D(th - M_PI / 2));
369  G4ThreeVector t0 = tfini * inter(n0, r0, n2, r2, nf, rf);
370  G4ThreeVector t1 = tfini * inter(n1, r1, n2, r2, nf, rf);
371  G4ThreeVector t2 = tfini * inter(n0, r0, n2, r3, nf, rf);
372  G4ThreeVector t3 = tfini * inter(n1, r1, n2, r3, nf, rf);
373  // G4cout<<t0<<" "<<t1<<" "<<t2<<" "<<t3<<" "<<G4RotateX3D(th-M_PI/2)*G4Point3D(0,404,zs[i].zup-zs[i].zlow)<<G4endl;
374 
375  G4ThreeVector thfin[8];
376  thfin[0] = G4ThreeVector(t0.x(), t0.y(), -dz);
377  thfin[1] = G4ThreeVector(t1.x(), t1.y(), -dz);
378  thfin[2] = G4ThreeVector(t3.x(), t3.y(), -dz);
379  thfin[3] = G4ThreeVector(t2.x(), t2.y(), -dz);
380  thfin[4] = G4ThreeVector(t0.x(), t0.y(), +dz);
381  thfin[5] = G4ThreeVector(t1.x(), t1.y(), +dz);
382  thfin[6] = G4ThreeVector(t3.x(), t3.y(), +dz);
383  thfin[7] = G4ThreeVector(t2.x(), t2.y(), +dz);
384 
385  G4VSolid* thfinsolid = new BelleCrystal("thfinsolid", 4, thfin);
386  G4LogicalVolume* thfinlogical = new G4LogicalVolume(thfinsolid, Materials::get("A5052"), "thfinlogical", 0, 0, 0);
387  thfinlogical->SetVisAttributes(att("iron2"));
388  new G4PVPlacement(G4Translate3D(0, 1251.6, zs[i].zlow)*tfin, thfinlogical, "thfinphysical", sectorlogical, false, 0, overlap);
389  }
390  }
391 
392  // cppcheck-suppress knownConditionTrueFalse
393  if (b_pipe) {
394  double zh = (1219. + 2288.) / 2, dz = (-1219. + 2288.) / 2;
395  G4VSolid* tubec = new G4Tubs("tubec", 0, 12. / 2, zh, 0, 2 * M_PI);
396  G4VSolid* tubew = new G4Tubs("tubew", 0, 10. / 2, zh, 0, 2 * M_PI);
397  G4LogicalVolume* tubeclogical = new G4LogicalVolume(tubec, Materials::get("C1220"), "tubec", 0, 0, 0);
398  G4LogicalVolume* tubewlogical = new G4LogicalVolume(tubew, Materials::get("G4_WATER"), "tubew", 0, 0, 0);
399  new G4PVPlacement(G4Translate3D(0, 0, 0), tubewlogical, "tubewphysical", tubeclogical, false, 1, overlap);
400  tubeclogical->SetVisAttributes(att("iron2"));
401  new G4PVPlacement(G4Translate3D(25, 1640, dz), tubeclogical, "tubecphysical", sectorlogical, false, 0, overlap);
402  new G4PVPlacement(G4Translate3D(-35, 1640, dz), tubeclogical, "tubecphysical", sectorlogical, false, 1, overlap);
403  }
404 
405  // cppcheck-suppress knownConditionTrueFalse
406  if (b_preamplifier) {
407  for (vector<cplacement_t>::const_iterator it = bp.begin(); it != bp.end(); ++it) {
408  G4Transform3D twc = get_transform(*it);
409  int indx = it - bp.begin();
410  auto pv = new G4PVPlacement(twc * G4TranslateZ3D(300 / 2 + 0.250 + get_pa_box_height() / 2), get_preamp(),
411  suf("phys_barrel_preamplifier", indx),
412  sectorlogical, false, 72 + 9 * (indx / 2) + (indx % 2), 0);
413  if (overlap)pv->CheckOverlaps(1000);
414  }
415  }
416 
417  // cppcheck-suppress knownConditionTrueFalse
418  if (b_crystal_holder) {
419  G4VSolid* holder0 = new G4Box("holder0", 19. / 2, 38. / 2, 6. / 2);
420  G4VSolid* holder1 = new G4Box("holder1", 7. / 2, 26. / 2, 7. / 2);
421  G4VSolid* holder = new G4SubtractionSolid("holder", holder0, holder1);
422  G4LogicalVolume* holderlogical = new G4LogicalVolume(holder, Materials::get("A5052"), "holderlogical", 0, 0, 0);
423  holderlogical->SetVisAttributes(att("holder"));
424  for (vector<cplacement_t>::const_iterator it = bp.begin(); it != bp.end(); ++it) {
425  G4Transform3D twc = get_transform(*it);
426  int indx = it - bp.begin();
427  auto pv = new G4PVPlacement(twc * G4TranslateZ3D(300 / 2 + 0.250 + get_pa_box_height() + 38. / 2)*G4RotateZ3D(M_PI / 2)*G4RotateX3D(
428  M_PI / 2), holderlogical, "holderphysical", sectorlogical, false, indx, 0);
429  if (overlap)pv->CheckOverlaps(1000);
430  }
431  }
432 
433  // cppcheck-suppress knownConditionTrueFalse
434  if (b_crystals) {
435  // vector<shape_t*> cryst = load_shapes("/ecl/data/crystal_shape_barrel.dat");
436  vector<shape_t*> cryst = load_shapes(m_sap, ECLParts::barrel);
437  vector<G4LogicalVolume*> wrapped_crystals;
438  for (auto it = cryst.begin(); it != cryst.end(); it++) {
439  shape_t* s = *it;
440  wrapped_crystals.push_back(wrapped_crystal(s, "barrel", 0.17 - 0.0325));
441  }
442 
443  for (vector<cplacement_t>::const_iterator it = bp.begin(); it != bp.end(); ++it) {
444  const cplacement_t& t = *it;
445  auto s = find_if(cryst.begin(), cryst.end(), [&t](const shape_t* shape) {return shape->nshape == t.nshape;});
446  if (s == cryst.end()) continue;
447 
448  G4Transform3D twc = get_transform(t);
449  int indx = it - bp.begin();
450  new G4PVPlacement(twc, wrapped_crystals[s - cryst.begin()], suf("ECLBarrelWrappedCrystal_Physical", indx), sectorlogical, false,
451  72 + 9 * (indx / 2) + (indx % 2), overlap);
452  }
453  }
454 
455  // cppcheck-suppress knownConditionTrueFalse
456  if (b_crystal_support) {
457  G4LogicalVolume* tl[8];
458  for (int i = 0; i < 8; i++) {
459  G4VSolid* t = get_crystal_support(i);
460  tl[i] = new G4LogicalVolume(t, Materials::get("SUS304"), "crystal_support_logical", 0, 0, 0);
461  tl[i]->SetVisAttributes(att("iron2"));
462  }
463 
464  double offset[] = { -945, -607, -305, -45, 45, 305, 607, 909, 1240, 1628, 2040};
465  int ns[] = {8, 5, 6, 7, 7, 6, 5, 4, 3, 2, 1};
466  double offsetb[] = {263.4 - 210 - 30, 27.4, 10.6, 16.1, 26.2, 21.7, 35, 22 + 240};
467  bool flip[] = {1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1};
468  for (int i = 0; i < 11; i++) {
469  G4Transform3D tflip = flip[i] ? G4RotateY3D(-M_PI / 2) : G4RotateY3D(M_PI / 2);
470  G4Transform3D tr0(G4RotateZ3D(2.5 / 180 * M_PI)*G4Translate3D(0, 1251.6,
471  0)*G4RotateZ3D(1.25 / 180 * M_PI)*G4Translate3D(25 + 3. / 2 + 0.5 / 2, 404 - 25 - 3.7 + 24.5 - 25.0 * tand(2.5) - 0.5,
472  offset[i])*tflip * G4Translate3D(-offsetb[ns[i] - 1], 0, 0)*G4RotateX3D(M_PI));
473  new G4PVPlacement(tr0, tl[ns[i] - 1], "crystal_support_physical", sectorlogical, false, 0, overlap);
474  G4Transform3D tr1(G4RotateZ3D(-2.5 / 180 * M_PI)*G4Translate3D(0, 1251.6,
475  0)*G4RotateZ3D(1.25 / 180 * M_PI)*G4Translate3D(-25 - 3. / 2 - 0.5 / 2, 404 - 25 - 3.7 + 24.5 - 25.0 * tand(2.5) - 0.5,
476  offset[i])*tflip * G4Translate3D(-offsetb[ns[i] - 1], 0, 0)*G4RotateX3D(M_PI));
477  new G4PVPlacement(tr1, tl[ns[i] - 1], "crystal_support_physical", sectorlogical, false, 1, overlap);
478  }
479  }
480 
481 }
Belle2::ECL::shape_t
Definition: shapes.h:39
Belle2::geometry
Common code concerning the geometry representation of the detector.
Definition: CreatorBase.h:28
Belle2::ECL::Point_t
Definition: BelleCrystal.h:40
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::GeoECLCreator::barrel
void barrel(G4LogicalVolume &)
Make the ECL barrel and then place elements inside it.
Definition: barrel.cc:56
Belle2::ECL::BelleCrystal
a Belle crystal in Geant4
Definition: BelleCrystal.h:45
Belle2::ECL::zr_t
Definition: BelleLathe.h:35
Belle2::ECL::BelleLathe
Definition: BelleLathe.h:60
Belle2::ECL::cplacement_t
Definition: shapes.h:52