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