Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4EllipticalTube.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4EllipticalTube implementation
 27 //
 28 // Author: David C. Williams (davidw@scipp.ucsc.edu)
 29 // Revision: Evgueni Tcherniaev (evgueni.tcherniaev@cern.ch), 23.12.2019
 30 // --------------------------------------------------------------------
 31 
 32 #include "G4EllipticalTube.hh"
 33 
 34 #if !(defined(G4GEOM_USE_UELLIPTICALTUBE) && defined(G4GEOM_USE_SYS_USOLIDS))
 35 
 36 #include "G4GeomTools.hh"
 37 #include "G4RandomTools.hh"
 38 #include "G4ClippablePolygon.hh"
 39 #include "G4AffineTransform.hh"
 40 #include "G4VoxelLimits.hh"
 41 #include "G4BoundingEnvelope.hh"
 42 
 43 #include "Randomize.hh"
 44 
 45 #include "G4VGraphicsScene.hh"
 46 #include "G4VisExtent.hh"
 47 
 48 #include "G4AutoLock.hh"
 49 
 50 namespace
 51 {
 52   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 53 }
 54 
 55 using namespace CLHEP;
 56 
 57 //////////////////////////////////////////////////////////////////////////
 58 //
 59 // Constructor
 60 
 61 G4EllipticalTube::G4EllipticalTube( const G4String &name,
 62                                           G4double Dx,
 63                                           G4double Dy,
 64                                           G4double Dz )
 65   : G4VSolid(name), fDx(Dx), fDy(Dy), fDz(Dz)
 66 {
 67   CheckParameters();
 68 }
 69 
 70 //////////////////////////////////////////////////////////////////////////
 71 //
 72 // Fake default constructor - sets only member data and allocates memory
 73 //                            for usage restricted to object persistency.
 74 
 75 G4EllipticalTube::G4EllipticalTube( __void__& a )
 76   : G4VSolid(a), halfTolerance(0.), fDx(0.), fDy(0.), fDz(0.),
 77     fRsph(0.), fDDx(0.), fDDy(0.), fSx(0.), fSy(0.), fR(0.),
 78     fQ1(0.), fQ2(0.), fScratch(0.)
 79 {
 80 }
 81 
 82 //////////////////////////////////////////////////////////////////////////
 83 //
 84 // Destructor
 85 
 86 G4EllipticalTube::~G4EllipticalTube()
 87 {
 88   delete fpPolyhedron; fpPolyhedron = nullptr;
 89 }
 90 
 91 //////////////////////////////////////////////////////////////////////////
 92 //
 93 // Copy constructor
 94 
 95 G4EllipticalTube::G4EllipticalTube(const G4EllipticalTube& rhs)
 96   : G4VSolid(rhs), halfTolerance(rhs.halfTolerance),
 97     fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
 98     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
 99     fRsph(rhs.fRsph), fDDx(rhs.fDDx), fDDy(rhs.fDDy),
