Belle II Software  release-05-02-19
BelleLathe.cc
1 #include "ecl/geometry/BelleLathe.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 std;
15 using namespace Belle2;
16 using namespace ECL;
17 
18 #define COMPARE 0
19 #define PERFCOUNTER 0
20 #if PERFCOUNTER==1
21 typedef int counter_t[6];
22 map<string, counter_t> counterl;
23 #define COUNTER(x) counterl[GetName()][x]++
24 //#define MATCHOUT(x) //if(GetName().find("sv_crystalcontainersolid")==0) cout<<x<<endl;
25 #else
26 #define COUNTER(x)
27 //
28 #endif
29 
30 //#define MATCHOUT(x) G4cout<<GetName()<<" "<<x<<G4endl;
31 #define MATCHOUT(x)
32 
33 struct Plane_t {
34  G4ThreeVector n;// Normal unit vector (x,y,z)
35  double d; // offset (d)
36  // => n.x*x + n.y*y + n.z*z + d = 0
37 };
38 
39 namespace Belle2 {
44  namespace ECL {
45  inline double dotxy(const G4ThreeVector& p, const G4ThreeVector& n)
46  {
47  return p.x() * n.x() + p.y() * n.y();
48  }
49  }
51 }
52 
53 ostream& operator <<(ostream& o, const zr_t& v)
54 {
55  return o << "{" << v.z << ", " << v.r << "}";
56 }
57 
58 struct curl_t {
59  G4ThreeVector v;
60  explicit curl_t(const G4ThreeVector& _v): v(_v) {}
61 };
62 
63 ostream& operator <<(ostream& o, const curl_t& c)
64 {
65  return o << "{" << c.v.x() << ", " << c.v.y() << ", " << c.v.z() << "}, ";
66 }
67 
68 
69 BelleLathe::BelleLathe(const G4String& pName, double phi0, double dphi, const vector<zr_t>& c)
70  : G4CSGSolid(pName)
71 {
72  Init(c, phi0, dphi);
73 }
74 
75 BelleLathe::BelleLathe(const G4String& pName, double phi0, double dphi, int n, double* z, double* rin, double* rout)
76  : G4CSGSolid(pName)
77 {
78  vector<zr_t> contour;
79  for (int i = 0; i < n; i++) {
80  zr_t t = {z[i], rin[i]};
81  contour.push_back(t);
82  }
83  for (int i = n - 1; i >= 0; i--) {
84  zr_t t = {z[i], rout[i]};
85  contour.push_back(t);
86  }
87 
88  Init(contour, phi0, dphi);
89 }
90 
91 void BelleLathe::Init(const vector<zr_t>& c, double phi0, double dphi)
92 {
93  vector<zr_t> contour = c;
94  // remove duplicated vertices
95  do {
96  vector<zr_t>::iterator it0 = contour.begin(), it1 = it0 + 1;
97  for (; it1 != contour.end();) {
98  const zr_t& s0 = *it0, &s1 = *it1;
99  if (abs(s0.z - s1.z) < kCarTolerance && abs(s0.r - s1.r) < kCarTolerance)
100  it1 = contour.erase(it1);
101  else {
102  ++it0; ++it1;
103  }
104  }
105  const zr_t& s0 = *it0, &s1 = contour[0];
106  if (abs(s0.z - s1.z) < kCarTolerance && abs(s0.r - s1.r) < kCarTolerance) contour.erase(it0);
107  } while (0);
108 
109  // remove vertices on the same line
110  do {
111  auto inc = [&contour](vector<zr_t>::iterator & it) -> void {
112  if (++it >= contour.end()) it = contour.begin();
113  };
114  auto dec = [&contour](vector<zr_t>::iterator & it) -> void {
115  if (--it < contour.begin()) it = (++contour.rbegin()).base();
116  };
117  vector<zr_t>::iterator it0 = contour.begin(), it1 = it0 + 1, it2 = it1 + 1;
118  for (; it0 != contour.end();) {
119  const zr_t& s0 = *it0, &s1 = *it1, &s2 = *it2;
120  double dr2 = s2.r - s0.r, dz2 = s2.z - s0.z;
121  double d = (s1.z - s0.z) * dr2 - (s1.r - s0.r) * dz2;
122  if (d * d < kCarTolerance * kCarTolerance * (dr2 * dr2 + dz2 * dz2)) {
123  it1 = contour.erase(it1); it2 = it1; inc(it2); it0 = it1; dec(it0);
124  } else {
125  ++it0; inc(it1); inc(it2);
126  }
127  }
128  } while (0);
129 
130  auto isClockwise = [&contour]() -> bool {
131  double sum = 0;
132  zr_t p0 = contour[0];
133  for (int i = 1, imax = contour.size(); i < imax; i++)
134  {
135  zr_t p1 = contour[i]; sum += (p1.z - p0.z) * (p1.r + p0.r);
136  p0 = p1;
137  }
138  zr_t p1 = contour[0]; sum += (p1.z - p0.z) * (p1.r + p0.r);
139  return sum > 0;
140  };
141 
142  if (isClockwise()) {
143  // std::ostringstream message;
144  // message << "Polygon is not in anti-clockwise order: " << GetName() << "\nReversing order...";
145  // G4Exception("BelleLathe::BelleLathe()", "BelleLathe", JustWarning, message);
146  std::reverse(contour.begin(), contour.end());
147  }
148  fcontour = contour;
149 
150  auto convexside = [this](cachezr_t& s, double eps) -> void {
151  s.isconvex = false;
152  if (s.dz > 0) return;
153  vector<zr_t>::const_iterator it = fcontour.begin();
154  double a = s.dz * s.is, b = s.dr * s.is, cc = b * s.z - a * s.r;
155  bool dp = false, dm = false;
156  s.isconvex = true;
157  do {
158  const zr_t& p = *it;
159  double d = a * p.r - b * p.z + cc; // distance to line
160  dm = dm || (d < -eps);
161  dp = dp || (d > eps);
162  if (dm && dp) {s.isconvex = false; return;}
163  } while (++it != fcontour.end());
164  };
165 
166  frmin = kInfinity;
167  frmax = -kInfinity;
168  fzmin = kInfinity;
169  fzmax = -kInfinity;
170  fcache.reserve(fcontour.size());
171  for (int i = 0, n = fcontour.size(); i < n; i++) {
172  const zr_t& s0 = fcontour[i], &s1 = fcontour[(i + 1) % n];
173  cachezr_t t;
174  t.z = s0.z;
175  t.r = s0.r;
176  t.dz = s1.z - s0.z;
177  t.dr = s1.r - s0.r;
178  t.s2 = t.dz * t.dz + t.dr * t.dr;
179  t.is2 = 1 / t.s2;
180  t.is = sqrt(t.is2);
181  t.zmin = min(s0.z, s1.z);
182  t.zmax = max(s0.z, s1.z);
183  t.r2min = pow(min(s0.r, s1.r), 2);
184  t.r2max = pow(max(s0.r, s1.r), 2);
185  t.ta = (s1.r - s0.r) / (s1.z - s0.z);
186  convexside(t, kCarTolerance);
187  fcache.push_back(t);
188 
189  frmax = max(frmax, s0.r);
190  frmin = min(frmin, s0.r);
191  fzmax = max(fzmax, s0.z);
192  fzmin = min(fzmin, s0.z);
193  }
194 
195  fphi = phi0;
196  fdphi = dphi;
197 
198  fdphi = std::min(2 * M_PI, fdphi);
199  fdphi = std::max(0.0, fdphi);
200  fc0 = cos(fphi);
201  fs0 = sin(fphi);
202  fc1 = cos(fphi + fdphi);
203  fs1 = sin(fphi + fdphi);
204 
205  fn0x = fs0;
206  fn0y = -fc0;
207  fn1x = -fs1;
208  fn1y = fc1;
209  fgtpi = fdphi > M_PI;
210  ftwopi = abs(fdphi - 2 * M_PI) < kCarTolerance;
211 
212  // cout << ftwopi << " " << fgtpi << " " << fn0y << " " << fn0x << " " << fn1y << " " << fn1x << endl;
213 
214  for (int i = 0, n = fcontour.size(); i < n; i++) {
215  const zr_t& s = fcontour[i];
216  fz.push_back(s.z);
217  }
218  sort(fz.begin(), fz.end());
219  fz.erase(std::unique(fz.begin(), fz.end()), fz.end());
220 
221  for (int i = 1, ni = fz.size(); i < ni; i++) {
222  double a = fz[i - 1], b = fz[i];
223  findx.push_back(fseg.size());
224  for (int j = 0, nj = fcache.size(); j < nj; j++) {
225  const cachezr_t& sj = fcache[j];
226  double cc = sj.zmin, d = sj.zmax;
227  if (cc != d and b > cc and d > a) { // overlap
228  fseg.push_back(j);
229  }
230  }
231  }
232  findx.push_back(fseg.size());
233 
234  getvolarea();
235 
236 #if COMPARE>0
237  auto getpolycone = [](const G4String & pName, double phi0, double dphi, const vector<zr_t>& c) -> G4GenericPolycone* {
238  vector<double> r, z;
239  r.reserve(c.size());
240  z.reserve(c.size());
241  for (int i = 0, imax = c.size(); i < imax; i++)
242  {
243  r.push_back(c[i].r);
244  z.push_back(c[i].z);
245  }
246  return new G4GenericPolycone(pName, phi0, dphi, c.size(), r.data(), z.data());
247  };
248  fshape = getpolycone(GetName(), phi0, dphi, fcontour);
249 #else
250  fshape = NULL;
251 #endif
252 // StreamInfo(G4cout);
253 }
254 
255 // Nominal constructor for BelleLathe whose parameters are to be set by
256 // a G4VParamaterisation later. Check and set half-widths as well as
257 // angles: final check of coplanarity
258 BelleLathe::BelleLathe(const G4String& pName)
259  : G4CSGSolid(pName)
260 {
261  vector<zr_t> a;
262  Init(a, 0, 2 * M_PI);
263 }
264 
265 // Fake default constructor - sets only member data and allocates memory
266 // for usage restricted to object persistency.
267 BelleLathe::BelleLathe(__void__& a)
268  : G4CSGSolid(a)
269 {
270  vector<zr_t> b;
271  Init(b, 0, 2 * M_PI);
272 }
273 
274 // Destructor
276 {
277 #if PERFCOUNTER==1
278  cout << GetName() << " ";
279  for (int i = 0; i < 6; i++) cout << counterl[GetName()][i] << " "; cout << endl;
280 #endif
281 }
282 
283 // Copy constructor
285  : G4CSGSolid(rhs), fcontour(rhs.fcontour), fcache(rhs.fcache), fz(rhs.fz),
286  findx(rhs.findx), fseg(rhs.fseg), farea(rhs.farea), ftlist(rhs.ftlist),
287  fphi(rhs.fphi), fdphi(rhs.fdphi), fs0(rhs.fs0), fc0(rhs.fc0), fs1(rhs.fs1),
288  fc1(rhs.fc1), fn0x(rhs.fn0x), fn0y(rhs.fn0y), fn1x(rhs.fn1x), fn1y(rhs.fn1y),
289  frmin(rhs.frmin), frmax(rhs.frmax), fzmin(rhs.fzmin), fzmax(rhs.fzmax),
290  fgtpi(rhs.fgtpi), ftwopi(rhs.ftwopi), fshape(rhs.fshape), fsurf(rhs.fsurf)
291 {
292 }
293 
294 // Assignment operator
296 {
297  // Check assignment to self
298  if (this == &rhs) { return *this; }
299 
300  // Copy base class data
301  G4CSGSolid::operator=(rhs);
302 
303  // Copy data
304  fcontour = rhs.fcontour;
305  fcache = rhs.fcache;
306  fz = rhs.fz;
307  findx = rhs.findx;
308  fseg = rhs.fseg;
309  farea = rhs.farea;
310  ftlist = rhs.ftlist;
311  fphi = rhs.fphi;
312  fdphi = rhs.fdphi;
313  fs0 = rhs.fs0;
314  fc0 = rhs.fc0;
315  fs1 = rhs.fs1;
316  fc1 = rhs.fc1;
317  fn0x = rhs.fn0x;
318  fn0y = rhs.fn0y;
319  fn1x = rhs.fn1x;
320  fn1y = rhs.fn1y;
321  frmin = rhs.frmin;
322  frmax = rhs.frmax;
323  fzmin = rhs.fzmin;
324  fzmax = rhs.fzmax;
325  fgtpi = rhs.fgtpi;
326  ftwopi = rhs.ftwopi;
327  fshape = rhs.fshape;
328  fsurf = rhs.fsurf;
329  return *this;
330 }
331 
332 
333 // Dispatch to parameterisation for replication mechanism dimension
334 // computation & modification.
335 void BelleLathe::ComputeDimensions(G4VPVParameterisation*,
336  const G4int,
337  const G4VPhysicalVolume*)
338 {
339  G4Exception("BelleLathe::ComputeDimensions()",
340  "GeomSolids0001", FatalException,
341  "BelleLathe does not support Parameterisation.");
342  // std::cout<<"ComputeDimensions"<<std::endl;
343  // p->ComputeDimensions(*this,n,pRep);
344 }
345 
346 vector<double> quadsolve(double a, double b, double c)
347 {
348  // solve equation a*t^2 + b*t + c = 0 taking care intermediate rounding errors
349  vector<double> t(2);
350  b *= 0.5;
351  double D = b * b - a * c;
352  if (D >= 0) {
353  double sD = sqrt(D);
354  double sum = b + ((b > 0) ? sD : -sD);
355  double t0 = -c / sum;
356  double t1 = -sum / a;
357  t[0] = t0;
358  t[1] = t1;
359  } else {
360  t.clear();
361  }
362 
363  return t;
364 }
365 
366 inline int quadsolve(double a, double b, double c, double& t0, double& t1)
367 {
368  // solve equation a*t^2 + b*t + c = 0 taking care intermediate rounding errors
369  b *= 0.5;
370  double D = b * b - a * c;
371  if (D >= 0) {
372  double sD = sqrt(D);
373  double sum = b + ((b > 0) ? sD : -sD);
374  t0 = -c / sum;
375  t1 = -sum / a;
376  return 2;
377  }
378  return 0;
379 }
380 
381 struct solution_t {double t, s;};
382 vector<solution_t> extremum(double A, double B, double C, double D, double E, double F)
383 {
384  // extremum of Fun(t,s) = A*t*t + B*t*s + C*s*s + D*t + E*s + F => dFun/ds = 0
385  vector<solution_t> res;
386  if (abs(B) < abs(A)) {
387  double a = 4 * A * C - B * B;
388  double b = 2 * (2 * A * E - B * D);
389  double c = 4 * A * F - D * D;
390  vector<double> ss = quadsolve(a, b, c);
391  for (auto s : ss) {
392  if (fpclassify(s) == FP_INFINITE) continue;
393  double t = -(s * B + D) / (2 * A);
394  solution_t r = {t, s};
395  res.push_back(r);
396  }
397  } else {
398  double B2 = B * B, CD = C * D, BE = B * E;
399  double a = A * (4 * A * C - B2);
400  double b = 2 * A * (2 * CD - BE);
401  double c = D * (CD - BE) + B2 * F;
402  vector<double> ts = quadsolve(a, b, c);
403  for (auto t : ts) {
404  if (fpclassify(t) == FP_INFINITE) continue;
405  double s = -(2 * t * A + D) / B;
406  solution_t r = {t, s};
407  res.push_back(r);
408  }
409  }
410  return res;
411 }
412 
413 // calculate all ray solid's surface intersection return ordered vector
414 vector<double> BelleLathe::linecross(const G4ThreeVector& p, const G4ThreeVector& n) const
415 {
416  auto hitside = [this, &p, &n](double t, double zmin, double zmax) -> bool {
417  double z = p.z() + n.z() * t;
418  bool k = zmin < z && z <= zmax;
419  if (k && !ftwopi)
420  {
421  double x = p.x() + n.x() * t;
422  double y = p.y() + n.y() * t;
423  k = k && insector(x, y);
424  }
425  return k;
426  };
427 
428  auto hitzside = [this, &p, &n](double t, double r2min, double r2max) -> bool {
429  double x = p.x() + n.x() * t;
430  double y = p.y() + n.y() * t;
431  double r2 = x * x + y * y;
432  bool k = r2min <= r2 && r2 < r2max;
433  if (k && !ftwopi)
434  {
435  k = k && insector(x, y);
436  }
437  return k;
438  };
439 
440  vector<double> tc;
441  double inz = 1 / n.z();
442  double nn = Belle2::ECL::dotxy(n, n), np = dotxy(n, p), pp = dotxy(p, p);
443  for (const cachezr_t& s : fcache) { // loop over sides
444  if (s.dz == 0.0) { // z-plane
445  double t = (s.z - p.z()) * inz;
446  if (hitzside(t, s.r2min, s.r2max)) { tc.push_back(t); }
447  } else {
448  double ta = s.ta;
449  double A, B, R2;
450  if (s.dr == 0.0) { // cylinder
451  double R = s.r;
452  R2 = R * R;
453 
454  A = -nn;
455  B = np;
456  } else { // cone
457  double taz = ta * (p.z() - s.z);
458  double R = taz + s.r;
459  R2 = R * R;
460 
461  double nzta = n.z() * ta;
462  A = nzta * nzta - nn;
463  B = np - nzta * R;
464  }
465  double D = B * B + (pp - R2) * A;
466  if (D > 0) {
467  double sD = sqrt(D), iA = 1 / A;
468  double t0 = (B + sD) * iA, t1 = (B - sD) * iA;
469  if (hitside(t0, s.zmin, s.zmax)) tc.push_back(t0);
470  if (hitside(t1, s.zmin, s.zmax)) tc.push_back(t1);
471  }
472  }
473  }
474 
475  if (!ftwopi) {
476  do { // side at phi0
477  double d = fn0x * p.x() + fn0y * p.y();
478  double vn = fn0x * n.x() + fn0y * n.y();
479  double t = -d / vn;
480  G4ThreeVector r = p + n * t;
481  zr_t zr = {r.z(), fc0 * r.x() + fs0 * r.y()};
482  if (vn != 0 && wn_poly(zr) == 2) tc.push_back(t);
483  } while (0);
484 
485  do { // side at phi0+dphi
486  double d = fn1x * p.x() + fn1y * p.y();
487  double vn = fn1x * n.x() + fn1y * n.y();
488  double t = -d / vn;
489  G4ThreeVector r = p + n * t;
490  zr_t zr = {r.z(), fc1 * r.x() + fs1 * r.y()};
491  if (vn != 0 && wn_poly(zr) == 2) tc.push_back(t);
492  } while (0);
493  }
494 
495  sort(tc.begin(), tc.end());
496  return tc;
497 }
498 
499 // Calculate extent under transform and specified limit
500 G4bool BelleLathe::CalculateExtent(const EAxis A,
501  const G4VoxelLimits& bb,
502  const G4AffineTransform& T,
503  G4double& pMin, G4double& pMax) const
504 {
505  auto maxdist = [this](const G4ThreeVector & n) -> G4ThreeVector {
506  G4ThreeVector r;
507  int i = 0, nsize = fcache.size();
508  if (ftwopi || insector(n.x(), n.y())) // n in sector
509  {
510  double nr = hypot(n.x(), n.y()), nz = n.z();
511  double dmax = -kInfinity, R = 0, Z = 0;
512  do {
513  const cachezr_t& s = fcache[i];
514  double d1 = nz * s.z + nr * s.r;
515  if (dmax < d1) { R = s.r; Z = s.z; dmax = d1;}
516  } while (++i < nsize);
517  if (nr > 0) {
518  R /= nr;
519  r.set(R * n.x(), R * n.y(), Z);
520  } else {
521  double phi = fphi + 0.5 * fdphi;
522  r.set(R * cos(phi), R * sin(phi), Z);
523  }
524  } else {
525  double dmax = -kInfinity;
526  do {
527  const cachezr_t& s = fcache[i];
528  // check both sides
529  G4ThreeVector rf(-fn0y * s.r, fn0x * s.r, s.z), rl(fn1y * s.r, -fn1x * s.r, s.z);
530  double d0 = rf * n, d1 = rl * n;
531  // cout<<rf<<" "<<rl<<endl;
532  if (dmax < d0) { r = rf; dmax = d0;}
533  if (dmax < d1) { r = rl; dmax = d1;}
534  } while (++i < nsize);
535  }
536  return r;
537  };
538 
539  struct seg_t {int i0, i1;};
540  auto clip = [](vector<G4ThreeVector>& vlist, vector<seg_t>& slist, const G4ThreeVector & n, double dist) {
541  vector<seg_t> snew;
542  vector<int> lone;
543 
544  vector<double> d;
545  for (const G4ThreeVector& v : vlist) d.push_back(v * n + dist);
546 
547  for (seg_t s : slist) {
548  double prod = d[s.i0] * d[s.i1];
549  // cout<<d[s.i0]<<" "<<d[s.i1]<<endl;
550  if (prod < 0) { // segment crosses plane - break it
551  G4ThreeVector rn = (vlist[s.i0] * d[s.i1] - vlist[s.i1] * d[s.i0]) * (1 / (d[s.i1] - d[s.i0]));
552  lone.push_back(vlist.size()); vlist.push_back(rn);
553  if (d[s.i0] < 0) {
554  s = {lone.back(), s.i1};
555  } else {
556  s = {s.i0, lone.back()};
557  }
558  } else if (prod == 0) { // segment end on plane
559  if (d[s.i0] == 0 && d[s.i1] > 0) {
560  lone.push_back(s.i0);
561  } else if (d[s.i0] > 0 && d[s.i1] == 0) {
562  lone.push_back(s.i1);
563  } else continue;
564  } else {
565  if (d[s.i0] < 0) continue; // segment below plane
566  }
567  snew.push_back(s);
568  }
569 
570  double dmax = -1e99;
571  int imax = -1, jmax = -1;
572  // search for the most distant points on the clipping plane
573  for (unsigned int i = 0; i < lone.size(); i++) {
574  for (unsigned int j = i + 1; j < lone.size(); j++) {
575  double d2 = (vlist[lone[i]] - vlist[lone[j]]).mag2();
576  if (d2 > dmax) { imax = lone[i]; jmax = lone[j];}
577  }
578  }
579 
580  // close the new polygon by creating new segments
581  if (imax >= 0) {
582  G4ThreeVector k = vlist[jmax] - vlist[imax];
583  sort(lone.begin(), lone.end(), [&k, &vlist, &imax](int i, int j) {return k * (vlist[i] - vlist[imax]) < k * (vlist[j] - vlist[imax]);});
584 
585  for (unsigned int i = 0; i < lone.size(); i += 2) {
586  seg_t t = {lone[i], lone[i + 1]};
587  for (const seg_t& s : snew) {
588  if (t.i1 == s.i0) { snew.push_back(t); break;}
589  if (t.i0 == s.i0) { swap(t.i0, t.i1); snew.push_back(t); break;}
590  }
591  }
592  }
593  swap(slist, snew);
594  };
595 
596  auto PhiCrossN = [this, clip](const vector<Plane_t>& planes) {
597  // unordered clipped phi-sides vertices within
598  // limiting planes
599  vector<G4ThreeVector> vlist; // vertex list
600  vector<seg_t> slist; // segment list
601  vector<G4ThreeVector> res;
602 
603  int nsize = fcache.size();
604  vlist.reserve(nsize);
605  slist.reserve(nsize);
606  for (int iphi = 0; iphi < 2; iphi++) {
607  vlist.clear();
608  slist.clear();
609  // phi-side directional vector is (kx,ky)
610  double kx = iphi ? -fn0y : fn1y, ky = iphi ? fn0x : -fn1x;
611  do {
612  int i = 0;
613  do {
614  const cachezr_t& s = fcache[i];
615  G4ThreeVector r(kx * s.r, ky * s.r, s.z);
616  vlist.push_back(r);
617  seg_t t = {i, i + 1};
618  slist.push_back(t);
619  } while (++i < nsize - 1);
620  const cachezr_t& s = fcache[nsize - 1];
621  G4ThreeVector r(kx * s.r, ky * s.r, s.z);
622  vlist.push_back(r);
623  seg_t t = {nsize - 1, 0};
624  slist.push_back(t);
625  } while (0);
626 
627  // clip phi-side polygon by limiting planes
628  for (const Plane_t& p : planes) {
629  // cout<<p.n<<" "<<p.d<<endl;
630  clip(vlist, slist, p.n, p.d);
631  // for(auto t:vlist) cout<<t<<" "; cout<<endl;
632  }
633  vector<bool> bv(vlist.size(), false);
634 
635  for (vector<seg_t>::const_iterator it = slist.begin(); it != slist.end(); ++it) {
636  bv[(*it).i0] = true;
637  bv[(*it).i1] = true;
638  }
639 
640  for (unsigned int i = 0; i < vlist.size(); i++) {
641  if (!bv[i]) continue;
642  res.push_back(vlist[i]);
643  }
644  }
645  return res;
646  };
647 
648  auto RCross = [this](const G4ThreeVector & op, const G4ThreeVector & k, const G4ThreeVector & u) {
649  // plane with origin at op and normal vector n = [k x u], k and u are orthogonal k*u = 0
650  // plane equation r = t*k + s*u + op
651  vector<solution_t> ts;
652  int nsize = fcache.size();
653  int i = 0;
654  do {
655  const cachezr_t& seg = fcache[i];
656  // r0 -- cone radius at z0, c -- cone axis
657  // cone equation is (r0 + tg * ((r-c0)*c))^2 = (r-c0)^2 - ((r-c0)*c)^2
658  double r0 = seg.r, z0 = seg.z, tg = seg.ta, tg2 = tg * tg;
659  double rtg = r0 * tg;
660 
661  G4ThreeVector o(op.x(), op.y(), op.z() - z0);
662 
663  double ko = k * o, uo = u * o, ck = k.z(), cu = u.z(), co = o.z();
664  double k2 = 1, u2 = 1, o2 = o * o;
665  double ck2 = ck * ck, cu2 = cu * cu, co2 = co * co;
666  double dr2 = r0 * r0 - o2;
667  if (seg.dz != 0.0) {
668  double q0 = 1 + tg2;
669  double q1 = co * q0 + rtg;
670 
671  double F00 = co2 * q0 + 2 * co * rtg + dr2;
672  double F10 = 2 * (ck * q1 - ko);
673  double F20 = ck2 * q0 - k2;
674  double F01 = 2 * (cu * q1 - uo);
675  double F11 = 2 * ck * cu * q0;
676  double F02 = cu2 * q0 - u2;
677 
678  vector<solution_t> res = extremum(F02, F11, F20, F01, F10, F00);
679  for (const solution_t& r : res) {
680  double t = r.s, s = r.t;
681  G4ThreeVector p = t * k + s * u + op;
682  if (seg.zmin < p.z() && p.z() < seg.zmax) {
683  solution_t e = {t, s};
684  if (ftwopi || insector(p.x(), p.y()))
685  ts.push_back(e);
686  }
687  }
688  }
689  double a = -(ck2 * u2 + cu2 * k2);
690  if (a != 0) {
691  if (abs(cu) > abs(ck)) {
692  double b = 2 * (ck * (cu * uo - co * u2) - cu2 * ko);
693  double c = co * (2 * cu * uo - co * u2) + cu2 * dr2;
694  vector<double> tv = quadsolve(a, b, c);
695  for (double t : tv) {
696  double s = -(co + ck * t) / cu;
697  G4ThreeVector p = t * k + s * u + op;
698  if (ftwopi || insector(p.x(), p.y())) {
699  solution_t e = {t, s};
700  ts.push_back(e);
701  }
702  }
703  } else {
704  double b = 2 * (cu * (ck * ko - co * k2) - ck2 * uo);
705  double c = co * (2 * ck * ko - co * k2) + ck2 * dr2;
706  vector<double> sv = quadsolve(a, b, c);
707  for (double s : sv) {
708  double t = -(co + cu * s) / ck;
709  G4ThreeVector p = t * k + s * u + op;
710  if (ftwopi || insector(p.x(), p.y())) {
711  solution_t e = {t, s};
712  ts.push_back(e);
713  }
714  }
715  }
716  }
717  } while (++i < nsize);
718  return ts;
719  };
720 
721  bool b1 = false, b2 = false;
722  G4ThreeVector n0, n1, n2;
723  switch (A) {
724  case kXAxis: n0.set(1, 0, 0); n1.set(0, 1, 0); n2.set(0, 0, 1); b1 = bb.IsYLimited(); b2 = bb.IsZLimited(); break;
725  case kYAxis: n0.set(0, 1, 0); n1.set(1, 0, 0); n2.set(0, 0, 1); b1 = bb.IsXLimited(); b2 = bb.IsZLimited(); break;
726  case kZAxis: n0.set(0, 0, 1); n1.set(1, 0, 0); n2.set(0, 1, 0); b1 = bb.IsXLimited(); b2 = bb.IsYLimited(); break;
727  default: break;
728  }
729 
730  double dmin1 = -kInfinity, dmax1 = kInfinity;
731  if (b1) {
732  switch (A) {
733  case kXAxis: dmin1 = bb.GetMinYExtent(); dmax1 = bb.GetMaxYExtent(); break;
734  case kYAxis: dmin1 = bb.GetMinXExtent(); dmax1 = bb.GetMaxXExtent(); break;
735  case kZAxis: dmin1 = bb.GetMinXExtent(); dmax1 = bb.GetMaxXExtent(); break;
736  default: break;
737  }
738  }
739 
740  double dmin2 = -kInfinity, dmax2 = kInfinity;
741  if (b2) {
742  switch (A) {
743  case kXAxis: dmin2 = bb.GetMinZExtent(); dmax2 = bb.GetMaxZExtent(); break;
744  case kYAxis: dmin2 = bb.GetMinZExtent(); dmax2 = bb.GetMaxZExtent(); break;
745  case kZAxis: dmin2 = bb.GetMinYExtent(); dmax2 = bb.GetMaxYExtent(); break;
746  default: break;
747  }
748  }
749 
750  G4AffineTransform iT = T.Inverse();
751  // axis to solid coordinates
752  G4ThreeVector n0t = iT.TransformAxis(n0);
753  G4ThreeVector smin = n0t * kInfinity, smax = (-kInfinity) * n0t; // extremum points in solid coordinate system
754  double pmin = kInfinity, pmax = -pmin;
755  if (b1 && b2) {
756  G4ThreeVector corners[] = {n1* dmin1 + n2 * dmin2, n1* dmax1 + n2 * dmin2, n1* dmax1 + n2 * dmax2, n1* dmin1 + n2 * dmax2};
757  for (G4ThreeVector& c : corners) iT.ApplyPointTransform(c); // to solid coordinates
758 
759  vector<Plane_t> planes;
760  for (int i = 0; i < 4; i++) {
761  const G4ThreeVector& c0 = corners[i], &c1 = corners[(i + 1) % 4];
762  vector<double> dists = linecross(c0, n0t);
763  // cout<<"c0 "<<c0<<endl;
764  for (double t : dists) {
765  G4ThreeVector p = n0t * t + c0;
766  double tt = t + c0 * n0t;
767  // cout<<p<<" "<<tt<<endl;
768  if (pmax < tt) { pmax = tt; smax = p;}
769  if (pmin > tt) { pmin = tt; smin = p;}
770  }
771 
772  G4ThreeVector u = c1 - c0, un = u.unit();
773  vector<solution_t> ts = RCross(c0, n0t, un);
774  double umax = u.mag();
775  for (solution_t r : ts) {
776  if (0 < r.s && r.s < umax) {
777  double tt = r.t + c0 * n0t;
778  G4ThreeVector p = n0t * r.t + un * r.s + c0;
779  // cout<<r.t<<" "<<r.s<<" "<<smax<<endl;
780  if (pmax < tt) { pmax = tt; smax = p;}
781  if (pmin > tt) { pmin = tt; smin = p;}
782  }
783  }
784  planes.push_back({ -un, un * c1});
785  }
786 
787  vector<G4ThreeVector> vside = PhiCrossN(planes);
788  for (G4ThreeVector& p : vside) {
789  // cout<<p<<endl;
790  double tt = n0t * p;
791  if (pmax < tt) { pmax = tt; smax = p;}
792  if (pmin > tt) { pmin = tt; smin = p;}
793  }
794 
795  } else if (b1 || b2) {
796  G4ThreeVector limits[2], u;
797  if (b1) {
798  limits[0] = n1 * dmin1;
799  limits[1] = n1 * dmax1;
800  u = iT.TransformAxis(n2);
801  } else {
802  limits[0] = n2 * dmin2;
803  limits[1] = n2 * dmax2;
804  u = iT.TransformAxis(n1);
805  }
806 
807  for (G4ThreeVector& c : limits) iT.ApplyPointTransform(c); // to solid coordinates
808  for (int i = 0; i < 2; i++) {
809  vector<solution_t> ts = RCross(limits[i], n0t, u);
810  for (solution_t r : ts) {
811  double tt = r.t + limits[i] * n0t;
812  G4ThreeVector p = n0t * r.t + u * r.s + limits[i];
813  // cout<<r.t<<" "<<r.s<<" "<<endl;
814  if (pmax < tt) { pmax = tt; smax = p;}
815  if (pmin > tt) { pmin = tt; smin = p;}
816  }
817  }
818 
819  vector<Plane_t> planes(2);
820  G4ThreeVector n;
821  if (b1) {
822  n = iT.TransformAxis(n1);
823  } else {
824  n = iT.TransformAxis(n2);
825  }
826  planes[0] = { n, -limits[0]* n};
827  planes[1] = { -n, limits[1]* n};
828  vector<G4ThreeVector> vside = PhiCrossN(planes);
829 
830  for (G4ThreeVector& p : vside) {
831  // double t = n0t*(p-limits[0]);
832  double tt = n0t * p;
833  // cout<<tt<<" "<<p<<" "<<endl;
834  if (pmax < tt) { pmax = tt; smax = p;}
835  if (pmin > tt) { pmin = tt; smin = p;}
836  }
837  }
838  // maximal distance in +- directions
839  G4ThreeVector rp = maxdist(n0t), rm = maxdist(-n0t);
840  if (bb.Inside(T.TransformPoint(rm))) {
841  double tt = rm * n0t;
842  if (pmin > tt) {pmin = tt; smin = rm;}
843  }
844  if (bb.Inside(T.TransformPoint(rp))) {
845  double tt = rp * n0t;
846  if (pmax < tt) {pmax = tt; smax = rp;}
847  }
848 
849  // to mother volume coordinate system
850  T.ApplyPointTransform(smin);
851  T.ApplyPointTransform(smax);
852  pmin = n0 * smin;
853  pmax = n0 * smax;
854 
855  pmin -= kCarTolerance;
856  pmax += kCarTolerance;
857  // bool hit = pmin > -kInfinity && pmax < kInfinity;
858  bool hit = pmin < pmax;
859 
860 #if COMPARE==10
861  auto surfhit = [this, &bb, &T, &n0, &n0t](double & pmin, double & pmax, bool print = false)->bool {
862  const int N = 1000 * 1000;
863  if (fsurf.size() == 0) for (int i = 0; i < N; i++) fsurf.push_back(GetPointOnSurface());
864 
865  int umin = -1, umax = -1;
866  double wmin = 1e99, wmax = -1e99;
867  for (int i = 0; i < N; i++)
868  {
869  if (bb.Inside(T.TransformPoint(fsurf[i]))) {
870  double w = n0t * fsurf[i];
871  if (wmin > w) {wmin = w; umin = i;}
872  if (wmax < w) {wmax = w; umax = i;}
873  }
874  }
875  if (print)cout << umin << " " << umax << " " << wmin << " " << wmax << endl;
876  if (umin >= 0 && umax >= 0)
877  {
878  G4ThreeVector qmin = fsurf[umin], qmax = fsurf[umax];
879  T.ApplyPointTransform(qmin);
880  T.ApplyPointTransform(qmax);
881  pmin = n0 * qmin, pmax = n0 * qmax;
882  return true;
883  }
884  return false;
885  };
886 
887  bool res = fshape->CalculateExtent(A, bb, T, pMin, pMax);
888  double srfmin = kInfinity, srfmax = -srfmin;
889  bool shit = surfhit(srfmin, srfmax);
890  double diff = kCarTolerance;
891  diff = 10;
892  // if (abs(pmin - pMin) > diff || abs(pmax - pMax) > diff || hit != res) {
893  if ((abs(pmin - srfmin) > diff || abs(pmax - srfmax) > diff) && shit) {
894  cout << "===================================\n";
895  cout << GetName() << " " << fcache.size() << " " << fphi << " " << fdphi << " " << ftwopi << "\n";
896  cout << hit << " " << res << " " << b1 << " " << b2 << "\n";
897  if (shit) {
898  cout << "ss " << srfmin << " " << srfmax << "\n";
899  } else {
900  cout << "ss : not in bounding box" << "\n";
901  }
902  cout << "my " << pmin << " " << pmax << "\n";
903  cout << "tc " << pMin << " " << pMax << "\n";
904  cout << "df " << pmin - pMin << " " << pmax - pMax << "\n";
905  G4ThreeVector bmin(bb.GetMinXExtent(), bb.GetMinYExtent(), bb.GetMinZExtent());
906  G4ThreeVector bmax(bb.GetMaxXExtent(), bb.GetMaxYExtent(), bb.GetMaxZExtent());
907  cout << "Axis=" << A << " " << bmin << " " << bmax << " " << T << "\n";
908  cout << rp << " " << rm << "\n";
909  cout << smin << " " << smax << "\n";
910  cout << flush;
911  // _exit(0);
912  }
913  // cout<<endl;
914 #endif
915  pMin = pmin;
916  pMax = pmax;
917 
918  return hit;
919 }
920 
921 // True if (x,y) is within the shape rotation
922 inline bool BelleLathe::insector(double x, double y) const
923 {
924  double d0 = fn0x * x + fn0y * y;
925  double d1 = fn1x * x + fn1y * y;
926  bool b0 = d0 < 0, b1 = d1 < 0;
927  return fgtpi ? b0 || b1 : b0 && b1;
928 }
929 
930 int BelleLathe::wn_poly(const zr_t& r) const
931 {
932  int wn = 0;
933  vector<double>::const_iterator it = upper_bound(fz.begin(), fz.end(), r.z);
934  // cout<<r<<" "<<fz.size()<<" "<<it-fz.begin()<<endl;
935  if (it != fz.begin() && it != fz.end()) {
936  int k = it - fz.begin();
937  for (int i = findx[k - 1]; i != findx[k]; i++) {
938  const cachezr_t& s = fcache[fseg[i]];
939  double dz = r.z - s.z, dr = r.r - s.r;
940  double crs = s.dr * dz - s.dz * dr;
941  wn -= (crs > 0) - (crs < 0);
942  }
943  }
944  return wn;
945 }
946 
947 double BelleLathe::mindist(const zr_t& r) const
948 {
949  double d = kInfinity;
950  int i = 0, n = fcache.size();
951  do {
952  const cachezr_t& s = fcache[i];
953  double dz = r.z - s.z, dr = r.r - s.r;
954  double dot = s.dz * dz + s.dr * dr; // projection of the point on the segement
955  if (dot < 0) {
956  d = min(d, dz * dz + dr * dr); // distance to the first point of the segement
957  } else if (dot <= s.s2) { // point should be within the segment
958  double crs = s.dr * dz - s.dz * dr;
959  d = min(d, crs * crs * s.is2);
960  }
961  } while (++i < n);
962  d = sqrt(d);
963  d = (wn_poly(r) == 2) ? -d : d;
964  return d;
965 }
966 
967 // Return whether point inside/outside/on surface, using tolerance
968 EInside BelleLathe::Inside(const G4ThreeVector& p) const
969 {
970  COUNTER(0);
971  const double delta = 0.5 * kCarTolerance;
972  EInside res = kInside;
973  if (!ftwopi) {
974  double d0 = fn0x * p.x() + fn0y * p.y();
975  double d1 = fn1x * p.x() + fn1y * p.y();
976  if (fgtpi) {
977  if (d0 > delta && d1 > delta) { res = kOutside;}
978  else if (d0 > -delta && d1 > -delta) { res = kSurface;}
979  } else {
980  if (d0 > delta || d1 > delta) { res = kOutside;}
981  else if (d0 > -delta || d1 > -delta) { res = kSurface;}
982  }
983  }
984  if (res != kOutside) {
985  zr_t r = {p.z(), p.perp()};
986  double d = mindist(r);
987  if (res == kSurface && d < delta) res = kSurface;
988  else if (d > delta) res = kOutside;
989  else if (d > -delta) res = kSurface;
990  else res = kInside;
991  }
992 
993 #if COMPARE==1
994  EInside dd = fshape->Inside(p);
995  if (1 || dd != res) {
996  double d0 = fn0x * p.x() + fn0y * p.y();
997  double d1 = fn1x * p.x() + fn1y * p.y();
998  // if (abs(d0) > kCarTolerance && abs(d1) > kCarTolerance) {
999  int oldprec = cout.precision(16);
1000  zr_t r = {p.z(), p.perp()};
1001  cout << GetName() << " Inside(p) " << p << " " << r << " my=" << res << " tc=" << dd <<
1002  " dist=" << mindist(r) << " " << d0 << " " << d1 << endl;
1003  cout.precision(oldprec);
1004  // }
1005  }
1006 #endif
1007  MATCHOUT("BelleLathe::Inside(p) " << p << " res= " << res);
1008  return res;
1009 }
1010 
1011 zr_t BelleLathe::normal(const zr_t& r, double& d2) const
1012 {
1013  double d = std::numeric_limits<double>::infinity(), t = 0;
1014  int iseg = -1;
1015  for (int i = 0, imax = fcache.size(); i < imax; i++) {
1016  const cachezr_t& s = fcache[i];
1017  double dz = r.z - s.z, dr = r.r - s.r;
1018  double dot = s.dz * dz + s.dr * dr; // projection of the point on the segement
1019  if (dot < 0) {
1020  double dist = dz * dz + dr * dr; // distance to the first point of the segement
1021  if (dist < d) { d = dist; t = dot * s.is2; iseg = i;}
1022  } else if (dot <= s.s2) { // point should be within the segment
1023  double crs = s.dr * dz - s.dz * dr;
1024  double dist = crs * crs * s.is2;
1025  if (dist < d) { d = dist; t = dot * s.is2; iseg = i;}
1026  }
1027  }
1028  d2 = d;
1029 
1030  auto getn = [this](int i)->zr_t{
1031  int imax = fcache.size();
1032  int i0 = i;
1033  if (i == -1) i0 = imax;
1034  const cachezr_t& s = fcache[i0];
1035  double is = sqrt(s.is2);
1036  return {s.dr * is, -s.dz * is};
1037  };
1038  return getn(iseg);
1039 
1040  if (t < 0.0) {
1041  const cachezr_t& s = fcache[iseg];
1042  zr_t dist = {r.z - s.z, r.r - s.r};
1043  double dist2 = dist.z * dist.z + dist.r * dist.r;
1044  if (dist2 > 1e-18) {
1045  double q = 1 / sqrt(dist2);
1046  if (wn_poly(r) == 2) q = -q;
1047  return {dist.z * q, dist.r * q};
1048  } else {
1049  zr_t n = getn(iseg), np = getn(iseg - 1);
1050  n.z += np.z; n.r += np.r;
1051  double n2 = n.z * n.z + n.r * n.r;
1052  double q = 1 / sqrt(n2);
1053  n.z *= q, n.r *= q;
1054  return n;
1055  }
1056  }
1057  return getn(iseg);
1058 }
1059 
1060 // Calculate side nearest to p, and return normal
1061 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
1062 G4ThreeVector BelleLathe::SurfaceNormal(const G4ThreeVector& p) const
1063 {
1064  COUNTER(1);
1065 
1066  auto side = [this](const zr_t & r, double d, int iside) {
1067  double nx = (iside) ? fn1x : fn0x, ny = (iside) ? fn1y : fn0y;
1068  if (wn_poly(r) == 2) return G4ThreeVector(nx, ny, 0);
1069  double cphi = (iside) ? fc1 : fc0, sphi = (iside) ? fs1 : fc0;
1070 
1071  double d2; zr_t n = normal(r, d2);
1072  double x = cphi * n.r, y = sphi * n.r;
1073  double u = sqrt(d2);
1074  d2 += d * d;
1075  G4ThreeVector res;
1076  if (d2 > 0) {
1077  double q = 1 / sqrt(d2);
1078  double cpsi = u * q, spsi = d * q;
1079  res.set(x * cpsi - y * spsi, x * spsi + y * cpsi, n.z);
1080  }
1081  res.set(x, y, n.z);
1082  return res;
1083  };
1084 
1085  G4ThreeVector res;
1086  zr_t r = {p.z(), p.perp()};
1087  double d2; zr_t n = normal(r, d2);
1088  double d = sqrt(d2);
1089  double pt = hypot(p.x(), p.y());
1090 
1091  if (pt > 0) {
1092  double ir = n.r / pt;
1093  res = G4ThreeVector(ir * p.x(), ir * p.y(), n.z);
1094  } else
1095  res = G4ThreeVector(n.r, 0, n.z);
1096 
1097  if (!ftwopi) {
1098  double d0 = fn0x * p.x() + fn0y * p.y();
1099  double d1 = fn1x * p.x() + fn1y * p.y();
1100  zr_t r0 = {p.z(), fc0 * p.x() + fs0 * p.y()}; // projection on plane phi
1101  zr_t r1 = {p.z(), fc1 * p.x() + fs1 * p.y()}; // projection on plane phi+dphi
1102  if (fgtpi) {
1103  if (d0 > 0 && d1 > 0) { // outside sector
1104  if (d0 < d1) {
1105  res = side(r0, d0, 0); goto exit;
1106  } else {
1107  res = side(r1, -d1, 1); goto exit;
1108  }
1109  } else {// inside sector
1110  if (wn_poly(r) == 2) { // point p inside the solid
1111  if (abs(d0) < d && abs(d0) < abs(d1)) { res = G4ThreeVector(fn0x, fn0y, 0); goto exit;}
1112  if (abs(d1) < d && abs(d1) < abs(d0)) { res = G4ThreeVector(fn1x, fn1y, 0); goto exit;}
1113  }
1114  }
1115  } else {
1116  if (d0 > 0 || d1 > 0) { // outside sector
1117  if (d0 < 0) {
1118  res = side(r1, -d1, 1); goto exit;
1119  } else {
1120  res = side(r0, d0, 0); goto exit;
1121  }
1122  } else {
1123  if (wn_poly(r) == 2) { // point p inside the solid
1124  if (abs(d0) < d && abs(d0) < abs(d1)) { res = G4ThreeVector(fn0x, fn0y, 0); goto exit;}
1125  if (abs(d1) < d && abs(d1) < abs(d0)) { res = G4ThreeVector(fn1x, fn1y, 0); goto exit;}
1126  }
1127  }
1128  }
1129  }
1130 exit:
1131 #if COMPARE==1
1132  G4ThreeVector dd = fshape->SurfaceNormal(p);
1133  if ((res - dd).mag() > 1e-11) {
1134  int oldprec = cout.precision(16);
1135  EInside inside = fshape->Inside(p);
1136  cout << GetName() << " SurfaceNormal(p) " << p << " " << res << " " << dd << " " << res - dd << " " << inside << endl;
1137  cout.precision(oldprec);
1138  // _exit(0);
1139  }
1140 #endif
1141  MATCHOUT("BelleLathe::SurfaceNormal(p,n) " << p << " res= " << res);
1142  return res;
1143 }
1144 
1145 // Calculate exact shortest distance to any boundary from outside
1146 // This is the best fast estimation of the shortest distance to trap
1147 // - Returns 0 is ThreeVector inside
1148 G4double BelleLathe::DistanceToIn(const G4ThreeVector& p) const
1149 {
1150  COUNTER(2);
1151  double d = 0;
1152  // int sector = 0, plane = 0;
1153  if (ftwopi) {
1154  zr_t r = {p.z(), p.perp()};
1155  d = max(mindist(r), 0.0);
1156  } else {
1157  double d0 = fn0x * p.x() + fn0y * p.y();
1158  double d1 = fn1x * p.x() + fn1y * p.y();
1159 
1160  if (fgtpi) {
1161  if (d0 > 0 && d1 > 0) { // outside sector
1162  if (d0 < d1) {
1163  zr_t r = {p.z(), -fn0y * p.x() + fn0x * p.y()}; // projection on plane
1164  d = sqrt(pow(max(mindist(r), 0.0), 2) + d0 * d0);
1165  } else {
1166  zr_t r = {p.z(), fn1y * p.x() - fn1x * p.y()}; // projection on plane
1167  d = sqrt(pow(max(mindist(r), 0.0), 2) + d1 * d1);
1168  }
1169  } else {
1170  zr_t r = {p.z(), p.perp()};
1171  d = max(mindist(r), 0.0);
1172  }
1173  } else {
1174  if (d0 > 0 || d1 > 0) { // outside sector
1175  if (d0 < 0) {
1176  zr_t r = {p.z(), fn1y * p.x() - fn1x * p.y()}; // projection on plane
1177  d = sqrt(pow(max(mindist(r), 0.0), 2) + d1 * d1);
1178  } else {
1179  zr_t r = {p.z(), -fn0y * p.x() + fn0x * p.y()}; // projection on plane
1180  d = sqrt(pow(max(mindist(r), 0.0), 2) + d0 * d0);
1181  }
1182  } else {
1183  zr_t r = {p.z(), p.perp()};
1184  d = max(mindist(r), 0.0);
1185  }
1186  }
1187  }
1188 #if COMPARE==1
1189  // double dd = fshape->Inside(p) == 2 ? 0.0 : fshape->DistanceToIn(p);
1190  double dd = fshape->DistanceToIn(p);
1191  // if (abs(d - dd) > kCarTolerance) {
1192  if (dd > d && abs(d - dd) > kCarTolerance) {
1193  int oldprec = cout.precision(16);
1194  EInside inside = fshape->Inside(p);
1195  zr_t r = {p.z(), p.perp()};
1196  cout << GetName() << " DistanceToIn(p) " << p << " " << r << " " << d << " " << dd << " " << d - dd << " " << inside << endl;
1197  cout.precision(oldprec);
1198  // exit(0);
1199  }
1200 #endif
1201  MATCHOUT("BelleLathe::DistanceToIn(p) " << p << " res= " << d);
1202  return d;
1203 }
1204 
1205 // Calculate exact shortest distance to any boundary from inside
1206 // - Returns 0 is ThreeVector outside
1207 G4double BelleLathe::DistanceToOut(const G4ThreeVector& p) const
1208 {
1209  // return ref->DistanceToOut(p);
1210  COUNTER(3);
1211  zr_t r = {p.z(), p.perp()};
1212  double d = mindist(r);
1213  if (!ftwopi) {
1214  double d0 = fn0x * p.x() + fn0y * p.y();
1215  double d1 = fn1x * p.x() + fn1y * p.y();
1216  if (fgtpi) {
1217  d = max(d, min(d0, d1));
1218  } else {
1219  d = max(d, max(d0, d1));
1220  }
1221  }
1222  d = max(-d, 0.0);
1223 #if COMPARE==1
1224  double dd = fshape->Inside(p) == 0 ? 0.0 : fshape->DistanceToOut(p);
1225  if (abs(d - dd) > kCarTolerance) {
1226  int oldprec = cout.precision(16);
1227  zr_t r = {p.z(), p.perp()};
1228  // cout<<r<<endl;
1229  cout << GetName() << " DistanceToOut(p) " << p << " " << r << " " << d << " " << dd << " " << d - dd << endl;
1230  cout.precision(oldprec);
1231  }
1232 #endif
1233  MATCHOUT("BelleLathe::DistanceToOut(p) " << p.x() << " " << p.y() << " " << p.z() << " res= " << d);
1234  return d;
1235 }
1236 
1237 // Calculate distance to shape from outside - return kInfinity if no intersection
1238 G4double BelleLathe::DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& n) const
1239 {
1240  // return fshape->DistanceToIn(p, n);
1241  auto getnormal = [this, &p, &n](int i, double t) ->G4ThreeVector{
1242  const int imax = fcache.size();
1243  G4ThreeVector o;
1244  if (i < 0)
1245  {
1246  } else if (i < imax)
1247  {
1248  const cachezr_t& s = fcache[i];
1249  if (s.dz == 0.0) {
1250  o.setZ(copysign(1, s.dr));
1251  } else {
1252  double x = p.x() + n.x() * t;
1253  double y = p.y() + n.y() * t;
1254  double sth = s.dr * s.is, cth = -s.dz * s.is;
1255  double ir = cth / sqrt(x * x + y * y);
1256  o.set(x * ir, y * ir, sth);
1257  }
1258  } else if (i == imax)
1259  {
1260  o.setX(fn0x), o.setY(fn0y);
1261  } else {
1262  o.setX(fn1x), o.setY(fn1y);
1263  }
1264  return o;
1265  };
1266 
1267  auto hitside = [this, &p, &n](double t, const cachezr_t& s) -> bool {
1268  double z = p.z() + n.z() * t;
1269  // cout<<t<<" "<<x<<" "<<y<<" "<<z<<endl;
1270  bool k = s.zmin < z && z <= s.zmax;
1271  if (k && !ftwopi)
1272  {
1273  double x = p.x() + n.x() * t;
1274  double y = p.y() + n.y() * t;
1275  k = k && insector(x, y);
1276  }
1277  return k;
1278  };
1279 
1280  auto hitzside = [this, &p, &n](double t, const cachezr_t& s) -> bool {
1281  double x = p.x() + n.x() * t;
1282  double y = p.y() + n.y() * t;
1283  double r2 = x * x + y * y;
1284  bool k = s.r2min <= r2 && r2 < s.r2max;
1285  if (k && !ftwopi)
1286  {
1287  k = k && insector(x, y);
1288  }
1289  return k;
1290  };
1291 
1292  auto hitphi0side = [this, &p, &n](double t) -> bool {
1293  double x = p.x() + n.x() * t;
1294  double y = p.y() + n.y() * t;
1295  double r = x * fc0 + y * fs0;
1296  if (r >= frmin)
1297  {
1298  double z = p.z() + n.z() * t;
1299  zr_t zr = {z, r};
1300  return wn_poly(zr) == 2;
1301  }
1302  return false;
1303  };
1304 
1305  auto hitphi1side = [this, &p, &n](double t) -> bool {
1306  double x = p.x() + n.x() * t;
1307  double y = p.y() + n.y() * t;
1308  double r = x * fc1 + y * fs1;
1309  if (r >= frmin)
1310  {
1311  double z = p.z() + n.z() * t;
1312  zr_t zr = {z, r};
1313  return wn_poly(zr) == 2;
1314  }
1315  return false;
1316  };
1317 
1318  double tmin = kInfinity;
1319  const int imax = fcache.size();
1320  int iseg = -1, isurface = -1;
1321  const double delta = 0.5 * kCarTolerance;
1322  double inz = 1 / n.z();
1323  double nn = dotxy(n, n), np = dotxy(n, p), pp = dotxy(p, p);
1324  double pz = p.z(), pr = sqrt(pp);
1325  for (int i = 0; i < imax; i++) { // loop over sides
1326  const cachezr_t& s = fcache[i];
1327  double dz = pz - s.z, dr = pr - s.r;
1328  double d = dz * s.dr - dr * s.dz;
1329  bool surface = false;
1330  if (abs(d * s.is) < delta) {
1331  double dot = dz * s.dz + dr * s.dr;
1332  if (dot >= 0 && dot <= s.s2) {
1333  surface = true;
1334  isurface = i;
1335  }
1336  }
1337  if (s.dz == 0.0) { // z-plane
1338  if (!surface) {
1339  double t = -dz * inz;
1340  if (0 < t && t < tmin && hitzside(t, s)) { tmin = t; iseg = i;}
1341  }
1342  } else {
1343  double A, B, R2;
1344  if (s.dr == 0.0) { // cylinder
1345  double R = s.r;
1346  R2 = R * R;
1347 
1348  A = -nn;
1349  B = np;
1350  } else { // cone
1351  double taz = s.ta * dz;
1352  double R = taz + s.r;
1353  R2 = R * R;
1354 
1355  double nzta = n.z() * s.ta;
1356  A = nzta * nzta - nn;
1357  B = np - nzta * R;
1358  }
1359  double C = pp - R2;
1360  double D = B * B + C * A;
1361  if (D > 0) {
1362  // double sD = sqrt(D), iA = 1 / A;
1363  // double t0 = (B + sD) * iA, t1 = (B - sD) * iA;
1364  double sD = sqrt(D), sum = B + copysign(sD, B);
1365  double t0 = -C / sum, t1 = sum / A;
1366  if (surface) { // exclude solution on surface
1367  if (abs(t0) > abs(t1)) {
1368  if (t0 > 0 && t0 < tmin && hitside(t0, s)) { tmin = t0; iseg = i;}
1369  } else {
1370  if (t1 > 0 && t1 < tmin && hitside(t1, s)) { tmin = t1; iseg = i;}
1371  }
1372  } else {
1373  if (t0 > 0 && t0 < tmin && hitside(t0, s)) { tmin = t0; iseg = i;}
1374  if (t1 > 0 && t1 < tmin && hitside(t1, s)) { tmin = t1; iseg = i;}
1375  }
1376  }
1377  }
1378  }
1379 
1380  if (!ftwopi) {
1381  do { // side at phi0
1382  double vn = fn0x * n.x() + fn0y * n.y();
1383  if (vn < 0) {
1384  double d = fn0x * p.x() + fn0y * p.y();
1385  double t = -d / vn;
1386  if (hitphi0side(t)) {
1387  bool surface = std::abs(d) < delta;
1388  if (surface) {
1389  tmin = 0; iseg = imax + 0;
1390  } else {
1391  if (0 < t && t < tmin) {tmin = t; iseg = imax + 0;}
1392 
1393  }
1394  }
1395  }
1396  } while (0);
1397 
1398  do { // side at phi0+dphi
1399  double vn = fn1x * n.x() + fn1y * n.y();
1400  if (vn < 0) {
1401  double d = fn1x * p.x() + fn1y * p.y();
1402  double t = -d / vn;
1403  if (hitphi1side(t)) {
1404  bool surface = std::abs(d) < delta;
1405  if (surface) {
1406  tmin = 0; iseg = imax + 1;
1407  } else {
1408  if (0 < t && t < tmin) { tmin = t; iseg = imax + 1;}
1409  }
1410  }
1411  }
1412  } while (0);
1413  }
1414 
1415  if (iseg >= 0) {
1416  if (getnormal(iseg, tmin)*n > 0) tmin = 0;
1417  // if (getnormal(iseg, tmin)*n > 0) tmin = kInfinity; // mimic genericpolycone
1418  }
1419 
1420  auto convex = [this, imax](int i) -> bool{
1421  if (i < imax)
1422  return fcache[i].isconvex;
1423  else
1424  return !fgtpi;
1425  };
1426 
1427  if (tmin >= 0 && tmin < kInfinity) {
1428  if (isurface >= 0) if (convex(isurface) && getnormal(isurface, 0)*n >= 0) tmin = kInfinity;
1429  } else {
1430  if (Inside(p) == kSurface) {
1431  if (isurface >= 0) {
1432  tmin = (getnormal(isurface, 0) * n >= 0) ? kInfinity : 0; // mimic genericpolycone
1433  }
1434  }
1435  }
1436 
1437 #if COMPARE==1
1438  // double dd = fshape->Inside(p) == 2 ? 0.0 : fshape->DistanceToIn(p, n);
1439  double dd = fshape->DistanceToIn(p, n);
1440  if (abs(tmin - dd) > 1e-10) {
1441  int oldprec = cout.precision(16);
1442  EInside inside = fshape->Inside(p);
1443  cout << GetName() << " DistanceToIn(p,v) " << p << " " << n << " " << tmin << " " << dd << " " << tmin - dd << " " << inside << " "
1444  << Inside(p) << " iseg = " << iseg << " " << isurface << endl;
1445  if (isurface >= 0) cout << getnormal(isurface, 0) << endl;
1446  cout.precision(oldprec);
1447  }
1448 #endif
1449  tmin = max(0.0, tmin);
1450  MATCHOUT("BelleLathe::DistanceToIn(p,n) " << p << " " << n << " res= " << tmin);
1451  return tmin;
1452 }
1453 
1454 // Calculate distance to surface of shape from inside
1455 G4double BelleLathe::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& n,
1456  const G4bool calcNorm, G4bool* IsValid, G4ThreeVector* _n) const
1457 {
1458  // return fshape->DistanceToOut(p, n, calcNorm, IsValid, _n);
1459  auto getnormal = [this, &p, &n](int i, double t)->G4ThreeVector{
1460  const int imax = fcache.size();
1461  G4ThreeVector o;
1462  if (i < 0)
1463  {
1464  } else if (i < imax)
1465  {
1466  const cachezr_t& s = fcache[i];
1467  if (s.dz == 0.0) {
1468  o.setZ(copysign(1, s.dr));
1469  } else {
1470  double x = p.x() + n.x() * t;
1471  double y = p.y() + n.y() * t;
1472  double sth = s.dr * s.is, cth = -s.dz * s.is;
1473  double ir = cth / sqrt(x * x + y * y);
1474  o.set(x * ir, y * ir, sth);
1475  }
1476  } else if (i == imax)
1477  {
1478  o.setX(fn0x), o.setY(fn0y);
1479  } else {
1480  o.setX(fn1x), o.setY(fn1y);
1481  }
1482  return o;
1483  };
1484 
1485  double nn = dotxy(n, n), np = dotxy(n, p), pp = dotxy(p, p);
1486  auto hitside = [this, &p, &n, nn, np, pp](double t, const cachezr_t& s) -> bool {
1487  double z = p.z() + n.z() * t;
1488  // cout<<t<<" "<<x<<" "<<y<<" "<<z<<endl;
1489  double dot = n.z() * s.dr * sqrt(pp + ((np + np) + nn * t) * t) - s.dz * (np + nn * t);
1490  bool k = s.zmin < z && z <= s.zmax && dot > 0;
1491  if (k && !ftwopi)
1492  {
1493  double x = p.x() + n.x() * t;
1494  double y = p.y() + n.y() * t;
1495 
1496  k = k && insector(x, y);
1497  }
1498  return k;
1499  };
1500 
1501  auto hitzside = [this, &p, &n](double t, const cachezr_t& s) -> bool {
1502  double x = p.x() + n.x() * t;
1503  double y = p.y() + n.y() * t;
1504  double r2 = x * x + y * y;
1505  bool k = s.dr * n.z() > 0 && s.r2min <= r2 && r2 < s.r2max;
1506  if (k && !ftwopi)
1507  {
1508  k = k && insector(x, y);
1509  }
1510  return k;
1511  };
1512 
1513  auto hitphi0side = [this, &p, &n](double t) -> bool {
1514  double x = p.x() + n.x() * t;
1515  double y = p.y() + n.y() * t;
1516  double r = x * fc0 + y * fs0;
1517  if (r >= frmin)
1518  {
1519  double z = p.z() + n.z() * t;
1520  zr_t zr = {z, r};
1521  return wn_poly(zr) == 2;
1522  }
1523  return false;
1524  };
1525 
1526  auto hitphi1side = [this, &p, &n](double t) -> bool {
1527  double x = p.x() + n.x() * t;
1528  double y = p.y() + n.y() * t;
1529  double r = x * fc1 + y * fs1;
1530  if (r >= frmin)
1531  {
1532  double z = p.z() + n.z() * t;
1533  zr_t zr = {z, r};
1534  return wn_poly(zr) == 2;
1535  }
1536  return false;
1537  };
1538 
1539  COUNTER(5);
1540  double tmin = kInfinity;
1541 
1542  const int imax = fcache.size();
1543  int iseg = -1, isurface = -1;
1544 
1545  const double delta = 0.5 * kCarTolerance;
1546  double inz = 1 / n.z();
1547  double pz = p.z(), pr = sqrt(pp);
1548 
1549  for (int i = 0; i < imax; i++) {
1550  const cachezr_t& s = fcache[i];
1551 
1552  double d = (pz - s.z) * s.dr - (pr - s.r) * s.dz;
1553  bool surface = abs(d * s.is) < delta;
1554  if (surface) isurface = i;
1555  if (s.dz == 0.0) {
1556  double t = (s.z - p.z()) * inz;
1557  if (surface) {
1558  if (hitzside(t, s)) {tmin = 0; iseg = i; break;}
1559  } else {
1560  if (0 < t && t < tmin && hitzside(t, s)) {tmin = t; iseg = i;}
1561  }
1562  } else {
1563  double A, B, R2;
1564  if (s.dr == 0.0) { // cylinder
1565  double R = s.r;
1566  R2 = R * R;
1567 
1568  A = -nn;
1569  B = np;
1570  } else { // cone
1571  double taz = s.ta * (p.z() - s.z);
1572  double R = taz + s.r;
1573  R2 = R * R;
1574 
1575  double nzta = n.z() * s.ta;
1576  A = nzta * nzta - nn;
1577  B = np - nzta * R;
1578  }
1579  double C = pp - R2;
1580  double D = B * B + C * A;
1581  if (D > 0) {
1582  double sD = sqrt(D);
1583  double sum = B + copysign(sD, B);
1584  double t0 = -C / sum, t1 = sum / A;
1585  // cout<<t0<<" "<<t1<<" "<<endl;
1586  if (surface) { //exclude solution on surface
1587  if (abs(t0) < abs(t1)) {
1588  if (hitside(t0, s)) { tmin = 0; iseg = i; break;}
1589  if (0 < t1 && t1 < tmin && hitside(t1, s)) { tmin = t1; iseg = i;}
1590  } else {
1591  if (hitside(t1, s)) { tmin = 0; iseg = i; break;}
1592  if (0 < t0 && t0 < tmin && hitside(t0, s)) { tmin = t0; iseg = i;}
1593  }
1594  } else { // check both solutions
1595  if (0 < t0 && t0 < tmin && hitside(t0, s)) { tmin = t0; iseg = i;}
1596  if (0 < t1 && t1 < tmin && hitside(t1, s)) { tmin = t1; iseg = i;}
1597  }
1598  }
1599  }
1600  // cout<<i<<" "<<iseg<<" "<<sqrt(d*d*s.is2)<<" "<<tmin<<" "<<s.zmin<<" "<<s.zmax<<endl;
1601  }
1602 
1603  if (!ftwopi) {
1604  do { // side at phi0
1605  double vn = fn0x * n.x() + fn0y * n.y();
1606  if (vn > 0) {
1607  double d = fn0x * p.x() + fn0y * p.y();
1608  double t = -d / vn;
1609  if (hitphi0side(t)) {
1610  bool surface = std::abs(d) < delta;
1611  if (surface) {
1612  tmin = 0; iseg = imax + 0;
1613  } else {
1614  if (0 < t && t < tmin) {tmin = t; iseg = imax + 0;}
1615 
1616  }
1617  }
1618  }
1619  } while (0);
1620 
1621  do { // side at phi0+dphi
1622  double vn = fn1x * n.x() + fn1y * n.y();
1623  if (vn > 0) {
1624  double d = fn1x * p.x() + fn1y * p.y();
1625  double t = -d / vn;
1626  if (hitphi1side(t)) {
1627  bool surface = std::abs(d) < delta;
1628  if (surface) {
1629  tmin = 0; iseg = imax + 1;
1630  } else {
1631  if (0 < t && t < tmin) { tmin = t; iseg = imax + 1;}
1632  }
1633  }
1634  }
1635  } while (0);
1636  }
1637 
1638  auto convex = [this, imax](int i) -> bool{
1639  if (i < imax)
1640  return fcache[i].isconvex;
1641  else
1642  return !fgtpi;
1643  };
1644 
1645  if (calcNorm) {
1646  if (tmin >= 0 && tmin < kInfinity) {
1647  *_n = getnormal(iseg, tmin);
1648  *IsValid = convex(iseg);
1649  } else {
1650  if (Inside(p) == kSurface) {
1651  if (isurface >= 0) {
1652  *_n = getnormal(isurface, tmin);
1653  *IsValid = convex(isurface);
1654  tmin = 0;
1655  }
1656  } else {
1657  *IsValid = false;
1658  }
1659  }
1660  }
1661  // cout<<"tmin "<<tmin<<endl;
1662 #if COMPARE==1
1663  do {
1664  // G4ThreeVector p0(1210.555046, -292.9578965, -36.71671492);
1665  // if ((p - p0).mag() > 1e-2)continue;
1666  bool isvalid;
1667  G4ThreeVector norm;
1668  double dd = fshape->DistanceToOut(p, n, calcNorm, &isvalid, &norm);
1669  if (abs(tmin - dd) > 1e-10 || (calcNorm && *IsValid != isvalid)) {
1670  int oldprec = cout.precision(16);
1671  cout << GetName() << " DistanceToOut(p,v) p,n =" << curl_t(p) << curl_t(n) << " calcNorm=" << calcNorm
1672  << " myInside=" << Inside(p) << " tmin=" << tmin << " dd=" << dd << " d=" << tmin - dd << " iseg=" << iseg << " isurf=" << isurface
1673  << " ";
1674  if (calcNorm) cout << "myIsValid = " << *IsValid << " tIsValid=" << isvalid << " myn=" << (*_n) << " tn=" << (norm);
1675  cout << endl;
1676  cout.precision(oldprec);
1677  // _exit(0);
1678  }
1679  } while (0);
1680 #endif
1681  MATCHOUT("BelleLathe::DistanceToOut(p,n) " << p << " " << n << " res= " << tmin);
1682  return tmin;
1683 }
1684 
1685 void BelleLathe::eartrim() const
1686 {
1687  ftlist.clear();
1688  unsigned int n = fcontour.size();
1689  vector<int> indx;
1690  for (unsigned int i = 0; i < n; i++) indx.push_back(i);
1691  int count = 0;
1692  while (indx.size() > 3 && ++count < 200) {
1693  unsigned int ni = indx.size();
1694  for (unsigned int i = 0; i < ni; i++) {
1695  int i0 = indx[i], i1 = indx[(i + 1) % ni], i2 = indx[(i + 2) % ni];
1696  double nx = fcontour[i2].z - fcontour[i0].z;
1697  double ny = fcontour[i2].r - fcontour[i0].r;
1698  double d1 = nx * (fcontour[i1].r - fcontour[i0].r) - ny * (fcontour[i1].z - fcontour[i0].z);
1699  bool ear = true;
1700  for (unsigned int j = 0; j < ni - 3; j++) {
1701  int k = indx[(i + 3 + j) % ni];
1702  double d = nx * (fcontour[k].r - fcontour[i0].r) - ny * (fcontour[k].z - fcontour[i0].z);
1703  if (d * d1 > 0) {
1704  ear = false;
1705  break;
1706  }
1707  }
1708  if (ear) {
1709  triangle_t t = {i0, i1, i2};
1710  ftlist.push_back(t);
1711  indx.erase(indx.begin() + (i + 1) % ni);
1712  break;
1713  }
1714  }
1715  }
1716  if (indx.size() == 3) {
1717  triangle_t t = {indx[0], indx[1], indx[2]};
1718  ftlist.push_back(t);
1719  }
1720 }
1721 
1722 void BelleLathe::getvolarea()
1723 {
1724  double vol = 0;
1725  for (const cachezr_t& s : fcache) vol += s.dz * ((3 * s.r) * (s.r + s.dr) + s.dr * s.dr);
1726  fCubicVolume = -fdphi * vol / 6;
1727 
1728  double totalarea = 0;
1729  if (!ftwopi) {
1730  eartrim();
1731  for (const triangle_t& t : ftlist) {
1732  const zr_t& p0 = fcontour[t.i0], &p1 = fcontour[t.i1], &p2 = fcontour[t.i2];
1733  double area = (p1.z - p0.z) * (p2.r - p0.r) - (p1.r - p0.r) * (p2.z - p0.z);
1734  totalarea += abs(area);
1735  farea.push_back(totalarea);
1736  }
1737  }
1738 
1739  for (const cachezr_t& s : fcache) {
1740  double area = fdphi * (s.r + 0.5 * s.dr) * sqrt(s.dz * s.dz + s.dr * s.dr);
1741  totalarea += area;
1742  farea.push_back(totalarea);
1743  }
1744  fSurfaceArea = farea.back();
1745 }
1746 
1747 G4ThreeVector BelleLathe::GetPointOnSurface() const
1748 {
1749  auto GetPointOnTriangle = [this](const triangle_t& t)-> G4ThreeVector{
1750  // barycentric coordinates
1751  double a1 = CLHEP::RandFlat::shoot(0., 1.), a2 = CLHEP::RandFlat::shoot(0., 1.);
1752  if (a1 + a2 > 1) { a1 = 1 - a1; a2 = 1 - a2;}
1753  double a0 = 1 - (a1 + a2);
1754  const zr_t& p0 = fcontour[t.i0], &p1 = fcontour[t.i1], &p2 = fcontour[t.i2];
1755  zr_t p = {p0.z* a0 + p1.z* a1 + p2.z * a2, p0.r* a0 + p1.r* a1 + p2.r * a2};
1756  double c, s;
1757  if (CLHEP::RandFlat::shoot(0., 1.) > 0.5) // select phi side
1758  {
1759  c = -fn0y; s = fn0x;
1760  } else {
1761  c = fn1y; s = -fn1x;
1762  }
1763  G4ThreeVector r1(p.r * c, p.r * s, p.z);
1764  return r1;
1765  };
1766 
1767  double rnd = CLHEP::RandFlat::shoot(0., farea.back());
1768  std::vector<double>::const_iterator it = std::lower_bound(farea.begin(), farea.end(), rnd);
1769  unsigned int i = it - farea.begin();
1770 
1771  if (!ftwopi) {
1772  if (i < ftlist.size()) {
1773  return GetPointOnTriangle(ftlist[i]);
1774  } else {
1775  i -= ftlist.size();
1776  }
1777  }
1778 
1779  const cachezr_t& s = fcache[i];
1780  double I = 2 * s.r + s.dr;
1781  double Iw = CLHEP::RandFlat::shoot(0., I);
1782  double q = sqrt(Iw * s.dr + s.r * s.r);
1783  double t = Iw / (q + s.r);
1784  double z = s.z + s.dz * t;
1785  double r = s.r + s.dr * t;
1786  double phi = CLHEP::RandFlat::shoot(fphi, fphi + fdphi);
1787  double x = r * cos(phi), y = r * sin(phi);
1788  return G4ThreeVector(x, y, z);
1789 }
1790 
1791 // GetEntityType
1792 G4GeometryType BelleLathe::GetEntityType() const
1793 {
1794  return G4String("BelleLathe");
1795 }
1796 
1797 // Make a clone of the object
1798 G4VSolid* BelleLathe::Clone() const
1799 {
1800  return new BelleLathe(*this);
1801 }
1802 
1803 // Stream object contents to an output stream
1804 std::ostream& BelleLathe::StreamInfo(std::ostream& os) const
1805 {
1806  G4int oldprc = os.precision(16);
1807  os << "-----------------------------------------------------------\n"
1808  << " *** Dump for solid - " << GetName() << " ***\n"
1809  << " ===================================================\n"
1810  << " Solid type: BelleLathe\n"
1811  << " Contour: " << fcontour.size() << " sides, {z, r} points \n";
1812  for (int i = 0, imax = fcontour.size(); i < imax; i++) {
1813  os << fcontour[i] << ", ";
1814  }
1815  os << "\n";
1816  for (int i = 0, imax = fcontour.size(); i < imax; i++) {
1817  os << fcache[i].isconvex << ", ";
1818  }
1819  os << "\n";
1820  os << "phi0 = " << fphi << ", dphi = " << fdphi << ", Full Circle = " << (ftwopi ? "yes" : "no") << "\n";
1821  double xmin = fzmin - 0.05 * (fzmax - fzmin), xmax = fzmax + 0.05 * (fzmax - fzmin);
1822  double ymin = frmin - 0.05 * (frmax - frmin), ymax = frmax + 0.05 * (frmax - frmin);
1823  os << " BB: " << xmin << ", " << xmax << ", " << ymin << ", " << ymax << endl;
1824  os << "-----------------------------------------------------------\n";
1825  os.precision(oldprc);
1826 
1827  return os;
1828 }
1829 
1830 // Methods for visualisation
1831 void BelleLathe::DescribeYourselfTo(G4VGraphicsScene& scene) const
1832 {
1833  scene.AddSolid(*this);
1834 }
1835 
1836 // Function to define the bounding box
1837 void BelleLathe::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
1838 {
1839  std::vector<vector_t> points;
1840  vector_t point;
1841 
1842  // Placeholder vectors
1843  const double inf = std::numeric_limits<double>::infinity();
1844  G4ThreeVector minimum(inf, inf, inf), maximum(-inf, -inf, -inf);
1845 
1846  // Outer vertices
1847  if (fdphi < 2 * M_PI) { // Only axis-crossings are relevant if the shape is a full circle
1848  point.x = frmax * cos(fphi); point.y = frmax * sin(fphi);
1849  points.push_back(point);
1850  point.x = frmax * cos(fphi + fdphi); point.y = frmax * sin(fphi + fdphi);
1851  points.push_back(point);
1852  }
1853 
1854  // Inner vertices
1855  point.x = frmin * cos(fphi); point.y = frmin * sin(fphi);
1856  points.push_back(point);
1857  if (frmin != 0) { // Avoid duplicate (0,0)
1858  point.x = frmin * cos(fphi + fdphi); point.y = frmin * sin(fphi + fdphi);
1859  points.push_back(point);
1860  }
1861 
1862  // Check if shape crosses an axis
1863  // Because the inside is concave, extremums will always be at vertices
1864  // Because the outside is convex, extremums may be at an axis-crossing
1865  if (insector(0, 1)) {
1866  point.x = 0; point.y = frmax;
1867  points.push_back(point);
1868  }
1869  if (insector(-1, 0)) {
1870  point.x = -frmax; point.y = 0;
1871  points.push_back(point);
1872  }
1873  if (insector(0, -1)) {
1874  point.x = 0; point.y = -frmax;
1875  points.push_back(point);
1876  }
1877  if (insector(1, 0)) {
1878  point.x = frmax; point.y = 0;
1879  points.push_back(point);
1880  }
1881 
1882  for (std::vector<vector_t>::size_type i = 0; i != points.size(); i++) {
1883  // global min in x
1884  if (points[i].x < minimum.x())
1885  minimum.setX(points[i].x);
1886 
1887  // global min in y
1888  if (points[i].y < minimum.y())
1889  minimum.setY(points[i].y);
1890 
1891  // global max in x
1892  if (points[i].x > maximum.x())
1893  maximum.setX(points[i].x);
1894 
1895  // global max in y
1896  if (points[i].y > maximum.y())
1897  maximum.setY(points[i].y);
1898  }
1899 
1900  // Set z extremum values
1901  minimum.setZ(fzmin);
1902  maximum.setZ(fzmax);
1903 
1904  // Assign
1905  pMin = minimum;
1906  pMax = maximum;
1907 }
1908 
1909 void takePolyhedron(const HepPolyhedron& p)
1910 {
1911  int i, nnode, iNodes[5], iVis[4], iFaces[4];
1912 
1913  for (int iface = 1; iface <= p.GetNoFacets(); iface++) {
1914  p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
1915  for (i = 0; i < nnode; i++) {
1916  if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) { //G.Barrand
1917  // processor_error = 1;
1918  G4cerr
1919  << "BooleanProcessor::takePolyhedron : problem 1."
1920  << G4endl;
1921  }
1922  if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) { //G.Barrand
1923  // processor_error = 1;
1924  G4cerr
1925  << "BooleanProcessor::takePolyhedron : problem 2. "
1926  << i << " " << iFaces[i] << " " << p.GetNoFacets() << G4endl;
1927  }
1928  }
1929  }
1930 }
1931 
1932 PolyhedronBelleLathe::PolyhedronBelleLathe(const std::vector<zr_t>& v, const std::vector<triangle_t>& t, double phi, double dphi)
1933 {
1934  int nphi = GetNumberOfRotationSteps();
1935  bool twopi = abs(dphi - 2 * M_PI) < 1e-6;
1936  int n = v.size();
1937  if (twopi) {
1938  int nv = n * nphi;
1939  int nf = nv;
1940  AllocateMemory(nv, nf);
1941 
1942  auto vnum = [nphi, n](int iphi, int ip) {
1943  return (iphi % nphi) * n + (ip % n) + 1;
1944  };
1945 
1946  int fcount = 1;
1947  double dfi = dphi / nphi;
1948  for (int i = 0; i < nphi; i++) {
1949  double fi = phi + i * dfi;
1950  double cf = cos(fi), sf = sin(fi);
1951  for (int j = 0; j < n; j++) pV[vnum(i, j)].set(v[j].r * cf, v[j].r * sf, v[j].z);
1952  for (int j = 0; j < n; j++) pF[fcount++ ] = G4Facet(vnum(i, j), 0, vnum(i, j + 1), 0, vnum(i + 1, j + 1), 0, vnum(i + 1, j), 0);
1953  }
1954  } else {
1955  // cout<<"NPHI = "<<nphi<<" "<<phi<<" "<<dphi<<endl;
1956  nphi = int(nphi * (dphi / (2 * M_PI)) + 0.5);
1957  nphi = nphi > 3 ? nphi : 3;
1958 
1959  // cout<<"NPHI = "<<nphi<<endl;
1960 
1961  int nv = n * nphi;
1962  int nf = n * (nphi - 1) + 2 * t.size();
1963  AllocateMemory(nv, nf);
1964 
1965  auto vnum = [n](int iphi, int ip) {
1966  return iphi * n + (ip % n) + 1;
1967  };
1968 
1969  int fcount = 1;
1970  double dfi = dphi / (nphi - 1);
1971  for (int i = 0; i < nphi; i++) {
1972  double fi = phi + i * dfi;
1973  double cf = cos(fi), sf = sin(fi);
1974  for (int j = 0; j < n; j++) pV[vnum(i, j)].set(v[j].r * cf, v[j].r * sf, v[j].z);
1975  if (i == nphi - 1) break;
1976  for (int j = 0; j < n; j++) pF[fcount++] = G4Facet(vnum(i, j), 0, vnum(i, j + 1), 0, vnum(i + 1, j + 1), 0, vnum(i + 1, j), 0);
1977  }
1978 
1979  for (const triangle_t& k : t) pF[fcount++] = G4Facet(vnum(0, k.i0), 0, vnum(0, k.i2), 0, vnum(0, k.i1), 0, 0, 0);
1980  int i = nphi - 1;
1981  for (const triangle_t& k : t) pF[fcount++] = G4Facet(vnum(i, k.i0), 0, vnum(i, k.i1), 0, vnum(i, k.i2), 0, 0, 0);
1982 
1983  }
1984  SetReferences();
1985  // takePolyhedron(*this);
1986 }
1987 
1989 
1990 G4Polyhedron* BelleLathe::CreatePolyhedron() const
1991 {
1992  eartrim();
1993  return new PolyhedronBelleLathe(fcontour, ftlist, fphi, fdphi);
1994 }
1995 
1996 #if 0
1997 #include <immintrin.h>
1998 double mindistsimd(const zr_t& r, const vector<cachezr_t>& contour)
1999 {
2000  double d = kInfinity;
2001  int i = 0, n = contour.size();
2002  __m128d zero = _mm_set_sd(0);
2003  __m128d one = _mm_set_sd(1);
2004  double wn = 0;
2005  do {
2006  const cachezr_t& s = contour[i];
2007  double dz = r.z - s.z, dr = r.r - s.r;
2008  double crs = s.dr * dz - s.dz * dr;
2009  double dot = s.dz * dz + s.dr * dr; // projection of the point on the segement
2010  // if(s.zmin<=r.z&&r.z<s.zmax) wn -= (crs>0) - (crs<0);
2011  __m128d crssd = _mm_set_sd(crs);
2012  __m128d maskgt = _mm_cmpgt_sd(crssd, zero);
2013  __m128d masklt = _mm_cmplt_sd(crssd, zero);
2014  __m128d left = _mm_sub_sd(_mm_and_pd(maskgt, one), _mm_and_pd(masklt, one));
2015  __m128d z = _mm_set_sd(s.z);
2016  __m128d mask = _mm_and_pd(_mm_cmple_sd(_mm_set_sd(s.zmin), z), _mm_cmplt_sd(z, _mm_set_sd(s.zmax)));
2017  left = _mm_and_pd(mask, left);
2018  double du = dz * dz + dr * dr;
2019  double dv = crs * crs * s.is2;
2020 
2021  masklt = _mm_cmplt_sd(_mm_set_sd(dot), zero);
2022  maskgt = _mm_cmpgt_sd(_mm_set_sd(dot), _mm_set_sd(s.s2));
2023 
2024  __m128d uu = _mm_or_pd(_mm_and_pd(maskgt, _mm_set_sd(kInfinity)), _mm_andnot_pd(maskgt, _mm_set_sd(dv)));
2025  __m128d vv = _mm_or_pd(_mm_and_pd(masklt, _mm_set_sd(min(d, du))), _mm_andnot_pd(masklt, _mm_set_sd(min(d, uu[0]))));
2026  wn -= left[0];
2027  d = vv[0];
2028  } while (++i < n);
2029  d = sqrt(d);
2030  d = (wn == 2) ? -d : d;
2031  // cout<<wn<<" "<<d<<endl;
2032  // cout<<sqrt(dp)<<" "<<sqrt(dm)<<endl;
2033  return d;
2034 }
2035 inline int left(const zr_t& r0, const zr_t& r1, const zr_t& r)
2036 {
2037  double d = (r1.z - r0.z) * (r.r - r0.r) - (r.z - r0.z) * (r1.r - r0.r);
2038  return (d > 0) - (d < 0);
2039 }
2040 
2041 inline int checkside(const zr_t& s0, const zr_t& s1, const zr_t& r)
2042 {
2043  double zmin = min(s0.z, s1.z), zmax = max(s0.z, s1.z);
2044  if (zmin <= r.z && r.z < zmax) return left(s0, s1, r);
2045  return 0;
2046 }
2047 
2048 int wn_poly(const zr_t& r, const vector<zr_t>& contour)
2049 {
2050  int wn = 0; // the winding number counter
2051  int i = 0, n = contour.size() - 1;
2052  do {
2053  wn += checkside(contour[i], contour[i + 1], r);
2054  } while (++i < n);
2055  wn += checkside(contour[n], contour[0], r);
2056  return wn;
2057 }
2058 
2059 double mindist(const zr_t& r, const vector<zr_t>& contour)
2060 {
2061  int wn = 0;
2062  double d = kInfinity;
2063  auto dist = [&contour, &d, &r, &wn](int i0, int i1)->void {
2064  const zr_t& s0 = contour[i0], &s1 = contour[i1];
2065  double zmin = min(s0.z, s1.z), zmax = max(s0.z, s1.z);
2066  double sz = s1.z - s0.z, sr = s1.r - s0.r;
2067  double dz = r.z - s0.z, dr = r.r - s0.r;
2068  double crs = dz * sr - sz * dr;
2069  if (zmin <= r.z && r.z < zmax) wn -= (crs > 0) - (crs < 0);
2070  double dot = sz * dz + sr * dr; // projection of the point on the segement
2071  double s2 = sz * sz + sr * sr;
2072  if (dot > s2) return; // point should be within the segment
2073  if (dot < 0)
2074  {
2075  double d2 = dz * dz + dr * dr; // distance to the first point of the segement
2076  d = min(d, d2);
2077  } else{
2078  d = min(d, crs* crs / s2);
2079  }
2080  // cout<<i0<<" "<<s0.z<<" "<<s0.r<<" "<<d<<" "<<wn<<endl;
2081  };
2082  int i = 0, n = contour.size() - 1;
2083  do {dist(i, i + 1);} while (++i < n);
2084  dist(n, 0);
2085  d = sqrt(d);
2086  d = (wn == 2) ? -d : d;
2087  // cout<<wn<<" "<<d<<endl;
2088  // cout<<sqrt(dp)<<" "<<sqrt(dm)<<endl;
2089  return d;
2090 }
2091 
2092 double mindist(const zr_t& r, const vector<cachezr_t>& contour)
2093 {
2094  double d = kInfinity;
2095  int wn = 0, i = 0, n = contour.size();
2096  do {
2097  const cachezr_t& s = contour[i];
2098  double dz = r.z - s.z, dr = r.r - s.r;
2099  double crs = s.dr * dz - s.dz * dr;
2100  double dot = s.dz * dz + s.dr * dr; // projection of the point on the segement
2101  if (s.zmin <= r.z && r.z < s.zmax) wn -= (crs > 0) - (crs < 0);
2102  if (dot > s.s2) continue; // point should be within the segment
2103  if (dot < 0) {
2104  d = min(d, dz * dz + dr * dr); // distance to the first point of the segement
2105  } else {
2106  d = min(d, crs * crs * s.is2);
2107  }
2108  // cout<<i<<" "<<s.z<<" "<<s.r<<" "<<d<<" "<<wn<<endl;
2109  // cout<<i<<" "<<d<<" "<<wn<<endl;
2110  } while (++i < n);
2111  d = sqrt(d);
2112  d = (wn == 2) ? -d : d;
2113  // cout<<wn<<" "<<d<<endl;
2114  return d;
2115 }
2116 
2117 
2118 #endif
Belle2::ECL::PolyhedronBelleLathe::~PolyhedronBelleLathe
virtual ~PolyhedronBelleLathe()
destructor
Definition: BelleLathe.cc:1988
solution_t
Definition: BelleLathe.cc:381
Belle2::operator<<
std::ostream & operator<<(std::ostream &output, const IntervalOfValidity &iov)
Definition: IntervalOfValidity.cc:196
prepareAsicCrosstalkSimDB.e
e
aux.
Definition: prepareAsicCrosstalkSimDB.py:53
Belle2::ECL::BelleLathe::operator=
BelleLathe & operator=(const BelleLathe &rhs)
assignment operator
Definition: BelleLathe.cc:295
Belle2::ECL::BelleLathe::~BelleLathe
virtual ~BelleLathe()
Destructor.
Definition: BelleLathe.cc:275
Belle2::ECL::PolyhedronBelleLathe::PolyhedronBelleLathe
PolyhedronBelleLathe(const std::vector< zr_t > &, const std::vector< triangle_t > &, double, double)
constructor
Definition: BelleLathe.cc:1932
Belle2::ECL::BelleLathe::BelleLathe
BelleLathe(const G4String &pName)
Constructor for "nominal" BelleLathe whose parameters are to be set by a G4VPVParamaterisation later.
Definition: BelleLathe.cc:258
curl_t
Definition: BelleLathe.cc:58
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::BelleLathe::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Two vectors define an axis-parallel bounding box for the shape.
Definition: BelleLathe.cc:1837
Belle2::ECL::zr_t
Definition: BelleLathe.h:35
Belle2::ECL::Plane_t
Definition: BelleCrystal.h:34
prepareAsicCrosstalkSimDB.u
u
merged u1 and u2
Definition: prepareAsicCrosstalkSimDB.py:46
Belle2::ECL::BelleLathe::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
calculate the extent of the volume
Definition: BelleLathe.cc:500
Belle2::ECL::BelleLathe
Definition: BelleLathe.h:60
Belle2::ECL::vector_t
struct for a three vector
Definition: BelleLathe.h:50
Belle2::ECL::triangle_t
struct for a triangle
Definition: BelleLathe.h:53
Belle2::rn
TString rn()
Get random string.
Definition: tools.h:41
Belle2::ECL::cachezr_t
Definition: BelleLathe.h:39
Belle2::ECL::BelleLathe::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
compute the dimensions
Definition: BelleLathe.cc:335
Belle2::ECL::PolyhedronBelleLathe
Belle lathe polyhedron.
Definition: BelleLathe.h:169