Belle II Software development
shapes.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// The shape structs below expose their parameters both as named members and as
10// the t[] array of an anonymous union; cppcheck does not see the named ones.
11// cppcheck-suppress-file unusedStructMember
12#include "ecl/geometry/shapes.h"
13#include <G4TessellatedSolid.hh>
14#include <G4TriangularFacet.hh>
15#include <G4QuadrangularFacet.hh>
16#include <G4ExtrudedSolid.hh>
17#include <G4Trap.hh>
18#include "ecl/geometry/BelleCrystal.h"
19#include <iostream>
20#include <fstream>
21#include <iterator>
22#include <CLHEP/Matrix/Vector.h>
23#include <CLHEP/Matrix/Matrix.h>
24#include <G4TwoVector.hh>
25#include <G4Vector3D.hh>
26#include <framework/utilities/FileSystem.h>
27#include <framework/logging/Logger.h>
28
29using namespace std;
30
31namespace Belle2 {
36 namespace ECL {
37 double cosd(double x) {return cos(x * (M_PI / 180));}
38 double sind(double x) {return sin(x * (M_PI / 180));}
39 double tand(double x) {return tan(x * (M_PI / 180));}
40
41 double volume(const G4ThreeVector& v0, const G4ThreeVector& v1, const G4ThreeVector& v2, const G4ThreeVector& v3)
42 {
43 G4ThreeVector x = v1 - v0, y = v2 - v0, z = v3 - v0;
44 return x.cross(y) * z;
45 }
46
47 G4ThreeVector newvertex(double d, const G4ThreeVector& v0, const G4ThreeVector& v1, const G4ThreeVector& v2,
48 const G4ThreeVector& v3)
49 {
50 G4ThreeVector x = v1 - v0, y = v2 - v0, z = v3 - v0;
51 G4ThreeVector nxy = x.cross(y).unit();
52 G4ThreeVector nyz = y.cross(z).unit();
53 G4ThreeVector nzx = z.cross(x).unit();
54
55 CLHEP::HepMatrix A(3, 3);
56 A[0][0] = nxy.x(), A[0][1] = nxy.y(), A[0][2] = nxy.z();
57 A[1][0] = nyz.x(), A[1][1] = nyz.y(), A[1][2] = nyz.z();
58 A[2][0] = nzx.x(), A[2][1] = nzx.y(), A[2][2] = nzx.z();
59
60 CLHEP::HepVector B(3);
61 B[0] = v0 * nxy - d;
62 B[1] = v0 * nyz - d;
63 B[2] = v0 * nzx - d;
64 // int ierr;
65 // cout<<A<<" "<<B<<" "<<A.inverse(ierr)<<endl;
66 CLHEP::HepVector r = A.inverse() * B; //CLHEP::solve(A, B);
67
68 return G4ThreeVector(r[0], r[1], r[2]);
69 }
70
71 G4ThreeVector moveto(double r, double phi)
72 {
73 return G4ThreeVector(r * cosd(phi), r * sind(phi), 0);
74 };
75
76
77 G4VSolid* shape_t::get_solid(const string& prefix, double wrapthick, G4Translate3D& shift) const
78 {
79 return get_bellecrystal(prefix, wrapthick, shift);
80 }
81
82 G4ThreeVector centerofgravity(const map<int, G4ThreeVector>& v, int i0, int n)
83 {
84 double cx = 0, cy = 0, A = 0;
85 for (int j = 0; j < n; j++) {
86 int j0 = j + i0, j1 = ((j + 1) % n) + i0;
87 const G4ThreeVector& v0 = v.find(j0)->second, &v1 = v.find(j1)->second;
88 double t = v0.x() * v1.y() - v0.y() * v1.x();
89 // cout<<j0<<" "<<j1<<" "<<t<<" "<<v0.x()<<" "<<v0.y()<<endl;
90 A += t;
91 cx += (v0.x() + v1.x()) * t;
92 cy += (v0.y() + v1.y()) * t;
93 }
94 cx /= 3 * A;
95 cy /= 3 * A;
96 return G4ThreeVector(cx, cy, v.find(i0)->second.z());
97 }
98
99 Point_t centerofgravity(const Point_t* b, const Point_t* e)
100 {
101 double cx = 0, cy = 0, A = 0;
102 int n = e - b;
103 for (int j = 0; j < n; j++) {
104 int j0 = j, j1 = ((j + 1) % n);
105 const Point_t& v0 = b[j0], &v1 = b[j1];
106 double t = v0.x * v1.y - v0.y * v1.x;
107 A += t;
108 cx += (v0.x + v1.x) * t;
109 cy += (v0.y + v1.y) * t;
110 }
111 cx /= 3 * A;
112 cy /= 3 * A;
113 Point_t t = {cx, cy};
114 return t;
115 }
116
118 struct quadrilateral_t: public shape_t {
119 quadrilateral_t() {}
120 virtual ~quadrilateral_t() override {}
122 virtual map<int, G4ThreeVector> make_verticies(double wrapthick) const = 0;
123
125 G4VSolid* get_tesselatedsolid(const string& prefix, double wrapthick, G4Translate3D& shift UNUSED) const override
126 {
127 map<int, G4ThreeVector> v = make_verticies(wrapthick);
128
129 string name(prefix);
130 name += to_string(nshape);
131 G4TessellatedSolid* s = new G4TessellatedSolid(name.c_str());
132
133 // Now add the facets to the solid
134 // top plane
135 s->AddFacet(new G4QuadrangularFacet(v[1], v[4], v[3], v[2], ABSOLUTE));
136 // bottom plane
137 s->AddFacet(new G4QuadrangularFacet(v[5], v[6], v[7], v[8], ABSOLUTE));
138 // lateral sides
139 s->AddFacet(new G4QuadrangularFacet(v[1], v[2], v[6], v[5], ABSOLUTE));
140 s->AddFacet(new G4QuadrangularFacet(v[2], v[3], v[7], v[6], ABSOLUTE));
141 s->AddFacet(new G4QuadrangularFacet(v[3], v[4], v[8], v[7], ABSOLUTE));
142 s->AddFacet(new G4QuadrangularFacet(v[4], v[1], v[5], v[8], ABSOLUTE));
143
144 // Finally declare the solid is complete
145 s->SetSolidClosed(true);
146 return s;
147 }
148
150 G4VSolid* get_extrudedsolid(const string& prefix, double wrapthick, G4Translate3D& shift UNUSED) const override
151 {
152 map<int, G4ThreeVector> v = make_verticies(wrapthick);
153
154 string name(prefix);
155 name += to_string(nshape);
156
157 std::vector<G4TwoVector> p1, p2;
158 p1.push_back(G4TwoVector(v[1].x(), v[1].y()));
159 p1.push_back(G4TwoVector(v[2].x(), v[2].y()));
160 p1.push_back(G4TwoVector(v[3].x(), v[3].y()));
161 p1.push_back(G4TwoVector(v[4].x(), v[4].y()));
162
163 p2.push_back(G4TwoVector(v[1 + 4].x(), v[1 + 4].y()));
164 p2.push_back(G4TwoVector(v[2 + 4].x(), v[2 + 4].y()));
165 p2.push_back(G4TwoVector(v[3 + 4].x(), v[3 + 4].y()));
166 p2.push_back(G4TwoVector(v[4 + 4].x(), v[4 + 4].y()));
167
168 double sum = 0, smin = 1e9, smax = -1e9;
169 for (int i = 0; i < 4; i++) {
170 for (int j = i + 1; j < 4; j++) {
171 double s2 = (p2[j] - p2[i]).mag2() / (p1[j] - p1[i]).mag2();
172 double s = sqrt(s2);
173 sum += s;
174 if (s > smax) smax = s;
175 if (s < smin) smin = s;
176 }
177 }
178 double ave = sum / 6;
179
180 double scale = ave;
181 G4TwoVector off1(0, 0), off2(scale * p1[0].x() - p2[0].x(), scale * p1[0].y() - p2[0].y());
182
183 return new G4ExtrudedSolid(name, p1, abs(v[1].z()), off1, 1, -off2, scale);
184 }
185
187 G4VSolid* get_trapezoid(const string& prefix, double wrapthick, G4Translate3D& shift) const override
188 {
189 map<int, G4ThreeVector> v = make_verticies(wrapthick);
190
191 // make sides @ +-Y parallel
192 G4ThreeVector b = v[4] - v[1];
193 G4ThreeVector d12 = v[2] - v[1], d43 = v[3] - v[4];
194 double b2 = b.mag2();
195 double h1 = b.cross(d12).mag2();
196 double h4 = b.cross(d43).mag2();
197 // if(abs(h1/h4-1)<0.02) cout<<"Shape = "<<nshape<<" "<<h1<<" "<<h4<<" "<<h1/h4-1<<endl;
198 if (h1 < h4) {
199 G4ThreeVector d24 = v[4] - v[2];
200 double d43b = d43 * b, s = (b2 * (d24 * d43) - (d24 * b) * d43b) / (d43b * d43b - b2 * d43.mag2());
201 v[3] = v[4] + s * d43;
202 v[7] = v[8] + s * (v[7] - v[8]);
203 } else {
204 G4ThreeVector d31 = v[1] - v[3];
205 double d12b = d12 * b, s = (b2 * (d31 * d12) - (d31 * b) * d12b) / (d12b * d12b - b2 * d12.mag2());
206 v[2] = v[1] + s * d12;
207 v[6] = v[5] + s * (v[6] - v[5]);
208 }
209
210 string name(prefix);
211 name += to_string(nshape);
212
213 G4ThreeVector pt[8];
214 pt[0] = v[4];
215 pt[1] = v[1];
216 pt[2] = v[3];
217 pt[3] = v[2];
218 pt[4] = v[8];
219 pt[5] = v[5];
220 pt[6] = v[7];
221 pt[7] = v[6];
222
223 auto alignz = [&](int i, int j) {pt[j].setZ(pt[i].z());};
224 auto aligny = [&](int i, int j) {pt[j].setY(pt[i].y());};
225
226 alignz(0, 1);
227 alignz(0, 2);
228 alignz(0, 3);
229
230 alignz(4, 1 + 4);
231 alignz(4, 2 + 4);
232 alignz(4, 3 + 4);
233
234 aligny(0, 1);
235 aligny(2, 3);
236
237 aligny(4, 5);
238 aligny(6, 7);
239
240 double dx = pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() + pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x();
241 dx /= 8;
242 double dy = pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y();
243 dy /= 4;
244 for (int j = 0; j < 8; j++) {
245 pt[j].setX(pt[j].x() - dx);
246 pt[j].setY(pt[j].y() - dy);
247 }
248
249 // int oprec = cout.precision(17);
250 // cout<<dx<<" "<<dy<<endl;
251 // cout.precision(oprec);
252
253 shift = G4Translate3D(dx, dy, 0);
254 G4VSolid* shape = new G4Trap(name.c_str(), pt);
255 // cout<<name<<" "<<shape->GetCubicVolume()*1e-3<<" "<<Volume<<endl;
256 // G4VSolid *shape = new BelleCrystal(name.c_str(), 4, pt);
257 return shape;
258 }
259
261 G4VSolid* get_bellecrystal(const string& prefix, double wrapthick, G4Translate3D& shift UNUSED) const override
262 {
263 map<int, G4ThreeVector> v = make_verticies(wrapthick);
264
265 string name(prefix);
266 name += to_string(nshape);
267
268 G4ThreeVector pt[8];
269 pt[0] = v[4];
270 pt[1] = v[1];
271 pt[2] = v[3];
272 pt[3] = v[2];
273 pt[4] = v[8];
274 pt[5] = v[5];
275 pt[6] = v[7];
276 pt[7] = v[6];
277
278 for (int i = 0; i < 8; i++) pt[i] = v[i + 1];
279 G4VSolid* shape = new BelleCrystal(name.c_str(), 4, pt);
280 // cout<<name<<" "<<shape->GetCubicVolume()*1e-3<<" "<<Volume<<endl;
281 // for(int i=0;i<100000;i++){
282 // G4ThreeVector a = shape->GetPointOnSurface();
283 // cout<<a.x()<<" "<<a.y()<<" "<<a.z()<<endl;
284 // }
285 // exit(0);
286 return shape;
287 }
288 };
289
291 struct quadrilateral_barrel_t: public quadrilateral_t {
292 union {
293 struct {
294 double A, B, H, a, b, h, alpha, beta, betap, gamma, Volume, Weight;
295 };
296 double t[12] = {};
297 };
298 quadrilateral_barrel_t() {}
299 virtual ~quadrilateral_barrel_t() override {}
300
302 bool istrap() const override
303 {
304 return true;
305 }
306
308 map<int, G4ThreeVector> make_verticies(double wrapthick) const override
309 {
310 map<int, G4ThreeVector> v;
311 // ensure sides to be parallel
312 double wh = h * h, wH = H * H, wnorm = wh + wH;
313 double tn = ((b - a) / 2 / h * wh + (B - A) / 2 / H * wH) / wnorm, d = h * tn, D = H * tn;
314 // double tn = tan(alpha*M_PI/180), d = h*tn, D = H*tn;
315 double m = (a + b) * 0.5, M = (A + B) * 0.5;
316
317 const double eps = 0.5e-3; // crystal sides are defined with 0.5 micron precision
318 if (fabs(a - (m - d)) > eps || fabs(b - (m + d)) > eps || fabs(A - (M - D)) > eps || fabs(B - (M + D)) > eps) {
319 double alfa = atan(tn) * 180 / M_PI;
320 B2WARNING("Cannot make parallel sides better than 0.5 mcm: alpha =" << alpha << " alpha from sides = " << alfa << " da = " << a -
321 (m - d) << " db = " << b - (m + d) << " dA = " << A - (M - D) << " dB = " << B - (M + D));
322 }
323
324 v[1] = G4ThreeVector(-(m - d) / 2, h / 2, -150);
325 v[2] = G4ThreeVector(-(m + d) / 2, -h / 2, -150);
326 v[3] = G4ThreeVector((m + d) / 2, -h / 2, -150);
327 v[4] = G4ThreeVector((m - d) / 2, h / 2, -150);
328 v[5] = G4ThreeVector(-(M - D) / 2, H / 2, 150);
329 v[6] = G4ThreeVector(-(M + D) / 2, -H / 2, 150);
330 v[7] = G4ThreeVector((M + D) / 2, -H / 2, 150);
331 v[8] = G4ThreeVector((M - D) / 2, H / 2, 150);
332
333 if (wrapthick != 0) {
334 map<int, G4ThreeVector> nv;
335 nv[1] = newvertex(wrapthick, v[1], v[5], v[2], v[4]);
336 nv[2] = newvertex(wrapthick, v[2], v[6], v[3], v[1]);
337 nv[3] = newvertex(wrapthick, v[3], v[7], v[4], v[2]);
338 nv[4] = newvertex(wrapthick, v[4], v[8], v[1], v[3]);
339 nv[5] = newvertex(wrapthick, v[5], v[1], v[8], v[6]);
340 nv[6] = newvertex(wrapthick, v[6], v[2], v[5], v[7]);
341 nv[7] = newvertex(wrapthick, v[7], v[3], v[6], v[8]);
342 nv[8] = newvertex(wrapthick, v[8], v[4], v[7], v[5]);
343 std::swap(nv, v);
344 }
345
346 // if(nshape==1){ for(int j=1;j<=8;j++) cout<<v[j]<<" "; cout<<endl;}
347
348 return v;
349 }
350 };
351
353 struct quadrilateral_endcap_t: public quadrilateral_t {
354 union {
355 struct {
356 double A, B, C, D, a, b, c, d, H_aA, H_dD, dg13, dg24, dg57, dg68, a1, a2, a3, a4, Volume, Weight;
357 };
358 double t[20] = {};
359 };
360 quadrilateral_endcap_t() {}
361 virtual ~quadrilateral_endcap_t() override {}
362
364 bool istrap() const override
365 {
366 double h1 = sind(a1) * D, h4 = sind(a4) * C;
367 return abs(h1 - h4) < 0.01 * h1;
368 }
369
371 map<int, G4ThreeVector> make_verticies(double wrapthick) const override
372 {
373 double minh = std::min(sind(a1) * D, sind(a4) * C);
374 double h2 = minh / 2;
375 double db = (tand(a1 - 90) - tand(a4 - 90)) * h2 / 2;
376
377 map<int, G4ThreeVector> v;
378 v[5] = G4ThreeVector(-A / 2 + db, -h2, -150);
379 v[6] = v[5] + moveto(D, a1);
380 v[8] = G4ThreeVector(A / 2 + db, -h2, -150);
381 v[7] = v[8] + moveto(C, 180 - a4);
382
383 // adjust position of v[2],v[3],v[7],v[6] to have all 4 points in a plane, other combination already in a plane by the construction
384 G4ThreeVector vB = v[7] - v[6], va(a, 0, 0), vd = moveto(d, a1), vc = moveto(c, 180 - a4);
385 double delta = vB.cross(va + vc - vd).z() / vB.cross(vc + vd).z(); // delta should be very small ~10^-6 or less
386
387 vd *= 1 + delta;
388 vc *= 1 - delta;
389
390 v[1] = v[5] + G4ThreeVector((H_aA * cosd(a1) + H_dD) / sind(a1), H_aA, 300);
391 v[2] = v[1] + vd;
392 v[4] = v[1] + va;
393 v[3] = v[4] + vc;
394
395 for (int j = 1; j <= 8; j++) v[j] = G4ThreeVector(v[j].x(), -v[j].y(), -v[j].z());
396
397 // G4ThreeVector c0 = centerofgravity(v, 1, 4);
398 // G4ThreeVector c1 = centerofgravity(v, 5, 4);
399 // G4ThreeVector cz = 0.5*(c0+c1);
400 // cout<<c0<<" "<<c1<<" "<<cz<<endl;
401
402 if (wrapthick != 0) {
403 map<int, G4ThreeVector> nv;
404 nv[1] = newvertex(wrapthick, v[1], v[5], v[2], v[4]);
405 nv[2] = newvertex(wrapthick, v[2], v[6], v[3], v[1]);
406 nv[3] = newvertex(wrapthick, v[3], v[7], v[4], v[2]);
407 nv[4] = newvertex(wrapthick, v[4], v[8], v[1], v[3]);
408 nv[5] = newvertex(wrapthick, v[5], v[1], v[8], v[6]);
409 nv[6] = newvertex(wrapthick, v[6], v[2], v[5], v[7]);
410 nv[7] = newvertex(wrapthick, v[7], v[3], v[6], v[8]);
411 nv[8] = newvertex(wrapthick, v[8], v[4], v[7], v[5]);
412 std::swap(nv, v);
413 }
414 // if(nshape==2){ for(int j=1;j<=8;j++) cout<<v[j]<<" "; cout<<endl;}
415 return v;
416 }
417 };
418
420 struct pent_t: public shape_t {
421 union {
422 struct {
423 double A, C, D, a, c, d, B, b, H_aA, H_dD, dg13, dg24, dg57, dg68, a1, a4, a2, a3, a9, Volume, Weight;
424 };
425 double t[21] = {};
426 };
427 bool _adjusted;
428 pent_t(): _adjusted(false) {}
429 virtual ~pent_t() override {}
430
432 void adjust()
433 {
434 if (!_adjusted) {
435 // adjust sizes to have flat sides
436 H_dD -= 0.00005551197484235;
437 B += 0.0011245236213532729;
438 b += -0.00044853029662963;
439 _adjusted = true;
440 }
441 }
442
444 bool istrap() const override { return false;}
445
447 map<int, G4ThreeVector> make_verticies(double wrapthick) const
448 {
449 assert(_adjusted);
450
451 double h2 = cosd(a1 - 90) * D / 2;
452 map<int, G4ThreeVector> v;
453 v[5] = G4ThreeVector(-A / 2, -h2, -150);
454 v[6] = v[5] + moveto(D, a1);
455 v[10] = v[6] + moveto(B, a1 + a2 - 180);
456 v[8] = G4ThreeVector(A / 2, -h2, -150);
457 v[7] = v[8] + moveto(D, 180 - a1);
458
459 v[1] = v[5] + G4ThreeVector((H_aA * cosd(a1) + H_dD) / sind(a1), H_aA, 300);
460 v[2] = v[1] + moveto(d, a1);
461 v[9] = v[2] + moveto(b, a1 + a2 - 180);
462 v[4] = v[1] + moveto(a, 0);
463 v[3] = v[4] + moveto(d, 180 - a1);
464
465 for (int j = 1; j <= 10; j++) v[j] = G4ThreeVector(v[j].x(), -v[j].y(), -v[j].z());
466 if (wrapthick != 0) {
467 map<int, G4ThreeVector> nv;
468 nv[ 1] = newvertex(wrapthick, v[1], v[5], v[2], v[4]);
469 nv[ 2] = newvertex(wrapthick, v[2], v[6], v[9], v[1]);
470 nv[ 9] = newvertex(wrapthick, v[9], v[10], v[3], v[2]);
471 nv[ 3] = newvertex(wrapthick, v[3], v[7], v[4], v[9]);
472 nv[ 4] = newvertex(wrapthick, v[4], v[8], v[1], v[3]);
473 nv[ 5] = newvertex(wrapthick, v[5], v[1], v[8], v[6]);
474 nv[ 6] = newvertex(wrapthick, v[6], v[2], v[5], v[10]);
475 nv[10] = newvertex(wrapthick, v[10], v[9], v[6], v[7]);
476 nv[ 7] = newvertex(wrapthick, v[7], v[3], v[10], v[8]);
477 nv[ 8] = newvertex(wrapthick, v[8], v[4], v[7], v[5]);
478 std::swap(nv, v);
479 }
480 return v;
481 }
482
484 G4VSolid* get_tesselatedsolid(const string& prefix, double wrapthick, G4Translate3D& shift UNUSED) const override
485 {
486 if (nshape != 36) return nullptr; // only one crystal has pentagon shape
487
488 map<int, G4ThreeVector> v = make_verticies(wrapthick);
489
490 string name(prefix);
491 name += to_string(nshape);
492 G4TessellatedSolid* s = new G4TessellatedSolid(name.c_str());
493
494 // Now add the facets to the solid
495 // top plane
496 s->AddFacet(new G4QuadrangularFacet(v[1], v[4], v[3], v[2], ABSOLUTE));
497 s->AddFacet(new G4TriangularFacet(v[2], v[3], v[9], ABSOLUTE));
498
499 //bottom plane
500 s->AddFacet(new G4QuadrangularFacet(v[5], v[6], v[7], v[8], ABSOLUTE));
501 s->AddFacet(new G4TriangularFacet(v[6], v[10], v[7], ABSOLUTE));
502
503 //sides
504 s->AddFacet(new G4QuadrangularFacet(v[1], v[2], v[6], v[5], ABSOLUTE));
505 s->AddFacet(new G4QuadrangularFacet(v[2], v[9], v[10], v[6], ABSOLUTE));
506 s->AddFacet(new G4QuadrangularFacet(v[9], v[3], v[7], v[10], ABSOLUTE));
507 s->AddFacet(new G4QuadrangularFacet(v[3], v[4], v[8], v[7], ABSOLUTE));
508 s->AddFacet(new G4QuadrangularFacet(v[4], v[1], v[5], v[8], ABSOLUTE));
509
510 // Finally declare the solid is complete
511 s->SetSolidClosed(true);
512 return s;
513 }
514
516 G4VSolid* get_extrudedsolid(const string& prefix, double wrapthick, G4Translate3D& shift UNUSED) const override
517 {
518 map<int, G4ThreeVector> v = make_verticies(wrapthick);
519
520 string name(prefix);
521 name += to_string(nshape);
522
523 std::vector<G4TwoVector> p1, p2;
524 p1.push_back(G4TwoVector(v[1].x(), v[1].y()));
525 p1.push_back(G4TwoVector(v[2].x(), v[2].y()));
526 p1.push_back(G4TwoVector(v[9].x(), v[9].y()));
527 p1.push_back(G4TwoVector(v[3].x(), v[3].y()));
528 p1.push_back(G4TwoVector(v[4].x(), v[4].y()));
529
530 p2.push_back(G4TwoVector(v[1 + 4].x(), v[1 + 4].y()));
531 p2.push_back(G4TwoVector(v[2 + 4].x(), v[2 + 4].y()));
532 p2.push_back(G4TwoVector(v[ 10].x(), v[ 10].y()));
533 p2.push_back(G4TwoVector(v[3 + 4].x(), v[3 + 4].y()));
534 p2.push_back(G4TwoVector(v[4 + 4].x(), v[4 + 4].y()));
535
536 double sum = 0, smin = 1e9, smax = -1e9;
537 int count = 0;
538 for (int i = 0; i < 5; i++) {
539 for (int j = i + 1; j < 5; j++) {
540 double s2 = (p2[j] - p2[i]).mag2() / (p1[j] - p1[i]).mag2();
541 double s = sqrt(s2);
542 sum += s;
543 if (s > smax) smax = s;
544 if (s < smin) smin = s;
545 count++;
546 }
547 }
548 double ave = sum / count;
549
550 double scale = ave;
551 G4TwoVector off1(0, 0), off2(scale * p1[0].x() - p2[0].x(), scale * p1[0].y() - p2[0].y());
552
553 return new G4ExtrudedSolid(name, p1, abs(v[1].z()), off1, 1, -off2, scale);
554 }
555
557 G4VSolid* get_trapezoid(const string& prefix, double wrapthick, G4Translate3D& shift) const override
558 {
559 if (nshape != 36) return nullptr; // only one crystal has pentagon shape
560
561 map<int, G4ThreeVector> v = make_verticies(wrapthick);
562
563 string name(prefix);
564 name += to_string(nshape);
565
566 G4ThreeVector pt[8];
567 pt[0] = v[4];
568 pt[1] = v[1];
569 pt[2] = v[3];
570 pt[3] = v[2];
571 pt[4] = v[8];
572 pt[5] = v[5];
573 pt[6] = v[7];
574 pt[7] = v[6];
575
576 auto alignz = [&](int i, int j) {pt[j].setZ(pt[i].z());};
577 auto aligny = [&](int i, int j) {pt[j].setY(pt[i].y());};
578
579 alignz(0, 1);
580 alignz(0, 2);
581 alignz(0, 3);
582
583 alignz(4, 1 + 4);
584 alignz(4, 2 + 4);
585 alignz(4, 3 + 4);
586
587 aligny(0, 1);
588 aligny(2, 3);
589
590 aligny(4, 5);
591 aligny(6, 7);
592
593
594 double dx = pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() + pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x();
595 dx /= 8;
596 double dy = pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y();
597 dy /= 4;
598 for (int j = 0; j < 8; j++) {
599 pt[j].setX(pt[j].x() - dx);
600 pt[j].setY(pt[j].y() - dy);
601 }
602 shift = G4Translate3D(dx, dy, 0);
603
604 // cout<<name<<" "<<dx<<" "<<dy<<endl;
605
606 G4VSolid* shape = new G4Trap(name.c_str(), pt);
607 // G4VSolid *shape = new BelleCrystal(name.c_str(), 4, pt);
608 return shape;
609 }
610
612 G4VSolid* get_bellecrystal(const string& prefix, double wrapthick, G4Translate3D& shift UNUSED) const override
613 {
614 map<int, G4ThreeVector> v = make_verticies(wrapthick);
615
616 string name(prefix);
617 name += to_string(nshape);
618
619 G4ThreeVector pt[10];
620 pt[0] = v[1];
621 pt[1] = v[2];
622 pt[2] = v[9];
623 pt[3] = v[3];
624 pt[4] = v[4];
625 pt[5] = v[5];
626 pt[6] = v[6];
627 pt[7] = v[10];
628 pt[8] = v[7];
629 pt[9] = v[8];
630
631 G4VSolid* shape = new BelleCrystal(name.c_str(), 5, pt);
632 // cout<<name<<" "<<shape->GetCubicVolume()*1e-3<<" "<<Volume<<endl;
633 // for(int i=0;i<100000;i++){
634 // G4ThreeVector a = shape->GetPointOnSurface();
635 // cout<<a.x()<<" "<<a.y()<<" "<<a.z()<<endl;
636 // }
637 // exit(0);
638 // if(prefix.find("crystal")!=string::npos){
639 // pentspeed(shape);
640 // exit(0);
641 // }
642 // if(prefix.find("wrap")!=string::npos){
643 // pentspeed2(shape);
644 // exit(0);
645 // }
646 return shape;
647 }
648 };
649
650 vector<shape_t*> load_shapes(const string& fname)
651 {
652 vector<shape_t*> shapes;
653 std::string fnamef = Belle2::FileSystem::findFile(fname);
654
655 ifstream IN(fnamef.c_str());
656 string tmp;
657 while (getline(IN, tmp)) {
658 size_t ic = tmp.find("#");
659 if (ic != string::npos) tmp.erase(ic);
660 istringstream iss(tmp);
661 vector<string> t;
662 copy(istream_iterator<string>(iss), istream_iterator<string>(), back_inserter(t));
663 if (t.size() > 0) {
664 shape_t* shape = nullptr;
665 if (t.size() == 21) {
666 shape = new quadrilateral_endcap_t();
667 quadrilateral_endcap_t& trap = static_cast<quadrilateral_endcap_t&>(*shape);
668
669 istringstream in(t[0]);
670 in >> trap.nshape;
671
672 for (size_t i = 1; i < t.size(); i++) {
673 in.str(t[i]); in.seekg(0, ios_base::beg);
674 in >> trap.t[i - 1];
675 }
676 } else if (t.size() == 13) {
677 shape = new quadrilateral_barrel_t();
678 quadrilateral_barrel_t& trap = static_cast<quadrilateral_barrel_t&>(*shape);
679
680 istringstream in(t[0]);
681 in >> trap.nshape;
682
683 for (size_t i = 1; i < t.size(); i++) {
684 in.str(t[i]); in.seekg(0, ios_base::beg);
685 in >> trap.t[i - 1];
686 }
687 } else if (t.size() == 22) {
688 shape = new pent_t();
689 pent_t& pent = static_cast<pent_t&>(*shape);
690
691 istringstream in(t[0]);
692 in >> pent.nshape;
693
694 for (size_t i = 1; i < t.size(); i++) {
695 in.str(t[i]); in.seekg(0, ios_base::beg);
696 in >> pent.t[i - 1];
697 }
698 pent.adjust();
699 }
700 shapes.push_back(shape);
701 }
702 }
703
704 return shapes;
705
706 // for(unsigned int i=0;i<shapes.size();i++){
707 // shape_t &t = *shapes[i];
708 // crystals[t.nshape] = t.get_solid(prefix, 0.250-0.006);
709 // }
710
711 // return crystals;
712 }
713
714 vector<cplacement_t> load_placements(const string& fname)
715 {
716 vector<cplacement_t> plcmnt;
717 std::string fnamef = Belle2::FileSystem::findFile(fname);
718
719 ifstream IN(fnamef.c_str());
720 string tmp;
721 while (getline(IN, tmp)) {
722 size_t ic = tmp.find("#");
723 if (ic != string::npos) tmp.erase(ic);
724 istringstream iss(tmp);
725 vector<string> t;
726 copy(istream_iterator<string>(iss), istream_iterator<string>(), back_inserter(t));
727 if (t.size() == 7) {
728 cplacement_t p;
729 istringstream in(t[0]);
730 in >> p.nshape;
731 in.str(t[1]); in.seekg(0, ios_base::beg);
732 in >> p.Rphi1;
733 in.str(t[2]); in.seekg(0, ios_base::beg);
734 in >> p.Rtheta;
735 in.str(t[3]); in.seekg(0, ios_base::beg);
736 in >> p.Rphi2;
737 in.str(t[4]); in.seekg(0, ios_base::beg);
738 in >> p.Pr;
739 in.str(t[5]); in.seekg(0, ios_base::beg);
740 in >> p.Ptheta;
741 in.str(t[6]); in.seekg(0, ios_base::beg);
742 in >> p.Pphi;
743 plcmnt.push_back(p);
744 }
745 }
746
747 return plcmnt;
748 }
749
750 G4Transform3D get_transform(const cplacement_t& t)
751 {
752 G4Transform3D r = G4Rotate3D(t.Rphi1, G4Vector3D(sin(t.Rtheta) * cos(t.Rphi2), sin(t.Rtheta) * sin(t.Rphi2), cos(t.Rtheta)));
753 G4Transform3D p = G4Translate3D(t.Pr * sin(t.Ptheta) * cos(t.Pphi), t.Pr * sin(t.Ptheta) * sin(t.Pphi), t.Pr * cos(t.Ptheta));
754 return p * r;
755 }
756
757 Belle2::ECLCrystalsShapeAndPosition loadCrystalsShapeAndPosition()
758 {
759 stringstream buffer;
760 auto fillbuffer = [&buffer](const string & fname) {
761 string path = Belle2::FileSystem::findFile(fname);
762 ifstream IN(path.c_str());
763 buffer.clear(); buffer.str("");
764 buffer << IN.rdbuf();
765 };
766
767 Belle2::ECLCrystalsShapeAndPosition a;
768 fillbuffer("/ecl/data/crystal_shape_forward.dat"); a.setShapeForward(buffer.str());
769 fillbuffer("/ecl/data/crystal_shape_barrel.dat"); a.setShapeBarrel(buffer.str());
770 fillbuffer("/ecl/data/crystal_shape_backward.dat"); a.setShapeBackward(buffer.str());
771 fillbuffer("/ecl/data/crystal_placement_forward.dat"); a.setPlacementForward(buffer.str());
772 fillbuffer("/ecl/data/crystal_placement_barrel.dat"); a.setPlacementBarrel(buffer.str());
773 fillbuffer("/ecl/data/crystal_placement_backward.dat"); a.setPlacementBackward(buffer.str());
774 return a;
775 // Belle2::IntervalOfValidity iov(0, 0, -1, -1); // IOV (0,0,-1,-1) is valid for all runs and experiments
776 // Belle2::Database::Instance().storeData<Belle2::ECLCrystalsShapeAndPosition>(&a, iov);
777 }
778
779 vector<shape_t*> load_shapes(stringstream& IN)
780 {
781 vector<shape_t*> shapes;
782 string tmp;
783 while (getline(IN, tmp)) {
784 size_t ic = tmp.find("#");
785 if (ic != string::npos) tmp.erase(ic);
786 istringstream iss(tmp);
787 vector<string> t;
788 copy(istream_iterator<string>(iss), istream_iterator<string>(), back_inserter(t));
789 if (t.size() > 0) {
790 shape_t* shape = nullptr;
791 if (t.size() == 21) {
792 shape = new quadrilateral_endcap_t();
793 quadrilateral_endcap_t& trap = static_cast<quadrilateral_endcap_t&>(*shape);
794
795 istringstream in(t[0]);
796 in >> trap.nshape;
797
798 for (size_t i = 1; i < t.size(); i++) {
799 in.str(t[i]); in.seekg(0, ios_base::beg);
800 in >> trap.t[i - 1];
801 }
802 } else if (t.size() == 13) {
803 shape = new quadrilateral_barrel_t();
804 quadrilateral_barrel_t& trap = static_cast<quadrilateral_barrel_t&>(*shape);
805
806 istringstream in(t[0]);
807 in >> trap.nshape;
808
809 for (size_t i = 1; i < t.size(); i++) {
810 in.str(t[i]); in.seekg(0, ios_base::beg);
811 in >> trap.t[i - 1];
812 }
813 } else if (t.size() == 22) {
814 shape = new pent_t();
815 pent_t& pent = static_cast<pent_t&>(*shape);
816
817 istringstream in(t[0]);
818 in >> pent.nshape;
819
820 for (size_t i = 1; i < t.size(); i++) {
821 in.str(t[i]); in.seekg(0, ios_base::beg);
822 in >> pent.t[i - 1];
823 }
824 pent.adjust();
825 }
826 shapes.push_back(shape);
827 }
828 }
829
830 return shapes;
831 }
832
833 vector<cplacement_t> load_placements(stringstream& IN)
834 {
835 vector<cplacement_t> plcmnt;
836 string tmp;
837 while (getline(IN, tmp)) {
838 size_t ic = tmp.find("#");
839 if (ic != string::npos) tmp.erase(ic);
840 istringstream iss(tmp);
841 vector<string> t;
842 copy(istream_iterator<string>(iss), istream_iterator<string>(), back_inserter(t));
843 if (t.size() == 7) {
844 cplacement_t p;
845 istringstream in(t[0]);
846 in >> p.nshape;
847 in.str(t[1]); in.seekg(0, ios_base::beg);
848 in >> p.Rphi1;
849 in.str(t[2]); in.seekg(0, ios_base::beg);
850 in >> p.Rtheta;
851 in.str(t[3]); in.seekg(0, ios_base::beg);
852 in >> p.Rphi2;
853 in.str(t[4]); in.seekg(0, ios_base::beg);
854 in >> p.Pr;
855 in.str(t[5]); in.seekg(0, ios_base::beg);
856 in >> p.Ptheta;
857 in.str(t[6]); in.seekg(0, ios_base::beg);
858 in >> p.Pphi;
859 plcmnt.push_back(p);
860 }
861 }
862
863 return plcmnt;
864 }
865
866 vector<cplacement_t> load_placements(const Belle2::ECLCrystalsShapeAndPosition* crystals, enum ECLParts part)
867 {
868 // Belle2::DBObjPtr<Belle2::ECLCrystalsShapeAndPosition> crystals;
869 // if (!crystals.isValid()) B2FATAL("No crystal's data in the database.");
870
871 // stringstream buffer;
872 // auto fillbuffer = [&buffer](const string &fname) {
873 // string path = Belle2::FileSystem::findFile(fname);
874 // ifstream IN(path.c_str());
875 // buffer.clear(); buffer.str("");
876 // buffer << IN.rdbuf();
877 // };
878
879 // Belle2::ECLCrystalsShapeAndPosition *crystals = new Belle2::ECLCrystalsShapeAndPosition();
880 // fillbuffer("/ecl/data/crystal_shape_forward.dat"); crystals->setShapeForward(buffer.str());
881 // fillbuffer("/ecl/data/crystal_shape_barrel.dat"); crystals->setShapeBarrel(buffer.str());
882 // fillbuffer("/ecl/data/crystal_shape_backward.dat"); crystals->setShapeBackward(buffer.str());
883 // fillbuffer("/ecl/data/crystal_placement_forward.dat"); crystals->setPlacementForward(buffer.str());
884 // fillbuffer("/ecl/data/crystal_placement_barrel.dat"); crystals->setPlacementBarrel(buffer.str());
885 // fillbuffer("/ecl/data/crystal_placement_backward.dat"); crystals->setPlacementBackward(buffer.str());
886
887 stringstream IN;
888 if (part == ECLParts::forward)
889 IN.str(crystals->getPlacementForward());
890 else if (part == ECLParts::barrel)
891 IN.str(crystals->getPlacementBarrel());
892 else if (part == ECLParts::backward)
893 IN.str(crystals->getPlacementBackward());
894 return load_placements(IN);
895 }
896
897 vector<shape_t*> load_shapes(const Belle2::ECLCrystalsShapeAndPosition* crystals, enum ECLParts part)
898 {
899 stringstream IN;
900 if (part == ECLParts::forward)
901 IN.str(crystals->getShapeForward());
902 else if (part == ECLParts::barrel)
903 IN.str(crystals->getShapeBarrel());
904 else if (part == ECLParts::backward)
905 IN.str(crystals->getShapeBackward());
906 return load_shapes(IN);
907 }
908
909 void testtest()
910 {
911 Belle2::ECLCrystalsShapeAndPosition crystals = loadCrystalsShapeAndPosition();
912 vector<cplacement_t> bp = load_placements(&crystals, ECLParts::forward);
913 for (vector<cplacement_t>::const_iterator it = bp.begin(); it != bp.end(); ++it) {
914 const cplacement_t& t = *it;
915 cout << t.nshape << " " << t.Rphi1 << " " << t.Rtheta << " " << t.Rphi2 << " " << t.Pr << " " << t.Ptheta << " " << t.Pphi << endl;
916 }
917 exit(0);
918 }
919 }
921}
Belle2::ECLCrystalsShapeAndPosition::getPlacementBarrel
const std::string & getPlacementBarrel() const
Return crystal placement in barrel.
Definition ECLCrystalsShapeAndPosition.h:78
Belle2::ECLCrystalsShapeAndPosition::getShapeBackward
const std::string & getShapeBackward() const
Return crystal shape in backward endcap.
Definition ECLCrystalsShapeAndPosition.h:58
Belle2::ECLCrystalsShapeAndPosition::getShapeForward
const std::string & getShapeForward() const
Return crystal shape in forward endcap.
Definition ECLCrystalsShapeAndPosition.h:38
Belle2::ECLCrystalsShapeAndPosition::getShapeBarrel
const std::string & getShapeBarrel() const
Return crystal shape in barrel.
Definition ECLCrystalsShapeAndPosition.h:48
Belle2::ECLCrystalsShapeAndPosition::getPlacementBackward
const std::string & getPlacementBackward() const
Return crystal placement in backward endcap.
Definition ECLCrystalsShapeAndPosition.h:88
Belle2::ECLCrystalsShapeAndPosition::getPlacementForward
const std::string & getPlacementForward() const
Return crystal placement in forward endcap.
Definition ECLCrystalsShapeAndPosition.h:68
Belle2::ECL::BelleCrystal
a Belle crystal in Geant4
Definition BelleCrystal.h:37
Belle2::FileSystem::findFile
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
Definition FileSystem.cc:159
Belle2::tan
double tan(double a)
tan for double
Definition beamHelpers.h:31
Belle2::atan
double atan(double a)
atan for double
Definition beamHelpers.h:34
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17
std
STL namespace.
Belle2::ECL::Point_t
struct for Point
Definition BelleCrystal.h:31
Belle2::ECL::cplacement_t
placement struct
Definition shapes.h:49
Belle2::ECL::pent_t
pentagon shape
Definition shapes.cc:420
Belle2::ECL::pent_t::get_tesselatedsolid
G4VSolid * get_tesselatedsolid(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get tessellated solid
Definition shapes.cc:484
Belle2::ECL::pent_t::istrap
bool istrap() const override
is trapped?
Definition shapes.cc:444
Belle2::ECL::pent_t::get_extrudedsolid
G4VSolid * get_extrudedsolid(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get extruded solid
Definition shapes.cc:516
Belle2::ECL::pent_t::get_bellecrystal
G4VSolid * get_bellecrystal(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get Belle crystal
Definition shapes.cc:612
Belle2::ECL::pent_t::adjust
void adjust()
adjust sizes to have flat sides
Definition shapes.cc:432
Belle2::ECL::pent_t::make_verticies
map< int, G4ThreeVector > make_verticies(double wrapthick) const
create map of vertices
Definition shapes.cc:447
Belle2::ECL::pent_t::get_trapezoid
G4VSolid * get_trapezoid(const string &prefix, double wrapthick, G4Translate3D &shift) const override
get trapezoid
Definition shapes.cc:557
Belle2::ECL::pent_t::_adjusted
bool _adjusted
are sizes adjusted?
Definition shapes.cc:427
Belle2::ECL::quadrilateral_barrel_t
quadrilateral struct for barrel
Definition shapes.cc:291
Belle2::ECL::quadrilateral_barrel_t::istrap
bool istrap() const override
is trapped
Definition shapes.cc:302
Belle2::ECL::quadrilateral_barrel_t::make_verticies
map< int, G4ThreeVector > make_verticies(double wrapthick) const override
create map of vertices
Definition shapes.cc:308
Belle2::ECL::quadrilateral_endcap_t
quadrilateral struct for end cap
Definition shapes.cc:353
Belle2::ECL::quadrilateral_endcap_t::istrap
bool istrap() const override
is trapped?
Definition shapes.cc:364
Belle2::ECL::quadrilateral_endcap_t::make_verticies
map< int, G4ThreeVector > make_verticies(double wrapthick) const override
create map of vertices
Definition shapes.cc:371
Belle2::ECL::quadrilateral_t::get_tesselatedsolid
G4VSolid * get_tesselatedsolid(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get tessellated solid
Definition shapes.cc:125
Belle2::ECL::quadrilateral_t::make_verticies
virtual map< int, G4ThreeVector > make_verticies(double wrapthick) const =0
create map of vertices
Belle2::ECL::quadrilateral_t::get_extrudedsolid
G4VSolid * get_extrudedsolid(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get extruded solid
Definition shapes.cc:150
Belle2::ECL::quadrilateral_t::get_bellecrystal
G4VSolid * get_bellecrystal(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get Belle crystal
Definition shapes.cc:261
Belle2::ECL::quadrilateral_t::get_trapezoid
G4VSolid * get_trapezoid(const string &prefix, double wrapthick, G4Translate3D &shift) const override
get trapezoid
Definition shapes.cc:187
Belle2::ECL::shape_t
shape
Definition shapes.h:29
Belle2::ECL::shape_t::nshape
int nshape
shapes
Definition shapes.h:30
Belle2::ECL::shape_t::get_solid
G4VSolid * get_solid(const std::string &prefix, double wrapthick, G4Translate3D &shift) const
get solid
Definition shapes.cc:77
Belle2::ECL::shape_t::get_bellecrystal
virtual G4VSolid * get_bellecrystal(const std::string &prefix, double wrapthick, G4Translate3D &shift) const =0
get Belle crystal