100     fSx(rhs.fSx), fSy(rhs.fSy), fR(rhs.fR),
101     fQ1(rhs.fQ1), fQ2(rhs.fQ2), fScratch(rhs.fScratch)
102 {
103 }
104 
105 //////////////////////////////////////////////////////////////////////////
106 //
107 // Assignment operator
108 
109 G4EllipticalTube& G4EllipticalTube::operator = (const G4EllipticalTube& rhs) 
110 {
111    // Check assignment to self
112    //
113    if (this == &rhs)  { return *this; }
114 
115    // Copy base class data
116    //
117    G4VSolid::operator=(rhs);
118 
119    // Copy data
120    //
121    halfTolerance = rhs.halfTolerance;
122    fDx = rhs.fDx;
123    fDy = rhs.fDy;
124    fDz = rhs.fDz;
125    fCubicVolume = rhs.fCubicVolume;
126    fSurfaceArea = rhs.fSurfaceArea;
127 
128    fRsph = rhs.fRsph;
129    fDDx  = rhs.fDDx;
130    fDDy  = rhs.fDDy;
131    fSx   = rhs.fSx;
132    fSy   = rhs.fSy;
133    fR    = rhs.fR;
134    fQ1   = rhs.fQ1;
135    fQ2   = rhs.fQ2;
136    fScratch = rhs.fScratch;
137 
138    fRebuildPolyhedron = false;
139    delete fpPolyhedron; fpPolyhedron = nullptr;
140 
141    return *this;
142 }
143 
144 //////////////////////////////////////////////////////////////////////////
145 //
146 // Check dimensions
147 
148 void G4EllipticalTube::CheckParameters()
149 {
150   // Check dimensions
151   //
152   halfTolerance = 0.5*kCarTolerance; // half tolerance
153   G4double dmin = 2*kCarTolerance;
154   if (fDx < dmin || fDy < dmin || fDz < dmin)
155   {
156     std::ostringstream message;
157     message << "Invalid (too small or negative) dimensions for Solid: "
158             << GetName()
159             << "\n  Dx = " << fDx
160             << "\n  Dy = " << fDy
161             << "\n  Dz = " << fDz;
162     G4Exception("G4EllipticalTube::CheckParameters()", "GeomSolids0002",
163           FatalException, message);
164   }
165 
166   // Set pre-calculatated values
167   //
168   halfTolerance = 0.5*kCarTolerance; // half tolerance
169   fRsph = std::sqrt(fDx * fDx + fDy * fDy + fDz * fDz); // radius of surrounding sphere
170   fDDx = fDx * fDx; // X semi-axis squared
171   fDDy = fDy * fDy; // Y semi-axis squared
172 
173   fR = std::min(fDx, fDy); // resulting radius, after scaling elipse to circle
174   fSx = fR / fDx; // X scale factor
175   fSy = fR / fDy; // Y scale factor
176 
177   fQ1 = 0.5 / fR; // distance approxiamtion dist = Q1 * (x^2 + y^2) - Q2
178   fQ2 = 0.5 * (fR + halfTolerance * halfTolerance / fR);
179   fScratch = 2. * fR * fR * DBL_EPSILON; // scratch within calculation error thickness
180   // fScratch = (B * B / A) * (2. + halfTolerance / A) * halfTolerance; // alternative
181 }
182 
183 //////////////////////////////////////////////////////////////////////////
184 //
185 // Get bounding box
186 
187 void G4EllipticalTube::BoundingLimits( G4ThreeVector& pMin,
188                                        G4ThreeVector& pMax ) const
189 {
190   pMin.set(-fDx,-fDy,-fDz);
191   pMax.set( fDx, fDy, fDz);
192 }
193 
194 //////////////////////////////////////////////////////////////////////////
195 //
196 // Calculate extent under transform and specified limit
197 
198 G4bool
199 G4EllipticalTube::CalculateExtent( const EAxis pAxis,
200                                    const G4VoxelLimits& pVoxelLimit,
201                                    const G4AffineTransform& pTransform,
202                                          G4double& pMin, G4double& pMax ) const
203 {
204   G4ThreeVector bmin, bmax;
205   G4bool exist;
206 
207   // Check bounding box (bbox)
208   //
209   BoundingLimits(bmin,bmax);
210   G4BoundingEnvelope bbox(bmin,bmax);
211 #ifdef G4BBOX_EXTENT
212   return bbox.CalculateExtent(pAxis,pVoxelLimit, pTransform, pMin, pMax);
213 #endif
214   if (bbox.BoundingBoxVsVoxelLimits(pAxis, pVoxelLimit, pTransform, pMin, pMax))
215   {
216     return exist = pMin < pMax;
217   }
218 
219   G4double dx = fDx;
220   G4double dy = fDy;
221   G4double dz = fDz;
222 
223   // Set bounding envelope (benv) and calculate extent
224   //
225   const G4int NSTEPS = 24; // number of steps for whole circle
226   G4double ang = twopi/NSTEPS;
227 
228   G4double sinHalf = std::sin(0.5*ang);
229   G4double cosHalf = std::cos(0.5*ang);
230   G4double sinStep = 2.*sinHalf*cosHalf;
231   G4double cosStep = 1. - 2.*sinHalf*sinHalf;
232   G4double sx = dx/cosHalf;
233   G4double sy = dy/cosHalf;
234 
235   G4double sinCur = sinHalf;
236   G4double cosCur = cosHalf;
237   G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
238   for (G4int k=0; k<NSTEPS; ++k)
239   {
240     baseA[k].set(sx*cosCur,sy*sinCur,-dz);
241     baseB[k].set(sx*cosCur,sy*sinCur, dz);
242 
243     G4double sinTmp = sinCur;
244     sinCur = sinCur*cosStep + cosCur*sinStep;
245     cosCur = cosCur*cosStep - sinTmp*sinStep;
246   }
247 
248   std::vector<const G4ThreeVectorList *> polygons(2);
249   polygons[0] = &baseA;
250   polygons[1] = &baseB;
251   G4BoundingEnvelope benv(bmin, bmax, polygons);
252   exist = benv.CalculateExtent(pAxis, pVoxelLimit, pTransform, pMin, pMax);
253   return exist;
254 }
255 
256 //////////////////////////////////////////////////////////////////////////
257 //
258 // Determine where is point: inside, outside or on surface
259 //
260 
261 EInside G4EllipticalTube::Inside( const G4ThreeVector& p ) const
262 {
263   G4double x = p.x() * fSx;
264   G4double y = p.y() * fSy;
265   G4double distR = fQ1 * (x * x + y * y) - fQ2;
266   G4double distZ = std::abs(p.z()) - fDz;
267   G4double dist = std::max(distR, distZ);
268 
269   if (dist > halfTolerance) return kOutside;
270   return (dist > -halfTolerance) ? kSurface : kInside;
271 }
272 
273 //////////////////////////////////////////////////////////////////////////
274 //
275 // Return unit normal at surface closest to p
276 
277 G4ThreeVector G4EllipticalTube::SurfaceNormal( const G4ThreeVector& p ) const
278 {
279   G4ThreeVector norm(0, 0, 0);
280   G4int nsurf = 0;
281 
282   // check lateral surface
283   G4double x = p.x() * fSx;
284   G4double y = p.y() * fSy;
285   G4double distR = fQ1 * (x * x + y * y) - fQ2;
286   if (std::abs(distR) <= halfTolerance)
287   {
288     norm = G4ThreeVector(p.x() * fDDy, p.y() * fDDx, 0.).unit();
289     ++nsurf;
290   }
291 
292   // check lateral bases
293   G4double distZ = std::abs(p.z()) - fDz;
294   if (std::abs(distZ) <= halfTolerance)
295   {
296     norm.setZ(p.z() < 0 ? -1. : 1.);
297     ++nsurf;
298   }
299 
300   // return normal
301   if (nsurf == 1) return norm;
302   else if (nsurf > 1) return norm.unit(); // edge
303   else
304   {
305     // Point is not on the surface
306     //
307 #ifdef G4SPECDEBUG
308     std::ostringstream message;
309     G4long oldprc = message.precision(16);
310     message << "Point p is not on surface (!?) of solid: "
311             << GetName() << G4endl;
312     message << "Position:\n";
313     message << "   p.x() = " << p.x()/mm << " mm\n";
314     message << "   p.y() = " << p.y()/mm << " mm\n";
315     message << "   p.z() = " << p.z()/mm << " mm";
316     G4cout.precision(oldprc);
317     G4Exception("G4EllipticalTube::SurfaceNormal(p)", "GeomSolids1002",
318                 JustWarning, message );
319     DumpInfo();
320 #endif
321     return ApproxSurfaceNormal(p);
322   }
323 }
324 
325 //////////////////////////////////////////////////////////////////////////
326 //
327 // Find surface nearest to point and return corresponding normal.
328 // The algorithm is similar to the algorithm used in Inside().
329 // This method normally should not be called.
330 
331 G4ThreeVector
332 G4EllipticalTube::ApproxSurfaceNormal( const G4ThreeVector& p ) const
333 {
334   G4double x = p.x() * fSx;
335   G4double y = p.y() * fSy;
336   G4double distR = fQ1 * (x * x + y * y) - fQ2;
337   G4double distZ = std::abs(p.z()) - fDz;
338   if (distR > distZ && (x * x + y * y) > 0)
339     return G4ThreeVector(p.x() * fDDy, p.y() * fDDx, 0.).unit();
340   else
341     return {0, 0, (p.z() < 0 ? -1. : 1.)};
342 }
343 
344 //////////////////////////////////////////////////////////////////////////
345 //
346 // Calculate distance to shape from outside, along normalised vector,
347 // return kInfinity if no intersection, or distance < halfTolerance
348 
349 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p,
350                                          const G4ThreeVector& v ) const
351 {
352   G4double offset = 0.;
353   G4ThreeVector pcur = p;
354 
355   // Check if point is flying away
356   //
357   G4double safex = std::abs(pcur.x()) - fDx;
358   G4double safey = std::abs(pcur.y()) - fDy;
359   G4double safez = std::abs(pcur.z()) - fDz;
360 
361   if (safez >= -halfTolerance && pcur.z() * v.z() >= 0.) return kInfinity;
362   if (safey >= -halfTolerance && pcur.y() * v.y() >= 0.) return kInfinity;
363   if (safex >= -halfTolerance && pcur.x() * v.x() >= 0.) return kInfinity;
364 
365   // Relocate point, if required
366   //
367   G4double Dmax = 32. * fRsph;
368   if (std::max(std::max(safex, safey), safez) > Dmax)
369   {
370     offset = (1. - 1.e-08) * pcur.mag() - 2. * fRsph;
371     pcur += offset * v;
372     G4double dist = DistanceToIn(pcur, v);
373     return (dist == kInfinity) ? kInfinity : dist + offset;
374   }
375 
376   // Scale elliptical tube to cylinder
377   //
378   G4double px = pcur.x() * fSx;
379   G4double py = pcur.y() * fSy;
380   G4double pz = pcur.z();
381   G4double vx = v.x() * fSx;
382   G4double vy = v.y() * fSy;
383   G4double vz = v.z();
384 
385   // Set coefficients of quadratic equation: A t^2 + 2B t + C = 0
386   //
387   G4double rr = px * px + py * py;
388   G4double A  = vx * vx + vy * vy;
389   G4double B  = px * vx + py * vy;
390   G4double C  = rr - fR * fR;
391   G4double D  = B * B - A * C;
392 
393   // Check if point is flying away relative to lateral surface
394   //
395   G4double distR  = fQ1 * rr - fQ2;
396   G4bool parallelToZ = (A < DBL_EPSILON || std::abs(vz) >= 1.);
397   if (distR >= -halfTolerance && (B >= 0. || parallelToZ)) return kInfinity;
398 
399   // Find intersection with Z planes
400   //
401   G4double invz  = (vz == 0) ? DBL_MAX : -1./vz;
402   G4double dz    = std::copysign(fDz, invz);
403   G4double tzmin = (pz - dz) * invz;
404   G4double tzmax = (pz + dz) * invz;
405 
406   // Solve qudratic equation. There are two cases special where D <= 0:
407   //   1) trajectory parallel to Z axis (A = 0, B = 0, C - any, D = 0)
408   //   2) touch (D = 0) or no intersection (D < 0) with lateral surface
409   //
410   if (parallelToZ) return (tzmin<halfTolerance) ? offset : tzmin + offset; // 1)
411   if (D <= A * A * fScratch) return kInfinity; // 2)
412 
413   // Find roots of quadratic equation
414   G4double tmp = -B - std::copysign(std::sqrt(D), B);
415   G4double t1 = tmp / A;
416   G4double t2 = C / tmp;
417   G4double trmin = std::min(t1, t2);
418   G4double trmax = std::max(t1, t2);
419 
420   // Return distance
421   G4double tin  = std::max(tzmin, trmin);
422   G4double tout = std::min(tzmax, trmax);
423 
424   if (tout <= tin + halfTolerance) return kInfinity; // touch or no hit
425   return (tin<halfTolerance) ? offset : tin + offset;
426 }
427 
428 //////////////////////////////////////////////////////////////////////////
429 //
430 // Estimate distance to the surface from outside,
431 // returns 0 if point is inside
432 
433 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p ) const
434 {
435   // safety distance to bounding box
436   G4double distX = std::abs(p.x()) - fDx;
437   G4double distY = std::abs(p.y()) - fDy;
438   G4double distZ = std::abs(p.z()) - fDz;
439   G4double distB = std::max(std::max(distX, distY), distZ);
440   // return (distB < 0) ? 0 : distB;
441 
442   // safety distance to lateral surface
443   G4double x = p.x() * fSx;
444   G4double y = p.y() * fSy;
445   G4double distR = std::sqrt(x * x + y * y) - fR;
446 
447   // return SafetyToIn
448   G4double dist = std::max(distB, distR);
449   return (dist < 0) ? 0 : dist;
450 }
451 
452 //////////////////////////////////////////////////////////////////////////
453 //
454 // Calculate distance to shape from inside and find normal
455 // at exit point, if required
456 // - when leaving the surface, return 0
457 
458 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p,
459                                           const G4ThreeVector& v,
460                                           const G4bool calcNorm,
461                                                 G4bool* validNorm,
462                                                 G4ThreeVector* n ) const
463 {
464   // Check if point flying away relative to Z planes
465   //
466   G4double pz = p.z();
467   G4double vz = v.z();
468   G4double distZ = std::abs(pz) - fDz;
469   if (distZ >= -halfTolerance && pz * vz > 0)
470   {
471     if (calcNorm)
472     {
473       *validNorm = true;
474       n->set(0, 0, (pz < 0) ? -1. : 1.);
475     }
476     return 0.;
477   }
478   G4double tzmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz, vz) - pz) / vz;
479 
480   // Scale elliptical tube to cylinder
481   //
482   G4double px = p.x() * fSx;
483   G4double py = p.y() * fSy;
484   G4double vx = v.x() * fSx;
485   G4double vy = v.y() * fSy;
486 
487   // Check if point is flying away relative to lateral surface
488   //
489   G4double rr = px * px + py * py;
490   G4double B  = px * vx + py * vy;
491   G4double distR  = fQ1 * rr - fQ2;
492   if (distR >= -halfTolerance && B > 0.)
493   {
494     if (calcNorm)
495     {
496       *validNorm = true;
497       *n = G4ThreeVector(px * fDDy, py * fDDx, 0.).unit();
498     }
499     return 0.;
500   }
501 
502   // Just in case check if point is outside, normally it should never be
503   //
504   if (std::max(distZ, distR) > halfTolerance)
505   {
506 #ifdef G4SPECDEBUG
507     std::ostringstream message;
508     G4long oldprc = message.precision(16);
509     message << "Point p is outside (!?) of solid: "
510             << GetName() << G4endl;
511     message << "Position:  " << p << G4endl;;
512     message << "Direction: " << v;
513     G4cout.precision(oldprc);
514     G4Exception("G4EllipticalTube::DistanceToOut(p,v)", "GeomSolids1002",
515                 JustWarning, message );
516     DumpInfo();
517 #endif
518     if (calcNorm)
519     {
520       *validNorm = true;
521       *n = ApproxSurfaceNormal(p);
522     }
523     return 0.;
524   }
525 
526   // Set coefficients of quadratic equation: A t^2 + 2B t + C = 0
527   //
528   G4double A  = vx * vx + vy * vy;
529   G4double C  = rr - fR * fR;
530   G4double D  = B * B - A * C;
531 
532   // Solve qudratic equation. There are two special cases where D <= 0:
533   //   1) trajectory parallel to Z axis (A = 0, B = 0, C - any, D = 0)
534   //   2) touch (D = 0) or no intersection (D < 0) with lateral surface
535   //
536   G4bool parallelToZ = (A < DBL_EPSILON || std::abs(vz) >= 1.);
537   if (parallelToZ) // 1)
538   {
539     if (calcNorm)
540     {
541       *validNorm = true;
542       n->set(0, 0, (vz < 0) ? -1. : 1.);
543     }
544     return tzmax;
545   }
546   if (D <= A * A * fScratch) // 2)
547   {
548     if (calcNorm)
549     {
550       *validNorm = true;
551       *n = G4ThreeVector(px * fDDy, py * fDDx, 0.).unit();
552     }
553     return 0.;
554   }
555 
556   // Find roots of quadratic equation
557   G4double tmp = -B - std::copysign(std::sqrt(D), B);
558   G4double t1 = tmp / A;
559   G4double t2 = C / tmp;
560   G4double trmax = std::max(t1, t2);
561 
562   // Return distance
563   G4double tmax = std::min(tzmax, trmax);
564 
565   // Set normal, if required, and return distance
566   //
567   if (calcNorm)
568   {
569     *validNorm = true;
570     G4ThreeVector pnew = p + tmax * v;
571     if (tmax == tzmax)
572       n->set(0, 0, (pnew.z() < 0) ? -1. : 1.);
573     else
574       *n = G4ThreeVector(pnew.x() * fDDy, pnew.y() * fDDx, 0.).unit();
575   }
576   return tmax;
577 }
578 
579 //////////////////////////////////////////////////////////////////////////
580 //
581 // Estimate distance to the surface from inside,
582 // returns 0 if point is outside
583 //
584 
585 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p ) const
586 {
587 #ifdef G4SPECDEBUG
588   if( Inside(p) == kOutside )
589   {
590     std::ostringstream message;
591     G4long oldprc = message.precision(16);
592     message << "Point p is outside (!?) of solid: " << GetName() << "\n"
593             << "Position:\n"
594             << "   p.x() = "  << p.x()/mm << " mm\n"
595             << "   p.y() = "  << p.y()/mm << " mm\n"
596             << "   p.z() = "  << p.z()/mm << " mm";
597     message.precision(oldprc) ;
598     G4Exception("G4ElliptocalTube::DistanceToOut(p)", "GeomSolids1002",
599                 JustWarning, message);
600     DumpInfo();
601   }
602 #endif
603   // safety distance to Z-bases
604   G4double distZ = fDz - std::abs(p.z());
605 
606   // safety distance lateral surface
607   G4double x = p.x() * fSx;
608   G4double y = p.y() * fSy;
609   G4double distR = fR - std::sqrt(x * x + y * y);
610 
611   // return SafetyToOut
612   G4double dist = std::min(distZ, distR);
613   return (dist < 0) ? 0 : dist;
614 }
615 
616 //////////////////////////////////////////////////////////////////////////
617 //
618 // GetEntityType
619 
620 G4GeometryType G4EllipticalTube::GetEntityType() const
621 {
622   return {"G4EllipticalTube"};
623 }
624 
625 //////////////////////////////////////////////////////////////////////////
626 //
627 // Make a clone of the object
628 
629 G4VSolid* G4EllipticalTube::Clone() const
630 {
631   return new G4EllipticalTube(*this);
632 }
633 
634 //////////////////////////////////////////////////////////////////////////
635 //
636 // Return volume
637 
638 G4double G4EllipticalTube::GetCubicVolume()
639 {
640   if (fCubicVolume == 0.)
641   {
642     fCubicVolume = twopi * fDx * fDy * fDz;
643   }
644   return fCubicVolume;
645 }
646 
647 //////////////////////////////////////////////////////////////////////////
648 //
649 // Return cached surface area
650 
651 G4double G4EllipticalTube::GetCachedSurfaceArea() const
652 {
653   G4ThreadLocalStatic G4double cached_Dx = 0;
654   G4ThreadLocalStatic G4double cached_Dy = 0;
655   G4ThreadLocalStatic G4double cached_Dz = 0;
656   G4ThreadLocalStatic G4double cached_area = 0;
657   if (cached_Dx != fDx || cached_Dy != fDy || cached_Dz != fDz)
658   {
659     cached_Dx = fDx;
660     cached_Dy = fDy;
661     cached_Dz = fDz;
662     cached_area = 2.*(pi*fDx*fDy + G4GeomTools::EllipsePerimeter(fDx, fDy)*fDz);
663   }
664   return cached_area;
665 }
666 
667 //////////////////////////////////////////////////////////////////////////
668 //
669 // Return surface area
670 
671 G4double G4EllipticalTube::GetSurfaceArea()
672 {
673   if(fSurfaceArea == 0.)
674   {
675     fSurfaceArea = GetCachedSurfaceArea();
676   }
677   return fSurfaceArea;
678 }
679 
680 //////////////////////////////////////////////////////////////////////////
681 //
682 // Stream object contents to output stream
683 
684 std::ostream& G4EllipticalTube::StreamInfo(std::ostream& os) const
685 {
686   G4long oldprc = os.precision(16);
687   os << "-----------------------------------------------------------\n"
688      << "    *** Dump for solid - " << GetName() << " ***\n"
689      << "    ===================================================\n"
690      << " Solid type: G4EllipticalTube\n"
691      << " Parameters: \n"
692      << "    length Z: " << fDz/mm << " mm \n"
693      << "    lateral surface equation: \n"
694      << "       (X / " << fDx << ")^2 + (Y / " << fDy << ")^2 = 1 \n"
695      << "-----------------------------------------------------------\n";
696   os.precision(oldprc);
697 
698   return os;
699 }
700 
701 
702 //////////////////////////////////////////////////////////////////////////
703 //
704 // Pick up a random point on the surface 
705 
706 G4ThreeVector G4EllipticalTube::GetPointOnSurface() const
707 {
708   // Select surface (0 - base at -Z, 1 - base at +Z, 2 - lateral surface)
709   //
710   G4double sbase = pi * fDx * fDy;
711   G4double ssurf = GetCachedSurfaceArea();
712   G4double select = ssurf * G4UniformRand();
713 
714   G4int k = 0;
715   if (select > sbase) k = 1;
716   if (select > 2. * sbase) k = 2;
717 
718   // Pick random point on selected surface (rejection sampling)
719   //
720   G4ThreeVector p;
721   switch (k) {
722     case 0: // base at -Z
723     {
724       G4TwoVector rho = G4RandomPointInEllipse(fDx, fDy);
725       p.set(rho.x(), rho.y(), -fDz);
726       break;
727     }
728     case 1: // base at +Z
729     {
730       G4TwoVector rho = G4RandomPointInEllipse(fDx, fDy);
731       p.set(rho.x(), rho.y(), fDz);
732       break;
733     }
734     case 2: // lateral surface
735     {
736       G4TwoVector rho = G4RandomPointOnEllipse(fDx, fDy);
737       p.set(rho.x(), rho.y(), (2. * G4UniformRand() - 1.) * fDz);
738       break;
739     }
740   }
741   return p;
742 }
743 
744 
745 //////////////////////////////////////////////////////////////////////////
746 //
747 // CreatePolyhedron
748 
749 G4Polyhedron* G4EllipticalTube::CreatePolyhedron() const
750 {
751   // create cylinder with radius=1...
752   //
753   G4Polyhedron* eTube = new G4PolyhedronTube(0., 1., fDz);
754 
755   // apply non-uniform scaling...
756   //
757   eTube->Transform(G4Scale3D(fDx, fDy, 1.));
758   return eTube;
759 }
760 
761 //////////////////////////////////////////////////////////////////////////
762 //
763 // GetPolyhedron
764 
765 G4Polyhedron* G4EllipticalTube::GetPolyhedron () const
766 {
767   if (fpPolyhedron == nullptr ||
768       fRebuildPolyhedron ||
769       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
770       fpPolyhedron->GetNumberOfRotationSteps())
771   {
772     G4AutoLock l(&polyhedronMutex);
773     delete fpPolyhedron;
774     fpPolyhedron = CreatePolyhedron();
775     fRebuildPolyhedron = false;
776     l.unlock();
777   }
778   return fpPolyhedron;
779 }
780 
781 //////////////////////////////////////////////////////////////////////////
782 //
783 // DescribeYourselfTo
784 
785 void G4EllipticalTube::DescribeYourselfTo( G4VGraphicsScene& scene ) const
786 {
787   scene.AddSolid (*this);
788 }
789 
790 //////////////////////////////////////////////////////////////////////////
791 //
792 // GetExtent
793 
794 G4VisExtent G4EllipticalTube::GetExtent() const
795 {
796   return { -fDx, fDx, -fDy, fDy, -fDz, fDz };
797 }
798 
799 #endif // !defined(G4GEOM_USE_UELLIPTICALTUBE) || !defined(G4GEOM_USE_SYS_USOLIDS)
800