Belle II Software  release-05-02-19
BelleCrystal.cc
1 #include "ecl/geometry/BelleCrystal.h"
2 
3 #include "globals.hh"
4 
5 #include "G4VoxelLimits.hh"
6 #include "G4AffineTransform.hh"
7 
8 #include "G4VPVParameterisation.hh"
9 
10 #include "G4VGraphicsScene.hh"
11 
12 #include "CLHEP/Random/RandFlat.h"
13 
14 using namespace Belle2;
15 using namespace std;
16 using namespace ECL;
17 
18 #define DO_QUOTE(X) #X
19 #define QUOTE(X) DO_QUOTE(X)
20 
21 #define PERFCOUNTER 0
22 //#define MATCHNAME sv_fwd_crystal2
23 
24 #if PERFCOUNTER==1
25 int counter[6];
26 # define COUNTER(x) counter[x]++
27 #else
28 # define COUNTER(x)
29 #endif
30 
31 #ifdef MATCHNAME
32 # define MATCHOUT(x) if(fmatch) cout<<x<<"\n";
33 #else
34 # define MATCHOUT(x) {}
35 #endif
36 
37 // Accuracy of coplanarity
38 const G4double kCoplanar_Tolerance = 1E-4;
39 
40 BelleCrystal::BelleCrystal(const G4String& pName, int n,
41  const G4ThreeVector* pt)
42  : G4CSGSolid(pName), nsides(n), fPlanes(nsides), fx(2 * nsides)
43 {
44 #if PERFCOUNTER==1
45  memset(fcounter, 0, sizeof(fcounter));
46 #endif
47 #ifdef MATCHNAME
48  fmatch = GetName() == QUOTE(MATCHNAME);
49  cout.precision(17);
50 #endif
51  // ref = new G4Trap("dummy",pt);
52  fDz = abs(pt[2 * nsides - 1].z());
53  for (unsigned int i = 0; i < 2 * nsides; i++) {
54  fx[i].x = pt[i].x();
55  fx[i].y = pt[i].y();
56  }
57  // cout<<pName<<" "<<fx.data()<<endl;
58  // for(int i=0; i<2*nsides + 2*(nsides-2); i++) ivertx(i);
59 
60  auto isConvex = [](std::vector<Point_t>::const_iterator begin, std::vector<Point_t>::const_iterator end) -> bool {
61  bool sign = false;
62  int np = end - begin;
63  for (int i0 = 0; i0 < np; i0++)
64  {
65  int i1 = (i0 + 1) % np, i2 = (i0 + 2) % np;
66  const Point_t& r2 = *(begin + i2), &r1 = *(begin + i1), &r0 = *(begin + i0);
67  double dx1 = r2.x - r1.x, dy1 = r2.y - r1.y;
68  double dx2 = r0.x - r1.x, dy2 = r0.y - r1.y;
69  double x = dx1 * dy2 - dy1 * dx2;
70  if (i0 == 0) sign = x > 0;
71  else {
72  if (sign != (x > 0)) return false;
73  }
74  }
75  return true;
76  };
77 
78  auto isClockwise = [](std::vector<Point_t>::const_iterator begin, std::vector<Point_t>::const_iterator end) -> bool {
79  std::vector<Point_t>::const_iterator it = begin;
80  Point_t r0 = *it++;
81  double sum = 0;
82  for (; it != end;)
83  {
84  Point_t r1 = *it++; sum += (r1.x - r0.x) * (r1.y + r0.y);
85  r0 = r1;
86  }
87  Point_t r1 = *begin; sum += (r1.x - r0.x) * (r1.y + r0.y);
88  return sum > 0;
89  };
90 
91  if (!isConvex(fx.begin(), fx.begin() + nsides)) {
92  std::ostringstream message;
93  message << "At -z polygon is not convex: " << GetName();
94  G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
95  }
96  if (!isConvex(fx.begin() + nsides, fx.end())) {
97  std::ostringstream message;
98  message << "At +z polygon is not convex: " << GetName();
99  G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
100  }
101 
102  if (isClockwise(fx.begin(), fx.begin() + nsides)) {
103  std::ostringstream message;
104  message << "At -z polygon is not in anti-clockwise order: " << GetName();
105  G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
106  }
107  if (isClockwise(fx.begin() + nsides, fx.end())) {
108  std::ostringstream message;
109  message << "At +z polygon is not in anti-clockwise order: " << GetName();
110  G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
111  }
112 
113  // sanity check
114  const G4ThreeVector* pm = pt, *pp = pt + nsides;
115  bool zm = pm[0].z() < 0, zp = pp[0].z() > 0;
116  for (unsigned int i = 1; i < nsides; i++) {
117  zm = zm && (pm[i - 1].z() - pm[i].z()) < kCarTolerance;
118  zp = zp && (pp[i - 1].z() - pp[i].z()) < kCarTolerance;
119  }
120  bool zpm = abs(pm[0].z() + pp[0].z()) < kCarTolerance;
121  if (!(zm && zp && zpm)) {
122  std::ostringstream message;
123  message << "Invalid vertice coordinates for Solid: " << GetName() << " " << zp << " " << zm << " " << zpm;
124  G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
125  }
126  if (nsides > 3) {
127  for (unsigned int i = 0; i < nsides; i++)
128  MakePlane(pt[i], pt[i + nsides], pt[((i + 1) % (nsides)) + nsides], pt[(i + 1) % nsides], fPlanes[i]);
129  } else {
130  std::ostringstream message;
131  message << "Wrong number of sides for Belle Crystal: " << GetName() << " nsides = " << nsides;
132  G4Exception("BelleCrystal::BelleCrystal()", "BelleCrystal", FatalException, message);
133  }
134 }
135 
136 // Nominal constructor for BelleCrystal
137 BelleCrystal::BelleCrystal(const G4String& pName)
138  : G4CSGSolid(pName), nsides(0), fDz(0)
139 {
140 }
141 
142 // Fake default constructor - sets only member data and allocates memory
143 // for usage restricted to object persistency.
144 BelleCrystal::BelleCrystal(__void__& a)
145  : G4CSGSolid(a), nsides(0), fDz(0)
146 {
147 }
148 
149 // Destructor
150 BelleCrystal::~BelleCrystal()
151 {
152 #if PERFCOUNTER==1
153  cout << GetName() << " ";
154  for (int i = 0; i < 6; i++) cout << counter[i] << " "; cout << endl;
155  exit(0);
156 #endif
157  // delete ref;
158 }
159 
160 // Copy constructor
161 BelleCrystal::BelleCrystal(const BelleCrystal& rhs)
162  : G4CSGSolid(rhs), nsides(rhs.nsides), fDz(rhs.fDz), fPlanes(rhs.fPlanes), fx(rhs.fx)
163 {
164 }
165 
166 // Assignment operator
167 BelleCrystal& BelleCrystal::operator = (const BelleCrystal& rhs)
168 {
169  // Check assignment to self
170  if (this == &rhs) { return *this; }
171 
172  // Copy base class data
173  G4CSGSolid::operator=(rhs);
174 
175  // Copy data
176  nsides = rhs.nsides;
177  fDz = rhs.fDz;
178  fx = rhs.fx;
179  fPlanes = rhs.fPlanes;
180 
181  return *this;
182 }
183 
184 // Calculate the coef's of the plane p1->p2->p3->p4->p1
185 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
186 // infront of the plane (i.e. from normal direction).
187 //
188 // Return true if the ThreeVectors are coplanar + set coef;s
189 // false if ThreeVectors are not coplanar
190 G4bool BelleCrystal::MakePlane(const G4ThreeVector& p1, const G4ThreeVector& p2,
191  const G4ThreeVector& p3, const G4ThreeVector& p4,
192  Plane_t& plane) const
193 {
194  G4ThreeVector v12 = p2 - p1;
195  G4ThreeVector v13 = p3 - p1;
196  G4ThreeVector v14 = p4 - p1;
197  G4ThreeVector Vcross = v12.cross(v13);
198 
199  if (std::fabs(Vcross.dot(v14) / (Vcross.mag()*v14.mag())) > kCoplanar_Tolerance) {
200  std::ostringstream message;
201  message << "Verticies are not in the same plane: " << GetName() << " volume =" << Vcross.dot(v14);
202  G4Exception("BelleCrystal::MakePlane()", "BelleCrystal", FatalException, message);
203  }
204 
205  double a = +(p4.y() - p2.y()) * (p3.z() - p1.z()) - (p3.y() - p1.y()) * (p4.z() - p2.z());
206  double b = -(p4.x() - p2.x()) * (p3.z() - p1.z()) + (p3.x() - p1.x()) * (p4.z() - p2.z());
207  double c = +(p4.x() - p2.x()) * (p3.y() - p1.y()) - (p3.x() - p1.x()) * (p4.y() - p2.y());
208  double sd = sqrt(a * a + b * b + c * c); // so now vector plane.(a,b,c) is unit
209  a /= sd;
210  b /= sd;
211  c /= sd;
212  plane.n = G4ThreeVector(a, b, c);
213  plane.d = -(a * p1.x() + b * p1.y() + c * p1.z());
214 
215  // plane.n = (p4-p2).cross(p3-p1).unit();
216  // plane.d =-(plane.n*p1);
217  return true;
218 }
219 
220 // Dispatch to parameterisation for replication mechanism dimension
221 // computation & modification.
222 void BelleCrystal::ComputeDimensions(G4VPVParameterisation*,
223  const G4int,
224  const G4VPhysicalVolume*)
225 {
226  G4Exception("BelleCrystal::ComputeDimensions()",
227  "GeomSolids0001", FatalException,
228  "BelleCrystal does not support Parameterisation.");
229  // std::cout<<"ComputeDimensions"<<std::endl;
230  // p->ComputeDimensions(*this,n,pRep);
231 }
232 
233 // by component min
234 G4ThreeVector min(const G4ThreeVector& a, const G4ThreeVector& b)
235 {
236  return G4ThreeVector(std::min(a.x(), b.x()), std::min(a.y(), b.y()), std::min(a.z(), b.z()));
237 }
238 
239 // by component max
240 G4ThreeVector max(const G4ThreeVector& a, const G4ThreeVector& b)
241 {
242  return G4ThreeVector(std::max(a.x(), b.x()), std::max(a.y(), b.y()), std::max(a.z(), b.z()));
243 }
244 
245 // Calculate extent under transform and specified limit
246 G4bool BelleCrystal::CalculateExtent(const EAxis pAxis,
247  const G4VoxelLimits& bb,
248  const G4AffineTransform& pTransform,
249  G4double& pMin, G4double& pMax) const
250 {
251  const double inf = std::numeric_limits<double>::infinity();
252  G4ThreeVector v0(inf, inf, inf), v1(-inf, -inf, -inf);
253  for (unsigned int i = 0; i < 2 * nsides; i++) {
254  G4ThreeVector v = pTransform.TransformPoint(vertex(i));
255  v0 = min(v0, v);
256  v1 = max(v1, v);
257  }
258  G4ThreeVector b0(bb.GetMinXExtent(), bb.GetMinYExtent(), bb.GetMinZExtent());
259  G4ThreeVector b1(bb.GetMaxXExtent(), bb.GetMaxYExtent(), bb.GetMaxZExtent());
260 
261  v0 = max(v0, b0);
262  v1 = min(v1, b1);
263 
264  switch (pAxis) {
265  case kXAxis: pMin = v0.x(); pMax = v1.x(); break;
266  case kYAxis: pMin = v0.y(); pMax = v1.y(); break;
267  case kZAxis: pMin = v0.z(); pMax = v1.z(); break;
268  default: break;
269  }
270 
271  G4bool flag = false;
272  if ((pMin < inf) || (pMax > -inf)) {
273  flag = true;
274  // Add tolerance to avoid precision troubles
275  pMin -= kCarTolerance ;
276  pMax += kCarTolerance ;
277  }
278  // do {
279  // double t_pMin, t_pMax;
280  // bool t_flag = ref->CalculateExtent(pAxis, bb, pTransform, t_pMin, t_pMax);
281  // if(flag!=t_flag||abs(t_pMin-pMin)>kCarTolerance||abs(t_pMax-pMax)>kCarTolerance){
282  // cout<<GetName()<<" "<<pAxis<<" "<<bb<<" "<<pTransform<<t_pMin-pMin<<" "<<t_pMax-pMax<<endl;
283  // exit(0);
284  // }
285  // } while(0);
286  return flag;
287 }
288 
289 // Return whether point inside/outside/on surface, using tolerance
290 EInside BelleCrystal::Inside(const G4ThreeVector& p) const
291 {
292  COUNTER(0);
293  MATCHOUT("Inside(p) " << p.x() << " " << p.y() << " " << p.z());
294  // return ref->Inside(p);
295  double d = std::fabs(p.z()) - fDz;
296  unsigned int i = 0;
297  do {
298  const Plane_t& t = fPlanes[i++];
299  d = max(t.n * p + t.d, d);
300  } while (i < nsides);
301  const G4double delta = 0.5 * kCarTolerance;
302  int in = 0;
303  in += d <= delta;
304  in += d <= -delta;
305 
306  EInside res = EInside(in);
307  // if(res != ref->Inside(p)){
308  // cout<<GetName()<<" "<<p<<endl;
309  // exit(0);
310  // }
311  return res;
312 }
313 
314 // Calculate side nearest to p, and return normal
315 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
316 G4ThreeVector BelleCrystal::SurfaceNormal(const G4ThreeVector& p) const
317 {
318  // return ref->SurfaceNormal(p);
319  COUNTER(1);
320  MATCHOUT("SurfaceNormal(p) " << p.x() << " " << p.y() << " " << p.z());
321  const G4double delta = 0.5 * kCarTolerance;
322  G4double safe = kInfinity;
323  unsigned int iside = 0, kside = 0;
324  vector<double> adist(nsides + 2);
325  auto dist = [delta, &safe, &iside, &kside, &adist](double d, unsigned int i) -> void {
326  d = std::abs(d);
327  adist[i] = d;
328  iside = (d < safe) ? i : iside;
329  safe = min(d, safe);
330  kside += d < delta;
331  };
332 
333  unsigned int i = 0;
334  do {
335  const Plane_t& t = fPlanes[i];
336  dist(t.n * p + t.d, i++);
337  } while (i < nsides);
338  do {
339  dist(p.z() - fDz, i++);
340  dist(p.z() + fDz, i++);
341  } while (0);
342 
343  G4ThreeVector res;
344  if (kside > 1) {
345  i = 0;
346  do {
347  if (adist[i++] < delta) res += fPlanes[i].n;
348  } while (i < nsides);
349  do {
350  if (adist[i++] < delta) res.setZ(res.z() + 1.0);
351  if (adist[i++] < delta) res.setZ(res.z() - 1.0);
352  } while (0);
353  res = res.unit();
354  } else {
355  if (iside < nsides) {
356  res = fPlanes[iside].n;
357  } else if (iside == nsides) {
358  res = G4ThreeVector(0, 0, 1);
359  } else {
360  res = G4ThreeVector(0, 0, -1);
361  }
362  }
363  // if((res- ref->SurfaceNormal(p)).mag()>kCarTolerance){
364  // cout<<GetName()<<" "<<p<<" "<<res-ref->SurfaceNormal(p)<<endl;
365  // exit(0);
366  // }
367 
368  return res;
369 }
370 
371 // Calculate exact shortest distance to any boundary from outside
372 // This is the best fast estimation of the shortest distance to trap
373 // - Returns 0 is ThreeVector inside
374 G4double BelleCrystal::DistanceToIn(const G4ThreeVector& p) const
375 {
376  // return ref->DistanceToIn(p);
377  COUNTER(2);
378  MATCHOUT("DistanceToIn(p) " << p.x() << " " << p.y() << " " << p.z());
379  G4double d = std::fabs(p.z()) - fDz;
380  unsigned int i = 0;
381  do {
382  const Plane_t& t = fPlanes[i++];
383  d = max(t.n * p + t.d, d);
384  } while (i < nsides);
385  d = max(d, 0.0);
386 
387  // if(abs(d - ref->DistanceToIn(p))>kCarTolerance){
388  // cout<<GetName()<<" "<<p<<" "<<d-ref->DistanceToIn(p)<<endl;
389  // exit(0);
390  // }
391  // if(fmatch&&d==0.0) cout<<"inside\n";
392  return d;
393 }
394 
395 // Calculate exact shortest distance to any boundary from inside
396 // - Returns 0 is ThreeVector outside
397 G4double BelleCrystal::DistanceToOut(const G4ThreeVector& p) const
398 {
399  // return ref->DistanceToOut(p);
400  COUNTER(3);
401  MATCHOUT("DistanceToOut(p) " << p.x() << " " << p.y() << " " << p.z());
402  G4double d = fDz - std::fabs(p.z());
403  unsigned int i = 0;
404  do {
405  const Plane_t& t = fPlanes[i++];
406  d = min(-(t.n * p + t.d), d);
407  } while (i < nsides);
408  d = max(d, 0.0);
409  // if(abs(d - ref->DistanceToOut(p))>0){
410  // cout<<GetName()<<" "<<p<<" "<<d<<" "<<ref->DistanceToOut(p)<<" "<<fDz-ref->GetZHalfLength()<<endl;
411  // exit(0);
412  // }
413  return d;
414 }
415 
416 // Calculate distance to shape from outside - return kInfinity if no intersection
417 //
418 // ALGORITHM:
419 // For each component, calculate pair of minimum and maximum intersection
420 // values for which the particle is in the extent of the shape
421 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
422 G4double BelleCrystal::DistanceToIn(const G4ThreeVector& p,
423  const G4ThreeVector& v) const
424 {
425  // return ref->DistanceToIn(p,v);
426  COUNTER(4);
427  MATCHOUT("DistanceToIn(p,v) " << p.x() << " " << p.y() << " " << p.z() << " " << v.x() << " " << v.y() << " " << v.z());
428  G4double delta = 0.5 * kCarTolerance;
429  double tin = 0, tout = kInfinity;
430  auto outside = [delta, &tin, &tout](double d, double nv) -> bool {
431  double t = -d / nv;
432  if (nv < 0)
433  tin = max(tin , t);
434  else {
435  if (d >= -delta) return true;
436  tout = min(tout, t);
437  }
438  return false;
439  };
440 
441  unsigned int i = 0;
442  do {
443  const Plane_t& t = fPlanes[i++];
444  if (outside(t.n * p + t.d, t.n * v)) {tin = kInfinity; goto exit;}
445  } while (i < nsides);
446  do {
447  if (outside(p.z() - fDz, v.z())) {tin = kInfinity; goto exit;}
448  if (outside(-p.z() - fDz, -v.z())) {tin = kInfinity; goto exit;}
449  } while (0);
450 exit:
451  double res = (tin > tout) ? kInfinity : tin;
452  // G4int oldprc = cout.precision(16);
453  // if(nsides==5) cout<<GetName()<<" "<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<v.x()<<" "<<v.y()<<" "<<v.z()<<" "<<res<<endl;
454  // cout.precision(oldprc);
455 
456  // if(abs(res - ref->DistanceToIn(p,v))>kCarTolerance){
457  // cout<<GetName()<<" "<<p<<" "<<v<<" "<<res-ref->DistanceToIn(p,v)<<endl;
458  // exit(0);
459  // }
460  return res;
461 }
462 
463 // Calculate distance to surface of shape from inside
464 // Calculate distance to x/y/z planes - smallest is exiting distance
465 G4double BelleCrystal::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
466  const G4bool calcNorm, G4bool* IsValid, G4ThreeVector* n) const
467 {
468  // return ref->DistanceToOut(p,v,calcNorm,IsValid,n);
469  COUNTER(5);
470  MATCHOUT("DistanceToOut(p,v) " << p.x() << " " << p.y() << " " << p.z() << " " << v.x() << " " << v.y() << " " << v.z() << " " <<
471  calcNorm);
472  const G4double delta = 0.5 * kCarTolerance;
473  unsigned int iside = 10;
474  G4double lmin = kInfinity;
475 
476  auto outside = [delta, &lmin, &iside](double d, double nv, unsigned int i) -> bool {
477  if (d > delta) // definitly outside
478  {
479  lmin = 0; iside = i; return true;
480  } else {
481  bool c = nv > 0;
482  if (d < -delta) // definitly inside
483  {
484  if (c) { // has to point to the same direction
485  d = -d;
486  if (d < nv * lmin) {lmin = d / nv; iside = i;}
487  }
488  } else { // surface
489  if (c) // points outside
490  {
491  lmin = 0; iside = i; return true;
492  }
493  }
494  }
495  return false;
496  };
497 
498  unsigned int i = 0;
499  do {
500  const Plane_t& t = fPlanes[i];
501  if (outside(t.n * p + t.d, t.n * v, i++)) goto exit;
502  } while (i < nsides);
503  do {
504  if (outside(p.z() - fDz, v.z(), i++)) goto exit;
505  if (outside(-p.z() - fDz, -v.z(), i++)) goto exit;
506  } while (0);
507 exit:
508  if (calcNorm) {
509  *IsValid = true;
510  if (iside < nsides) {
511  *n = fPlanes[iside].n;
512  } else if (iside == nsides) {
513  *n = G4ThreeVector(0, 0, 1);
514  } else {
515  *n = G4ThreeVector(0, 0, -1);
516  }
517  }
518  // G4int oldprc = cout.precision(16);
519  // if(nsides==5) cout<<GetName()<<" "<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<v.x()<<" "<<v.y()<<" "<<v.z()<<" "<<calcNorm<<" "<<lmin<<endl;
520  // cout.precision(oldprc);
521  // do {
522  // G4bool t_IsValid; G4ThreeVector t_n;
523  // double t_lmin = ref->DistanceToOut(p,v,calcNorm,&t_IsValid,&t_n);
524  // if(abs(lmin - t_lmin)>kCarTolerance||*IsValid!=t_IsValid||(*n - t_n).mag()>kCarTolerance){
525  // cout<<GetName()<<" "<<p<<" "<<v<<" "<<lmin-t_lmin<<(*n-t_n)<<endl;
526  // exit(0);
527  // }
528  // } while(0);
529  return lmin;
530 }
531 
532 // GetEntityType
533 G4GeometryType BelleCrystal::GetEntityType() const
534 {
535  return G4String("BelleCrystal");
536 }
537 
538 // Make a clone of the object
539 G4VSolid* BelleCrystal::Clone() const
540 {
541  return new BelleCrystal(*this);
542 }
543 
544 // Stream object contents to an output stream
545 std::ostream& BelleCrystal::StreamInfo(std::ostream& os) const
546 {
547  G4int oldprc = os.precision(16);
548  os << "-----------------------------------------------------------\n"
549  << " *** Dump for solid - " << GetName() << " ***\n"
550  << " ===================================================\n"
551  << " Solid type: BelleCrystal\n"
552  << " Parameters: \n"
553  << " crystal side plane equations:\n"
554  << "-----------------------------------------------------------\n";
555  os.precision(oldprc);
556 
557  return os;
558 }
559 
560 // triangle vertices are ordered in anti-clockwise order looking from outside the solid
561 unsigned int ivertx4[][3] = {{0, 5, 4}, {0, 1, 5}, {1, 7, 5}, {1, 3, 7}, {3, 6, 7}, {3, 2, 6}, {2, 4, 6}, {2, 0, 4}, {0, 3, 1}, {0, 2, 3}, {4, 5, 7}, {4, 7, 6}};
562 unsigned int ivertx5[][3] = {{0, 6, 5}, {0, 1, 6}, {1, 7, 6}, {1, 2, 7}, {2, 8, 7}, {2, 3, 8}, {3, 9, 8}, {3, 4, 9}, {4, 5, 9}, {4, 0, 5}, {0, 2, 1}, {0, 3, 2}, {0, 4, 3}, {5, 6, 7}, {5, 7, 8}, {5, 8, 9}};
563 
564 const unsigned int* BelleCrystal::ivertx(unsigned int it) const
565 {
566  static unsigned int vert[3];
567  // if(nsides==4) return ivertx4[it];
568  if (nsides == 0) return ivertx4[it];
569  // else if(nsides==5) { vert[0] = ivertx5[it][0], vert[1] = ivertx5[it][1], vert[2] = ivertx5[it][2];}
570  else {
571  if (it < 2 * nsides) {
572  int j = it / 2;
573  if ((it & 1) == 0) {
574  vert[0] = j, vert[2] = j + nsides, vert[1] = ((j + 1) % nsides) + nsides;
575  } else {
576  vert[0] = j, vert[2] = ((j + 1) % nsides) + nsides, vert[1] = (j + 1) % nsides;
577  }
578  } else if (it < 2 * nsides + (nsides - 2)) {
579  int j = it - 2 * nsides;
580  vert[0] = 0, vert[2] = 1 + j, vert[1] = 2 + j;
581  } else {
582  int j = it - (2 * nsides + (nsides - 2));
583  vert[0] = nsides, vert[2] = nsides + 2 + j, vert[1] = nsides + 1 + j;
584  }
585  }
586  // cout<<it<<" {"<<vert[0]<<", "<<vert[1]<<", "<<vert[2]<<"}"<<endl;
587  return vert;
588 }
589 
590 G4ThreeVector crossplanes(const Plane_t& p0, const Plane_t& p1, double Dz)
591 {
592  const G4ThreeVector& n0 = p0.n;
593  const G4ThreeVector& n1 = p1.n;
594  double A = n0.z() * Dz + p0.d;
595  double B = n1.z() * Dz + p1.d;
596  double iD = 1 / (n0.y() * n1.x() - n0.x() * n1.y());
597  double x = (A * n1.y() - B * n0.y()) * iD;
598  double y = (B * n0.x() - A * n1.x()) * iD;
599  double z = Dz;
600  return G4ThreeVector(x, y, z);
601 }
602 
603 #define PACK(k0,k1,w,n) ((((k0)<<(0))+((k1)<<(w)))<<((2*w)*n))
604 G4ThreeVector BelleCrystal::vertex(unsigned int i) const
605 {
606  double Dz = (i < nsides) ? -fDz : fDz;
607 
608  // if(nsides==4) {
609  // const unsigned char w = 2, mask = ((1<<w)-1);
610  // unsigned int map = PACK(0,2,w,0)+PACK(0,3,w,1)+PACK(1,2,w,2)+PACK(1,3,w,3); // 20 bits
611  // map >>= (2*w)*(i%4);
612  // int k0 = map&mask, k1 = (map>>2)&mask;
613  // return crossplanes(fPlanes[k0],fPlanes[k1],Dz);
614  // } else if(nsides==5){
615  // const unsigned char w = 3, mask = ((1<<w)-1);
616  // unsigned int map = PACK(0,4,w,0)+PACK(1,0,w,1)+PACK(2,1,w,2)+PACK(3,2,w,3)+PACK(4,3,w,4); // 30 bits
617  // map >>= (2*w)*(i%5);
618  // int k0 = map&mask, k1 = (map>>w)&mask;
619  // return crossplanes(fPlanes[k0],fPlanes[k1],Dz);
620  // }
621  return G4ThreeVector(fx[i].x, fx[i].y, Dz);
622 }
623 #undef PACK
624 
625 // uniformly sampled point over triangle
626 G4ThreeVector BelleCrystal::GetPointOnTriangle(int it) const
627 {
628  // barycentric coordinates
629  double a1 = CLHEP::RandFlat::shoot(0., 1.), a2 = CLHEP::RandFlat::shoot(0., 1.);
630  if (a1 + a2 > 1) { a1 = 1 - a1; a2 = 1 - a2;}
631  double a0 = 1 - (a1 + a2);
632  const unsigned int* iv = ivertx(it);
633  G4ThreeVector p0 = vertex(iv[0]), p1 = vertex(iv[1]), p2 = vertex(iv[2]);
634  G4ThreeVector r1(p0.x()*a0 + p1.x()*a1 + p2.x()*a2,
635  p0.y()*a0 + p1.y()*a1 + p2.y()*a2,
636  p0.z()*a0 + p1.z()*a1 + p2.z()*a2);
637  return r1;
638 }
639 
640 double BelleCrystal::area(int it, double& vol) const
641 {
642  const unsigned int* iv = ivertx(it);
643  G4ThreeVector p0 = vertex(iv[0]), p1 = vertex(iv[1]), p2 = vertex(iv[2]);
644  // std::cout<<it<<" "<<p0<<" "<<p1<<" "<<p2<<std::endl;
645  G4ThreeVector n = (p1 - p0).cross(p2 - p0);
646  vol = n * p0;
647  return 0.5 * n.mag();
648 }
649 
650 double BelleCrystal::getvolarea() const
651 {
652  fareas.clear();
653  int nt = 2 * nsides + 2 * (nsides - 2); // total number of triangles
654  fareas.reserve(nt);
655  double totarea = 0, totvol = 0;
656  for (int i = 0; i < nt; i++) {
657  double v, s = area(i, v);
658  totarea += s;
659  totvol += v;
660  fareas.push_back(totarea);
661  // cout<<i<<" "<<s<<" "<<v<<endl;
662  }
663  return totvol / 6;
664 }
665 
666 G4double BelleCrystal::GetCubicVolume()
667 {
668  // cout<<GetName()<<" GetCubicVolume "<<fCubicVolume<<endl;
669  if (fCubicVolume == 0.0) {
670  fCubicVolume = getvolarea();
671  fSurfaceArea = fareas.back();
672  }
673  return fCubicVolume;
674 }
675 
676 G4double BelleCrystal::GetSurfaceArea()
677 {
678  if (fSurfaceArea == 0.0) {
679  fCubicVolume = getvolarea();
680  fSurfaceArea = fareas.back();
681  }
682  return fSurfaceArea;
683 }
684 
685 // GetPointOnSurface
686 G4ThreeVector BelleCrystal::GetPointOnSurface() const
687 {
688  if (fareas.size() == 0) getvolarea();
689  double r = CLHEP::RandFlat::shoot(0., fareas.back());
690  std::vector<double>::const_iterator it = std::lower_bound(fareas.begin(), fareas.end(), r);
691  return (it != fareas.end()) ? GetPointOnTriangle(it - fareas.begin()) : GetPointOnSurface();
692 }
693 
694 // Methods for visualisation
695 void BelleCrystal::DescribeYourselfTo(G4VGraphicsScene& scene) const
696 {
697  scene.AddSolid(*this);
698 }
699 
700 PolyhedronBelleCrystal::PolyhedronBelleCrystal(int n, const G4ThreeVector* pt)
701 {
702  int nsides = n / 2;
703  AllocateMemory(n, nsides + 2 * (nsides - 2));
704  for (int i = 0; i < n; i++) pV[i + 1] = pt[i];
705 
706  int count = 1;
707  for (int j = 0; j < nsides; j++)
708  pF[count++] = G4Facet(1 + j, 0, 1 + (j + 1) % nsides, 0, 1 + ((j + 1) % nsides) + nsides, 0, 1 + j + nsides, 0);
709  for (int j = 0; j < nsides - 2; j++) pF[count++] = G4Facet(1 + nsides, 0, 2 + j + nsides, 0, 3 + j + nsides, 0, 0, 0);
710  for (int j = 0; j < nsides - 2; j++) pF[count++] = G4Facet(1, 0, 3 + j, 0, 2 + j, 0, 0, 0);
711  SetReferences();
712 }
713 
714 PolyhedronBelleCrystal::~PolyhedronBelleCrystal() {}
715 
716 G4Polyhedron* BelleCrystal::CreatePolyhedron() const
717 {
718  int np = 2 * nsides;
719  vector<G4ThreeVector> pt(np);
720  for (int i = 0; i < np; i++) pt[i] = vertex(i);
721  return new PolyhedronBelleCrystal(np, pt.data());
722 }
723 
724 // Function to define the bounding box
725 void BelleCrystal::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
726 {
727  // All extremums will be at vertices
728  int np = 2 * nsides;
729  G4ThreeVector pt;
730 
731  // Placeholder vectors
732  const double inf = std::numeric_limits<double>::infinity();
733  G4ThreeVector minimum(inf, inf, inf), maximum(-inf, -inf, -inf);
734 
735  for (int i = 0; i < np; i++) {
736  pt = vertex(i);
737 
738  // Assign new minimum and new maximum
739  minimum = min(minimum, pt);
740  maximum = max(maximum, pt);
741  }
742 
743  // Assign
744  pMin = minimum;
745  pMax = maximum;
746 }
Belle2::ECL::BelleCrystal::GetPointOnTriangle
G4ThreeVector GetPointOnTriangle(int) const
Returns a random point on the surface of one of the faces.
Definition: BelleCrystal.cc:626
Belle2::ECL::BelleCrystal::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Two vectors define an axis-parallel bounding box for the shape.
Definition: BelleCrystal.cc:725
Belle2::ECL::BelleCrystal::ivertx
const unsigned int * ivertx(unsigned int i) const
get the ith vertex
Definition: BelleCrystal.cc:564
Belle2::ECL::BelleCrystal::BelleCrystal
BelleCrystal(const G4String &pName)
Constructor for "nominal" BelleCrystal.
Definition: BelleCrystal.cc:137
Belle2::ECL::BelleCrystal::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
calculate the extent of the volume
Definition: BelleCrystal.cc:246
Belle2::ECL::PolyhedronBelleCrystal::PolyhedronBelleCrystal
PolyhedronBelleCrystal(int, const G4ThreeVector *)
constructor
Definition: BelleCrystal.cc:700
Belle2::ECL::Point_t
Definition: BelleCrystal.h:40
Belle2::ECL::BelleCrystal::GetSurfaceArea
G4double GetSurfaceArea()
get the surface area
Definition: BelleCrystal.cc:676
Belle2::ECL::BelleCrystal::~BelleCrystal
virtual ~BelleCrystal()
Destructor.
Definition: BelleCrystal.cc:150
Belle2::ECL::BelleCrystal::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
compute the dimensions
Definition: BelleCrystal.cc:222
Belle2::ECL::BelleCrystal::nsides
unsigned int nsides
the number of sides
Definition: BelleCrystal.h:133
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::ECL::Plane_t
Definition: BelleCrystal.h:34
Belle2::ECL::PolyhedronBelleCrystal
Belle crystal in polyhedron.
Definition: BelleCrystal.h:142
Belle2::ECL::BelleCrystal::GetCubicVolume
G4double GetCubicVolume()
get the Cubic volume
Definition: BelleCrystal.cc:666
Belle2::ECL::BelleCrystal::area
double area(int, double &) const
triangle area
Definition: BelleCrystal.cc:640
Belle2::ECL::PolyhedronBelleCrystal::~PolyhedronBelleCrystal
virtual ~PolyhedronBelleCrystal()
destructor
Definition: BelleCrystal.cc:714
Belle2::ECL::BelleCrystal::operator=
BelleCrystal & operator=(const BelleCrystal &rhs)
assignment operator
Definition: BelleCrystal.cc:167