Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4Paraboloid.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 // Implementation for G4Paraboloid class
 27 //
 28 // Author : Lukas Lindroos (CERN), July 2007
 29 // Revised: Tatiana Nikitina (CERN)
 30 // --------------------------------------------------------------------
 31 
 32 #include "globals.hh"
 33 
 34 #include "G4Paraboloid.hh"
 35 
 36 #if !(defined(G4GEOM_USE_UPARABOLOID) && defined(G4GEOM_USE_SYS_USOLIDS))
 37 
 38 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"
 40 #include "G4BoundingEnvelope.hh"
 41 
 42 #include "meshdefs.hh"
 43 
 44 #include "Randomize.hh"
 45 
 46 #include "G4VGraphicsScene.hh"
 47 #include "G4VisExtent.hh"
 48 
 49 #include "G4AutoLock.hh"
 50 
 51 namespace
 52 {
 53   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 54 }
 55 
 56 using namespace CLHEP;
 57 
 58 ///////////////////////////////////////////////////////////////////////////////
 59 //
 60 // constructor - check parameters
 61 //
 62 G4Paraboloid::G4Paraboloid(const G4String& pName,
 63                                  G4double pDz,
 64                                  G4double pR1,
 65                                  G4double pR2)
 66  : G4VSolid(pName)
 67 {
 68   if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) )
 69   {
 70     std::ostringstream message;
 71     message << "Invalid dimensions. Negative Input Values or R1>=R2 - "
 72             << GetName();
 73     G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002", 
 74                 FatalErrorInArgument, message,
 75                 "Z half-length must be larger than zero or R1>=R2.");
 76   }
 77 
 78   r1 = pR1;
 79   r2 = pR2;
 80   dz = pDz;
 81 
 82   // r1^2 = k1 * (-dz) + k2
 83   // r2^2 = k1 * ( dz) + k2
 84   // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
 85   // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
 86 
 87   k1 = (r2 * r2 - r1 * r1) / 2 / dz;
 88   k2 = (r2 * r2 + r1 * r1) / 2;
 89 }
 90 
 91 ///////////////////////////////////////////////////////////////////////////////
 92 //
 93 // Fake default constructor - sets only member data and allocates memory
 94 //                            for usage restricted to object persistency.
 95 //
 96 G4Paraboloid::G4Paraboloid( __void__& a )
 97   : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
 98 {
 99 }
