| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4GenericTrap implementation
27 //
28 // Authors:
29 // Tatiana Nikitina, CERN; Ivana Hrivnacova,
30 // Adapted from Root Arb8 implementation by
31 //
32 // 27.05.2024 - Evgueni Tcherniaev, complete
33 // -------------------------------------------
34
35 #include "G4GenericTrap.hh"
36
37 #if !defined(G4GEOM_USE_UGENERICTRAP)
38
39 #include <iomanip>
40
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4VoxelLimits.hh"
44 #include "G4AffineTransform.hh"
45 #include "G4BoundingEnvelope.hh"
46 #include "G4QuickRand.hh"
47
48 #include "G4GeomTools.hh"
49
50 #include "G4VGraphicsScene.hh"
51 #include "G4Polyhedron.hh"
52 #include "G4VisExtent.hh"
53
54 #include "G4AutoLock.hh"
55
56 namespace
57 {
58 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZ
59 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZ
60 }
61
62 //////////////////////////////////////////////
63 //
64 // Constructor
65 //
66 G4GenericTrap::G4GenericTrap(const G4String& n
67 const std::vector
68 : G4VSolid(name)
69 {
70 halfTolerance = 0.5*kCarTolerance;
71 CheckParameters(halfZ, vertices);
72 ComputeLateralSurfaces();
73 ComputeBoundingBox();
74 ComputeScratchLength();
75 }
76
77 //////////////////////////////////////////////
78 //
79 // Fake default constructor - sets only member
80 // for usage restri
81 G4GenericTrap::G4GenericTrap(__void__& a)
82 : G4VSolid(a)
83 {
84 }
85
86 //////////////////////////////////////////////
87 //
88 // Destructor
89
90 G4GenericTrap::~G4GenericTrap()
91 {
92 }
93
94 //////////////////////////////////////////////
95 //
96 // Copy constructor
97 //
98 G4GenericTrap::G4GenericTrap(const G4GenericTr
99 : G4VSolid(rhs),
100 halfTolerance(rhs.halfTolerance), fScratch
101 fDz(rhs.fDz), fVertices(rhs.fVertices), fI
102 fMinBBox(rhs.fMinBBox), fMaxBBox(rhs.fMaxB
103 fVisSubdivisions(rhs.fVisSubdivisions),
104 fSurfaceArea(rhs.fSurfaceArea), fCubicVolu
105 {
106 for (auto i = 0; i < 5; ++i) { fTwist[i] = r
107 for (auto i = 0; i < 4; ++i) { fDelta[i] = r
108 for (auto i = 0; i < 8; ++i) { fPlane[i] = r
109 for (auto i = 0; i < 4; ++i) { fSurf[i] = rh
110 for (auto i = 0; i < 4; ++i) { fArea[i] = rh
111 }
112
113 //////////////////////////////////////////////
114 //
115 // Assignment
116 //
117 G4GenericTrap& G4GenericTrap::operator=(const
118 {
119 // Check assignment to self
120 if (this == &rhs) { return *this; }
121
122 // Copy base class data
123 G4VSolid::operator=(rhs);
124
125 // Copy data
126 halfTolerance = rhs.halfTolerance; fScratch
127 fDz = rhs.fDz; fVertices = rhs.fVertices; fI
128 fMinBBox = rhs.fMinBBox; fMaxBBox = rhs.fMax
129 fVisSubdivisions = rhs.fVisSubdivisions;
130 fSurfaceArea = rhs.fSurfaceArea; fCubicVolum
131
132 for (auto i = 0; i < 8; ++i) { fVertices[i]
133 for (auto i = 0; i < 5; ++i) { fTwist[i] = r
134 for (auto i = 0; i < 4; ++i) { fDelta[i] = r
135 for (auto i = 0; i < 8; ++i) { fPlane[i] = r
136 for (auto i = 0; i < 4; ++i) { fSurf[i] = rh
137 for (auto i = 0; i < 4; ++i) { fArea[i] = rh
138
139 fRebuildPolyhedron = false;
140 delete fpPolyhedron; fpPolyhedron = nullptr;
141
142 return *this;
143 }
144
145 //////////////////////////////////////////////
146 //
147 // Returns position of the point (inside/outsi
148 //
149 EInside G4GenericTrap::Inside(const G4ThreeVec
150 {
151 G4double px = p.x(), py = p.y(), pz = p.z();
152
153 // intersect edges by z = pz plane
154 G4TwoVector pp[4];
155 G4double z = (pz + fDz);
156 for (auto i = 0; i < 4; ++i) { pp[i] = fVert
157
158 // estimate distance to the solid
159 G4double dist = std::abs(pz) - fDz;
160 for (auto i = 0; i < 4; ++i)
161 {
162 if (fTwist[i] == 0.)
163 {
164 G4double dd = fSurf[i].D*px + fSurf[i].E
165 dist = std::max(dd, dist);
166 }
167 else
168 {
169 auto j = (i + 1)%4;
170 G4TwoVector a = pp[i];
171 G4TwoVector b = pp[j];
172 G4double dx = b.x() - a.x();
173 G4double dy = b.y() - a.y();
174 G4double leng = std::sqrt(dx*dx + dy*dy)
175 G4double dd = (dx*(py - a.y()) - dy*(px
176 dist = std::max(dd, dist);
177 }
178 }
179 return (dist > halfTolerance) ? kOutside :
180 ((dist > -halfTolerance) ? kSurface : kIns
181 }
182
183 //////////////////////////////////////////////
184 //
185 // Return unit normal to surface at p
186 //
187 G4ThreeVector G4GenericTrap::SurfaceNormal( co
188 {
189 G4double halfToleranceSquared = halfToleranc
190 G4double px = p.x(), py = p.y(), pz = p.z();
191 G4double nx = 0, ny = 0, nz = 0;
192
193 // intersect edges by z = pz plane
194 G4TwoVector pp[4];
195 G4double tz = (pz + fDz);
196 for (auto i = 0; i < 4; ++i) { pp[i] = fVert
197
198 // check bottom and top faces
199 G4double dz = std::abs(pz) - fDz;
200 nz = std::copysign(G4double(std::abs(dz) <=
201
202 // check lateral faces
203 for (auto i = 0; i < 4; ++i)
204 {
205 if (fTwist[i] == 0.)
206 {
207 G4double dd = fSurf[i].D*px + fSurf[i].E
208 if (std::abs(dd) <= halfTolerance)
209 {
210 nx += fSurf[i].D;
211 ny += fSurf[i].E;
212 nz += fSurf[i].F;
213 }
214 }
215 else
216 {
217 auto j = (i + 1)%4;
218 G4TwoVector a = pp[i];
219 G4TwoVector b = pp[j];
220 G4double dx = b.x() - a.x();
221 G4double dy = b.y() - a.y();
222 G4double ll = dx*dx + dy*dy;
223 G4double dd = dx*(py - a.y()) - dy*(px -
224 if (dd*dd <= halfToleranceSquared*ll)
225 {
226 G4double x = fSurf[i].A*pz + fSurf[i].
227 G4double y = fSurf[i].B*pz + fSurf[i].
228 G4double z = fSurf[i].A*px + fSurf[i].
229 G4double inv = 1./std::sqrt(x*x + y*y
230 nx += x*inv;
231 ny += y*inv;
232 nz += z*inv;
233 }
234 }
235 }
236
237 // return normal
238 G4double mag2 = nx*nx + ny*ny + nz*nz;
239 if (mag2 == 0.) return ApproxSurfaceNormal(p
240 G4double mag = std::sqrt(mag2);
241 G4double inv = (mag == 1.) ? 1. : 1./mag;
242 return { nx*inv, ny*inv, nz*inv };
243 }
244
245 //////////////////////////////////////////////
246 //
247 // Find surface nearest to the point and retur
248 // Normally this method should not be called
249 //
250 G4ThreeVector G4GenericTrap::ApproxSurfaceNorm
251 {
252 #ifdef G4SPECSDEBUG
253 std::ostringstream message;
254 message.precision(16);
255 message << "Point p is not on surface of sol
256 << "Position: " << p << " is "
257 << ((Inside(p) == kInside) ? "inside
258 StreamInfo(message);
259 G4Exception("G4GenericTrap::ApproxSurfaceNor
260 JustWarning, message );
261 #endif
262 G4double px = p.x(), py = p.y(), pz = p.z();
263 G4double dist = -DBL_MAX;
264 auto iside = 0;
265
266 // intersect edges by z = pz plane
267 G4TwoVector pp[4];
268 G4double tz = (pz + fDz);
269 for (auto i = 0; i < 4; ++i) { pp[i] = fVert
270
271 // check lateral faces
272 for (auto i = 0; i < 4; ++i)
273 {
274 if (fTwist[i] == 0.)
275 {
276 G4double d = fSurf[i].D*px + fSurf[i].E*
277 if (d > dist) { dist = d; iside = i; }
278 }
279 else
280 {
281 auto j = (i + 1)%4;
282 G4TwoVector a = pp[i];
283 G4TwoVector b = pp[j];
284 G4double dx = b.x() - a.x();
285 G4double dy = b.y() - a.y();
286 G4double leng = std::sqrt(dx*dx + dy*dy)
287 G4double d = (dx*(py - a.y()) - dy*(px -
288 if (d > dist) { dist = d; iside = i; }
289 }
290 }
291 // check bottom and top faces
292 G4double distz = std::abs(pz) - fDz;
293 if (distz >= dist) return { 0., 0., std::cop
294
295 G4double x = fSurf[iside].A*pz + fSurf[iside
296 G4double y = fSurf[iside].B*pz + fSurf[iside
297 G4double z = fSurf[iside].A*px + fSurf[iside
298 G4double inv = 1./std::sqrt(x*x + y*y + z*z)
299 return { x*inv, y*inv, z*inv };
300 }
301
302 //////////////////////////////////////////////
303 //
304 // Compute distance to the surface from outsid
305 // return kInfinity if no intersection
306 //
307 G4double G4GenericTrap::DistanceToIn(const G4T
308 const G4T
309 {
310 G4double px = p.x(), py = p.y(), pz = p.z();
311 G4double vx = v.x(), vy = v.y(), vz = v.z();
312
313 // Find intersections with the bounding poly
314 //
315 if (std::abs(pz) - fDz >= -halfTolerance &&
316 G4double invz = (vz == 0) ? kInfinity : -1./
317 G4double dz = std::copysign(fDz,invz);
318 G4double xin = (pz - dz)*invz;
319 G4double xout = (pz + dz)*invz;
320
321 // Check plane faces
322 for (auto k = 0; k < 4; ++k)
323 {
324 if (fTwist[k] != 0) continue; // skip twis
325 const G4GenericTrapPlane& surf = fPlane[2*
326 G4double cosa = surf.A*vx + surf.B*vy + su
327 G4double dist = surf.A*px + surf.B*py + su
328 if (dist >= -halfTolerance)
329 {
330 if (cosa >= 0.) { return kInfinity; } //
331 G4double tmp = -dist/cosa;
332 xin = std::max(tmp, xin);
333 }
334 else
335 {
336 if (cosa > 0) { xout = std::min(-dist/co
337 if (cosa < 0) { xin = std::max(-dist/cos
338 }
339 }
340
341 // Check planes around twisted faces, each t
342 G4double tin = xin;
343 G4double tout = xout;
344 for (auto i = 0; i < 4; ++i)
345 {
346 if (fTwist[i] == 0) continue; // skip plan
347
348 // check intersection with 1st bounding pl
349 const G4GenericTrapPlane& surf1 = fPlane[2
350 G4double cosa = surf1.A*vx + surf1.B*vy +
351 G4double dist = surf1.A*px + surf1.B*py +
352 if (dist >= -halfTolerance)
353 {
354 if (cosa >= 0.) { return kInfinity; } //
355 G4double tmp = -dist/cosa;
356 tin = std::max(tmp, tin);
357 }
358 else
359 {
360 if (cosa > 0) { tout = std::min(-dist/co
361 if (cosa < 0) { tin = std::max(-dist/cos
362 }
363
364 // check intersection with 2nd bounding pl
365 const G4GenericTrapPlane& surf2 = fPlane[2
366 cosa = surf2.A*vx + surf2.B*vy + surf2.C*v
367 dist = surf2.A*px + surf2.B*py + surf2.C*p
368 if (dist >= -halfTolerance)
369 {
370 if (cosa >= 0.) { return kInfinity; } //
371 G4double tmp = -dist/cosa;
372 tin = std::max(tmp, tin);
373 }
374 else
375 {
376 if (cosa > 0) { tout = std::min(-dist/co
377 if (cosa < 0) { tin = std::max(-dist/cos
378 }
379 }
380 if (tout - tin <= halfTolerance) { return kI
381
382 // if point is outside of the bounding box
383 // then move it to the surface of the boundi
384 G4double t0 = 0., x0 = px, y0 = py, z0 = pz;
385 if (x0 < fMinBBox.x() - halfTolerance ||
386 y0 < fMinBBox.y() - halfTolerance ||
387 z0 < fMinBBox.z() - halfTolerance ||
388 x0 > fMaxBBox.x() + halfTolerance ||
389 y0 > fMaxBBox.y() + halfTolerance ||
390 z0 > fMaxBBox.z() + halfTolerance)
391 {
392 t0 = tin;
393 x0 += vx*t0;
394 y0 += vy*t0;
395 z0 += vz*t0;
396 }
397
398 // Check intersections with twisted faces
399 //
400 G4double ttin[2] = { DBL_MAX, DBL_MAX };
401 G4double ttout[2] = { tout, tout };
402
403 if (tin - xin < halfTolerance) ttin[0] = xin
404 if (xout - tout < halfTolerance) ttout[0] =
405 G4double tminimal = tin - halfTolerance;
406 G4double tmaximal = tout + halfTolerance;
407
408 constexpr G4double EPSILON = 1000.*DBL_EPSIL
409 for (auto i = 0; i < 4; ++i)
410 {
411 if (fTwist[i] == 0) continue; // skip plan
412
413 // twisted face, solve quadratic equation
414 G4double ABC = fSurf[i].A*vx + fSurf[i].B
415 G4double ABCF = fSurf[i].A*x0 + fSurf[i].B
416 G4double A = ABC*vz;
417 G4double B = 0.5*(fSurf[i].D*vx + fSurf[i]
418 G4double C = fSurf[i].D*x0 + fSurf[i].E*y0
419 if (std::abs(A) <= EPSILON)
420 {
421 // case 1: track is parallel to the surf
422 if (std::abs(B) <= EPSILON)
423 {
424 // check position of the track relativ
425 // for the reason of precision it's be
426 auto j = (i + 1)%4;
427 G4double z = (z0 + fDz);
428 G4TwoVector a = fVertices[i] + fDelta[
429 G4TwoVector b = fVertices[j] + fDelta[
430 G4double dx = b.x() - a.x();
431 G4double dy = b.y() - a.y();
432 G4double leng = std::sqrt(dx*dx + dy*d
433 G4double dist = dx*(y0 - a.y()) - dy*(
434 if (dist >= -halfTolerance*leng) { ret
435 continue;
436 }
437
438 // case 2: single root
439 G4double tmp = t0 - 0.5*C/B;
440 // compute normal at the intersection po
441 G4double x = px + vx*tmp;
442 G4double y = py + vy*tmp;
443 G4double z = pz + vz*tmp;
444 const G4GenericTrapSurface& surf = fSurf
445 G4double nx = surf.A*z + surf.D;
446 G4double ny = surf.B*z + surf.E;
447 G4double nz = surf.A*x + surf.B*y + 2.*s
448
449 if (nx*vx + ny*vy + nz*vz >= 0.) // poin
450 {
451 auto k = (i == 3) ? 0 : i + 1;
452 G4double tz = (pz + fDz);
453 G4TwoVector a = fVertices[i] + fDelta[
454 G4TwoVector b = fVertices[k] + fDelta[
455 G4double dx = b.x() - a.x();
456 G4double dy = b.y() - a.y();
457 G4double leng = std::sqrt(dx*dx + dy*d
458 G4double dist = dx*(py - a.y()) - dy*(
459 if (dist >= -halfTolerance*leng) { ret
460
461 if (tmp < tminimal || tmp > tmaximal)
462 if (std::abs(tmp - ttout[0]) < halfTol
463 if (tmp < ttout[0])
464 {
465 ttout[1] = ttout[0];
466 ttout[0] = tmp;
467 }
468 else { ttout[1] = std::min(tmp, ttout[
469 }
470 else // point is flying to inside
471 {
472 if (tmp < tminimal || tmp > tmaximal)
473 if (std::abs(tmp - ttin[0]) < halfTole
474 if (tmp < ttin[0])
475 {
476 ttin[1] = ttin[0];
477 ttin[0] = tmp;
478 }
479 else { ttin[1] = std::min(tmp, ttin[1]
480 }
481 continue;
482 }
483
484 // case 3: scratch or no intersection
485 G4double D = B*B - A*C;
486 if (D < 0.25*fScratch*fScratch*A*A)
487 {
488 if (A > 0) return kInfinity;
489 continue;
490 }
491
492 // case 4: two intersection points
493 G4double tmp = -B - std::copysign(std::sqr
494 G4double t1 = tmp/A + t0;
495 G4double t2 = C/tmp + t0;
496 G4double tsurfin = std::min(t1, t2);
497 G4double tsurfout = std::max(t1, t2);
498 if (A < 0) std::swap(tsurfin, tsurfout);
499
500 if (tsurfin >= tminimal && tsurfin <= tmax
501 {
502 if (std::abs(tsurfin - ttin[0]) >= halfT
503 {
504 if (tsurfin < ttin[0])
505 {
506 ttin[1] = ttin[0];
507 ttin[0] = tsurfin;
508 }
509 else { ttin[1] = std::min(tsurfin, tti
510 }
511 }
512 if (tsurfout >= tminimal && tsurfout <= tm
513 {
514 if (std::abs(tsurfout - ttout[0]) >= hal
515 {
516 if (tsurfout < ttout[0])
517 {
518 ttout[1] = ttout[0];
519 ttout[0] = tsurfout;
520 }
521 else { ttout[1] = std::min(tsurfout, t
522 }
523 }
524 }
525
526 // Compute distance to In
527 //
528 if (ttin[0] == DBL_MAX) { return kInfinity;
529
530 // single entry point
531 if (ttin[1] == DBL_MAX)
532 {
533 G4double distin = ttin[0];
534 G4double distout = (distin >= ttout[0] - h
535 G4double dist = (distout <= halfTolerance
536 return (dist < halfTolerance) ? 0. : dist;
537 }
538
539 // two entry points
540 if (ttin[1] < ttout[0])
541 {
542 G4double distin = ttin[1], distout = ttout
543 G4double dist = (distout <= halfTolerance
544 return (dist < halfTolerance) ? 0. : dist;
545 }
546
547 // check 1st pair of in-out points
548 G4double distin1 = ttin[0], distout1 = ttout
549 G4double dist1 = (distout1 <= halfTolerance
550 if (dist1 != kInfinity) { return (dist1 < ha
551
552 // check 2nd pair of in-out points
553 G4double distin2 = ttin[1], distout2 = ttout
554 G4double dist2 = (distout2 <= halfTolerance
555 return (dist2 < halfTolerance) ? 0. : dist2;
556 }
557
558 //////////////////////////////////////////////
559 //
560 // Estimate safety distance to the surface fro
561 //
562 G4double G4GenericTrap::DistanceToIn(const G4T
563 {
564 G4double px = p.x(), py = p.y(), pz = p.z();
565
566 // intersect edges by z = pz plane
567 G4TwoVector pp[4];
568 G4double z = (pz + fDz);
569 for (auto i = 0; i < 4; ++i) { pp[i] = fVert
570
571 // estimate distance to the solid
572 G4double dist = std::abs(pz) - fDz;
573 for (auto i = 0; i < 4; ++i)
574 {
575 if (fTwist[i] == 0.)
576 {
577 G4double dd = fSurf[i].D*px + fSurf[i].E
578 dist = std::max(dd, dist);
579 }
580 else
581 {
582 // comptute distance to the wedge (two p
583 auto j = (i + 1)%4;
584 G4double cx = pp[j].x() - pp[i].x();
585 G4double cy = pp[j].y() - pp[i].y();
586 G4double d = (pp[i].x() - px)*cy + (py -
587 G4ThreeVector na(-cy, cx, fDelta[i].x()*
588 G4ThreeVector nb(-cy, cx, fDelta[j].x()*
589 G4double amag2 = na.mag2();
590 G4double bmag2 = nb.mag2();
591 G4double distab = (amag2 > bmag2) ? d/st
592 dist = std::max(distab, dist);
593 }
594 }
595 return (dist > 0.) ? dist : 0.; // return sa
596 }
597
598 //////////////////////////////////////////////
599 //
600 // Calculate distance to surface from inside
601 //
602 G4double G4GenericTrap::DistanceToOut(const G4
603 const G4
604 const G4
605 G4
606 G4
607 {
608 G4double px = p.x(), py = p.y(), pz = p.z();
609 G4double vx = v.x(), vy = v.y(), vz = v.z();
610
611 // Check intersections with plane faces
612 //
613 if ((std::abs(pz) - fDz) >= -halfTolerance &
614 {
615 if (calcNorm)
616 {
617 *validNorm = true;
618 n->set(0., 0., std::copysign(1., pz));
619 }
620 return 0.;
621 }
622 G4double tout = (vz == 0) ? DBL_MAX : (std::
623 G4int iface = (vz < 0) ? -4 : -2; // little
624
625 for (auto i = 0; i < 4; ++i)
626 {
627 if (fTwist[i] != 0) continue; // skip twis
628 const G4GenericTrapPlane& surf = fPlane[2*
629 G4double cosa = surf.A*vx + surf.B*vy + su
630 if (cosa > 0)
631 {
632 G4double dist = surf.A*px + surf.B*py +
633 if (dist >= -halfTolerance)
634 {
635 if (calcNorm)
636 {
637 *validNorm = true;
638 n->set(surf.A, surf.B, surf.C);
639 }
640 return 0.;
641 }
642 G4double tmp = -dist/cosa;
643 if (tout > tmp) { tout = tmp; iface = i;
644 }
645 }
646
647 // Check intersections with twisted faces
648 //
649 constexpr G4double EPSILON = 1000.*DBL_EPSIL
650 for (auto i = 0; i < 4; ++i)
651 {
652 if (fTwist[i] == 0) continue; // skip plan
653
654 // twisted face, solve quadratic equation
655 const G4GenericTrapSurface& surf = fSurf[i
656 G4double ABC = surf.A*vx + surf.B*vy + su
657 G4double ABCF = surf.A*px + surf.B*py + su
658 G4double A = ABC*vz;
659 G4double B = 0.5*(surf.D*vx + surf.E*vy +
660 G4double C = surf.D*px + surf.E*py + ABCF*
661
662 if (std::abs(A) <= EPSILON)
663 {
664 // case 1: track is parallel to the surf
665 if (std::abs(B) <= EPSILON) { continue;
666
667 // case 2: single root
668 G4double tmp = -0.5*C/B;
669 // compute normal at intersection point
670 G4double x = px + vx*tmp;
671 G4double y = py + vy*tmp;
672 G4double z = pz + vz*tmp;
673 G4double nx = surf.A*z + surf.D;
674 G4double ny = surf.B*z + surf.E;
675 G4double nz = surf.A*x + surf.B*y + 2.*s
676
677 if (nx*vx + ny*vy + nz*vz > 0.) // point
678 {
679 auto k = (i + 1)%4;
680 G4double tz = (pz + fDz);
681 G4TwoVector a = fVertices[i] + fDelta[
682 G4TwoVector b = fVertices[k] + fDelta[
683 G4double dx = b.x() - a.x();
684 G4double dy = b.y() - a.y();
685 G4double leng = std::sqrt(dx*dx + dy*d
686 G4double dist = dx*(py - a.y()) - dy*(
687 if (dist >= -halfTolerance*leng) // po
688 {
689 if (calcNorm)
690 {
691 *validNorm = false;
692 G4double normx = surf.A*pz + surf.
693 G4double normy = surf.B*pz + surf.
694 G4double normz = surf.A*px + surf.
695 G4double inv = 1./std::sqrt(normx*
696 n->set(normx*inv, normy*inv, normz
697 }
698 return 0.;
699 }
700 if (tout > tmp) { tout = tmp; iface =
701 }
702 continue;
703 }
704
705 // case 3: scratch or no intersection
706 G4double D = B*B - A*C;
707 if (D < 0.25*fScratch*fScratch*A*A)
708 {
709 // check position of the point
710 auto j = (i + 1)%4;
711 G4double tz = pz + fDz;
712 G4TwoVector a = fVertices[i] + fDelta[i]
713 G4TwoVector b = fVertices[j] + fDelta[j]
714 G4double dx = b.x() - a.x();
715 G4double dy = b.y() - a.y();
716 G4double leng = std::sqrt(dx*dx + dy*dy)
717 G4double dist = dx*(py - a.y()) - dy*(px
718 if (dist <= -halfTolerance*leng) { cont
719 if (A > 0 || dist > halfTolerance*leng)
720 {
721 if (calcNorm)
722 {
723 *validNorm = false;
724 G4double nx = surf.A*pz + surf.D;
725 G4double ny = surf.B*pz + surf.E;
726 G4double nz = surf.A*px + surf.B*py
727 G4double inv = 1./std::sqrt(nx*nx +
728 n->set(nx*inv, ny*inv, nz*inv);
729 }
730 return 0.;
731 }
732 continue;
733 }
734
735 // case 4: two intersection points
736 G4double tmp = -B - std::copysign(std::sqr
737 G4double t1 = tmp / A;
738 G4double t2 = C / tmp;
739 G4double tmin = std::min(t1, t2);
740 G4double tmax = std::max(t1, t2);
741
742 if (A < 0) // concave profile: tmin(out) -
743 {
744 if (std::abs(tmax) < std::abs(tmin)) con
745 if (tmin <= halfTolerance) // point is o
746 {
747 G4double t = 0.5*(tmin + tmax);
748 G4double x = px + vx*t;
749 G4double y = py + vy*t;
750 G4double z = pz + vz*t;
751
752 auto j = (i + 1)%4;
753 G4double tz = z + fDz;
754 G4TwoVector a = fVertices[i] + fDelta[
755 G4TwoVector b = fVertices[j] + fDelta[
756 G4double dx = b.x() - a.x();
757 G4double dy = b.y() - a.y();
758 G4double leng = std::sqrt(dx*dx + dy*d
759 G4double dist = dx*(y - a.y()) - dy*(x
760 if (dist <= -halfTolerance*leng) cont
761 if (calcNorm)
762 {
763 *validNorm = false;
764 G4double nx = surf.A*pz + surf.D;
765 G4double ny = surf.B*pz + surf.E;
766 G4double nz = surf.A*px + surf.B*py
767 G4double inv = 1./std::sqrt(nx*nx +
768 n->set(nx*inv, ny*inv, nz*inv);
769 }
770 return 0.;
771 }
772 if (tout > tmin) { tout = tmin; iface =
773 }
774 else // convex profile: tmin(in) -> tmax(o
775 {
776 if (tmax < halfTolerance) // point is on
777 {
778 if (calcNorm)
779 {
780 *validNorm = false;
781 G4double nx = surf.A*pz + surf.D;
782 G4double ny = surf.B*pz + surf.E;
783 G4double nz = surf.A*px + surf.B*py
784 G4double inv = 1./std::sqrt(nx*nx +
785 n->set(nx*inv, ny*inv, nz*inv);
786 }
787 return 0.;
788 }
789 if (tout > tmax) { tout = tmax; iface =
790 }
791 }
792
793 // Compute normal, if required, and return d
794 //
795 if (tout < halfTolerance) tout = 0.;
796 if (calcNorm)
797 {
798 if (iface < 0)
799 {
800 *validNorm = true;
801 n->set(0, 0, iface + 3); // little trick
802 }
803 else
804 {
805 const G4GenericTrapSurface& surf = fSurf
806 if (fTwist[iface] == 0)
807 {
808 *validNorm = true;
809 n->set(surf.D, surf.E, surf.F);
810 }
811 else
812 {
813 *validNorm = false;
814 G4double x = px + vx*tout;
815 G4double y = py + vy*tout;
816 G4double z = pz + vz*tout;
817 G4double nx = surf.A*z + surf.D;
818 G4double ny = surf.B*z + surf.E;
819 G4double nz = surf.A*x + surf.B*y + 2.
820 G4double inv = 1./std::sqrt(nx*nx + ny
821 n->set(nx*inv, ny*inv, nz*inv);
822 }
823 }
824 }
825 return tout;
826 }
827
828 //////////////////////////////////////////////
829 //
830 // Estimate safety distance to the surface fro
831 //
832 G4double G4GenericTrap::DistanceToOut(const G4
833 {
834 G4double px = p.x(), py = p.y(), pz = p.z();
835
836 // intersect edges by z = pz plane
837 G4TwoVector pp[4];
838 G4double z = (pz + fDz);
839 for (auto i = 0; i < 4; ++i) { pp[i] = fVert
840
841 // estimate distance to the solid
842 G4double dist = std::abs(pz) - fDz;
843 for (auto i = 0; i < 4; ++i)
844 {
845 if (fTwist[i] == 0.)
846 {
847 G4double dd = fSurf[i].D*px + fSurf[i].E
848 dist = std::max(dd, dist);
849 }
850 else
851 {
852 // comptute distance to the wedge (two p
853 auto j = (i + 1)%4;
854 G4double cx = pp[j].x() - pp[i].x();
855 G4double cy = pp[j].y() - pp[i].y();
856 G4double d = (pp[i].x() - px)*cy + (py -
857 G4ThreeVector na(-cy, cx, fDelta[i].x()*
858 G4ThreeVector nb(-cy, cx, fDelta[j].x()*
859 G4double amag2 = na.mag2();
860 G4double bmag2 = nb.mag2();
861 G4double distab = (amag2 > bmag2) ? d/st
862 dist = std::max(distab, dist);
863 }
864 }
865 return (dist < 0.) ? -dist : 0.; // return s
866 }
867
868 //////////////////////////////////////////////
869 //
870 // Return bounding box
871 //
872 void G4GenericTrap::BoundingLimits(G4ThreeVect
873 G4ThreeVect
874 {
875 pMin = fMinBBox;
876 pMax = fMaxBBox;
877 }
878
879 //////////////////////////////////////////////
880 //
881 // Calculate extent under transform and specif
882 //
883 G4bool
884 G4GenericTrap::CalculateExtent(const EAxis pAx
885 const G4VoxelLi
886 const G4AffineT
887 G4double&
888 {
889 G4ThreeVector bmin, bmax;
890 G4bool exist;
891
892 // Check bounding box (bbox)
893 //
894 BoundingLimits(bmin,bmax);
895 G4BoundingEnvelope bbox(bmin,bmax);
896 #ifdef G4BBOX_EXTENT
897 return bbox.CalculateExtent(pAxis,pVoxelLimi
898 #endif
899 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
900 {
901 return exist = pMin < pMax;
902 }
903
904 // Set bounding envelope (benv) and calculat
905 //
906 // To build the bounding envelope with plane
907 // the trapezoid is subdivided in two triang
908 // duplication of vertices in the bases in a
909 // a convex polyhedron (some faces of the en
910 //
911 G4double dz = GetZHalfLength();
912 G4ThreeVectorList baseA(8), baseB(8);
913 for (G4int i = 0; i < 4; ++i)
914 {
915 G4TwoVector va = GetVertex(i);
916 G4TwoVector vb = GetVertex(i+4);
917 baseA[2*i].set(va.x(), va.y(),-dz);
918 baseB[2*i].set(vb.x(), vb.y(), dz);
919 }
920 for (G4int i = 0; i < 4; ++i)
921 {
922 G4int k1 = 2*i, k2 = (2*i + 2)%8;
923 G4double ax = (baseA[k2].x() - baseA[k1].x
924 G4double ay = (baseA[k2].y() - baseA[k1].y
925 G4double bx = (baseB[k2].x() - baseB[k1].x
926 G4double by = (baseB[k2].y() - baseB[k1].y
927 G4double znorm = ax*by - ay*bx;
928 baseA[k1+1] = (znorm < 0.0) ? baseA[k2] :
929 baseB[k1+1] = (znorm < 0.0) ? baseB[k1] :
930 }
931
932 std::vector<const G4ThreeVectorList *> polyg
933 polygons[0] = &baseA;
934 polygons[1] = &baseB;
935
936 G4BoundingEnvelope benv(bmin, bmax, polygons
937 exist = benv.CalculateExtent(pAxis,pVoxelLim
938 return exist;
939 }
940
941 //////////////////////////////////////////////
942 //
943 // Return type of the solid
944 //
945 G4GeometryType G4GenericTrap::GetEntityType()
946 {
947 return { "G4GenericTrap" };
948 }
949
950 //////////////////////////////////////////////
951 //
952 // Return true if not twisted, false otherwise
953 //
954 G4bool G4GenericTrap::IsFaceted() const
955 {
956 return (!fIsTwisted);
957 }
958
959 //////////////////////////////////////////////
960 //
961 // Make a clone of the solid
962 //
963 G4VSolid* G4GenericTrap::Clone() const
964 {
965 return new G4GenericTrap(*this);
966 }
967
968 //////////////////////////////////////////////
969 //
970 // Write parameters of the solid to the specif
971 //
972 std::ostream& G4GenericTrap::StreamInfo(std::o
973 {
974 G4long oldprc = os.precision(16);
975 os << "-------------------------------------
976 << " *** Dump for solid - " << GetName
977 << " =================================
978 << "Solid geometry type: " << GetEntityTy
979 << " half length Z: " << fDz/mm << "\n"
980 << " list of vertices:\n";
981 for (G4int i = 0; i < 8; ++i)
982 {
983 os << " #" << i << " " << fVertices[i]
984 }
985 os << "-------------------------------------
986 os.precision(oldprc);
987 return os;
988 }
989
990 //////////////////////////////////////////////
991 //
992 // Pick up a random point on the surface
993 //
994 G4ThreeVector G4GenericTrap::GetPointOnSurface
995 {
996 if (fArea[0] + fArea[1] + fArea[2] + fArea[3
997 {
998 G4AutoLock l(&lateralareaMutex);
999 fArea[0] = GetLateralFaceArea(0);
1000 fArea[1] = GetLateralFaceArea(1);
1001 fArea[2] = GetLateralFaceArea(2);
1002 fArea[3] = GetLateralFaceArea(3);
1003 l.unlock();
1004 }
1005
1006 constexpr G4int iface[6][4] =
1007 { {0,1,2,3}, {0,4,5,1}, {1,5,6,2}, {2,6,7
1008
1009 G4bool isTwisted[6] = {false};
1010 for (auto i = 0; i < 4; ++i) { isTwisted[i
1011
1012 G4double ssurf[6];
1013 G4TwoVector A = fVertices[3] - fVertices[1]
1014 G4TwoVector B = fVertices[2] - fVertices[0]
1015 G4TwoVector C = fVertices[7] - fVertices[5]
1016 G4TwoVector D = fVertices[6] - fVertices[4]
1017 ssurf[0] = (A.x()*B.y() - A.y()*B.x())*0.5;
1018 ssurf[1] = ssurf[0] + fArea[0];
1019 ssurf[2] = ssurf[1] + fArea[1];
1020 ssurf[3] = ssurf[2] + fArea[2];
1021 ssurf[4] = ssurf[3] + fArea[3];
1022 ssurf[5] = ssurf[4] + (C.x()*D.y() - C.y()*
1023
1024 G4double select = ssurf[5]*G4QuickRand();
1025 G4int k;
1026 for (k = 0; k < 5; ++k) { if (select <= ssu
1027
1028 G4int i0 = iface[k][0];
1029 G4int i1 = iface[k][1];
1030 G4int i2 = iface[k][2];
1031 G4int i3 = iface[k][3];
1032 G4ThreeVector pp[4];
1033 pp[0].set(fVertices[i0].x(), fVertices[i0].
1034 pp[1].set(fVertices[i1].x(), fVertices[i1].
1035 pp[2].set(fVertices[i2].x(), fVertices[i2].
1036 pp[3].set(fVertices[i3].x(), fVertices[i3].
1037
1038 G4ThreeVector point;
1039 if (isTwisted[k])
1040 { // twisted face, rejection sampling
1041 G4double maxlength = std::max((pp[2] - pp
1042 point = (pp[0] + pp[1] + pp[2] + pp[3])*0
1043 for (auto i = 0; i < 10000; ++i)
1044 {
1045 G4double u = G4QuickRand();
1046 G4ThreeVector p0 = pp[0] + (pp[1] - pp[
1047 G4ThreeVector p1 = pp[3] + (pp[2] - pp[
1048 G4double v = G4QuickRand()*(maxlength/(
1049 if (v >= 1.) continue;
1050 point = p0 + (p1 - p0)*v;
1051 break;
1052 }
1053 }
1054 else
1055 { // plane face
1056 G4double u = G4QuickRand();
1057 G4double v = G4QuickRand();
1058 if (u + v > 1.) { u = 1. - u; v = 1. - v;
1059 G4double ss = (((pp[2] - pp[0]).cross(pp[
1060 G4int j = (select > ssurf[k] - ss) ? 3 :
1061 point = pp[0] + (pp[2] - pp[0])*u + (pp[j
1062 }
1063 return point;
1064 }
1065
1066 /////////////////////////////////////////////
1067 //
1068 // Return volume of the solid
1069 //
1070 G4double G4GenericTrap::GetCubicVolume()
1071 {
1072 if (fCubicVolume == 0.0)
1073 {
1074 // diagonals
1075 G4TwoVector A = fVertices[3] - fVertices[
1076 G4TwoVector B = fVertices[2] - fVertices[
1077 G4TwoVector C = fVertices[7] - fVertices[
1078 G4TwoVector D = fVertices[6] - fVertices[
1079
1080 // kross products
1081 G4double AB = A.x()*B.y() - A.y()*B.x();
1082 G4double CD = C.x()*D.y() - C.y()*D.x();
1083 G4double AD = A.x()*D.y() - A.y()*D.x();
1084 G4double CB = C.x()*B.y() - C.y()*B.x();
1085
1086 fCubicVolume = ((AB + CD)/3. + (AD + CB)/
1087 }
1088 return fCubicVolume;
1089 }
1090
1091 /////////////////////////////////////////////
1092 //
1093 // Compute lateral face area
1094 //
1095 G4double G4GenericTrap::GetLateralFaceArea(G4
1096 {
1097 G4int i1 = iface, i2 = (i1 + 1)%4, i3 = i1
1098
1099 // plane face
1100 if (fTwist[iface] == 0.0)
1101 {
1102 G4ThreeVector p1(fVertices[i1].x(), fVert
1103 G4ThreeVector p2(fVertices[i2].x(), fVert
1104 G4ThreeVector p3(fVertices[i3].x(), fVert
1105 G4ThreeVector p4(fVertices[i4].x(), fVert
1106 return ((p4 - p1).cross(p3 - p2)).mag()*0
1107 }
1108
1109 // twisted face
1110 constexpr G4int NSTEP = 250;
1111 constexpr G4double dt = 1./NSTEP;
1112
1113 G4double x21 = fVertices[i2].x() - fVertice
1114 G4double y21 = fVertices[i2].y() - fVertice
1115 G4double x31 = fVertices[i3].x() - fVertice
1116 G4double y31 = fVertices[i3].y() - fVertice
1117 G4double x42 = fVertices[i4].x() - fVertice
1118 G4double y42 = fVertices[i4].y() - fVertice
1119 G4double x43 = fVertices[i4].x() - fVertice
1120 G4double y43 = fVertices[i4].y() - fVertice
1121
1122 G4double A = x21*y43 - y21*x43;
1123 G4double B0 = x21*y31 - y21*x31;
1124 G4double B1 = x42*y31 - y42*x31;
1125 G4double HH = 4*fDz*fDz;
1126 G4double invAA = 1./(A*A);
1127 G4double sqrtAA = 2.*std::abs(A);
1128 G4double invSqrtAA = 1./sqrtAA;
1129
1130 G4double area = 0.;
1131 for (G4int i = 0; i < NSTEP; ++i)
1132 {
1133 G4double t = (i + 0.5)*dt;
1134 G4double I = y21 + (y43 - y21)*t;
1135 G4double J = x21 + (x43 - x21)*t;
1136 G4double IIJJ = HH*(I*I + J*J);
1137 G4double B = B1*t + B0;
1138
1139 G4double aa = A*A;
1140 G4double bb = 2.*A*B;
1141 G4double cc = IIJJ + B*B;
1142
1143 G4double R1 = std::sqrt(aa + bb + cc);
1144 G4double R0 = std::sqrt(cc);
1145 G4double log1 = std::log(std::abs(sqrtAA*
1146 G4double log0 = std::log(std::abs(sqrtAA*
1147
1148 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0)
1149 }
1150 return area*dt;
1151 }
1152
1153 /////////////////////////////////////////////
1154 //
1155 // Return surface area of the solid
1156 //
1157 G4double G4GenericTrap::GetSurfaceArea()
1158 {
1159 if (fSurfaceArea == 0.0)
1160 {
1161 G4TwoVector A = fVertices[3] - fVertices[
1162 G4TwoVector B = fVertices[2] - fVertices[
1163 G4TwoVector C = fVertices[7] - fVertices[
1164 G4TwoVector D = fVertices[6] - fVertices[
1165 G4double S_bot = (A.x()*B.y() - A.y()*B.x
1166 G4double S_top = (C.x()*D.y() - C.y()*D.x
1167 fArea[0] = GetLateralFaceArea(0);
1168 fArea[1] = GetLateralFaceArea(1);
1169 fArea[2] = GetLateralFaceArea(2);
1170 fArea[3] = GetLateralFaceArea(3);
1171 fSurfaceArea = S_bot + S_top + fArea[0] +
1172 }
1173 return fSurfaceArea;
1174 }
1175
1176 /////////////////////////////////////////////
1177 //
1178 // Check parameters of the solid
1179 //
1180 void
1181 G4GenericTrap::CheckParameters(G4double halfZ
1182 const std::vec
1183 {
1184 // Check number of vertices
1185 //
1186 if (vertices.size() != 8)
1187 {
1188 std::ostringstream message;
1189 message << "Number of vertices is " << ve
1190 << " instead of 8 for Solid: " <<
1191 G4Exception("G4GenericTrap::CheckParamete
1192 FatalException, message);
1193 }
1194
1195 // Check dZ
1196 //
1197 if ((fDz = halfZ) < 2.*kCarTolerance)
1198 {
1199 std::ostringstream message;
1200 message << "Z-dimension is too small or n
1201 << " for Solid: " << GetName() <<
1202 G4Exception("G4GenericTrap::CheckParamete
1203 FatalException, message);
1204 }
1205
1206 // Check order of vertices
1207 //
1208 for (auto i = 0; i < 8; ++i) { fVertices.at
1209
1210 // Bottom polygon area
1211 G4TwoVector diag1 = fVertices[2] - fVertice
1212 G4TwoVector diag2 = fVertices[3] - fVertice
1213 G4double ldiagbot = std::max(diag1.mag(), d
1214 G4double zbot = diag1.x()*diag2.y() - diag1
1215 if (std::abs(zbot) < ldiagbot*kCarTolerance
1216
1217 // Top polygon area
1218 G4TwoVector diag3 = fVertices[6] - fVertice
1219 G4TwoVector diag4 = fVertices[7] - fVertice
1220 G4double ldiagtop = std::max(diag3.mag(), d
1221 G4double ztop = diag3.x()*diag4.y() - diag3
1222 if (std::abs(ztop) < ldiagtop*kCarTolerance
1223
1224 if (zbot*ztop < 0.)
1225 {
1226 std::ostringstream message;
1227 message << "Vertices of the bottom and to
1228 StreamInfo(message);
1229 G4Exception("G4GenericTrap::CheckParamete
1230 FatalException, message);
1231 }
1232 if ((zbot > 0.) || (ztop > 0.))
1233 {
1234 std::swap(fVertices[1], fVertices[3]);
1235 std::swap(fVertices[5], fVertices[7]);
1236 std::ostringstream message;
1237 message << "Vertices re-ordered in Solid:
1238 << "Vertices of the bottom and to
1239 G4Exception("G4GenericTrap::CheckParamete
1240 JustWarning, message);
1241 }
1242
1243 // Check for degeneracy
1244 //
1245 G4int ndegenerated = 0;
1246 for (auto i = 0; i < 4; ++i)
1247 {
1248 auto j = (i + 1)%4;
1249 G4double lbot = (fVertices[j] - fVertices
1250 G4double ltop = (fVertices[j + 4] - fVert
1251 ndegenerated += (std::max(lbot, ltop) < k
1252 }
1253 if (ndegenerated > 1 ||
1254 GetCubicVolume() < fDz*std::max(ldiagbo
1255 {
1256 std::ostringstream message;
1257 message << "Degenerated solid !\n";
1258 StreamInfo(message);
1259 G4Exception("G4GenericTrap::CheckParamete
1260 FatalException, message);
1261 }
1262
1263 // Check that the polygons are convex
1264 //
1265 G4bool isConvex = true;
1266 for (auto i = 0; i < 4; ++i)
1267 {
1268 auto j = (i + 1)%4;
1269 auto k = (j + 1)%4;
1270 G4TwoVector edge1 = fVertices[j] - fVerti
1271 G4TwoVector edge2 = fVertices[k] - fVerti
1272 isConvex = ((edge1.x()*edge2.y() - edge1.
1273 if (!isConvex) break;
1274 G4TwoVector edge3 = fVertices[j + 4] - fV
1275 G4TwoVector edge4 = fVertices[k + 4] - fV
1276 isConvex = ((edge3.x()*edge4.y() - edge3.
1277 if (!isConvex) break;
1278 }
1279 if (!isConvex)
1280 {
1281 std::ostringstream message;
1282 message << "The bottom and top faces must
1283 StreamInfo(message);
1284 G4Exception("G4GenericTrap::CheckParamete
1285 FatalException, message);
1286 }
1287 }
1288
1289 /////////////////////////////////////////////
1290 //
1291 // Compute surface equations and twist angles
1292 //
1293 void G4GenericTrap::ComputeLateralSurfaces()
1294 {
1295 for (auto i = 0; i < 4; ++i)
1296 {
1297 auto j = (i + 1)%4;
1298 G4ThreeVector p1(fVertices[j].x(), fVerti
1299 G4ThreeVector p2(fVertices[i].x(), fVerti
1300 G4ThreeVector p3(fVertices[j + 4].x(), fV
1301 G4ThreeVector p4(fVertices[i + 4].x(), fV
1302 G4ThreeVector ebot = p2 - p1;
1303 G4ThreeVector etop = p4 - p3;
1304 G4double lbot = ebot.mag();
1305 G4double ltop = etop.mag();
1306 G4double zcross = ebot.x()*etop.y() - ebo
1307 G4double eps = kCarTolerance*std::max(lbo
1308 if (std::min(lbot, ltop) < kCarTolerance
1309 { // plane surface: Dx + Ey + Fz + G = 0
1310 G4ThreeVector normal;
1311 if (std::max(lbot, ltop) < kCarToleranc
1312 {
1313 auto k = (j + 1)%4;
1314 auto l = (k + 1)%4;
1315 G4TwoVector vl = fVertices[l] + fVert
1316 G4TwoVector vi = fVertices[i] + fVert
1317 G4TwoVector vj = fVertices[j] + fVert
1318 G4TwoVector vk = fVertices[k] + fVert
1319 G4TwoVector vij = (vi - vl).unit() +
1320 G4ThreeVector epar = (p4 + p3 - p2 -
1321 G4ThreeVector eort = epar.cross(G4Thr
1322 normal = (eort.cross(epar)).unit();
1323 }
1324 else
1325 {
1326 normal = ((p4 - p1).cross(p3 - p2)).u
1327 }
1328 fSurf[i].D = fPlane[2*i].A = fPlane[2*i
1329 fSurf[i].E = fPlane[2*i].B = fPlane[2*i
1330 fSurf[i].F = fPlane[2*i].C = fPlane[2*i
1331 fSurf[i].G = fPlane[2*i].D = fPlane[2*i
1332 }
1333 else
1334 { // hyperbolic paraboloid: Axz + Byz + C
1335 fIsTwisted = true;
1336 G4double angle = std::acos(ebot.dot(eto
1337 if (angle > CLHEP::halfpi)
1338 {
1339 std::ostringstream message;
1340 message << "Twist on " << angle/CLHEP
1341 << " degrees, should not be m
1342 StreamInfo(message);
1343 G4Exception("G4GenericTrap::ComputeLa
1344 FatalException, message);
1345 }
1346 fTwist[i] = std::copysign(angle, zcross
1347 // set equation of twisted surface (hyp
1348 fSurf[i].A = 2.*fDz*(p4.y() - p3.y() -
1349 fSurf[i].B =-2.*fDz*(p4.x() - p3.x() -
1350 fSurf[i].C = ((p4.x() - p2.x())*(p3.y()
1351 fSurf[i].D = 2.*fDz*fDz*(p4.y() - p3.y(
1352 fSurf[i].E =-2.*fDz*fDz*(p4.x() - p3.x(
1353 fSurf[i].F = 2.*fDz*(p4.x()*p3.y() - p3
1354 fSurf[i].G = fDz*fDz*((p4.x() + p2.x())
1355 G4double magnitude = G4ThreeVector(fSu
1356 if (magnitude < kCarTolerance) continue
1357 fSurf[i].A /= magnitude;
1358 fSurf[i].B /= magnitude;
1359 fSurf[i].C /= magnitude;
1360 fSurf[i].D /= magnitude;
1361 fSurf[i].E /= magnitude;
1362 fSurf[i].F /= magnitude;
1363 fSurf[i].G /= magnitude;
1364 // set planes of bounding polyhedron
1365 G4ThreeVector normal1, normal2;
1366 G4ThreeVector c1, c2;
1367 if (fTwist[i] < 0.)
1368 {
1369 normal1 = ((p2 - p1).cross(p4 - p1)).
1370 normal2 = ((p3 - p4).cross(p1 - p4)).
1371 c1 = p1;
1372 c2 = p4;
1373 }
1374 else
1375 {
1376 normal1 = ((p3 - p2).cross(p1 - p2)).
1377 normal2 = ((p2 - p3).cross(p4 - p3)).
1378 c1 = p2;
1379 c2 = p3;
1380 }
1381 fPlane[2*i].A = normal1.x();
1382 fPlane[2*i].B = normal1.y();
1383 fPlane[2*i].C = normal1.z();
1384 fPlane[2*i].D = -normal1.dot(c1);
1385 fPlane[2*i + 1].A = normal2.x();
1386 fPlane[2*i + 1].B = normal2.y();
1387 fPlane[2*i + 1].C = normal2.z();
1388 fPlane[2*i + 1].D = -normal2.dot(c2);
1389 }
1390 fDelta[i] = (fVertices[i + 4] - fVertices
1391 }
1392 }
1393
1394 /////////////////////////////////////////////
1395 //
1396 // Set bounding box
1397 //
1398 void G4GenericTrap::ComputeBoundingBox()
1399 {
1400 G4double minX, maxX, minY, maxY;
1401 minX = maxX = fVertices[0].x();
1402 minY = maxY = fVertices[0].y();
1403 for (auto i = 1; i < 8; ++i)
1404 {
1405 minX = std::min(minX, fVertices[i].x());
1406 maxX = std::max(maxX, fVertices[i].x());
1407 minY = std::min(minY, fVertices[i].y());
1408 maxY = std::max(maxY, fVertices[i].y());
1409 }
1410 fMinBBox = G4ThreeVector(minX, minY,-fDz);
1411 fMaxBBox = G4ThreeVector(maxX, maxY, fDz);
1412
1413 // Check correctness of the bounding box
1414 //
1415 if (minX >= maxX || minY >= maxY || -fDz >=
1416 {
1417 std::ostringstream message;
1418 message << "Bad bounding box (min >= max)
1419 << GetName() << " !"
1420 << "\npMin = " << fMinBBox
1421 << "\npMax = " << fMaxBBox;
1422 G4Exception("G4GenericTrap::ComputeBoundi
1423 JustWarning, message);
1424 DumpInfo();
1425 }
1426 }
1427
1428 /////////////////////////////////////////////
1429 //
1430 // Set max length of a scratch
1431 //
1432 void G4GenericTrap::ComputeScratchLength()
1433 {
1434 G4double scratch = kInfinity;
1435 for (auto i = 0; i < 4; ++i)
1436 {
1437 if (fTwist[i] == 0.) continue; // skip pl
1438 auto k = (i + 1)%4;
1439 G4ThreeVector p1(fVertices[i].x(), fVerti
1440 G4ThreeVector p2(fVertices[k].x(), fVerti
1441 G4ThreeVector p3(fVertices[i + 4].x(), fV
1442 G4ThreeVector p4(fVertices[k + 4].x(), fV
1443 G4ThreeVector p0 = (p1 + p2 + p3 + p4)*0.
1444 G4ThreeVector norm = SurfaceNormal(p0);
1445 G4ThreeVector pp[2]; // points inside and
1446 pp[0] = p0 - norm * halfTolerance;
1447 pp[1] = p0 + norm * halfTolerance;
1448 G4ThreeVector vv[2]; // unit vectors alon
1449 vv[0] = (p4 - p1).unit();
1450 vv[1] = (p3 - p2).unit();
1451 // find intersection points and compute t
1452 for (auto ip = 0; ip < 2; ++ip)
1453 {
1454 G4double px = pp[ip].x();
1455 G4double py = pp[ip].y();
1456 G4double pz = pp[ip].z();
1457 for (auto iv = 0; iv < 2; ++iv)
1458 {
1459 G4double vx = vv[iv].x();
1460 G4double vy = vv[iv].y();
1461 G4double vz = vv[iv].z();
1462 // solve quadratic equation
1463 G4double ABC = fSurf[i].A*vx + fSurf
1464 G4double ABCF = fSurf[i].A*px + fSurf
1465 G4double A = ABC*vz;
1466 G4double B = 0.5*(fSurf[i].D*vx + fSu
1467 G4double C = fSurf[i].D*px + fSurf[i]
1468 G4double D = B*B - A*C;
1469 if (D < 0) continue;
1470 G4double leng = 2.*std::sqrt(D)/std::
1471 scratch = std::min(leng, scratch);
1472 }
1473 }
1474 }
1475 fScratch = std::max(kCarTolerance, scratch)
1476 }
1477
1478 /////////////////////////////////////////////
1479 //
1480 // GetPolyhedron
1481 //
1482 G4Polyhedron* G4GenericTrap::GetPolyhedron ()
1483 {
1484 if ( (fpPolyhedron == nullptr)
1485 || fRebuildPolyhedron
1486 || (fpPolyhedron->GetNumberOfRotationStep
1487 fpPolyhedron->GetNumberOfRotationStep
1488 {
1489 G4AutoLock l(&polyhedronMutex);
1490 delete fpPolyhedron;
1491 fpPolyhedron = CreatePolyhedron();
1492 fRebuildPolyhedron = false;
1493 l.unlock();
1494 }
1495 return fpPolyhedron;
1496 }
1497
1498 /////////////////////////////////////////////
1499 //
1500 // Method for visualisation
1501 //
1502 void G4GenericTrap::DescribeYourselfTo(G4VGra
1503 {
1504 scene.AddSolid(*this);
1505 }
1506
1507 /////////////////////////////////////////////
1508 //
1509 // Return VisExtent
1510 //
1511 G4VisExtent G4GenericTrap::GetExtent() const
1512 {
1513 return { fMinBBox.x(), fMaxBBox.x(),
1514 fMinBBox.y(), fMaxBBox.y(),
1515 fMinBBox.z(), fMaxBBox.z() };
1516 }
1517
1518 // ------------------------------------------
1519
1520 G4Polyhedron* G4GenericTrap::CreatePolyhedron
1521 {
1522 // Approximation of Twisted Side
1523 // Construct extra Points, if Twisted Side
1524 //
1525 G4Polyhedron* polyhedron;
1526 G4int nVertices, nFacets;
1527
1528 G4int subdivisions = 0;
1529 if (fIsTwisted)
1530 {
1531 if (GetVisSubdivisions() != 0)
1532 {
1533 subdivisions = GetVisSubdivisions();
1534 }
1535 else
1536 {
1537 // Estimation of Number of Subdivisions
1538 //
1539 G4double maxTwist = 0.;
1540 for(G4int i = 0; i < 4; ++i)
1541 {
1542 if (GetTwistAngle(i) > maxTwist) { ma
1543 }
1544
1545 // Computes bounding vectors for the sh
1546 //
1547 G4double Dx, Dy;
1548 Dx = 0.5*(fMaxBBox.x() - fMinBBox.y());
1549 Dy = 0.5*(fMaxBBox.y() - fMinBBox.y());
1550 if (Dy > Dx) { Dx = Dy; }
1551
1552 subdivisions = 8*G4int(maxTwist/(Dx*Dx*
1553 if (subdivisions < 4) { subdivisions =
1554 if (subdivisions > 30) { subdivisions =
1555 }
1556 }
1557 G4int sub4 = 4*subdivisions;
1558 nVertices = 8 + subdivisions*4;
1559 nFacets = 6 + subdivisions*4;
1560 G4double cf = 1./(subdivisions + 1);
1561 polyhedron = new G4Polyhedron(nVertices, nF
1562
1563 // Set vertices
1564 //
1565 G4int icur = 0;
1566 for (G4int i = 0; i < 4; ++i)
1567 {
1568 G4ThreeVector v(fVertices[i].x(),fVertice
1569 polyhedron->SetVertex(++icur, v);
1570 }
1571 for (G4int i = 0; i < subdivisions; ++i)
1572 {
1573 for (G4int j = 0; j < 4; ++j)
1574 {
1575 G4TwoVector u = fVertices[j]+cf*(i+1)*(
1576 G4ThreeVector v(u.x(),u.y(),-fDz+cf*2*f
1577 polyhedron->SetVertex(++icur, v);
1578 }
1579 }
1580 for (G4int i = 4; i < 8; ++i)
1581 {
1582 G4ThreeVector v(fVertices[i].x(),fVertice
1583 polyhedron->SetVertex(++icur, v);
1584 }
1585
1586 // Set facets
1587 //
1588 icur = 0;
1589 polyhedron->SetFacet(++icur, 1, 4, 3, 2); /
1590 for (G4int i = 0; i < subdivisions + 1; ++i
1591 {
1592 G4int is = i*4;
1593 polyhedron->SetFacet(++icur, 5+is, 8+is,
1594 polyhedron->SetFacet(++icur, 8+is, 7+is,
1595 polyhedron->SetFacet(++icur, 7+is, 6+is,
1596 polyhedron->SetFacet(++icur, 6+is, 5+is,
1597 }
1598 polyhedron->SetFacet(++icur, 5+sub4, 6+sub4
1599
1600 polyhedron->SetReferences();
1601 polyhedron->InvertFacets();
1602
1603 return polyhedron;
1604 }
1605
1606 /////////////////////////////////////////////
1607 //
1608 // Print out a warning if A has an unexpected
1609 //
1610 void G4GenericTrap::WarningSignA(const G4Stri
1611 const G4Thre
1612 {
1613 std::ostringstream message;
1614 message.precision(16);
1615 message << icase << " in " << GetName() <<
1616 << " p" << p << "\n"
1617 << " v" << v << "\n"
1618 << " A = " << A << " (is "
1619 << ((A < 0) ? "negative, instead of
1620 StreamInfo(message);
1621 const G4String function = "G4GenericTrap::D
1622 G4Exception(function, "GeomSolids1002", Jus
1623 }
1624
1625 /////////////////////////////////////////////
1626 //
1627 // Print out a warning if B has an unexpected
1628 //
1629 void G4GenericTrap::WarningSignB(const G4Stri
1630 G4double f,
1631 const G4Thre
1632 {
1633 std::ostringstream message;
1634 message.precision(16);
1635 message << icase << " in " << GetName() <<
1636 << " p" << p << "\n"
1637 << " v" << v << "\n"
1638 << " f = " << f << " B = " << B <
1639 << ((B < 0) ? "negative, instead of
1640 StreamInfo(message);
1641 const G4String function = "G4GenericTrap::D
1642 G4Exception(function, "GeomSolids1002", Jus
1643 }
1644
1645 /////////////////////////////////////////////
1646 //
1647 // Print out a warning in DistanceToIn(p,v)
1648 //
1649 void G4GenericTrap::WarningDistanceToIn(G4int
1650 G4dou
1651 const
1652 {
1653 G4String check = "";
1654 if (ttin[1] != DBL_MAX)
1655 {
1656 G4double tcheck = 0.5*(ttout[0] + ttin[1]
1657 if (Inside(p + v*tcheck) != kOutside) che
1658 }
1659
1660 auto position = Inside(p);
1661 std::ostringstream message;
1662 message.precision(16);
1663 message << k << "_Unexpected sequence of in
1664 << " position = " << ((position =
1665 << " p" << p << "\n"
1666 << " v" << v << "\n"
1667 << " range : [" << tmin << ",
1668 << " ttin[2] : "
1669 << ((ttin[0] == DBL_MAX) ? kInfinit
1670 << ((ttin[1] == DBL_MAX) ? kInfinit
1671 << " ttout[2] : "
1672 << ((ttout[0] == DBL_MAX) ? kInfini
1673 << ((ttout[1] == DBL_MAX) ? kInfini
1674 StreamInfo(message);
1675 G4Exception("G4GenericTrap::DistanceToIn(p,
1676 JustWarning, message );
1677 }
1678
1679 /////////////////////////////////////////////
1680 //
1681 // Print out a warning in DistanceToOut(p,v)
1682 //
1683 void G4GenericTrap::WarningDistanceToOut(cons
1684 cons
1685 G4do
1686 {
1687 auto position = Inside(p);
1688 std::ostringstream message;
1689 message.precision(16);
1690 message << "Unexpected final tout = " << to
1691 << " position = " << ((position =
1692 << " p" << p << "\n"
1693 << " v" << v << "\n";
1694 StreamInfo(message);
1695 G4Exception("G4GenericTrap::DistanceToOut(p
1696 JustWarning, message );
1697 }
1698
1699 #endif
1700