1 #include "ecl/geometry/BelleLathe.h"
5 #include "G4VoxelLimits.hh"
6 #include "G4AffineTransform.hh"
8 #include "G4VPVParameterisation.hh"
10 #include "G4VGraphicsScene.hh"
12 #include "CLHEP/Random/RandFlat.h"
21 typedef int counter_t[6];
22 map<string, counter_t> counterl;
23 #define COUNTER(x) counterl[GetName()][x]++
45 inline double dotxy(
const G4ThreeVector& p,
const G4ThreeVector& n)
47 return p.x() * n.x() + p.y() * n.y();
55 return o <<
"{" << v.z <<
", " << v.r <<
"}";
60 explicit curl_t(
const G4ThreeVector& _v): v(_v) {}
65 return o <<
"{" << c.v.x() <<
", " << c.v.y() <<
", " << c.v.z() <<
"}, ";
69 BelleLathe::BelleLathe(
const G4String& pName,
double phi0,
double dphi,
const vector<zr_t>& c)
75 BelleLathe::BelleLathe(
const G4String& pName,
double phi0,
double dphi,
int n,
double* z,
double* rin,
double* rout)
79 for (
int i = 0; i < n; i++) {
80 zr_t t = {z[i], rin[i]};
83 for (
int i = n - 1; i >= 0; i--) {
84 zr_t t = {z[i], rout[i]};
88 Init(contour, phi0, dphi);
91 void BelleLathe::Init(
const vector<zr_t>& c,
double phi0,
double dphi)
93 vector<zr_t> contour = c;
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);
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);
111 auto inc = [&contour](vector<zr_t>::iterator & it) ->
void {
112 if (++it >= contour.end()) it = contour.begin();
114 auto dec = [&contour](vector<zr_t>::iterator & it) ->
void {
115 if (--it < contour.begin()) it = (++contour.rbegin()).base();
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);
125 ++it0; inc(it1); inc(it2);
130 auto isClockwise = [&contour]() ->
bool {
132 zr_t p0 = contour[0];
133 for (
int i = 1, imax = contour.size(); i < imax; i++)
135 zr_t p1 = contour[i]; sum += (p1.z - p0.z) * (p1.r + p0.r);
138 zr_t p1 = contour[0]; sum += (p1.z - p0.z) * (p1.r + p0.r);
146 std::reverse(contour.begin(), contour.end());
150 auto convexside = [
this](
cachezr_t& s,
double eps) ->
void {
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;
159 double d = a * p.r - b * p.z + cc;
160 dm = dm || (d < -eps);
161 dp = dp || (d > eps);
162 if (dm && dp) {s.isconvex =
false;
return;}
163 }
while (++it != fcontour.end());
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];
178 t.s2 = t.dz * t.dz + t.dr * t.dr;
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);
189 frmax = max(frmax, s0.r);
190 frmin = min(frmin, s0.r);
191 fzmax = max(fzmax, s0.z);
192 fzmin = min(fzmin, s0.z);
198 fdphi = std::min(2 * M_PI, fdphi);
199 fdphi = std::max(0.0, fdphi);
202 fc1 = cos(fphi + fdphi);
203 fs1 = sin(fphi + fdphi);
209 fgtpi = fdphi > M_PI;
210 ftwopi = abs(fdphi - 2 * M_PI) < kCarTolerance;
214 for (
int i = 0, n = fcontour.size(); i < n; i++) {
215 const zr_t& s = fcontour[i];
218 sort(fz.begin(), fz.end());
219 fz.erase(std::unique(fz.begin(), fz.end()), fz.end());
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++) {
226 double cc = sj.zmin, d = sj.zmax;
227 if (cc != d and b > cc and d > a) {
232 findx.push_back(fseg.size());
237 auto getpolycone = [](
const G4String & pName,
double phi0,
double dphi,
const vector<zr_t>& c) -> G4GenericPolycone* {
241 for (
int i = 0, imax = c.size(); i < imax; i++)
246 return new G4GenericPolycone(pName, phi0, dphi, c.size(), r.data(), z.data());
248 fshape = getpolycone(GetName(), phi0, dphi, fcontour);
262 Init(a, 0, 2 * M_PI);
271 Init(b, 0, 2 * M_PI);
278 cout << GetName() <<
" ";
279 for (
int i = 0; i < 6; i++) cout << counterl[GetName()][i] <<
" "; cout << endl;
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)
298 if (
this == &rhs) {
return *
this; }
301 G4CSGSolid::operator=(rhs);
304 fcontour = rhs.fcontour;
337 const G4VPhysicalVolume*)
339 G4Exception(
"BelleLathe::ComputeDimensions()",
340 "GeomSolids0001", FatalException,
341 "BelleLathe does not support Parameterisation.");
346 vector<double> quadsolve(
double a,
double b,
double c)
351 double D = b * b - a * c;
354 double sum = b + ((b > 0) ? sD : -sD);
355 double t0 = -c / sum;
356 double t1 = -sum / a;
366 inline int quadsolve(
double a,
double b,
double c,
double& t0,
double& t1)
370 double D = b * b - a * c;
373 double sum = b + ((b > 0) ? sD : -sD);
382 vector<solution_t> extremum(
double A,
double B,
double C,
double D,
double E,
double F)
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);
392 if (fpclassify(s) == FP_INFINITE)
continue;
393 double t = -(s * B + D) / (2 * A);
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);
404 if (fpclassify(t) == FP_INFINITE)
continue;
405 double s = -(2 * t * A + D) / B;
414 vector<double> BelleLathe::linecross(
const G4ThreeVector& p,
const G4ThreeVector& n)
const
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;
421 double x = p.x() + n.x() * t;
422 double y = p.y() + n.y() * t;
423 k = k && insector(x, y);
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;
435 k = k && insector(x, y);
441 double inz = 1 / n.z();
442 double nn = Belle2::ECL::dotxy(n, n), np = dotxy(n, p), pp = dotxy(p, p);
445 double t = (s.z - p.z()) * inz;
446 if (hitzside(t, s.r2min, s.r2max)) { tc.push_back(t); }
457 double taz = ta * (p.z() - s.z);
458 double R = taz + s.r;
461 double nzta = n.z() * ta;
462 A = nzta * nzta - nn;
465 double D = B * B + (pp - R2) * A;
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);
477 double d = fn0x * p.x() + fn0y * p.y();
478 double vn = fn0x * n.x() + fn0y * n.y();
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);
486 double d = fn1x * p.x() + fn1y * p.y();
487 double vn = fn1x * n.x() + fn1y * n.y();
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);
495 sort(tc.begin(), tc.end());
501 const G4VoxelLimits& bb,
502 const G4AffineTransform& T,
503 G4double& pMin, G4double& pMax)
const
505 auto maxdist = [
this](
const G4ThreeVector & n) -> G4ThreeVector {
507 int i = 0, nsize = fcache.size();
508 if (ftwopi || insector(n.x(), n.y()))
510 double nr = hypot(n.x(), n.y()), nz = n.z();
511 double dmax = -kInfinity, R = 0, Z = 0;
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);
519 r.set(R * n.x(), R * n.y(), Z);
521 double phi = fphi + 0.5 * fdphi;
522 r.set(R * cos(phi), R * sin(phi), Z);
525 double dmax = -kInfinity;
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;
532 if (dmax < d0) { r = rf; dmax = d0;}
533 if (dmax < d1) { r = rl; dmax = d1;}
534 }
while (++i < nsize);
539 struct seg_t {
int i0, i1;};
540 auto clip = [](vector<G4ThreeVector>& vlist, vector<seg_t>& slist,
const G4ThreeVector & n,
double dist) {
545 for (
const G4ThreeVector& v : vlist) d.push_back(v * n + dist);
547 for (seg_t s : slist) {
548 double prod = d[s.i0] * d[s.i1];
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);
554 s = {lone.back(), s.i1};
556 s = {s.i0, lone.back()};
558 }
else if (prod == 0) {
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);
565 if (d[s.i0] < 0)
continue;
571 int imax = -1, jmax = -1;
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];}
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]);});
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;}
596 auto PhiCrossN = [
this, clip](
const vector<Plane_t>& planes) {
599 vector<G4ThreeVector> vlist;
601 vector<G4ThreeVector> res;
603 int nsize = fcache.size();
604 vlist.reserve(nsize);
605 slist.reserve(nsize);
606 for (
int iphi = 0; iphi < 2; iphi++) {
610 double kx = iphi ? -fn0y : fn1y, ky = iphi ? fn0x : -fn1x;
615 G4ThreeVector r(kx * s.r, ky * s.r, s.z);
617 seg_t t = {i, i + 1};
619 }
while (++i < nsize - 1);
621 G4ThreeVector r(kx * s.r, ky * s.r, s.z);
623 seg_t t = {nsize - 1, 0};
628 for (
const Plane_t& p : planes) {
630 clip(vlist, slist, p.n, p.d);
633 vector<bool> bv(vlist.size(),
false);
635 for (vector<seg_t>::const_iterator it = slist.begin(); it != slist.end(); ++it) {
640 for (
unsigned int i = 0; i < vlist.size(); i++) {
641 if (!bv[i])
continue;
642 res.push_back(vlist[i]);
648 auto RCross = [
this](
const G4ThreeVector & op,
const G4ThreeVector & k,
const G4ThreeVector & u) {
651 vector<solution_t> ts;
652 int nsize = fcache.size();
658 double r0 = seg.r, z0 = seg.z, tg = seg.ta, tg2 = tg * tg;
659 double rtg = r0 * tg;
661 G4ThreeVector o(op.x(), op.y(), op.z() - z0);
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;
669 double q1 = co * q0 + rtg;
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;
678 vector<solution_t> res = extremum(F02, F11, F20, F01, F10, F00);
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) {
684 if (ftwopi || insector(p.x(), p.y()))
689 double a = -(ck2 * u2 + cu2 * k2);
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())) {
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())) {
717 }
while (++i < nsize);
721 bool b1 =
false, b2 =
false;
722 G4ThreeVector n0, n1, n2;
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;
730 double dmin1 = -kInfinity, dmax1 = kInfinity;
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;
740 double dmin2 = -kInfinity, dmax2 = kInfinity;
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;
750 G4AffineTransform iT = T.Inverse();
752 G4ThreeVector n0t = iT.TransformAxis(n0);
753 G4ThreeVector smin = n0t * kInfinity, smax = (-kInfinity) * n0t;
754 double pmin = kInfinity, pmax = -pmin;
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);
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);
764 for (
double t : dists) {
765 G4ThreeVector p = n0t * t + c0;
766 double tt = t + c0 * n0t;
768 if (pmax < tt) { pmax = tt; smax = p;}
769 if (pmin > tt) { pmin = tt; smin = p;}
772 G4ThreeVector u = c1 - c0, un = u.unit();
773 vector<solution_t> ts = RCross(c0, n0t, un);
774 double umax = u.mag();
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;
780 if (pmax < tt) { pmax = tt; smax = p;}
781 if (pmin > tt) { pmin = tt; smin = p;}
784 planes.push_back({ -un, un * c1});
787 vector<G4ThreeVector> vside = PhiCrossN(planes);
788 for (G4ThreeVector& p : vside) {
791 if (pmax < tt) { pmax = tt; smax = p;}
792 if (pmin > tt) { pmin = tt; smin = p;}
795 }
else if (b1 || b2) {
796 G4ThreeVector limits[2], u;
798 limits[0] = n1 * dmin1;
799 limits[1] = n1 * dmax1;
800 u = iT.TransformAxis(n2);
802 limits[0] = n2 * dmin2;
803 limits[1] = n2 * dmax2;
804 u = iT.TransformAxis(n1);
807 for (G4ThreeVector& c : limits) iT.ApplyPointTransform(c);
808 for (
int i = 0; i < 2; i++) {
809 vector<solution_t> ts = RCross(limits[i], n0t, u);
811 double tt = r.t + limits[i] * n0t;
812 G4ThreeVector p = n0t * r.t + u * r.s + limits[i];
814 if (pmax < tt) { pmax = tt; smax = p;}
815 if (pmin > tt) { pmin = tt; smin = p;}
819 vector<Plane_t> planes(2);
822 n = iT.TransformAxis(n1);
824 n = iT.TransformAxis(n2);
826 planes[0] = { n, -limits[0]* n};
827 planes[1] = { -n, limits[1]* n};
828 vector<G4ThreeVector> vside = PhiCrossN(planes);
830 for (G4ThreeVector& p : vside) {
834 if (pmax < tt) { pmax = tt; smax = p;}
835 if (pmin > tt) { pmin = tt; smin = p;}
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;}
844 if (bb.Inside(T.TransformPoint(rp))) {
845 double tt = rp * n0t;
846 if (pmax < tt) {pmax = tt; smax = rp;}
850 T.ApplyPointTransform(smin);
851 T.ApplyPointTransform(smax);
855 pmin -= kCarTolerance;
856 pmax += kCarTolerance;
858 bool hit = pmin < pmax;
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());
865 int umin = -1, umax = -1;
866 double wmin = 1e99, wmax = -1e99;
867 for (
int i = 0; i < N; i++)
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;}
875 if (print)cout << umin <<
" " << umax <<
" " << wmin <<
" " << wmax << endl;
876 if (umin >= 0 && umax >= 0)
878 G4ThreeVector qmin = fsurf[umin], qmax = fsurf[umax];
879 T.ApplyPointTransform(qmin);
880 T.ApplyPointTransform(qmax);
881 pmin = n0 * qmin, pmax = n0 * qmax;
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;
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";
898 cout <<
"ss " << srfmin <<
" " << srfmax <<
"\n";
900 cout <<
"ss : not in bounding box" <<
"\n";
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";
922 inline bool BelleLathe::insector(
double x,
double y)
const
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;
930 int BelleLathe::wn_poly(
const zr_t& r)
const
933 vector<double>::const_iterator it = upper_bound(fz.begin(), fz.end(), r.z);
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++) {
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);
947 double BelleLathe::mindist(
const zr_t& r)
const
949 double d = kInfinity;
950 int i = 0, n = fcache.size();
953 double dz = r.z - s.z, dr = r.r - s.r;
954 double dot = s.dz * dz + s.dr * dr;
956 d = min(d, dz * dz + dr * dr);
957 }
else if (dot <= s.s2) {
958 double crs = s.dr * dz - s.dz * dr;
959 d = min(d, crs * crs * s.is2);
963 d = (wn_poly(r) == 2) ? -d : d;
968 EInside BelleLathe::Inside(
const G4ThreeVector& p)
const
971 const double delta = 0.5 * kCarTolerance;
972 EInside res = kInside;
974 double d0 = fn0x * p.x() + fn0y * p.y();
975 double d1 = fn1x * p.x() + fn1y * p.y();
977 if (d0 > delta && d1 > delta) { res = kOutside;}
978 else if (d0 > -delta && d1 > -delta) { res = kSurface;}
980 if (d0 > delta || d1 > delta) { res = kOutside;}
981 else if (d0 > -delta || d1 > -delta) { res = kSurface;}
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;
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();
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);
1007 MATCHOUT(
"BelleLathe::Inside(p) " << p <<
" res= " << res);
1011 zr_t BelleLathe::normal(
const zr_t& r,
double& d2)
const
1013 double d = std::numeric_limits<double>::infinity(), t = 0;
1015 for (
int i = 0, imax = fcache.size(); i < imax; i++) {
1017 double dz = r.z - s.z, dr = r.r - s.r;
1018 double dot = s.dz * dz + s.dr * dr;
1020 double dist = dz * dz + dr * dr;
1021 if (dist < d) { d = dist; t = dot * s.is2; iseg = i;}
1022 }
else if (dot <= s.s2) {
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;}
1030 auto getn = [
this](
int i)->
zr_t{
1031 int imax = fcache.size();
1033 if (i == -1) i0 = imax;
1035 double is = sqrt(s.is2);
1036 return {s.dr * is, -s.dz * is};
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};
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);
1062 G4ThreeVector BelleLathe::SurfaceNormal(
const G4ThreeVector& p)
const
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;
1071 double d2;
zr_t n = normal(r, d2);
1072 double x = cphi * n.r, y = sphi * n.r;
1073 double u = sqrt(d2);
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);
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());
1092 double ir = n.r / pt;
1093 res = G4ThreeVector(ir * p.x(), ir * p.y(), n.z);
1095 res = G4ThreeVector(n.r, 0, n.z);
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()};
1101 zr_t r1 = {p.z(), fc1 * p.x() + fs1 * p.y()};
1103 if (d0 > 0 && d1 > 0) {
1105 res = side(r0, d0, 0);
goto exit;
1107 res = side(r1, -d1, 1);
goto exit;
1110 if (wn_poly(r) == 2) {
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;}
1116 if (d0 > 0 || d1 > 0) {
1118 res = side(r1, -d1, 1);
goto exit;
1120 res = side(r0, d0, 0);
goto exit;
1123 if (wn_poly(r) == 2) {
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;}
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);
1141 MATCHOUT(
"BelleLathe::SurfaceNormal(p,n) " << p <<
" res= " << res);
1148 G4double BelleLathe::DistanceToIn(
const G4ThreeVector& p)
const
1154 zr_t r = {p.z(), p.perp()};
1155 d = max(mindist(r), 0.0);
1157 double d0 = fn0x * p.x() + fn0y * p.y();
1158 double d1 = fn1x * p.x() + fn1y * p.y();
1161 if (d0 > 0 && d1 > 0) {
1163 zr_t r = {p.z(), -fn0y * p.x() + fn0x * p.y()};
1164 d = sqrt(pow(max(mindist(r), 0.0), 2) + d0 * d0);
1166 zr_t r = {p.z(), fn1y * p.x() - fn1x * p.y()};
1167 d = sqrt(pow(max(mindist(r), 0.0), 2) + d1 * d1);
1170 zr_t r = {p.z(), p.perp()};
1171 d = max(mindist(r), 0.0);
1174 if (d0 > 0 || d1 > 0) {
1176 zr_t r = {p.z(), fn1y * p.x() - fn1x * p.y()};
1177 d = sqrt(pow(max(mindist(r), 0.0), 2) + d1 * d1);
1179 zr_t r = {p.z(), -fn0y * p.x() + fn0x * p.y()};
1180 d = sqrt(pow(max(mindist(r), 0.0), 2) + d0 * d0);
1183 zr_t r = {p.z(), p.perp()};
1184 d = max(mindist(r), 0.0);
1190 double dd = fshape->DistanceToIn(p);
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);
1201 MATCHOUT(
"BelleLathe::DistanceToIn(p) " << p <<
" res= " << d);
1207 G4double BelleLathe::DistanceToOut(
const G4ThreeVector& p)
const
1211 zr_t r = {p.z(), p.perp()};
1212 double d = mindist(r);
1214 double d0 = fn0x * p.x() + fn0y * p.y();
1215 double d1 = fn1x * p.x() + fn1y * p.y();
1217 d = max(d, min(d0, d1));
1219 d = max(d, max(d0, d1));
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()};
1229 cout << GetName() <<
" DistanceToOut(p) " << p <<
" " << r <<
" " << d <<
" " << dd <<
" " << d - dd << endl;
1230 cout.precision(oldprec);
1233 MATCHOUT(
"BelleLathe::DistanceToOut(p) " << p.x() <<
" " << p.y() <<
" " << p.z() <<
" res= " << d);
1238 G4double BelleLathe::DistanceToIn(
const G4ThreeVector& p,
const G4ThreeVector& n)
const
1241 auto getnormal = [
this, &p, &n](
int i,
double t) ->G4ThreeVector{
1242 const int imax = fcache.size();
1246 }
else if (i < imax)
1250 o.setZ(copysign(1, s.dr));
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);
1258 }
else if (i == imax)
1260 o.setX(fn0x), o.setY(fn0y);
1262 o.setX(fn1x), o.setY(fn1y);
1267 auto hitside = [
this, &p, &n](
double t,
const cachezr_t& s) ->
bool {
1268 double z = p.z() + n.z() * t;
1270 bool k = s.zmin < z && z <= s.zmax;
1273 double x = p.x() + n.x() * t;
1274 double y = p.y() + n.y() * t;
1275 k = k && insector(x, y);
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;
1287 k = k && insector(x, y);
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;
1298 double z = p.z() + n.z() * t;
1300 return wn_poly(zr) == 2;
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;
1311 double z = p.z() + n.z() * t;
1313 return wn_poly(zr) == 2;
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++) {
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) {
1339 double t = -dz * inz;
1340 if (0 < t && t < tmin && hitzside(t, s)) { tmin = t; iseg = i;}
1351 double taz = s.ta * dz;
1352 double R = taz + s.r;
1355 double nzta = n.z() * s.ta;
1356 A = nzta * nzta - nn;
1360 double D = B * B + C * A;
1364 double sD = sqrt(D), sum = B + copysign(sD, B);
1365 double t0 = -C / sum, t1 = sum / A;
1367 if (abs(t0) > abs(t1)) {
1368 if (t0 > 0 && t0 < tmin && hitside(t0, s)) { tmin = t0; iseg = i;}
1370 if (t1 > 0 && t1 < tmin && hitside(t1, s)) { tmin = t1; iseg = i;}
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;}
1382 double vn = fn0x * n.x() + fn0y * n.y();
1384 double d = fn0x * p.x() + fn0y * p.y();
1386 if (hitphi0side(t)) {
1387 bool surface = std::abs(d) < delta;
1389 tmin = 0; iseg = imax + 0;
1391 if (0 < t && t < tmin) {tmin = t; iseg = imax + 0;}
1399 double vn = fn1x * n.x() + fn1y * n.y();
1401 double d = fn1x * p.x() + fn1y * p.y();
1403 if (hitphi1side(t)) {
1404 bool surface = std::abs(d) < delta;
1406 tmin = 0; iseg = imax + 1;
1408 if (0 < t && t < tmin) { tmin = t; iseg = imax + 1;}
1416 if (getnormal(iseg, tmin)*n > 0) tmin = 0;
1420 auto convex = [
this, imax](
int i) ->
bool{
1422 return fcache[i].isconvex;
1427 if (tmin >= 0 && tmin < kInfinity) {
1428 if (isurface >= 0)
if (convex(isurface) && getnormal(isurface, 0)*n >= 0) tmin = kInfinity;
1430 if (Inside(p) == kSurface) {
1431 if (isurface >= 0) {
1432 tmin = (getnormal(isurface, 0) * n >= 0) ? kInfinity : 0;
1439 double dd = fshape->DistanceToIn(p, n);
1440 if (abs(tmin - dd) > 1
e-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);
1449 tmin = max(0.0, tmin);
1450 MATCHOUT(
"BelleLathe::DistanceToIn(p,n) " << p <<
" " << n <<
" res= " << tmin);
1455 G4double BelleLathe::DistanceToOut(
const G4ThreeVector& p,
const G4ThreeVector& n,
1456 const G4bool calcNorm, G4bool* IsValid, G4ThreeVector* _n)
const
1459 auto getnormal = [
this, &p, &n](
int i,
double t)->G4ThreeVector{
1460 const int imax = fcache.size();
1464 }
else if (i < imax)
1468 o.setZ(copysign(1, s.dr));
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);
1476 }
else if (i == imax)
1478 o.setX(fn0x), o.setY(fn0y);
1480 o.setX(fn1x), o.setY(fn1y);
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;
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;
1493 double x = p.x() + n.x() * t;
1494 double y = p.y() + n.y() * t;
1496 k = k && insector(x, y);
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;
1508 k = k && insector(x, y);
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;
1519 double z = p.z() + n.z() * t;
1521 return wn_poly(zr) == 2;
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;
1532 double z = p.z() + n.z() * t;
1534 return wn_poly(zr) == 2;
1540 double tmin = kInfinity;
1542 const int imax = fcache.size();
1543 int iseg = -1, isurface = -1;
1545 const double delta = 0.5 * kCarTolerance;
1546 double inz = 1 / n.z();
1547 double pz = p.z(), pr = sqrt(pp);
1549 for (
int i = 0; i < imax; i++) {
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;
1556 double t = (s.z - p.z()) * inz;
1558 if (hitzside(t, s)) {tmin = 0; iseg = i;
break;}
1560 if (0 < t && t < tmin && hitzside(t, s)) {tmin = t; iseg = i;}
1571 double taz = s.ta * (p.z() - s.z);
1572 double R = taz + s.r;
1575 double nzta = n.z() * s.ta;
1576 A = nzta * nzta - nn;
1580 double D = B * B + C * A;
1582 double sD = sqrt(D);
1583 double sum = B + copysign(sD, B);
1584 double t0 = -C / sum, t1 = sum / A;
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;}
1591 if (hitside(t1, s)) { tmin = 0; iseg = i;
break;}
1592 if (0 < t0 && t0 < tmin && hitside(t0, s)) { tmin = t0; iseg = i;}
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;}
1605 double vn = fn0x * n.x() + fn0y * n.y();
1607 double d = fn0x * p.x() + fn0y * p.y();
1609 if (hitphi0side(t)) {
1610 bool surface = std::abs(d) < delta;
1612 tmin = 0; iseg = imax + 0;
1614 if (0 < t && t < tmin) {tmin = t; iseg = imax + 0;}
1622 double vn = fn1x * n.x() + fn1y * n.y();
1624 double d = fn1x * p.x() + fn1y * p.y();
1626 if (hitphi1side(t)) {
1627 bool surface = std::abs(d) < delta;
1629 tmin = 0; iseg = imax + 1;
1631 if (0 < t && t < tmin) { tmin = t; iseg = imax + 1;}
1638 auto convex = [
this, imax](
int i) ->
bool{
1640 return fcache[i].isconvex;
1646 if (tmin >= 0 && tmin < kInfinity) {
1647 *_n = getnormal(iseg, tmin);
1648 *IsValid = convex(iseg);
1650 if (Inside(p) == kSurface) {
1651 if (isurface >= 0) {
1652 *_n = getnormal(isurface, tmin);
1653 *IsValid = convex(isurface);
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
1674 if (calcNorm) cout <<
"myIsValid = " << *IsValid <<
" tIsValid=" << isvalid <<
" myn=" << (*_n) <<
" tn=" << (norm);
1676 cout.precision(oldprec);
1681 MATCHOUT(
"BelleLathe::DistanceToOut(p,n) " << p <<
" " << n <<
" res= " << tmin);
1685 void BelleLathe::eartrim()
const
1688 unsigned int n = fcontour.size();
1690 for (
unsigned int i = 0; i < n; i++) indx.push_back(i);
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);
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);
1710 ftlist.push_back(t);
1711 indx.erase(indx.begin() + (i + 1) % ni);
1716 if (indx.size() == 3) {
1718 ftlist.push_back(t);
1722 void BelleLathe::getvolarea()
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;
1728 double totalarea = 0;
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);
1740 double area = fdphi * (s.r + 0.5 * s.dr) * sqrt(s.dz * s.dz + s.dr * s.dr);
1742 farea.push_back(totalarea);
1744 fSurfaceArea = farea.back();
1747 G4ThreeVector BelleLathe::GetPointOnSurface()
const
1749 auto GetPointOnTriangle = [
this](
const triangle_t& t)-> G4ThreeVector{
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};
1757 if (CLHEP::RandFlat::shoot(0., 1.) > 0.5)
1759 c = -fn0y; s = fn0x;
1761 c = fn1y; s = -fn1x;
1763 G4ThreeVector r1(p.r * c, p.r * s, p.z);
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();
1772 if (i < ftlist.size()) {
1773 return GetPointOnTriangle(ftlist[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);
1792 G4GeometryType BelleLathe::GetEntityType()
const
1794 return G4String(
"BelleLathe");
1798 G4VSolid* BelleLathe::Clone()
const
1804 std::ostream& BelleLathe::StreamInfo(std::ostream& os)
const
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] <<
", ";
1816 for (
int i = 0, imax = fcontour.size(); i < imax; i++) {
1817 os << fcache[i].isconvex <<
", ";
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);
1831 void BelleLathe::DescribeYourselfTo(G4VGraphicsScene& scene)
const
1833 scene.AddSolid(*
this);
1839 std::vector<vector_t> points;
1843 const double inf = std::numeric_limits<double>::infinity();
1844 G4ThreeVector minimum(inf, inf, inf), maximum(-inf, -inf, -inf);
1847 if (fdphi < 2 * M_PI) {
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);
1855 point.x = frmin * cos(fphi); point.y = frmin * sin(fphi);
1856 points.push_back(point);
1858 point.x = frmin * cos(fphi + fdphi); point.y = frmin * sin(fphi + fdphi);
1859 points.push_back(point);
1865 if (insector(0, 1)) {
1866 point.x = 0; point.y = frmax;
1867 points.push_back(point);
1869 if (insector(-1, 0)) {
1870 point.x = -frmax; point.y = 0;
1871 points.push_back(point);
1873 if (insector(0, -1)) {
1874 point.x = 0; point.y = -frmax;
1875 points.push_back(point);
1877 if (insector(1, 0)) {
1878 point.x = frmax; point.y = 0;
1879 points.push_back(point);
1882 for (std::vector<vector_t>::size_type i = 0; i != points.size(); i++) {
1884 if (points[i].x < minimum.x())
1885 minimum.setX(points[i].x);
1888 if (points[i].y < minimum.y())
1889 minimum.setY(points[i].y);
1892 if (points[i].x > maximum.x())
1893 maximum.setX(points[i].x);
1896 if (points[i].y > maximum.y())
1897 maximum.setY(points[i].y);
1901 minimum.setZ(fzmin);
1902 maximum.setZ(fzmax);
1909 void takePolyhedron(
const HepPolyhedron& p)
1911 int i, nnode, iNodes[5], iVis[4], iFaces[4];
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()) {
1919 <<
"BooleanProcessor::takePolyhedron : problem 1."
1922 if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) {
1925 <<
"BooleanProcessor::takePolyhedron : problem 2. "
1926 << i <<
" " << iFaces[i] <<
" " << p.GetNoFacets() << G4endl;
1934 int nphi = GetNumberOfRotationSteps();
1935 bool twopi = abs(dphi - 2 * M_PI) < 1e-6;
1940 AllocateMemory(nv, nf);
1942 auto vnum = [nphi, n](
int iphi,
int ip) {
1943 return (iphi % nphi) * n + (ip % n) + 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);
1956 nphi = int(nphi * (dphi / (2 * M_PI)) + 0.5);
1957 nphi = nphi > 3 ? nphi : 3;
1962 int nf = n * (nphi - 1) + 2 * t.size();
1963 AllocateMemory(nv, nf);
1965 auto vnum = [n](
int iphi,
int ip) {
1966 return iphi * n + (ip % n) + 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);
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);
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);
1990 G4Polyhedron* BelleLathe::CreatePolyhedron()
const
1997 #include <immintrin.h>
1998 double mindistsimd(
const zr_t& r,
const vector<cachezr_t>& contour)
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);
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;
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;
2021 masklt = _mm_cmplt_sd(_mm_set_sd(dot), zero);
2022 maskgt = _mm_cmpgt_sd(_mm_set_sd(dot), _mm_set_sd(s.s2));
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]))));
2030 d = (wn == 2) ? -d : d;
2035 inline int left(
const zr_t& r0,
const zr_t& r1,
const zr_t& r)
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);
2041 inline int checkside(
const zr_t& s0,
const zr_t& s1,
const zr_t& r)
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);
2048 int wn_poly(
const zr_t& r,
const vector<zr_t>& contour)
2051 int i = 0, n = contour.size() - 1;
2053 wn += checkside(contour[i], contour[i + 1], r);
2055 wn += checkside(contour[n], contour[0], r);
2059 double mindist(
const zr_t& r,
const vector<zr_t>& contour)
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;
2071 double s2 = sz * sz + sr * sr;
2072 if (dot > s2)
return;
2075 double d2 = dz * dz + dr * dr;
2078 d = min(d, crs* crs / s2);
2082 int i = 0, n = contour.size() - 1;
2083 do {dist(i, i + 1);}
while (++i < n);
2086 d = (wn == 2) ? -d : d;
2092 double mindist(
const zr_t& r,
const vector<cachezr_t>& contour)
2094 double d = kInfinity;
2095 int wn = 0, i = 0, n = contour.size();
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;
2101 if (s.zmin <= r.z && r.z < s.zmax) wn -= (crs > 0) - (crs < 0);
2102 if (dot > s.s2)
continue;
2104 d = min(d, dz * dz + dr * dr);
2106 d = min(d, crs * crs * s.is2);
2112 d = (wn == 2) ? -d : d;