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