100 
101 ///////////////////////////////////////////////////////////////////////////////
102 //
103 // Destructor
104 //
105 G4Paraboloid::~G4Paraboloid()
106 {
107   delete fpPolyhedron; fpPolyhedron = nullptr;
108 }
109 
110 ///////////////////////////////////////////////////////////////////////////////
111 //
112 // Copy constructor
113 //
114 G4Paraboloid::G4Paraboloid(const G4Paraboloid& rhs)
115   : G4VSolid(rhs),
116     fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume),
117     dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2)
118 {
119 }
120 
121 ///////////////////////////////////////////////////////////////////////////////
122 //
123 // Assignment operator
124 //
125 G4Paraboloid& G4Paraboloid::operator = (const G4Paraboloid& rhs) 
126 {
127    // Check assignment to self
128    //
129    if (this == &rhs)  { return *this; }
130 
131    // Copy base class data
132    //
133    G4VSolid::operator=(rhs);
134 
135    // Copy data
136    //
137    fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
138    dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2;
139    fRebuildPolyhedron = false;
140    delete fpPolyhedron; fpPolyhedron = nullptr;
141 
142    return *this;
143 }
144 
145 ///////////////////////////////////////////////////////////////////////////////
146 //
147 // Get bounding box
148 //
149 void G4Paraboloid::BoundingLimits(G4ThreeVector& pMin,
150                                   G4ThreeVector& pMax) const
151 {
152   pMin.set(-r2,-r2,-dz);
153   pMax.set( r2, r2, dz);
154 
155   // Check correctness of the bounding box
156   //
157   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
158   {
159     std::ostringstream message;
160     message << "Bad bounding box (min >= max) for solid: "
161             << GetName() << " !"
162             << "\npMin = " << pMin
163             << "\npMax = " << pMax;
164     G4Exception("G4Paraboloid::BoundingLimits()", "GeomMgt0001",
165                 JustWarning, message);
166     DumpInfo();
167   }
168 }
169 
170 ///////////////////////////////////////////////////////////////////////////////
171 //
172 // Calculate extent under transform and specified limit
173 //
174 G4bool
175 G4Paraboloid::CalculateExtent(const EAxis pAxis,
176                               const G4VoxelLimits& pVoxelLimit,
177                               const G4AffineTransform& pTransform,
178                                     G4double& pMin, G4double& pMax) const
179 {
180   G4ThreeVector bmin, bmax;
181 
182   // Get bounding box
183   BoundingLimits(bmin,bmax);
184 
185   // Find extent
186   G4BoundingEnvelope bbox(bmin,bmax);
187   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
188 }
189 
190 ///////////////////////////////////////////////////////////////////////////////
191 //
192 // Return whether point inside/outside/on surface
193 //
194 EInside G4Paraboloid::Inside(const G4ThreeVector& p) const
195 {
196   // First check is  the point is above or below the solid.
197   //
198   if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; }
199 
200   G4double rho2 = p.perp2(),
201            rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance),
202            A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
203  
204   if(A < 0 && sqr(A) > rhoSurfTimesTol2)
205   {
206     // Actually checking rho < radius of paraboloid at z = p.z().
207     // We're either inside or in lower/upper cutoff area.
208    
209     if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
210     {
211       // We're in the upper/lower cutoff area, sides have a paraboloid shape
212       // maybe further checks should be made to make these nicer
213 
214       return kSurface;
215     }
216     else
217     {
218       return kInside;
219     }
220   }
221   else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
222   {
223     // We're in the parabolic surface.
224 
225     return kSurface;
226   }
227   else
228   {
229     return kOutside;
230   }
231 }
232 
233 ///////////////////////////////////////////////////////////////////////////////
234 //
235 // SurfaceNormal
236 //
237 G4ThreeVector G4Paraboloid::SurfaceNormal( const G4ThreeVector& p) const
238 {
239   G4ThreeVector n(0, 0, 0);
240   if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance)
241   {
242     // If above or below just return normal vector for the cutoff plane.
243 
244     n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z()));
245   }
246   else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
247   {
248     // This means we're somewhere in the plane z = dz or z = -dz.
249     // (As far as the program is concerned anyway.
250 
251     if(p.z() < 0) // Are we in upper or lower plane?
252     {
253       if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance))
254       {
255         n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit();
256       }
257       else if(r1 < 0.5 * kCarTolerance
258            || p.perp2() > sqr(r1 - 0.5 * kCarTolerance))
259       {
260         n = G4ThreeVector(p.x(), p.y(), 0.).unit()
261           + G4ThreeVector(0., 0., -1.).unit();
262         n = n.unit();
263       }
264       else
265       {
266         n = G4ThreeVector(0., 0., -1.);
267       }
268     }
269     else
270     {
271       if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance))
272       {
273         n = G4ThreeVector(p.x(), p.y(), 0.).unit();
274       }
275       else if(r2 < 0.5 * kCarTolerance
276            || p.perp2() > sqr(r2 - 0.5 * kCarTolerance))
277       {
278         n = G4ThreeVector(p.x(), p.y(), 0.).unit()
279           + G4ThreeVector(0., 0., 1.).unit();
280         n = n.unit();
281       }
282       else
283       {
284         n = G4ThreeVector(0., 0., 1.);
285       }
286     }
287   }
288   else
289   {
290     G4double rho2 = p.perp2();
291     G4double rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance);
292     G4double A = rho2 - ((k1 *p.z() + k2)
293                + 0.25 * kCarTolerance * kCarTolerance);
294 
295     if(A < 0 && sqr(A) > rhoSurfTimesTol2)
296     {
297       // Actually checking rho < radius of paraboloid at z = p.z().
298       // We're inside.
299 
300       if(p.mag2() != 0) { n = p.unit(); }
301     }
302     else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
303     {
304       // We're in the parabolic surface.
305 
306       n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
307     }
308     else
309     {
310       n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
311     }
312   }
313 
314   if(n.mag2() == 0)
315   {
316     std::ostringstream message;
317     message << "No normal defined for this point p." << G4endl
318             << "          p = " << 1 / mm * p << " mm";
319     G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002",
320                 JustWarning, message);
321   }
322   return n;
323 }
324 
325 ///////////////////////////////////////////////////////////////////////////////
326 //
327 // Calculate distance to shape from outside, along normalised vector
328 // - return kInfinity if no intersection
329 //
330 //
331 G4double G4Paraboloid::DistanceToIn( const G4ThreeVector& p,
332                                      const G4ThreeVector& v  ) const
333 {
334   G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
335   G4double tol2 = kCarTolerance*kCarTolerance;
336   G4double tolh = 0.5*kCarTolerance;
337 
338   if((r2 != 0.0) && p.z() > - tolh + dz) 
339   {
340     // If the points is above check for intersection with upper edge.
341 
342     if(v.z() < 0)
343     {
344       G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz.
345       if(sqr(p.x() + v.x()*intersection)
346        + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance))
347       {
348         if(p.z() < tolh + dz)
349           { return 0; }
350         else
351           { return intersection; }
352       }
353     }
354     else  // Direction away, no possibility of intersection
355     {
356       return kInfinity;
357     }
358   }
359   else if((r1 != 0.0) && p.z() < tolh - dz)
360   {
361     // If the points is belove check for intersection with lower edge.
362 
363     if(v.z() > 0)
364     {
365       G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz.
366       if(sqr(p.x() + v.x()*intersection)
367        + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance))
368       {
369         if(p.z() > -tolh - dz)
370         {
371           return 0;
372         }
373         else
374         {
375           return intersection;
376         }
377       }
378     }
379     else  // Direction away, no possibility of intersection
380     {
381       return kInfinity;
382     }
383   }
384 
385   G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(),
386            vRho2 = v.perp2(), intersection,
387            B = (k1 * p.z() + k2 - rho2) * vRho2;
388 
389   if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) )
390     || (p.z() < - dz+kCarTolerance)
391     || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside.
392   {
393     // Is there a problem with squaring rho twice?
394 
395     if(vRho2<tol2) // Needs to be treated seperately.
396     {
397       intersection = ((rho2 - k2)/k1 - p.z())/v.z();
398       if(intersection < 0) { return kInfinity; }
399       else if(std::fabs(p.z() + v.z() * intersection) <= dz)
400       {
401         return intersection;
402       }
403       else
404       {
405         return kInfinity;
406       }
407     }
408     else if(A*A + B < 0) // No real intersections.
409     {
410       return kInfinity;
411     }
412     else
413     {
414       intersection = (A - std::sqrt(B + sqr(A))) / vRho2;
415       if(intersection < 0)
416       {
417         return kInfinity;
418       }
419       else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh)
420       {
421         return intersection;
422       }
423       else
424       {
425         return kInfinity;
426       }
427     }
428   }
429   else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2)
430   {
431     // If this is true we're somewhere in the border.
432 
433     G4ThreeVector normal(p.x(), p.y(), -k1/2);
434     if(normal.dot(v) <= 0)
435       { return 0; }
436     else
437       { return kInfinity; }
438   }
439   else
440   {
441     std::ostringstream message;
442     if(Inside(p) == kInside)
443     {
444       message << "Point p is inside! - " << GetName() << G4endl;
445     }
446     else
447     {
448       message << "Likely a problem in this function, for solid: " << GetName()
449               << G4endl;
450     }
451     message << "          p = " << p * (1/mm) << " mm" << G4endl
452             << "          v = " << v * (1/mm) << " mm";
453     G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002",
454                 JustWarning, message);
455     return 0;
456   }
457 }
458 
459 ///////////////////////////////////////////////////////////////////////////////
460 //
461 // Calculate distance (<= actual) to closest surface of shape from outside
462 // - Return zero if point inside
463 //
464 G4double G4Paraboloid::DistanceToIn(const G4ThreeVector& p) const
465 {
466   G4double safz = -dz+std::fabs(p.z());
467   if(safz<0.) { safz=0.; }
468   G4double safr = kInfinity;
469 
470   G4double rho = p.x()*p.x()+p.y()*p.y();
471   G4double paraRho = (p.z()-k2)/k1;
472   G4double sqrho = std::sqrt(rho);
473 
474   if(paraRho<0.)
475   {
476     safr=sqrho-r2;
477     if(safr>safz) { safz=safr; }
478     return safz;
479   }
480 
481   G4double sqprho = std::sqrt(paraRho);
482   G4double dRho = sqrho-sqprho;
483   if(dRho<0.) { return safz; }
484 
485   G4double talf = -2.*k1*sqprho;
486   G4double tmp  = 1+talf*talf;
487   if(tmp<0.) { return safz; }
488 
489   G4double salf = talf/std::sqrt(tmp);
490   safr = std::fabs(dRho*salf);
491   if(safr>safz) { safz=safr; }
492 
493   return safz;
494 }
495 
496 ///////////////////////////////////////////////////////////////////////////////
497 //
498 // Calculate distance to surface of shape from 'inside'
499 //
500 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p,
501                                     const G4ThreeVector& v,
502                                     const G4bool calcNorm,
503                                           G4bool* validNorm,
504                                           G4ThreeVector* n ) const
505 {
506   G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
507   G4double vRho2 = v.perp2(), intersection;
508   G4double tol2 = kCarTolerance*kCarTolerance;
509   G4double tolh = 0.5*kCarTolerance;
510 
511   if(calcNorm) { *validNorm = false; }
512 
513   // We have that the particle p follows the line x = p + s * v
514   // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and
515   // z = p.z() + s * v.z()
516   // The equation for all points on the surface (surface expanded for
517   // to include all z) x^2 + y^2 = k1 * z + k2 => .. =>
518   // => s = (A +- std::sqrt(A^2 + B)) / vRho2
519   // where:
520   //
521   G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
522   //
523   // and:
524   //
525   G4double B = (-rho2 + paraRho2) * vRho2;
526 
527   if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2
528     && std::fabs(p.z()) < dz - kCarTolerance)
529   {
530     // Make sure it's safely inside.
531 
532     if(v.z() > 0)
533     {
534       // It's heading upwards, check where it colides with the plane z = dz.
535       // When it does, is that in the surface of the paraboloid.
536       // z = p.z() + variable * v.z() for all points the particle can go.
537       // => variable = (z - p.z()) / v.z() so intersection must be:
538 
539       intersection = (dz - p.z()) / v.z();
540       G4ThreeVector ip = p + intersection * v; // Point of intersection.
541 
542       if(ip.perp2() < sqr(r2 + kCarTolerance))
543       {
544         if(calcNorm)
545         {
546           *n = G4ThreeVector(0, 0, 1);
547           if(r2 < tolh || ip.perp2() > sqr(r2 - tolh))
548           {
549             *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
550             *n = n->unit();
551           }
552           *validNorm = true;
553         }
554         return intersection;
555       }
556     }
557     else if(v.z() < 0)
558     {
559       // It's heading downwards, check were it colides with the plane z = -dz.
560       // When it does, is that in the surface of the paraboloid.
561       // z = p.z() + variable * v.z() for all points the particle can go.
562       // => variable = (z - p.z()) / v.z() so intersection must be:
563 
564       intersection = (-dz - p.z()) / v.z();
565       G4ThreeVector ip = p + intersection * v; // Point of intersection.
566 
567       if(ip.perp2() < sqr(r1 + tolh))
568       {
569         if(calcNorm)
570         {
571           *n = G4ThreeVector(0, 0, -1);
572           if(r1 < tolh || ip.perp2() > sqr(r1 - tolh))
573           {
574             *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
575             *n = n->unit();
576           }
577           *validNorm = true;
578         }
579         return intersection;
580       }
581     }
582 
583     // Now check for collisions with paraboloid surface.
584 
585     if(vRho2 == 0) // Needs to be treated seperately.
586     {
587       intersection = ((rho2 - k2)/k1 - p.z())/v.z();
588       if(calcNorm)
589       {
590         G4ThreeVector intersectionP = p + v * intersection;
591         *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
592         *n = n->unit();
593 
594         *validNorm = true;
595       }
596       return intersection;
597     }
598     else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
599     {
600       // intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
601       // The above calculation has a precision problem:
602       // known problem of solving quadratic equation with small A  
603 
604       A = A/vRho2;
605       B = (k1 * p.z() + k2 - rho2)/vRho2;
606       intersection = B/(-A + std::sqrt(B + sqr(A)));
607       if(calcNorm)
608       {
609         G4ThreeVector intersectionP = p + v * intersection;
610         *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
611         *n = n->unit();
612         *validNorm = true;
613       }
614       return intersection;
615     }
616     std::ostringstream message;
617     message << "There is no intersection between given line and solid!"
618             << G4endl
619             << "          p = " << p << G4endl
620             << "          v = " << v;
621     G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
622                 JustWarning, message);
623 
624     return kInfinity;
625   }
626   else if ( (rho2 < paraRho2 + kCarTolerance
627          || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 )
628          && std::fabs(p.z()) < dz + tolh) 
629   {
630     // If this is true we're somewhere in the border.
631     
632     G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
633 
634     if(std::fabs(p.z()) > dz - tolh)
635     {
636       // We're in the lower or upper edge
637       //
638       if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
639       {             // If we're heading out of the object that is treated here
640         if(calcNorm)
641         {
642           *validNorm = true;
643           if(p.z() > 0)
644             { *n = G4ThreeVector(0, 0, 1); }
645           else
646             { *n = G4ThreeVector(0, 0, -1); }
647         }
648         return 0;
649       }
650 
651       if(v.z() == 0)
652       {
653         // Case where we're moving inside the surface needs to be
654         // treated separately.
655         // Distance until it goes out through a side is returned.
656 
657         G4double r = (p.z() > 0)? r2 : r1;
658         G4double pDotV = p.dot(v);
659         A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y()));
660         intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2;
661 
662         if(calcNorm)
663         {
664           *validNorm = true;
665 
666           *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z()))
667               + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y()
668               * intersection, -k1/2).unit()).unit();
669         }
670         return intersection;
671       }
672     }
673     //
674     // Problem in the Logic :: Following condition for point on upper surface 
675     //                         and Vz<0  will return 0 (Problem #1015), but
676     //                         it has to return intersection with parabolic
677     //                         surface or with lower plane surface (z = -dz)
678     // The logic has to be  :: If not found intersection until now,
679     // do not exit but continue to search for possible intersection.
680     // Only for point situated on both borders (Z and parabolic)
681     // this condition has to be taken into account and done later
682     //
683     //
684     // else if(normal.dot(v) >= 0)
685     // {
686     //   if(calcNorm)
687     //   {
688     //     *validNorm = true;
689     //     *n = normal.unit();
690     //   }
691     //   return 0;
692     // }
693 
694     if(v.z() > 0)
695     {
696       // Check for collision with upper edge.
697 
698       intersection = (dz - p.z()) / v.z();
699       G4ThreeVector ip = p + intersection * v;
700 
701       if(ip.perp2() < sqr(r2 - tolh))
702       {
703         if(calcNorm)
704         {
705           *validNorm = true;
706           *n = G4ThreeVector(0, 0, 1);
707         }
708         return intersection;
709       }
710       else if(ip.perp2() < sqr(r2 + tolh))
711       {
712         if(calcNorm)
713         {
714           *validNorm = true;
715           *n = G4ThreeVector(0, 0, 1)
716              + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
717           *n = n->unit();
718         }
719         return intersection;
720       }
721     }
722     if( v.z() < 0)
723     {
724       // Check for collision with lower edge.
725 
726       intersection = (-dz - p.z()) / v.z();
727       G4ThreeVector ip = p + intersection * v;
728 
729       if(ip.perp2() < sqr(r1 - tolh))
730       {
731         if(calcNorm)
732         {
733           *validNorm = true;
734           *n = G4ThreeVector(0, 0, -1);
735         }
736         return intersection;
737       }
738       else if(ip.perp2() < sqr(r1 + tolh))
739       {
740         if(calcNorm)
741         {
742           *validNorm = true;
743           *n = G4ThreeVector(0, 0, -1)
744              + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
745           *n = n->unit();
746         }
747         return intersection;
748       }
749     }
750 
751     // Note: comparison with zero below would not be correct !
752     //
753     if(std::fabs(vRho2) > tol2) // precision error in the calculation of
754     {                           // intersection = (A+std::sqrt(B+sqr(A)))/vRho2
755       A = A/vRho2;
756       B = (k1 * p.z() + k2 - rho2);
757       if(std::fabs(B)>kCarTolerance)
758       {
759         B = (B)/vRho2;
760         intersection = B/(-A + std::sqrt(B + sqr(A)));
761       }
762       else                      // Point is On both borders: Z and parabolic
763       {                         // solution depends on normal.dot(v) sign
764         if(normal.dot(v) >= 0)
765         {
766           if(calcNorm)
767           {
768             *validNorm = true;
769             *n = normal.unit();
770           }
771           return 0;
772         }
773         intersection = 2.*A;
774       }
775     }
776     else
777     {
778       intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
779     }
780 
781     if(calcNorm)
782     {
783       *validNorm = true;
784       *n = G4ThreeVector(p.x() + intersection * v.x(), p.y()
785          + intersection * v.y(), - k1 / 2);
786       *n = n->unit();
787     }
788     return intersection;
789   }
790   else
791   {
792 #ifdef G4SPECSDEBUG
793     if(kOutside == Inside(p))
794     {
795       G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
796                   JustWarning, "Point p is outside!");
797     }
798     else
799       G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
800                   JustWarning, "There's an error in this functions code.");
801 #endif
802     return kInfinity;
803   }
804   return 0;
805 } 
806 
807 ///////////////////////////////////////////////////////////////////////////////
808 //
809 // Calculate distance (<=actual) to closest surface of shape from inside
810 //
811 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p) const
812 {
813   G4double safe=0.0,rho,safeR,safeZ ;
814   G4double tanRMax,secRMax,pRMax ;
815 
816 #ifdef G4SPECSDEBUG
817   if( Inside(p) == kOutside )
818   {
819      G4cout << G4endl ;
820      DumpInfo();
821      std::ostringstream message;
822      G4long oldprc = message.precision(16);
823      message << "Point p is outside !?" << G4endl
824              << "Position:" << G4endl
825              << "   p.x() = "   << p.x()/mm << " mm" << G4endl
826              << "   p.y() = "   << p.y()/mm << " mm" << G4endl
827              << "   p.z() = "   << p.z()/mm << " mm";
828      message.precision(oldprc) ;
829      G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002",
830                  JustWarning, message);
831   }
832 #endif
833 
834   rho = p.perp();
835   safeZ = dz - std::fabs(p.z()) ;
836 
837   tanRMax = (r2 - r1)*0.5/dz ;
838   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
839   pRMax   = tanRMax*p.z() + (r1+r2)*0.5 ;
840   safeR  = (pRMax - rho)/secRMax ;
841 
842   if (safeZ < safeR) { safe = safeZ; }
843   else { safe = safeR; }
844   if ( safe < 0.5 * kCarTolerance ) { safe = 0; }
845   return safe ;
846 }
847 
848 //////////////////////////////////////////////////////////////////////////
849 //
850 // G4EntityType
851 //
852 G4GeometryType G4Paraboloid::GetEntityType() const
853 {
854   return {"G4Paraboloid"};
855 }
856 
857 //////////////////////////////////////////////////////////////////////////
858 //
859 // Make a clone of the object
860 //
861 G4VSolid* G4Paraboloid::Clone() const
862 {
863   return new G4Paraboloid(*this);
864 }
865 
866 //////////////////////////////////////////////////////////////////////////
867 //
868 // Stream object contents to an output stream
869 //
870 std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const
871 {
872   G4long oldprc = os.precision(16);
873   os << "-----------------------------------------------------------\n"
874      << "    *** Dump for solid - " << GetName() << " ***\n"
875      << "    ===================================================\n"
876      << " Solid type: G4Paraboloid\n"
877      << " Parameters: \n"
878      << "    z half-axis:   " << dz/mm << " mm \n"
879      << "    radius at -dz: " << r1/mm << " mm \n"
880      << "    radius at dz:  " << r2/mm << " mm \n"
881      << "-----------------------------------------------------------\n";
882   os.precision(oldprc);
883 
884   return os;
885 }
886 
887 ////////////////////////////////////////////////////////////////////
888 //
889 // GetPointOnSurface
890 //
891 G4ThreeVector G4Paraboloid::GetPointOnSurface() const
892 {
893   G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea;
894   G4double z = G4RandFlat::shoot(0.,1.);
895   G4double phi = G4RandFlat::shoot(0., twopi);
896   if(pi*(sqr(r1) + sqr(r2))/A >= z)
897   {
898     G4double rho;
899     if(pi * sqr(r1) / A > z)
900     {
901       rho = r1 * std::sqrt(G4RandFlat::shoot(0., 1.));
902       return { rho * std::cos(phi), rho * std::sin(phi), -dz };
903     }
904     else
905     {
906       rho = r2 * std::sqrt(G4RandFlat::shoot(0., 1));
907       return { rho * std::cos(phi), rho * std::sin(phi), dz };
908     }
909   }
910   else
911   {
912     z = G4RandFlat::shoot(0., 1.)*2*dz - dz;
913     return { std::sqrt(z*k1 + k2)*std::cos(phi),
914              std::sqrt(z*k1 + k2)*std::sin(phi), z};
915   }
916 }
917 
918 /////////////////////////////////////////////////////////////////////////////
919 //
920 // Methods for visualisation
921 //
922 void G4Paraboloid::DescribeYourselfTo (G4VGraphicsScene& scene) const
923 {
924   scene.AddSolid(*this);
925 }
926 
927 G4Polyhedron* G4Paraboloid::CreatePolyhedron () const
928 {
929   return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
930 }
931 
932 G4Polyhedron* G4Paraboloid::GetPolyhedron () const
933 {
934   if (fpPolyhedron == nullptr ||
935       fRebuildPolyhedron ||
936       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
937       fpPolyhedron->GetNumberOfRotationSteps())
938   {
939     G4AutoLock l(&polyhedronMutex);
940     delete fpPolyhedron;
941     fpPolyhedron = CreatePolyhedron();
942     fRebuildPolyhedron = false;
943     l.unlock();
944   }
945   return fpPolyhedron;
946 }
947 
948 #endif // !defined(G4GEOM_USE_UPARABOLOID) || !defined(G4GEOM_USE_SYS_USOLIDS)
949