ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Paraboloid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Paraboloid.cc
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 
59 //
60 // constructor - check parameters
61 //
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 
92 //
93 // Fake default constructor - sets only member data and allocates memory
94 // for usage restricted to object persistency.
95 //
97  : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
98 {
99 }
100 
102 //
103 // Destructor
104 //
106 {
107  delete fpPolyhedron; fpPolyhedron = nullptr;
108 }
109 
111 //
112 // Copy constructor
113 //
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 
122 //
123 // Assignment operator
124 //
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  //
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 
146 //
147 // Get bounding box
148 //
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 
171 //
172 // Calculate extent under transform and specified limit
173 //
174 G4bool
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 
191 //
192 // Return whether point inside/outside/on surface
193 //
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 
234 //
235 // SurfaceNormal
236 //
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 
326 //
327 // Calculate distance to shape from outside, along normalised vector
328 // - return kInfinity if no intersection
329 //
330 //
332  const G4ThreeVector& v ) const
333 {
334  G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
336  G4double tolh = 0.5*kCarTolerance;
337 
338  if(r2 && 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 && 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 
460 //
461 // Calculate distance (<= actual) to closest surface of shape from outside
462 // - Return zero if point inside
463 //
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 
497 //
498 // Calculate distance to surface of shape from 'inside'
499 //
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;
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 
808 //
809 // Calculate distance (<=actual) to closest surface of shape from inside
810 //
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  G4int 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 
849 //
850 // G4EntityType
851 //
853 {
854  return G4String("G4Paraboloid");
855 }
856 
858 //
859 // Make a clone of the object
860 //
862 {
863  return new G4Paraboloid(*this);
864 }
865 
867 //
868 // Stream object contents to an output stream
869 //
870 std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const
871 {
872  G4int 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 
888 //
889 // GetPointOnSurface
890 //
892 {
894  G4double z = G4RandFlat::shoot(0.,1.);
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 G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz);
903  }
904  else
905  {
906  rho = r2 * std::sqrt(G4RandFlat::shoot(0., 1));
907  return G4ThreeVector(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 G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi),
914  std::sqrt(z*k1 + k2)*std::sin(phi), z);
915  }
916 }
917 
919 //
920 // Methods for visualisation
921 //
923 {
924  scene.AddSolid(*this);
925 }
926 
928 {
929  return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
930 }
931 
933 {
934  if (fpPolyhedron == nullptr ||
938  {
939  G4AutoLock l(&polyhedronMutex);
940  delete fpPolyhedron;
942  fRebuildPolyhedron = false;
943  l.unlock();
944  }
945  return fpPolyhedron;
946 }
947 
948 #endif // !defined(G4GEOM_USE_UPARABOLOID) || !defined(G4GEOM_USE_SYS_USOLIDS)