| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4GenericTrap implementation
27 //
28 // Authors:
29 // Tatiana Nikitina, CERN; Ivana Hrivnacova, IPN Orsay
30 // Adapted from Root Arb8 implementation by Andrei Gheata, CERN
31 //
32 // 27.05.2024 - Evgueni Tcherniaev, complete revision, speed up
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_INITIALIZER;
59 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZER;
60 }
61
62 ////////////////////////////////////////////////////////////////////////
63 //
64 // Constructor
65 //
66 G4GenericTrap::G4GenericTrap(const G4String& name, G4double halfZ,
67 const std::vector<G4TwoVector>& vertices)
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 data and allocates memory
80 // for usage restricted to object persistency.
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 G4GenericTrap& rhs)
99 : G4VSolid(rhs),
100 halfTolerance(rhs.halfTolerance), fScratch(rhs.fScratch),
101 fDz(rhs.fDz), fVertices(rhs.fVertices), fIsTwisted(rhs.fIsTwisted),
102 fMinBBox(rhs.fMinBBox), fMaxBBox(rhs.fMaxBBox),
103 fVisSubdivisions(rhs.fVisSubdivisions),
104 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
105 {
106 for (auto i = 0; i < 5; ++i) { fTwist[i] = rhs.fTwist[i]; }
107 for (auto i = 0; i < 4; ++i) { fDelta[i] = rhs.fDelta[i]; }
108 for (auto i = 0; i < 8; ++i) { fPlane[i] = rhs.fPlane[i]; }
109 for (auto i = 0; i < 4; ++i) { fSurf[i] = rhs.fSurf[i]; }
110 for (auto i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
111 }
112
113 ////////////////////////////////////////////////////////////////////////
114 //
115 // Assignment
116 //
117 G4GenericTrap& G4GenericTrap::operator=(const G4GenericTrap& rhs)
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 = rhs.fScratch;
127 fDz = rhs.fDz; fVertices = rhs.fVertices; fIsTwisted = rhs.fIsTwisted;
128 fMinBBox = rhs.fMinBBox; fMaxBBox = rhs.fMaxBBox;
129 fVisSubdivisions = rhs.fVisSubdivisions;
130 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
131
132 for (auto i = 0; i < 8; ++i) { fVertices[i] = rhs.fVertices[i]; }
133 for (auto i = 0; i < 5; ++i) { fTwist[i] = rhs.fTwist[i]; }
134 for (auto i = 0; i < 4; ++i) { fDelta[i] = rhs.fDelta[i]; }
135 for (auto i = 0; i < 8; ++i) { fPlane[i] = rhs.fPlane[i]; }
136 for (auto i = 0; i < 4; ++i) { fSurf[i] = rhs.fSurf[i]; }
137 for (auto i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
138
139 fRebuildPolyhedron = false;
140 delete fpPolyhedron; fpPolyhedron = nullptr;
141
142 return *this;
143 }
144
145 ////////////////////////////////////////////////////////////////////////
146 //
147 // Returns position of the point (inside/outside/surface)
148 //
149 EInside G4GenericTrap::Inside(const G4ThreeVector& p) const
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] = fVertices[i] + fDelta[i]*z; }
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*py + fSurf[i].F*pz + fSurf[i].G;
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 - a.x()))/leng;
176 dist = std::max(dd, dist);
177 }
178 }
179 return (dist > halfTolerance) ? kOutside :
180 ((dist > -halfTolerance) ? kSurface : kInside);
181 }
182
183 ////////////////////////////////////////////////////////////////////////
184 //
185 // Return unit normal to surface at p
186 //
187 G4ThreeVector G4GenericTrap::SurfaceNormal( const G4ThreeVector& p ) const
188 {
189 G4double halfToleranceSquared = halfTolerance*halfTolerance;
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] = fVertices[i] + fDelta[i]*tz; }
197
198 // check bottom and top faces
199 G4double dz = std::abs(pz) - fDz;
200 nz = std::copysign(G4double(std::abs(dz) <= halfTolerance), pz);
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*py + fSurf[i].F*pz + fSurf[i].G;
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 - a.x());
224 if (dd*dd <= halfToleranceSquared*ll)
225 {
226 G4double x = fSurf[i].A*pz + fSurf[i].D;
227 G4double y = fSurf[i].B*pz + fSurf[i].E;
228 G4double z = fSurf[i].A*px + fSurf[i].B*py + 2.*fSurf[i].C*pz + fSurf[i].F;
229 G4double inv = 1./std::sqrt(x*x + y*y + z*z);
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); // point is not on the surface
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 return corresponding normal
248 // Normally this method should not be called
249 //
250 G4ThreeVector G4GenericTrap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
251 {
252 #ifdef G4SPECSDEBUG
253 std::ostringstream message;
254 message.precision(16);
255 message << "Point p is not on surface of solid: " << GetName() << " !?\n"
256 << "Position: " << p << " is "
257 << ((Inside(p) == kInside) ? "inside" : "outside") << "\n";
258 StreamInfo(message);
259 G4Exception("G4GenericTrap::ApproxSurfaceNormal(p)", "GeomSolids1002",
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] = fVertices[i] + fDelta[i]*tz; }
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*py + fSurf[i].F*pz + fSurf[i].G;
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 - a.x()))/leng;
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::copysign(1., pz) };
294
295 G4double x = fSurf[iside].A*pz + fSurf[iside].D;
296 G4double y = fSurf[iside].B*pz + fSurf[iside].E;
297 G4double z = fSurf[iside].A*px + fSurf[iside].B*py + 2.*fSurf[iside].C*pz + fSurf[iside].F;
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 outside,
305 // return kInfinity if no intersection
306 //
307 G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p,
308 const G4ThreeVector& v) const
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 polyhedron
314 //
315 if (std::abs(pz) - fDz >= -halfTolerance && pz*vz >= 0) { return kInfinity; }
316 G4double invz = (vz == 0) ? kInfinity : -1./vz;
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 twisted faces
325 const G4GenericTrapPlane& surf = fPlane[2*k];
326 G4double cosa = surf.A*vx + surf.B*vy + surf.C*vz;
327 G4double dist = surf.A*px + surf.B*py + surf.C*pz + surf.D;
328 if (dist >= -halfTolerance)
329 {
330 if (cosa >= 0.) { return kInfinity; } // point flies away
331 G4double tmp = -dist/cosa;
332 xin = std::max(tmp, xin);
333 }
334 else
335 {
336 if (cosa > 0) { xout = std::min(-dist/cosa, xout); }
337 if (cosa < 0) { xin = std::max(-dist/cosa, xin); }
338 }
339 }
340
341 // Check planes around twisted faces, each twisted face is bounded by two planes
342 G4double tin = xin;
343 G4double tout = xout;
344 for (auto i = 0; i < 4; ++i)
345 {
346 if (fTwist[i] == 0) continue; // skip plane faces
347
348 // check intersection with 1st bounding plane
349 const G4GenericTrapPlane& surf1 = fPlane[2*i];
350 G4double cosa = surf1.A*vx + surf1.B*vy + surf1.C*vz;
351 G4double dist = surf1.A*px + surf1.B*py + surf1.C*pz + surf1.D;
352 if (dist >= -halfTolerance)
353 {
354 if (cosa >= 0.) { return kInfinity; } // point flies away
355 G4double tmp = -dist/cosa;
356 tin = std::max(tmp, tin);
357 }
358 else
359 {
360 if (cosa > 0) { tout = std::min(-dist/cosa, tout); }
361 if (cosa < 0) { tin = std::max(-dist/cosa, tin); }
362 }
363
364 // check intersection with 2nd bounding plane
365 const G4GenericTrapPlane& surf2 = fPlane[2*i + 1];
366 cosa = surf2.A*vx + surf2.B*vy + surf2.C*vz;
367 dist = surf2.A*px + surf2.B*py + surf2.C*pz + surf2.D;
368 if (dist >= -halfTolerance)
369 {
370 if (cosa >= 0.) { return kInfinity; } // point flies away
371 G4double tmp = -dist/cosa;
372 tin = std::max(tmp, tin);
373 }
374 else
375 {
376 if (cosa > 0) { tout = std::min(-dist/cosa, tout); }
377 if (cosa < 0) { tin = std::max(-dist/cosa, tin); }
378 }
379 }
380 if (tout - tin <= halfTolerance) { return kInfinity; } // touch or no hit
381
382 // if point is outside of the bounding box
383 // then move it to the surface of the bounding polyhedron
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] = xout;
405 G4double tminimal = tin - halfTolerance;
406 G4double tmaximal = tout + halfTolerance;
407
408 constexpr G4double EPSILON = 1000.*DBL_EPSILON;
409 for (auto i = 0; i < 4; ++i)
410 {
411 if (fTwist[i] == 0) continue; // skip plane faces
412
413 // twisted face, solve quadratic equation
414 G4double ABC = fSurf[i].A*vx + fSurf[i].B*vy + fSurf[i].C*vz;
415 G4double ABCF = fSurf[i].A*x0 + fSurf[i].B*y0 + fSurf[i].C*z0 + fSurf[i].F;
416 G4double A = ABC*vz;
417 G4double B = 0.5*(fSurf[i].D*vx + fSurf[i].E*vy + ABCF*vz + ABC*z0);
418 G4double C = fSurf[i].D*x0 + fSurf[i].E*y0 + ABCF*z0 + fSurf[i].G;
419 if (std::abs(A) <= EPSILON)
420 {
421 // case 1: track is parallel to the surface
422 if (std::abs(B) <= EPSILON)
423 {
424 // check position of the track relative to the surface,
425 // for the reason of precision it's better to use (x0,y0,z0) instead of (px,py,pz)
426 auto j = (i + 1)%4;
427 G4double z = (z0 + fDz);
428 G4TwoVector a = fVertices[i] + fDelta[i]*z;
429 G4TwoVector b = fVertices[j] + fDelta[j]*z;
430 G4double dx = b.x() - a.x();
431 G4double dy = b.y() - a.y();
432 G4double leng = std::sqrt(dx*dx + dy*dy);
433 G4double dist = dx*(y0 - a.y()) - dy*(x0 - a.x());
434 if (dist >= -halfTolerance*leng) { return kInfinity; }
435 continue;
436 }
437
438 // case 2: single root
439 G4double tmp = t0 - 0.5*C/B;
440 // compute normal at the intersection point and check direction
441 G4double x = px + vx*tmp;
442 G4double y = py + vy*tmp;
443 G4double z = pz + vz*tmp;
444 const G4GenericTrapSurface& surf = fSurf[i];
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.*surf.C*z + surf.F;
448
449 if (nx*vx + ny*vy + nz*vz >= 0.) // point is flying to outside
450 {
451 auto k = (i == 3) ? 0 : i + 1;
452 G4double tz = (pz + fDz);
453 G4TwoVector a = fVertices[i] + fDelta[i]*tz;
454 G4TwoVector b = fVertices[k] + fDelta[k]*tz;
455 G4double dx = b.x() - a.x();
456 G4double dy = b.y() - a.y();
457 G4double leng = std::sqrt(dx*dx + dy*dy);
458 G4double dist = dx*(py - a.y()) - dy*(px - a.x());
459 if (dist >= -halfTolerance*leng) { return kInfinity; } // point is flies away
460
461 if (tmp < tminimal || tmp > tmaximal) continue;
462 if (std::abs(tmp - ttout[0]) < halfTolerance) continue;
463 if (tmp < ttout[0])
464 {
465 ttout[1] = ttout[0];
466 ttout[0] = tmp;
467 }
468 else { ttout[1] = std::min(tmp, ttout[1]); }
469 }
470 else // point is flying to inside
471 {
472 if (tmp < tminimal || tmp > tmaximal) continue;
473 if (std::abs(tmp - ttin[0]) < halfTolerance) continue;
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::sqrt(D), B);
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 <= tmaximal)
501 {
502 if (std::abs(tsurfin - ttin[0]) >= halfTolerance)
503 {
504 if (tsurfin < ttin[0])
505 {
506 ttin[1] = ttin[0];
507 ttin[0] = tsurfin;
508 }
509 else { ttin[1] = std::min(tsurfin, ttin[1]); }
510 }
511 }
512 if (tsurfout >= tminimal && tsurfout <= tmaximal)
513 {
514 if (std::abs(tsurfout - ttout[0]) >= halfTolerance)
515 {
516 if (tsurfout < ttout[0])
517 {
518 ttout[1] = ttout[0];
519 ttout[0] = tsurfout;
520 }
521 else { ttout[1] = std::min(tsurfout, ttout[1]); }
522 }
523 }
524 }
525
526 // Compute distance to In
527 //
528 if (ttin[0] == DBL_MAX) { return kInfinity; } // no entry point
529
530 // single entry point
531 if (ttin[1] == DBL_MAX)
532 {
533 G4double distin = ttin[0];
534 G4double distout = (distin >= ttout[0] - halfTolerance) ? ttout[1] : ttout[0];
535 G4double dist = (distout <= halfTolerance || distout - distin <= halfTolerance) ? kInfinity : distin;
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[0];
543 G4double dist = (distout <= halfTolerance || distout - distin <= halfTolerance) ? kInfinity : distin;
544 return (dist < halfTolerance) ? 0. : dist;
545 }
546
547 // check 1st pair of in-out points
548 G4double distin1 = ttin[0], distout1 = ttout[0];
549 G4double dist1 = (distout1 <= halfTolerance || distout1 - distin1 <= halfTolerance) ? kInfinity : distin1;
550 if (dist1 != kInfinity) { return (dist1 < halfTolerance) ? 0. : dist1; }
551
552 // check 2nd pair of in-out points
553 G4double distin2 = ttin[1], distout2 = ttout[1];
554 G4double dist2 = (distout2 <= halfTolerance || distout2 - distin2 <= halfTolerance) ? kInfinity : distin2;
555 return (dist2 < halfTolerance) ? 0. : dist2;
556 }
557
558 ////////////////////////////////////////////////////////////////////////
559 //
560 // Estimate safety distance to the surface from outside
561 //
562 G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p) const
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] = fVertices[i] + fDelta[i]*z; }
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*py + fSurf[i].F*pz + fSurf[i].G;
578 dist = std::max(dd, dist);
579 }
580 else
581 {
582 // comptute distance to the wedge (two planes) in front of the surface
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 - pp[i].y())*cx;
587 G4ThreeVector na(-cy, cx, fDelta[i].x()*cy - fDelta[i].y()*cx);
588 G4ThreeVector nb(-cy, cx, fDelta[j].x()*cy - fDelta[j].y()*cx);
589 G4double amag2 = na.mag2();
590 G4double bmag2 = nb.mag2();
591 G4double distab = (amag2 > bmag2) ? d/std::sqrt(amag2) : d/std::sqrt(bmag2); // d > 0
592 dist = std::max(distab, dist);
593 }
594 }
595 return (dist > 0.) ? dist : 0.; // return safety distance
596 }
597
598 ////////////////////////////////////////////////////////////////////////
599 //
600 // Calculate distance to surface from inside
601 //
602 G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p,
603 const G4ThreeVector& v,
604 const G4bool calcNorm,
605 G4bool* validNorm,
606 G4ThreeVector* n) const
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 && pz*vz > 0.)
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::copysign(fDz, vz) - pz)/vz;
623 G4int iface = (vz < 0) ? -4 : -2; // little trick for z-normal: (-4+3)=-1, (-2+3)=+1
624
625 for (auto i = 0; i < 4; ++i)
626 {
627 if (fTwist[i] != 0) continue; // skip twisted faces
628 const G4GenericTrapPlane& surf = fPlane[2*i];
629 G4double cosa = surf.A*vx + surf.B*vy + surf.C*vz;
630 if (cosa > 0)
631 {
632 G4double dist = surf.A*px + surf.B*py + surf.C*pz + surf.D;
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_EPSILON;
650 for (auto i = 0; i < 4; ++i)
651 {
652 if (fTwist[i] == 0) continue; // skip plane faces
653
654 // twisted face, solve quadratic equation
655 const G4GenericTrapSurface& surf = fSurf[i];
656 G4double ABC = surf.A*vx + surf.B*vy + surf.C*vz;
657 G4double ABCF = surf.A*px + surf.B*py + surf.C*pz + surf.F;
658 G4double A = ABC*vz;
659 G4double B = 0.5*(surf.D*vx + surf.E*vy + ABCF*vz + ABC*pz);
660 G4double C = surf.D*px + surf.E*py + ABCF*pz + surf.G;
661
662 if (std::abs(A) <= EPSILON)
663 {
664 // case 1: track is parallel to the surface
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 and check direction
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.*surf.C*z + surf.F;
676
677 if (nx*vx + ny*vy + nz*vz > 0.) // point is flying to outside
678 {
679 auto k = (i + 1)%4;
680 G4double tz = (pz + fDz);
681 G4TwoVector a = fVertices[i] + fDelta[i]*tz;
682 G4TwoVector b = fVertices[k] + fDelta[k]*tz;
683 G4double dx = b.x() - a.x();
684 G4double dy = b.y() - a.y();
685 G4double leng = std::sqrt(dx*dx + dy*dy);
686 G4double dist = dx*(py - a.y()) - dy*(px - a.x());
687 if (dist >= -halfTolerance*leng) // point is on the surface
688 {
689 if (calcNorm)
690 {
691 *validNorm = false;
692 G4double normx = surf.A*pz + surf.D;
693 G4double normy = surf.B*pz + surf.E;
694 G4double normz = surf.A*px + surf.B*py + 2.*surf.C*pz + surf.F;
695 G4double inv = 1./std::sqrt(normx*normx + normy*normy + normz*normz);
696 n->set(normx*inv, normy*inv, normz*inv);
697 }
698 return 0.;
699 }
700 if (tout > tmp) { tout = tmp; iface = i; }
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]*tz;
713 G4TwoVector b = fVertices[j] + fDelta[j]*tz;
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 - a.x());
718 if (dist <= -halfTolerance*leng) { continue; } // point is inside
719 if (A > 0 || dist > halfTolerance*leng) // convex surface (or point is outside)
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 + 2.*surf.C*pz + surf.F;
727 G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
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::sqrt(D), B);
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) -> tmax(in)
743 {
744 if (std::abs(tmax) < std::abs(tmin)) continue; // point flies inside
745 if (tmin <= halfTolerance) // point is on external side of the surface
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[i]*tz;
755 G4TwoVector b = fVertices[j] + fDelta[j]*tz;
756 G4double dx = b.x() - a.x();
757 G4double dy = b.y() - a.y();
758 G4double leng = std::sqrt(dx*dx + dy*dy);
759 G4double dist = dx*(y - a.y()) - dy*(x - a.x());
760 if (dist <= -halfTolerance*leng) continue; // scratch
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 + 2.*surf.C*pz + surf.F;
767 G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
768 n->set(nx*inv, ny*inv, nz*inv);
769 }
770 return 0.;
771 }
772 if (tout > tmin) { tout = tmin; iface = i; }
773 }
774 else // convex profile: tmin(in) -> tmax(out)
775 {
776 if (tmax < halfTolerance) // point is on the surface
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 + 2.*surf.C*pz + surf.F;
784 G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
785 n->set(nx*inv, ny*inv, nz*inv);
786 }
787 return 0.;
788 }
789 if (tout > tmax) { tout = tmax; iface = i; }
790 }
791 }
792
793 // Compute normal, if required, and return distance to out
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: (-4+3)=-1, (-2+3)=+1
802 }
803 else
804 {
805 const G4GenericTrapSurface& surf = fSurf[iface];
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.*surf.C*z + surf.F;
820 G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
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 from inside
831 //
832 G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p) const
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] = fVertices[i] + fDelta[i]*z; }
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*py + fSurf[i].F*pz + fSurf[i].G;
848 dist = std::max(dd, dist);
849 }
850 else
851 {
852 // comptute distance to the wedge (two planes) in front of the surface
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 - pp[i].y())*cx;
857 G4ThreeVector na(-cy, cx, fDelta[i].x()*cy - fDelta[i].y()*cx);
858 G4ThreeVector nb(-cy, cx, fDelta[j].x()*cy - fDelta[j].y()*cx);
859 G4double amag2 = na.mag2();
860 G4double bmag2 = nb.mag2();
861 G4double distab = (amag2 > bmag2) ? d/std::sqrt(amag2) : d/std::sqrt(bmag2); // d < 0
862 dist = std::max(distab, dist);
863 }
864 }
865 return (dist < 0.) ? -dist : 0.; // return safety distance
866 }
867
868 ////////////////////////////////////////////////////////////////////////
869 //
870 // Return bounding box
871 //
872 void G4GenericTrap::BoundingLimits(G4ThreeVector& pMin,
873 G4ThreeVector& pMax) const
874 {
875 pMin = fMinBBox;
876 pMax = fMaxBBox;
877 }
878
879 ////////////////////////////////////////////////////////////////////////
880 //
881 // Calculate extent under transform and specified limits
882 //
883 G4bool
884 G4GenericTrap::CalculateExtent(const EAxis pAxis,
885 const G4VoxelLimits& pVoxelLimit,
886 const G4AffineTransform& pTransform,
887 G4double& pMin, G4double& pMax) const
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,pVoxelLimit,pTransform,pMin,pMax);
898 #endif
899 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
900 {
901 return exist = pMin < pMax;
902 }
903
904 // Set bounding envelope (benv) and calculate extent
905 //
906 // To build the bounding envelope with plane faces, each lateral face of
907 // the trapezoid is subdivided in two triangles. Subdivision is done by
908 // duplication of vertices in the bases in a way that the envelope be
909 // a convex polyhedron (some faces of the envelope can be degenerated)
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] : baseA[k1];
929 baseB[k1+1] = (znorm < 0.0) ? baseB[k1] : baseB[k2];
930 }
931
932 std::vector<const G4ThreeVectorList *> polygons(2);
933 polygons[0] = &baseA;
934 polygons[1] = &baseB;
935
936 G4BoundingEnvelope benv(bmin, bmax, polygons);
937 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
938 return exist;
939 }
940
941 ////////////////////////////////////////////////////////////////////////
942 //
943 // Return type of the solid
944 //
945 G4GeometryType G4GenericTrap::GetEntityType() const
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 specified output stream
971 //
972 std::ostream& G4GenericTrap::StreamInfo(std::ostream& os) const
973 {
974 G4long oldprc = os.precision(16);
975 os << "-----------------------------------------------------------\n"
976 << " *** Dump for solid - " << GetName() << " ***\n"
977 << " ===================================================\n"
978 << "Solid geometry type: " << GetEntityType() << "\n"
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] << "\n";
984 }
985 os << "-----------------------------------------------------------\n";
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() const
995 {
996 if (fArea[0] + fArea[1] + fArea[2] + fArea[3] == 0.)
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,3}, {3,7,4,0}, {4,5,6,7} };
1008
1009 G4bool isTwisted[6] = {false};
1010 for (auto i = 0; i < 4; ++i) { isTwisted[i + 1] = (fTwist[i] != 0.0); }
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; // -fDz face
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()*D.x())*.5; // +fDz face
1023
1024 G4double select = ssurf[5]*G4QuickRand();
1025 G4int k;
1026 for (k = 0; k < 5; ++k) { if (select <= ssurf[k]) break; }
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].y(), ((k == 5) ? fDz : -fDz));
1034 pp[1].set(fVertices[i1].x(), fVertices[i1].y(), ((k == 0) ? -fDz : fDz));
1035 pp[2].set(fVertices[i2].x(), fVertices[i2].y(), ((k == 0) ? -fDz : fDz));
1036 pp[3].set(fVertices[i3].x(), fVertices[i3].y(), ((k == 5) ? fDz : -fDz));
1037
1038 G4ThreeVector point;
1039 if (isTwisted[k])
1040 { // twisted face, rejection sampling
1041 G4double maxlength = std::max((pp[2] - pp[1]).mag(), (pp[3] - pp[0]).mag());
1042 point = (pp[0] + pp[1] + pp[2] + pp[3])*0.25;
1043 for (auto i = 0; i < 10000; ++i)
1044 {
1045 G4double u = G4QuickRand();
1046 G4ThreeVector p0 = pp[0] + (pp[1] - pp[0])*u;
1047 G4ThreeVector p1 = pp[3] + (pp[2] - pp[3])*u;
1048 G4double v = G4QuickRand()*(maxlength/(p1 - p0).mag());
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[3] - pp[0])).mag())*0.5;
1060 G4int j = (select > ssurf[k] - ss) ? 3 : 1;
1061 point = pp[0] + (pp[2] - pp[0])*u + (pp[j] - pp[0])*v;
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[1];
1076 G4TwoVector B = fVertices[2] - fVertices[0];
1077 G4TwoVector C = fVertices[7] - fVertices[5];
1078 G4TwoVector D = fVertices[6] - fVertices[4];
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)/6.)*fDz;
1087 }
1088 return fCubicVolume;
1089 }
1090
1091 ////////////////////////////////////////////////////////////////////////
1092 //
1093 // Compute lateral face area
1094 //
1095 G4double G4GenericTrap::GetLateralFaceArea(G4int iface) const
1096 {
1097 G4int i1 = iface, i2 = (i1 + 1)%4, i3 = i1 + 4, i4 = i2 + 4;
1098
1099 // plane face
1100 if (fTwist[iface] == 0.0)
1101 {
1102 G4ThreeVector p1(fVertices[i1].x(), fVertices[i1].y(),-fDz);
1103 G4ThreeVector p2(fVertices[i2].x(), fVertices[i2].y(),-fDz);
1104 G4ThreeVector p3(fVertices[i3].x(), fVertices[i3].y(), fDz);
1105 G4ThreeVector p4(fVertices[i4].x(), fVertices[i4].y(), fDz);
1106 return ((p4 - p1).cross(p3 - p2)).mag()*0.5;
1107 }
1108
1109 // twisted face
1110 constexpr G4int NSTEP = 250;
1111 constexpr G4double dt = 1./NSTEP;
1112
1113 G4double x21 = fVertices[i2].x() - fVertices[i1].x();
1114 G4double y21 = fVertices[i2].y() - fVertices[i1].y();
1115 G4double x31 = fVertices[i3].x() - fVertices[i1].x();
1116 G4double y31 = fVertices[i3].y() - fVertices[i1].y();
1117 G4double x42 = fVertices[i4].x() - fVertices[i2].x();
1118 G4double y42 = fVertices[i4].y() - fVertices[i2].y();
1119 G4double x43 = fVertices[i4].x() - fVertices[i3].x();
1120 G4double y43 = fVertices[i4].y() - fVertices[i3].y();
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*R1 + 2.*aa + bb));
1146 G4double log0 = std::log(std::abs(sqrtAA*R0 + bb));
1147
1148 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0) + IIJJ*invSqrtAA*(log1 - log0);
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[1];
1162 G4TwoVector B = fVertices[2] - fVertices[0];
1163 G4TwoVector C = fVertices[7] - fVertices[5];
1164 G4TwoVector D = fVertices[6] - fVertices[4];
1165 G4double S_bot = (A.x()*B.y() - A.y()*B.x())*0.5;
1166 G4double S_top = (C.x()*D.y() - C.y()*D.x())*0.5;
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] + fArea[1] + fArea[2] + fArea[3];
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::vector<G4TwoVector>& vertices)
1183 {
1184 // Check number of vertices
1185 //
1186 if (vertices.size() != 8)
1187 {
1188 std::ostringstream message;
1189 message << "Number of vertices is " << vertices.size()
1190 << " instead of 8 for Solid: " << GetName() << " !";
1191 G4Exception("G4GenericTrap::CheckParameters()", "GeomSolids0002",
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 negative (dZ = " << fDz << " mm)"
1201 << " for Solid: " << GetName() << " !";
1202 G4Exception("G4GenericTrap::CheckParameters()", "GeomSolids0002",
1203 FatalException, message);
1204 }
1205
1206 // Check order of vertices
1207 //
1208 for (auto i = 0; i < 8; ++i) { fVertices.at(i) = vertices[i]; }
1209
1210 // Bottom polygon area
1211 G4TwoVector diag1 = fVertices[2] - fVertices[0];
1212 G4TwoVector diag2 = fVertices[3] - fVertices[1];
1213 G4double ldiagbot = std::max(diag1.mag(), diag2.mag());
1214 G4double zbot = diag1.x()*diag2.y() - diag1.y()*diag2.x();
1215 if (std::abs(zbot) < ldiagbot*kCarTolerance) zbot = 0.;
1216
1217 // Top polygon area
1218 G4TwoVector diag3 = fVertices[6] - fVertices[4];
1219 G4TwoVector diag4 = fVertices[7] - fVertices[5];
1220 G4double ldiagtop = std::max(diag3.mag(), diag4.mag());
1221 G4double ztop = diag3.x()*diag4.y() - diag3.y()*diag4.x();
1222 if (std::abs(ztop) < ldiagtop*kCarTolerance) ztop = 0.;
1223
1224 if (zbot*ztop < 0.)
1225 {
1226 std::ostringstream message;
1227 message << "Vertices of the bottom and top polygons are defined in opposite directions !\n";
1228 StreamInfo(message);
1229 G4Exception("G4GenericTrap::CheckParameters()", "GeomSolids0002",
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: " << GetName() << " !\n"
1238 << "Vertices of the bottom and top polygons must be defined in a clockwise direction.";
1239 G4Exception("G4GenericTrap::CheckParameters()", "GeomSolids1001",
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[i]).mag();
1250 G4double ltop = (fVertices[j + 4] - fVertices[i + 4]).mag();
1251 ndegenerated += (std::max(lbot, ltop) < kCarTolerance);
1252 }
1253 if (ndegenerated > 1 ||
1254 GetCubicVolume() < fDz*std::max(ldiagbot, ldiagtop)*kCarTolerance)
1255 {
1256 std::ostringstream message;
1257 message << "Degenerated solid !\n";
1258 StreamInfo(message);
1259 G4Exception("G4GenericTrap::CheckParameters()", "GeomSolids0002",
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] - fVertices[i];
1271 G4TwoVector edge2 = fVertices[k] - fVertices[j];
1272 isConvex = ((edge1.x()*edge2.y() - edge1.y()*edge2.x()) < kCarTolerance);
1273 if (!isConvex) break;
1274 G4TwoVector edge3 = fVertices[j + 4] - fVertices[i + 4];
1275 G4TwoVector edge4 = fVertices[k + 4] - fVertices[j + 4];
1276 isConvex = ((edge3.x()*edge4.y() - edge3.y()*edge4.x()) < kCarTolerance);
1277 if (!isConvex) break;
1278 }
1279 if (!isConvex)
1280 {
1281 std::ostringstream message;
1282 message << "The bottom and top faces must be convex polygons !\n";
1283 StreamInfo(message);
1284 G4Exception("G4GenericTrap::CheckParameters()", "GeomSolids0002",
1285 FatalException, message);
1286 }
1287 }
1288
1289 ////////////////////////////////////////////////////////////////////////
1290 //
1291 // Compute surface equations and twist angles of lateral faces
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(), fVertices[j].y(), -fDz);
1299 G4ThreeVector p2(fVertices[i].x(), fVertices[i].y(), -fDz);
1300 G4ThreeVector p3(fVertices[j + 4].x(), fVertices[j + 4].y(), fDz);
1301 G4ThreeVector p4(fVertices[i + 4].x(), fVertices[i + 4].y(), fDz);
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() - ebot.y()*etop.x();
1307 G4double eps = kCarTolerance*std::max(lbot,ltop);
1308 if (std::min(lbot, ltop) < kCarTolerance || std::abs(zcross) < eps)
1309 { // plane surface: Dx + Ey + Fz + G = 0
1310 G4ThreeVector normal;
1311 if (std::max(lbot, ltop) < kCarTolerance) // degenerated face
1312 {
1313 auto k = (j + 1)%4; // N
1314 auto l = (k + 1)%4; // i | j
1315 G4TwoVector vl = fVertices[l] + fVertices[l + 4]; // +---+
1316 G4TwoVector vi = fVertices[i] + fVertices[i + 4]; // l / \ k
1317 G4TwoVector vj = fVertices[j] + fVertices[j + 4]; // +-------+
1318 G4TwoVector vk = fVertices[k] + fVertices[k + 4]; //
1319 G4TwoVector vij = (vi - vl).unit() + (vj - vk).unit();
1320 G4ThreeVector epar = (p4 + p3 - p2 - p1);
1321 G4ThreeVector eort = epar.cross(G4ThreeVector(vij.x(), vij.y(), 0.0));
1322 normal = (eort.cross(epar)).unit();
1323 }
1324 else
1325 {
1326 normal = ((p4 - p1).cross(p3 - p2)).unit();
1327 }
1328 fSurf[i].D = fPlane[2*i].A = fPlane[2*i + 1].A = normal.x();
1329 fSurf[i].E = fPlane[2*i].B = fPlane[2*i + 1].B = normal.y();
1330 fSurf[i].F = fPlane[2*i].C = fPlane[2*i + 1].C = normal.z();
1331 fSurf[i].G = fPlane[2*i].D = fPlane[2*i + 1].D = -normal.dot((p1 + p2 + p3 + p4)/4.);
1332 }
1333 else
1334 { // hyperbolic paraboloid: Axz + Byz + Czz + Dx + Ey + Fz + G = 0
1335 fIsTwisted = true;
1336 G4double angle = std::acos(ebot.dot(etop)/(lbot*ltop));
1337 if (angle > CLHEP::halfpi)
1338 {
1339 std::ostringstream message;
1340 message << "Twist on " << angle/CLHEP::deg
1341 << " degrees, should not be more than 90 degrees !";
1342 StreamInfo(message);
1343 G4Exception("G4GenericTrap::ComputeLateralSurfaces()", "GeomSolids0002",
1344 FatalException, message);
1345 }
1346 fTwist[i] = std::copysign(angle, zcross);
1347 // set equation of twisted surface (hyperbolic paraboloid)
1348 fSurf[i].A = 2.*fDz*(p4.y() - p3.y() - p2.y() + p1.y());
1349 fSurf[i].B =-2.*fDz*(p4.x() - p3.x() - p2.x() + p1.x());
1350 fSurf[i].C = ((p4.x() - p2.x())*(p3.y() - p1.y()) - (p4.y() - p2.y())*(p3.x() - p1.x()));
1351 fSurf[i].D = 2.*fDz*fDz*(p4.y() - p3.y() + p2.y() - p1.y());
1352 fSurf[i].E =-2.*fDz*fDz*(p4.x() - p3.x() + p2.x() - p1.x());
1353 fSurf[i].F = 2.*fDz*(p4.x()*p3.y() - p3.x()*p4.y() - p2.x()*p1.y() + p1.x()*p2.y());
1354 fSurf[i].G = fDz*fDz*((p4.x() + p2.x())*(p3.y() + p1.y()) - (p3.x() + p1.x())*(p4.y() + p2.y()));
1355 G4double magnitude = G4ThreeVector(fSurf[i].D, fSurf[i].E, fSurf[i].F).mag();
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)).unit();
1370 normal2 = ((p3 - p4).cross(p1 - p4)).unit();
1371 c1 = p1;
1372 c2 = p4;
1373 }
1374 else
1375 {
1376 normal1 = ((p3 - p2).cross(p1 - p2)).unit();
1377 normal2 = ((p2 - p3).cross(p4 - p3)).unit();
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[i])/(2*fDz);
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 >= fDz)
1416 {
1417 std::ostringstream message;
1418 message << "Bad bounding box (min >= max) for solid: "
1419 << GetName() << " !"
1420 << "\npMin = " << fMinBBox
1421 << "\npMax = " << fMaxBBox;
1422 G4Exception("G4GenericTrap::ComputeBoundingBox()", "GeomSolids1001",
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 plane face
1438 auto k = (i + 1)%4;
1439 G4ThreeVector p1(fVertices[i].x(), fVertices[i].y(), -fDz);
1440 G4ThreeVector p2(fVertices[k].x(), fVertices[k].y(), -fDz);
1441 G4ThreeVector p3(fVertices[i + 4].x(), fVertices[i + 4].y(), fDz);
1442 G4ThreeVector p4(fVertices[k + 4].x(), fVertices[k + 4].y(), fDz);
1443 G4ThreeVector p0 = (p1 + p2 + p3 + p4)*0.25; // center of the face
1444 G4ThreeVector norm = SurfaceNormal(p0);
1445 G4ThreeVector pp[2]; // points inside and outside the surface
1446 pp[0] = p0 - norm * halfTolerance;
1447 pp[1] = p0 + norm * halfTolerance;
1448 G4ThreeVector vv[2]; // unit vectors along the diagonals
1449 vv[0] = (p4 - p1).unit();
1450 vv[1] = (p3 - p2).unit();
1451 // find intersection points and compute the scratch
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[i].B*vy + fSurf[i].C*vz;
1464 G4double ABCF = fSurf[i].A*px + fSurf[i].B*py + fSurf[i].C*pz + fSurf[i].F;
1465 G4double A = ABC*vz;
1466 G4double B = 0.5*(fSurf[i].D*vx + fSurf[i].E*vy + ABCF*vz + ABC*pz);
1467 G4double C = fSurf[i].D*px + fSurf[i].E*py + ABCF*pz + fSurf[i].G;
1468 G4double D = B*B - A*C;
1469 if (D < 0) continue;
1470 G4double leng = 2.*std::sqrt(D)/std::abs(A);
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 () const
1483 {
1484 if ( (fpPolyhedron == nullptr)
1485 || fRebuildPolyhedron
1486 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1487 fpPolyhedron->GetNumberOfRotationSteps()) )
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(G4VGraphicsScene& scene) const
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() const
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 for smooth visualisation
1538 //
1539 G4double maxTwist = 0.;
1540 for(G4int i = 0; i < 4; ++i)
1541 {
1542 if (GetTwistAngle(i) > maxTwist) { maxTwist = GetTwistAngle(i); }
1543 }
1544
1545 // Computes bounding vectors for the shape
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*Dx)*fDz);
1553 if (subdivisions < 4) { subdivisions = 4; }
1554 if (subdivisions > 30) { subdivisions = 30; }
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, nFacets);
1562
1563 // Set vertices
1564 //
1565 G4int icur = 0;
1566 for (G4int i = 0; i < 4; ++i)
1567 {
1568 G4ThreeVector v(fVertices[i].x(),fVertices[i].y(),-fDz);
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)*(fVertices[j+4]-fVertices[j]);
1576 G4ThreeVector v(u.x(),u.y(),-fDz+cf*2*fDz*(i+1));
1577 polyhedron->SetVertex(++icur, v);
1578 }
1579 }
1580 for (G4int i = 4; i < 8; ++i)
1581 {
1582 G4ThreeVector v(fVertices[i].x(),fVertices[i].y(),fDz);
1583 polyhedron->SetVertex(++icur, v);
1584 }
1585
1586 // Set facets
1587 //
1588 icur = 0;
1589 polyhedron->SetFacet(++icur, 1, 4, 3, 2); // Z-plane
1590 for (G4int i = 0; i < subdivisions + 1; ++i)
1591 {
1592 G4int is = i*4;
1593 polyhedron->SetFacet(++icur, 5+is, 8+is, 4+is, 1+is);
1594 polyhedron->SetFacet(++icur, 8+is, 7+is, 3+is, 4+is);
1595 polyhedron->SetFacet(++icur, 7+is, 6+is, 2+is, 3+is);
1596 polyhedron->SetFacet(++icur, 6+is, 5+is, 1+is, 2+is);
1597 }
1598 polyhedron->SetFacet(++icur, 5+sub4, 6+sub4, 7+sub4, 8+sub4); // Z-plane
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 sign
1609 //
1610 void G4GenericTrap::WarningSignA(const G4String& method, const G4String& icase, G4double A,
1611 const G4ThreeVector& p, const G4ThreeVector& v) const
1612 {
1613 std::ostringstream message;
1614 message.precision(16);
1615 message << icase << " in " << GetName() << "\n"
1616 << " p" << p << "\n"
1617 << " v" << v << "\n"
1618 << " A = " << A << " (is "
1619 << ((A < 0) ? "negative, instead of positive)" : "positive, instead of negative)") << " !?\n";
1620 StreamInfo(message);
1621 const G4String function = "G4GenericTrap::DistanceTo" + method + "(p,v)";
1622 G4Exception(function, "GeomSolids1002", JustWarning, message );
1623 }
1624
1625 ////////////////////////////////////////////////////////////////////////
1626 //
1627 // Print out a warning if B has an unexpected sign
1628 //
1629 void G4GenericTrap::WarningSignB(const G4String& method, const G4String& icase,
1630 G4double f, G4double B,
1631 const G4ThreeVector& p, const G4ThreeVector& v) const
1632 {
1633 std::ostringstream message;
1634 message.precision(16);
1635 message << icase << " in " << GetName() << "\n"
1636 << " p" << p << "\n"
1637 << " v" << v << "\n"
1638 << " f = " << f << " B = " << B << " (is "
1639 << ((B < 0) ? "negative, instead of positive)" : "positive, instead of negative)") << " !?\n";
1640 StreamInfo(message);
1641 const G4String function = "G4GenericTrap::DistanceTo" + method + "(p,v)";
1642 G4Exception(function, "GeomSolids1002", JustWarning, message );
1643 }
1644
1645 ////////////////////////////////////////////////////////////////////////
1646 //
1647 // Print out a warning in DistanceToIn(p,v)
1648 //
1649 void G4GenericTrap::WarningDistanceToIn(G4int k, const G4ThreeVector& p, const G4ThreeVector& v,
1650 G4double tmin, G4double tmax,
1651 const G4double ttin[2], const G4double ttout[2]) 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) check = "\n !!! check point 0.5*(ttout[0] + ttin[1]) is NOT outside !!!";
1658 }
1659
1660 auto position = Inside(p);
1661 std::ostringstream message;
1662 message.precision(16);
1663 message << k << "_Unexpected sequence of intersections in solid: " << GetName() << " !?\n"
1664 << " position = " << ((position == kInside) ? "kInside" : ((position == kOutside) ? "kOutside" : "kSurface")) << "\n"
1665 << " p" << p << "\n"
1666 << " v" << v << "\n"
1667 << " range : [" << tmin << ", " << tmax << "]\n"
1668 << " ttin[2] : "
1669 << ((ttin[0] == DBL_MAX) ? kInfinity : ttin[0]) << ", "
1670 << ((ttin[1] == DBL_MAX) ? kInfinity : ttin[1]) << "\n"
1671 << " ttout[2] : "
1672 << ((ttout[0] == DBL_MAX) ? kInfinity : ttout[0]) << ", "
1673 << ((ttout[1] == DBL_MAX) ? kInfinity : ttout[1]) << check << "\n";
1674 StreamInfo(message);
1675 G4Exception("G4GenericTrap::DistanceToIn(p,v)", "GeomSolids1002",
1676 JustWarning, message );
1677 }
1678
1679 ////////////////////////////////////////////////////////////////////////
1680 //
1681 // Print out a warning in DistanceToOut(p,v)
1682 //
1683 void G4GenericTrap::WarningDistanceToOut(const G4ThreeVector& p,
1684 const G4ThreeVector& v,
1685 G4double tout) const
1686 {
1687 auto position = Inside(p);
1688 std::ostringstream message;
1689 message.precision(16);
1690 message << "Unexpected final tout = " << tout << " in solid: " << GetName() << " !?\n"
1691 << " position = " << ((position == kInside) ? "kInside" : ((position == kOutside) ? "kOutside" : "kSurface")) << "\n"
1692 << " p" << p << "\n"
1693 << " v" << v << "\n";
1694 StreamInfo(message);
1695 G4Exception("G4GenericTrap::DistanceToOut(p,v)", "GeomSolids1002",
1696 JustWarning, message );
1697 }
1698
1699 #endif
1700