Belle II Software  release-08-01-10
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, 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::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: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:148
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::atan
double atan(double a)
atan for double
Definition: beamHelpers.h:34
Belle2::tan
double tan(double a)
tan for double
Definition: beamHelpers.h:31
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
pentagon shape
Definition: shapes.cc:416
Belle2::ECL::pent_t::istrap
bool istrap() const override
is trapped?
Definition: shapes.cc:440
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::get_tesselatedsolid
G4VSolid * get_tesselatedsolid(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get tesselated solid
Definition: shapes.cc:480
Belle2::ECL::pent_t::adjust
void adjust()
adjust sizes to have flat sides
Definition: shapes.cc:428
Belle2::ECL::pent_t::_adjusted
bool _adjusted
are sizes adjusted?
Definition: shapes.cc:423
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::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::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_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_tesselatedsolid
G4VSolid * get_tesselatedsolid(const string &prefix, double wrapthick, G4Translate3D &shift UNUSED) const override
get tesselated 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_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