1 #include "ecl/geometry/BelleCrystal.h"
5 #include "G4VoxelLimits.hh"
6 #include "G4AffineTransform.hh"
8 #include "G4VPVParameterisation.hh"
10 #include "G4VGraphicsScene.hh"
12 #include "CLHEP/Random/RandFlat.h"
18 #define DO_QUOTE(X) #X
19 #define QUOTE(X) DO_QUOTE(X)
26 # define COUNTER(x) counter[x]++
32 # define MATCHOUT(x) if(fmatch) cout<<x<<"\n";
34 # define MATCHOUT(x) {}
38 const G4double kCoplanar_Tolerance = 1E-4;
41 const G4ThreeVector* pt)
42 : G4CSGSolid(pName), nsides(n), fPlanes(nsides), fx(2 * nsides)
45 memset(fcounter, 0,
sizeof(fcounter));
48 fmatch = GetName() == QUOTE(MATCHNAME);
52 fDz = abs(pt[2 * nsides - 1].z());
53 for (
unsigned int i = 0; i < 2 * nsides; i++) {
60 auto isConvex = [](std::vector<Point_t>::const_iterator begin, std::vector<Point_t>::const_iterator end) ->
bool {
63 for (
int i0 = 0; i0 < np; i0++)
65 int i1 = (i0 + 1) % np, i2 = (i0 + 2) % np;
66 const Point_t& r2 = *(begin + i2), &r1 = *(begin + i1), &r0 = *(begin + i0);
67 double dx1 = r2.x - r1.x, dy1 = r2.y - r1.y;
68 double dx2 = r0.x - r1.x, dy2 = r0.y - r1.y;
69 double x = dx1 * dy2 - dy1 * dx2;
70 if (i0 == 0) sign = x > 0;
72 if (sign != (x > 0))
return false;
78 auto isClockwise = [](std::vector<Point_t>::const_iterator begin, std::vector<Point_t>::const_iterator end) ->
bool {
79 std::vector<Point_t>::const_iterator it = begin;
84 Point_t r1 = *it++; sum += (r1.x - r0.x) * (r1.y + r0.y);
87 Point_t r1 = *begin; sum += (r1.x - r0.x) * (r1.y + r0.y);
91 if (!isConvex(fx.begin(), fx.begin() + nsides)) {
92 std::ostringstream message;
93 message <<
"At -z polygon is not convex: " << GetName();
94 G4Exception(
"BelleCrystal::BelleCrystal()",
"BelleCrystal", FatalException, message);
96 if (!isConvex(fx.begin() + nsides, fx.end())) {
97 std::ostringstream message;
98 message <<
"At +z polygon is not convex: " << GetName();
99 G4Exception(
"BelleCrystal::BelleCrystal()",
"BelleCrystal", FatalException, message);
102 if (isClockwise(fx.begin(), fx.begin() + nsides)) {
103 std::ostringstream message;
104 message <<
"At -z polygon is not in anti-clockwise order: " << GetName();
105 G4Exception(
"BelleCrystal::BelleCrystal()",
"BelleCrystal", FatalException, message);
107 if (isClockwise(fx.begin() + nsides, fx.end())) {
108 std::ostringstream message;
109 message <<
"At +z polygon is not in anti-clockwise order: " << GetName();
110 G4Exception(
"BelleCrystal::BelleCrystal()",
"BelleCrystal", FatalException, message);
114 const G4ThreeVector* pm = pt, *pp = pt + nsides;
115 bool zm = pm[0].z() < 0, zp = pp[0].z() > 0;
116 for (
unsigned int i = 1; i < nsides; i++) {
117 zm = zm && (pm[i - 1].z() - pm[i].z()) < kCarTolerance;
118 zp = zp && (pp[i - 1].z() - pp[i].z()) < kCarTolerance;
120 bool zpm = abs(pm[0].z() + pp[0].z()) < kCarTolerance;
121 if (!(zm && zp && zpm)) {
122 std::ostringstream message;
123 message <<
"Invalid vertice coordinates for Solid: " << GetName() <<
" " << zp <<
" " << zm <<
" " << zpm;
124 G4Exception(
"BelleCrystal::BelleCrystal()",
"BelleCrystal", FatalException, message);
127 for (
unsigned int i = 0; i < nsides; i++)
128 MakePlane(pt[i], pt[i + nsides], pt[((i + 1) % (nsides)) + nsides], pt[(i + 1) % nsides], fPlanes[i]);
130 std::ostringstream message;
131 message <<
"Wrong number of sides for Belle Crystal: " << GetName() <<
" nsides = " << nsides;
132 G4Exception(
"BelleCrystal::BelleCrystal()",
"BelleCrystal", FatalException, message);
138 : G4CSGSolid(pName), nsides(0), fDz(0)
145 : G4CSGSolid(a), nsides(0), fDz(0)
153 cout << GetName() <<
" ";
154 for (
int i = 0; i < 6; i++) cout << counter[i] <<
" "; cout << endl;
162 : G4CSGSolid(rhs), nsides(rhs.nsides), fDz(rhs.fDz), fPlanes(rhs.fPlanes), fx(rhs.fx)
170 if (
this == &rhs) {
return *
this; }
173 G4CSGSolid::operator=(rhs);
179 fPlanes = rhs.fPlanes;
190 G4bool BelleCrystal::MakePlane(
const G4ThreeVector& p1,
const G4ThreeVector& p2,
191 const G4ThreeVector& p3,
const G4ThreeVector& p4,
194 G4ThreeVector v12 = p2 - p1;
195 G4ThreeVector v13 = p3 - p1;
196 G4ThreeVector v14 = p4 - p1;
197 G4ThreeVector Vcross = v12.cross(v13);
199 if (std::fabs(Vcross.dot(v14) / (Vcross.mag()*v14.mag())) > kCoplanar_Tolerance) {
200 std::ostringstream message;
201 message <<
"Verticies are not in the same plane: " << GetName() <<
" volume =" << Vcross.dot(v14);
202 G4Exception(
"BelleCrystal::MakePlane()",
"BelleCrystal", FatalException, message);
205 double a = +(p4.y() - p2.y()) * (p3.z() - p1.z()) - (p3.y() - p1.y()) * (p4.z() - p2.z());
206 double b = -(p4.x() - p2.x()) * (p3.z() - p1.z()) + (p3.x() - p1.x()) * (p4.z() - p2.z());
207 double c = +(p4.x() - p2.x()) * (p3.y() - p1.y()) - (p3.x() - p1.x()) * (p4.y() - p2.y());
208 double sd = sqrt(a * a + b * b + c * c);
212 plane.n = G4ThreeVector(a, b, c);
213 plane.d = -(a * p1.x() + b * p1.y() + c * p1.z());
224 const G4VPhysicalVolume*)
226 G4Exception(
"BelleCrystal::ComputeDimensions()",
227 "GeomSolids0001", FatalException,
228 "BelleCrystal does not support Parameterisation.");
234 G4ThreeVector min(
const G4ThreeVector& a,
const G4ThreeVector& b)
236 return G4ThreeVector(std::min(a.x(), b.x()), std::min(a.y(), b.y()), std::min(a.z(), b.z()));
240 G4ThreeVector max(
const G4ThreeVector& a,
const G4ThreeVector& b)
242 return G4ThreeVector(std::max(a.x(), b.x()), std::max(a.y(), b.y()), std::max(a.z(), b.z()));
247 const G4VoxelLimits& bb,
248 const G4AffineTransform& pTransform,
249 G4double& pMin, G4double& pMax)
const
251 const double inf = std::numeric_limits<double>::infinity();
252 G4ThreeVector v0(inf, inf, inf), v1(-inf, -inf, -inf);
253 for (
unsigned int i = 0; i < 2 *
nsides; i++) {
254 G4ThreeVector v = pTransform.TransformPoint(vertex(i));
258 G4ThreeVector b0(bb.GetMinXExtent(), bb.GetMinYExtent(), bb.GetMinZExtent());
259 G4ThreeVector b1(bb.GetMaxXExtent(), bb.GetMaxYExtent(), bb.GetMaxZExtent());
265 case kXAxis: pMin = v0.x(); pMax = v1.x();
break;
266 case kYAxis: pMin = v0.y(); pMax = v1.y();
break;
267 case kZAxis: pMin = v0.z(); pMax = v1.z();
break;
272 if ((pMin < inf) || (pMax > -inf)) {
275 pMin -= kCarTolerance ;
276 pMax += kCarTolerance ;
290 EInside BelleCrystal::Inside(
const G4ThreeVector& p)
const
293 MATCHOUT(
"Inside(p) " << p.x() <<
" " << p.y() <<
" " << p.z());
295 double d = std::fabs(p.z()) - fDz;
298 const Plane_t& t = fPlanes[i++];
299 d = max(t.n * p + t.d, d);
301 const G4double delta = 0.5 * kCarTolerance;
306 EInside res = EInside(in);
316 G4ThreeVector BelleCrystal::SurfaceNormal(
const G4ThreeVector& p)
const
320 MATCHOUT(
"SurfaceNormal(p) " << p.x() <<
" " << p.y() <<
" " << p.z());
321 const G4double delta = 0.5 * kCarTolerance;
322 G4double safe = kInfinity;
323 unsigned int iside = 0, kside = 0;
324 vector<double> adist(
nsides + 2);
325 auto dist = [delta, &safe, &iside, &kside, &adist](
double d,
unsigned int i) ->
void {
328 iside = (d < safe) ? i : iside;
336 dist(t.n * p + t.d, i++);
339 dist(p.z() - fDz, i++);
340 dist(p.z() + fDz, i++);
347 if (adist[i++] < delta) res += fPlanes[i].n;
350 if (adist[i++] < delta) res.setZ(res.z() + 1.0);
351 if (adist[i++] < delta) res.setZ(res.z() - 1.0);
356 res = fPlanes[iside].n;
357 }
else if (iside ==
nsides) {
358 res = G4ThreeVector(0, 0, 1);
360 res = G4ThreeVector(0, 0, -1);
374 G4double BelleCrystal::DistanceToIn(
const G4ThreeVector& p)
const
378 MATCHOUT(
"DistanceToIn(p) " << p.x() <<
" " << p.y() <<
" " << p.z());
379 G4double d = std::fabs(p.z()) - fDz;
382 const Plane_t& t = fPlanes[i++];
383 d = max(t.n * p + t.d, d);
397 G4double BelleCrystal::DistanceToOut(
const G4ThreeVector& p)
const
401 MATCHOUT(
"DistanceToOut(p) " << p.x() <<
" " << p.y() <<
" " << p.z());
402 G4double d = fDz - std::fabs(p.z());
405 const Plane_t& t = fPlanes[i++];
406 d = min(-(t.n * p + t.d), d);
422 G4double BelleCrystal::DistanceToIn(
const G4ThreeVector& p,
423 const G4ThreeVector& v)
const
427 MATCHOUT(
"DistanceToIn(p,v) " << p.x() <<
" " << p.y() <<
" " << p.z() <<
" " << v.x() <<
" " << v.y() <<
" " << v.z());
428 G4double delta = 0.5 * kCarTolerance;
429 double tin = 0, tout = kInfinity;
430 auto outside = [delta, &tin, &tout](
double d,
double nv) ->
bool {
435 if (d >= -delta)
return true;
443 const Plane_t& t = fPlanes[i++];
444 if (outside(t.n * p + t.d, t.n * v)) {tin = kInfinity;
goto exit;}
447 if (outside(p.z() - fDz, v.z())) {tin = kInfinity;
goto exit;}
448 if (outside(-p.z() - fDz, -v.z())) {tin = kInfinity;
goto exit;}
451 double res = (tin > tout) ? kInfinity : tin;
465 G4double BelleCrystal::DistanceToOut(
const G4ThreeVector& p,
const G4ThreeVector& v,
466 const G4bool calcNorm, G4bool* IsValid, G4ThreeVector* n)
const
470 MATCHOUT(
"DistanceToOut(p,v) " << p.x() <<
" " << p.y() <<
" " << p.z() <<
" " << v.x() <<
" " << v.y() <<
" " << v.z() <<
" " <<
472 const G4double delta = 0.5 * kCarTolerance;
473 unsigned int iside = 10;
474 G4double lmin = kInfinity;
476 auto outside = [delta, &lmin, &iside](
double d,
double nv,
unsigned int i) ->
bool {
479 lmin = 0; iside = i;
return true;
486 if (d < nv * lmin) {lmin = d / nv; iside = i;}
491 lmin = 0; iside = i;
return true;
501 if (outside(t.n * p + t.d, t.n * v, i++))
goto exit;
504 if (outside(p.z() - fDz, v.z(), i++))
goto exit;
505 if (outside(-p.z() - fDz, -v.z(), i++))
goto exit;
511 *n = fPlanes[iside].n;
512 }
else if (iside ==
nsides) {
513 *n = G4ThreeVector(0, 0, 1);
515 *n = G4ThreeVector(0, 0, -1);
533 G4GeometryType BelleCrystal::GetEntityType()
const
535 return G4String(
"BelleCrystal");
539 G4VSolid* BelleCrystal::Clone()
const
545 std::ostream& BelleCrystal::StreamInfo(std::ostream& os)
const
547 G4int oldprc = os.precision(16);
548 os <<
"-----------------------------------------------------------\n"
549 <<
" *** Dump for solid - " << GetName() <<
" ***\n"
550 <<
" ===================================================\n"
551 <<
" Solid type: BelleCrystal\n"
553 <<
" crystal side plane equations:\n"
554 <<
"-----------------------------------------------------------\n";
555 os.precision(oldprc);
561 unsigned int ivertx4[][3] = {{0, 5, 4}, {0, 1, 5}, {1, 7, 5}, {1, 3, 7}, {3, 6, 7}, {3, 2, 6}, {2, 4, 6}, {2, 0, 4}, {0, 3, 1}, {0, 2, 3}, {4, 5, 7}, {4, 7, 6}};
562 unsigned int ivertx5[][3] = {{0, 6, 5}, {0, 1, 6}, {1, 7, 6}, {1, 2, 7}, {2, 8, 7}, {2, 3, 8}, {3, 9, 8}, {3, 4, 9}, {4, 5, 9}, {4, 0, 5}, {0, 2, 1}, {0, 3, 2}, {0, 4, 3}, {5, 6, 7}, {5, 7, 8}, {5, 8, 9}};
566 static unsigned int vert[3];
568 if (
nsides == 0)
return ivertx4[it];
580 vert[0] = 0, vert[2] = 1 + j, vert[1] = 2 + j;
590 G4ThreeVector crossplanes(
const Plane_t& p0,
const Plane_t& p1,
double Dz)
592 const G4ThreeVector& n0 = p0.n;
593 const G4ThreeVector& n1 = p1.n;
594 double A = n0.z() * Dz + p0.d;
595 double B = n1.z() * Dz + p1.d;
596 double iD = 1 / (n0.y() * n1.x() - n0.x() * n1.y());
597 double x = (A * n1.y() - B * n0.y()) * iD;
598 double y = (B * n0.x() - A * n1.x()) * iD;
600 return G4ThreeVector(x, y, z);
603 #define PACK(k0,k1,w,n) ((((k0)<<(0))+((k1)<<(w)))<<((2*w)*n))
604 G4ThreeVector BelleCrystal::vertex(
unsigned int i)
const
606 double Dz = (i <
nsides) ? -fDz : fDz;
621 return G4ThreeVector(fx[i].x, fx[i].y, Dz);
629 double a1 = CLHEP::RandFlat::shoot(0., 1.), a2 = CLHEP::RandFlat::shoot(0., 1.);
630 if (a1 + a2 > 1) { a1 = 1 - a1; a2 = 1 - a2;}
631 double a0 = 1 - (a1 + a2);
632 const unsigned int* iv =
ivertx(it);
633 G4ThreeVector p0 = vertex(iv[0]), p1 = vertex(iv[1]), p2 = vertex(iv[2]);
634 G4ThreeVector r1(p0.x()*a0 + p1.x()*a1 + p2.x()*a2,
635 p0.y()*a0 + p1.y()*a1 + p2.y()*a2,
636 p0.z()*a0 + p1.z()*a1 + p2.z()*a2);
642 const unsigned int* iv =
ivertx(it);
643 G4ThreeVector p0 = vertex(iv[0]), p1 = vertex(iv[1]), p2 = vertex(iv[2]);
645 G4ThreeVector n = (p1 - p0).cross(p2 - p0);
647 return 0.5 * n.mag();
650 double BelleCrystal::getvolarea()
const
655 double totarea = 0, totvol = 0;
656 for (
int i = 0; i < nt; i++) {
657 double v, s =
area(i, v);
660 fareas.push_back(totarea);
669 if (fCubicVolume == 0.0) {
670 fCubicVolume = getvolarea();
671 fSurfaceArea = fareas.back();
678 if (fSurfaceArea == 0.0) {
679 fCubicVolume = getvolarea();
680 fSurfaceArea = fareas.back();
686 G4ThreeVector BelleCrystal::GetPointOnSurface()
const
688 if (fareas.size() == 0) getvolarea();
689 double r = CLHEP::RandFlat::shoot(0., fareas.back());
690 std::vector<double>::const_iterator it = std::lower_bound(fareas.begin(), fareas.end(), r);
691 return (it != fareas.end()) ?
GetPointOnTriangle(it - fareas.begin()) : GetPointOnSurface();
695 void BelleCrystal::DescribeYourselfTo(G4VGraphicsScene& scene)
const
697 scene.AddSolid(*
this);
703 AllocateMemory(n, nsides + 2 * (nsides - 2));
704 for (
int i = 0; i < n; i++) pV[i + 1] = pt[i];
707 for (
int j = 0; j < nsides; j++)
708 pF[count++] = G4Facet(1 + j, 0, 1 + (j + 1) % nsides, 0, 1 + ((j + 1) % nsides) + nsides, 0, 1 + j + nsides, 0);
709 for (
int j = 0; j < nsides - 2; j++) pF[count++] = G4Facet(1 + nsides, 0, 2 + j + nsides, 0, 3 + j + nsides, 0, 0, 0);
710 for (
int j = 0; j < nsides - 2; j++) pF[count++] = G4Facet(1, 0, 3 + j, 0, 2 + j, 0, 0, 0);
716 G4Polyhedron* BelleCrystal::CreatePolyhedron()
const
719 vector<G4ThreeVector> pt(np);
720 for (
int i = 0; i < np; i++) pt[i] = vertex(i);
732 const double inf = std::numeric_limits<double>::infinity();
733 G4ThreeVector minimum(inf, inf, inf), maximum(-inf, -inf, -inf);
735 for (
int i = 0; i < np; i++) {
739 minimum = min(minimum, pt);
740 maximum = max(maximum, pt);