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