Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Paraboloid Class Reference

#include <G4Paraboloid.hh>

Inheritance diagram for G4Paraboloid:
Collaboration diagram for G4Paraboloid:

Public Member Functions

 G4Paraboloid (const G4String &pName, G4double pDz, G4double pR1, G4double pR2)
 
virtual ~G4Paraboloid ()
 
G4double GetZHalfLength () const
 
G4double GetRadiusMinusZ () const
 
G4double GetRadiusPlusZ () const
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
G4double CalculateSurfaceArea () const
 
void SetZHalfLength (G4double dz)
 
void SetRadiusMinusZ (G4double R1)
 
void SetRadiusPlusZ (G4double R2)
 
void Extent (G4ThreeVector &pMin, G4ThreeVector &pMax) const
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
G4GeometryType GetEntityType () const
 
G4VSolidClone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
G4ThreeVector GetPointOnSurface () const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4PolyhedronCreatePolyhedron () const
 
G4PolyhedronGetPolyhedron () const
 
 G4Paraboloid (__void__ &)
 
 G4Paraboloid (const G4Paraboloid &rhs)
 
G4Paraboloidoperator= (const G4Paraboloid &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
G4bool operator== (const G4VSolid &s) const
 
G4String GetName () const
 
void SetName (const G4String &name)
 
G4double GetTolerance () const
 
virtual void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
void DumpInfo () const
 
virtual G4VisExtent GetExtent () const
 
virtual const G4VSolidGetConstituentSolid (G4int no) const
 
virtual G4VSolidGetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolidGetDisplacedSolidPtr () const
 
virtual G4DisplacedSolidGetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolidoperator= (const G4VSolid &rhs)
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Protected Attributes

G4bool fRebuildPolyhedron
 
G4PolyhedronfpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Additional Inherited Members

- Protected Member Functions inherited from G4VSolid
void CalculateClippedPolygonExtent (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipCrossSection (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipBetweenSections (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipPolygon (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis) const
 

Detailed Description

Definition at line 75 of file G4Paraboloid.hh.

Constructor & Destructor Documentation

G4Paraboloid::G4Paraboloid ( const G4String pName,
G4double  pDz,
G4double  pR1,
G4double  pR2 
)

Definition at line 66 of file G4Paraboloid.cc.

70  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
71  fSurfaceArea(0.), fCubicVolume(0.)
72 
73 {
74  if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) )
75  {
76  std::ostringstream message;
77  message << "Invalid dimensions. Negative Input Values or R1>=R2 - "
78  << GetName();
79  G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002",
80  FatalErrorInArgument, message,
81  "Z half-length must be larger than zero or R1>=R2.");
82  }
83 
84  r1 = pR1;
85  r2 = pR2;
86  dz = pDz;
87 
88  // r1^2 = k1 * (-dz) + k2
89  // r2^2 = k1 * ( dz) + k2
90  // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
91  // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
92 
93  k1 = (r2 * r2 - r1 * r1) / 2 / dz;
94  k2 = (r2 * r2 + r1 * r1) / 2;
95 }
G4String GetName() const
G4Polyhedron * fpPolyhedron
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4bool fRebuildPolyhedron

Here is the call graph for this function:

Here is the caller graph for this function:

G4Paraboloid::~G4Paraboloid ( )
virtual

Definition at line 113 of file G4Paraboloid.cc.

114 {
115  delete fpPolyhedron; fpPolyhedron = 0;
116 }
G4Polyhedron * fpPolyhedron
G4Paraboloid::G4Paraboloid ( __void__ &  a)

Definition at line 102 of file G4Paraboloid.cc.

103  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
104  fSurfaceArea(0.), fCubicVolume(0.),
105  dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
106 {
107 }
G4Polyhedron * fpPolyhedron
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4bool fRebuildPolyhedron
G4Paraboloid::G4Paraboloid ( const G4Paraboloid rhs)

Definition at line 122 of file G4Paraboloid.cc.

123  : G4VSolid(rhs), fRebuildPolyhedron(false), fpPolyhedron(0),
124  fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume),
125  dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2)
126 {
127 }
G4Polyhedron * fpPolyhedron
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4bool fRebuildPolyhedron

Member Function Documentation

G4bool G4Paraboloid::CalculateExtent ( const EAxis  pAxis,
const G4VoxelLimits pVoxelLimit,
const G4AffineTransform pTransform,
G4double pmin,
G4double pmax 
) const
virtual

Implements G4VSolid.

Definition at line 194 of file G4Paraboloid.cc.

198 {
199  G4ThreeVector bmin, bmax;
200 
201  // Get bounding box
202  Extent(bmin,bmax);
203 
204  // Find extent
205  G4BoundingEnvelope bbox(bmin,bmax);
206  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
207 }
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const

Here is the call graph for this function:

G4double G4Paraboloid::CalculateSurfaceArea ( ) const
inline

Here is the caller graph for this function:

G4VSolid * G4Paraboloid::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 879 of file G4Paraboloid.cc.

880 {
881  return new G4Paraboloid(*this);
882 }
G4Paraboloid(const G4String &pName, G4double pDz, G4double pR1, G4double pR2)
Definition: G4Paraboloid.cc:66

Here is the call graph for this function:

G4Polyhedron * G4Paraboloid::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 945 of file G4Paraboloid.cc.

946 {
947  return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
948 }
static constexpr double twopi
Definition: G4SIunits.hh:76

Here is the caller graph for this function:

void G4Paraboloid::DescribeYourselfTo ( G4VGraphicsScene scene) const
virtual

Implements G4VSolid.

Definition at line 940 of file G4Paraboloid.cc.

941 {
942  scene.AddSolid(*this);
943 }
virtual void AddSolid(const G4Box &)=0

Here is the call graph for this function:

G4double G4Paraboloid::DistanceToIn ( const G4ThreeVector p,
const G4ThreeVector v 
) const
virtual

Implements G4VSolid.

Definition at line 349 of file G4Paraboloid.cc.

351 {
352  G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
354  G4double tolh = 0.5*kCarTolerance;
355 
356  if(r2 && p.z() > - tolh + dz)
357  {
358  // If the points is above check for intersection with upper edge.
359 
360  if(v.z() < 0)
361  {
362  G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz.
363  if(sqr(p.x() + v.x()*intersection)
364  + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance))
365  {
366  if(p.z() < tolh + dz)
367  { return 0; }
368  else
369  { return intersection; }
370  }
371  }
372  else // Direction away, no possibility of intersection
373  {
374  return kInfinity;
375  }
376  }
377  else if(r1 && p.z() < tolh - dz)
378  {
379  // If the points is belove check for intersection with lower edge.
380 
381  if(v.z() > 0)
382  {
383  G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz.
384  if(sqr(p.x() + v.x()*intersection)
385  + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance))
386  {
387  if(p.z() > -tolh - dz)
388  {
389  return 0;
390  }
391  else
392  {
393  return intersection;
394  }
395  }
396  }
397  else // Direction away, no possibility of intersection
398  {
399  return kInfinity;
400  }
401  }
402 
403  G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(),
404  vRho2 = v.perp2(), intersection,
405  B = (k1 * p.z() + k2 - rho2) * vRho2;
406 
407  if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) )
408  || (p.z() < - dz+kCarTolerance)
409  || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside.
410  {
411  // Is there a problem with squaring rho twice?
412 
413  if(vRho2<tol2) // Needs to be treated seperately.
414  {
415  intersection = ((rho2 - k2)/k1 - p.z())/v.z();
416  if(intersection < 0) { return kInfinity; }
417  else if(std::fabs(p.z() + v.z() * intersection) <= dz)
418  {
419  return intersection;
420  }
421  else
422  {
423  return kInfinity;
424  }
425  }
426  else if(A*A + B < 0) // No real intersections.
427  {
428  return kInfinity;
429  }
430  else
431  {
432  intersection = (A - std::sqrt(B + sqr(A))) / vRho2;
433  if(intersection < 0)
434  {
435  return kInfinity;
436  }
437  else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh)
438  {
439  return intersection;
440  }
441  else
442  {
443  return kInfinity;
444  }
445  }
446  }
447  else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2)
448  {
449  // If this is true we're somewhere in the border.
450 
451  G4ThreeVector normal(p.x(), p.y(), -k1/2);
452  if(normal.dot(v) <= 0)
453  { return 0; }
454  else
455  { return kInfinity; }
456  }
457  else
458  {
459  std::ostringstream message;
460  if(Inside(p) == kInside)
461  {
462  message << "Point p is inside! - " << GetName() << G4endl;
463  }
464  else
465  {
466  message << "Likely a problem in this function, for solid: " << GetName()
467  << G4endl;
468  }
469  message << " p = " << p * (1/mm) << " mm" << G4endl
470  << " v = " << v * (1/mm) << " mm";
471  G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002",
472  JustWarning, message);
473  return 0;
474  }
475 }
G4String GetName() const
static constexpr double mm
Definition: G4SIunits.hh:115
double perp2() const
static const G4double kInfinity
Definition: geomdefs.hh:42
double x() const
double B(double temperature)
double z() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
double A(double temperature)
EInside Inside(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4Paraboloid::DistanceToIn ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 482 of file G4Paraboloid.cc.

483 {
484  G4double safz = -dz+std::fabs(p.z());
485  if(safz<0) { safz=0; }
486  G4double safr = kInfinity;
487 
488  G4double rho = p.x()*p.x()+p.y()*p.y();
489  G4double paraRho = (p.z()-k2)/k1;
490  G4double sqrho = std::sqrt(rho);
491 
492  if(paraRho<0)
493  {
494  safr=sqrho-r2;
495  if(safr>safz) { safz=safr; }
496  return safz;
497  }
498 
499  G4double sqprho = std::sqrt(paraRho);
500  G4double dRho = sqrho-sqprho;
501  if(dRho<0) { return safz; }
502 
503  G4double talf = -2.*k1*sqprho;
504  G4double tmp = 1+talf*talf;
505  if(tmp<0.) { return safz; }
506 
507  G4double salf = talf/std::sqrt(tmp);
508  safr = std::fabs(dRho*salf);
509  if(safr>safz) { safz=safr; }
510 
511  return safz;
512 }
static const G4double kInfinity
Definition: geomdefs.hh:42
double x() const
double z() const
double y() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4Paraboloid::DistanceToOut ( const G4ThreeVector p,
const G4ThreeVector v,
const G4bool  calcNorm = G4bool(false),
G4bool validNorm = 0,
G4ThreeVector n = 0 
) const
virtual

Implements G4VSolid.

Definition at line 518 of file G4Paraboloid.cc.

523 {
524  G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
525  G4double vRho2 = v.perp2(), intersection;
527  G4double tolh = 0.5*kCarTolerance;
528 
529  if(calcNorm) { *validNorm = false; }
530 
531  // We have that the particle p follows the line x = p + s * v
532  // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and
533  // z = p.z() + s * v.z()
534  // The equation for all points on the surface (surface expanded for
535  // to include all z) x^2 + y^2 = k1 * z + k2 => .. =>
536  // => s = (A +- std::sqrt(A^2 + B)) / vRho2
537  // where:
538  //
539  G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
540  //
541  // and:
542  //
543  G4double B = (-rho2 + paraRho2) * vRho2;
544 
545  if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2
546  && std::fabs(p.z()) < dz - kCarTolerance)
547  {
548  // Make sure it's safely inside.
549 
550  if(v.z() > 0)
551  {
552  // It's heading upwards, check where it colides with the plane z = dz.
553  // When it does, is that in the surface of the paraboloid.
554  // z = p.z() + variable * v.z() for all points the particle can go.
555  // => variable = (z - p.z()) / v.z() so intersection must be:
556 
557  intersection = (dz - p.z()) / v.z();
558  G4ThreeVector ip = p + intersection * v; // Point of intersection.
559 
560  if(ip.perp2() < sqr(r2 + kCarTolerance))
561  {
562  if(calcNorm)
563  {
564  *n = G4ThreeVector(0, 0, 1);
565  if(r2 < tolh || ip.perp2() > sqr(r2 - tolh))
566  {
567  *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
568  *n = n->unit();
569  }
570  *validNorm = true;
571  }
572  return intersection;
573  }
574  }
575  else if(v.z() < 0)
576  {
577  // It's heading downwards, check were it colides with the plane z = -dz.
578  // When it does, is that in the surface of the paraboloid.
579  // z = p.z() + variable * v.z() for all points the particle can go.
580  // => variable = (z - p.z()) / v.z() so intersection must be:
581 
582  intersection = (-dz - p.z()) / v.z();
583  G4ThreeVector ip = p + intersection * v; // Point of intersection.
584 
585  if(ip.perp2() < sqr(r1 + tolh))
586  {
587  if(calcNorm)
588  {
589  *n = G4ThreeVector(0, 0, -1);
590  if(r1 < tolh || ip.perp2() > sqr(r1 - tolh))
591  {
592  *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
593  *n = n->unit();
594  }
595  *validNorm = true;
596  }
597  return intersection;
598  }
599  }
600 
601  // Now check for collisions with paraboloid surface.
602 
603  if(vRho2 == 0) // Needs to be treated seperately.
604  {
605  intersection = ((rho2 - k2)/k1 - p.z())/v.z();
606  if(calcNorm)
607  {
608  G4ThreeVector intersectionP = p + v * intersection;
609  *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
610  *n = n->unit();
611 
612  *validNorm = true;
613  }
614  return intersection;
615  }
616  else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
617  {
618  // intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
619  // The above calculation has a precision problem:
620  // known problem of solving quadratic equation with small A
621 
622  A = A/vRho2;
623  B = (k1 * p.z() + k2 - rho2)/vRho2;
624  intersection = B/(-A + std::sqrt(B + sqr(A)));
625  if(calcNorm)
626  {
627  G4ThreeVector intersectionP = p + v * intersection;
628  *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
629  *n = n->unit();
630  *validNorm = true;
631  }
632  return intersection;
633  }
634  std::ostringstream message;
635  message << "There is no intersection between given line and solid!"
636  << G4endl
637  << " p = " << p << G4endl
638  << " v = " << v;
639  G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
640  JustWarning, message);
641 
642  return kInfinity;
643  }
644  else if ( (rho2 < paraRho2 + kCarTolerance
645  || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 )
646  && std::fabs(p.z()) < dz + tolh)
647  {
648  // If this is true we're somewhere in the border.
649 
650  G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
651 
652  if(std::fabs(p.z()) > dz - tolh)
653  {
654  // We're in the lower or upper edge
655  //
656  if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
657  { // If we're heading out of the object that is treated here
658  if(calcNorm)
659  {
660  *validNorm = true;
661  if(p.z() > 0)
662  { *n = G4ThreeVector(0, 0, 1); }
663  else
664  { *n = G4ThreeVector(0, 0, -1); }
665  }
666  return 0;
667  }
668 
669  if(v.z() == 0)
670  {
671  // Case where we're moving inside the surface needs to be
672  // treated separately.
673  // Distance until it goes out through a side is returned.
674 
675  G4double r = (p.z() > 0)? r2 : r1;
676  G4double pDotV = p.dot(v);
677  A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y()));
678  intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2;
679 
680  if(calcNorm)
681  {
682  *validNorm = true;
683 
684  *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z()))
685  + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y()
686  * intersection, -k1/2).unit()).unit();
687  }
688  return intersection;
689  }
690  }
691  //
692  // Problem in the Logic :: Following condition for point on upper surface
693  // and Vz<0 will return 0 (Problem #1015), but
694  // it has to return intersection with parabolic
695  // surface or with lower plane surface (z = -dz)
696  // The logic has to be :: If not found intersection until now,
697  // do not exit but continue to search for possible intersection.
698  // Only for point situated on both borders (Z and parabolic)
699  // this condition has to be taken into account and done later
700  //
701  //
702  // else if(normal.dot(v) >= 0)
703  // {
704  // if(calcNorm)
705  // {
706  // *validNorm = true;
707  // *n = normal.unit();
708  // }
709  // return 0;
710  // }
711 
712  if(v.z() > 0)
713  {
714  // Check for collision with upper edge.
715 
716  intersection = (dz - p.z()) / v.z();
717  G4ThreeVector ip = p + intersection * v;
718 
719  if(ip.perp2() < sqr(r2 - tolh))
720  {
721  if(calcNorm)
722  {
723  *validNorm = true;
724  *n = G4ThreeVector(0, 0, 1);
725  }
726  return intersection;
727  }
728  else if(ip.perp2() < sqr(r2 + tolh))
729  {
730  if(calcNorm)
731  {
732  *validNorm = true;
733  *n = G4ThreeVector(0, 0, 1)
734  + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
735  *n = n->unit();
736  }
737  return intersection;
738  }
739  }
740  if( v.z() < 0)
741  {
742  // Check for collision with lower edge.
743 
744  intersection = (-dz - p.z()) / v.z();
745  G4ThreeVector ip = p + intersection * v;
746 
747  if(ip.perp2() < sqr(r1 - tolh))
748  {
749  if(calcNorm)
750  {
751  *validNorm = true;
752  *n = G4ThreeVector(0, 0, -1);
753  }
754  return intersection;
755  }
756  else if(ip.perp2() < sqr(r1 + tolh))
757  {
758  if(calcNorm)
759  {
760  *validNorm = true;
761  *n = G4ThreeVector(0, 0, -1)
762  + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
763  *n = n->unit();
764  }
765  return intersection;
766  }
767  }
768 
769  // Note: comparison with zero below would not be correct !
770  //
771  if(std::fabs(vRho2) > tol2) // precision error in the calculation of
772  { // intersection = (A+std::sqrt(B+sqr(A)))/vRho2
773  A = A/vRho2;
774  B = (k1 * p.z() + k2 - rho2);
775  if(std::fabs(B)>kCarTolerance)
776  {
777  B = (B)/vRho2;
778  intersection = B/(-A + std::sqrt(B + sqr(A)));
779  }
780  else // Point is On both borders: Z and parabolic
781  { // solution depends on normal.dot(v) sign
782  if(normal.dot(v) >= 0)
783  {
784  if(calcNorm)
785  {
786  *validNorm = true;
787  *n = normal.unit();
788  }
789  return 0;
790  }
791  intersection = 2.*A;
792  }
793  }
794  else
795  {
796  intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
797  }
798 
799  if(calcNorm)
800  {
801  *validNorm = true;
802  *n = G4ThreeVector(p.x() + intersection * v.x(), p.y()
803  + intersection * v.y(), - k1 / 2);
804  *n = n->unit();
805  }
806  return intersection;
807  }
808  else
809  {
810 #ifdef G4SPECSDEBUG
811  if(kOutside == Inside(p))
812  {
813  G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
814  JustWarning, "Point p is outside!");
815  }
816  else
817  G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
818  JustWarning, "There's an error in this functions code.");
819 #endif
820  return kInfinity;
821  }
822  return 0;
823 }
double perp2() const
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
double dot(const Hep3Vector &) const
double B(double temperature)
double z() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
double A(double temperature)
tuple v
Definition: test.py:18
EInside Inside(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
Hep3Vector unit() const
double y() const
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4Paraboloid::DistanceToOut ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 829 of file G4Paraboloid.cc.

830 {
831  G4double safe=0.0,rho,safeR,safeZ ;
832  G4double tanRMax,secRMax,pRMax ;
833 
834 #ifdef G4SPECSDEBUG
835  if( Inside(p) == kOutside )
836  {
837  G4cout << G4endl ;
838  DumpInfo();
839  std::ostringstream message;
840  G4int oldprc = message.precision(16);
841  message << "Point p is outside !?" << G4endl
842  << "Position:" << G4endl
843  << " p.x() = " << p.x()/mm << " mm" << G4endl
844  << " p.y() = " << p.y()/mm << " mm" << G4endl
845  << " p.z() = " << p.z()/mm << " mm";
846  message.precision(oldprc) ;
847  G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002",
848  JustWarning, message);
849  }
850 #endif
851 
852  rho = p.perp();
853  safeZ = dz - std::fabs(p.z()) ;
854 
855  tanRMax = (r2 - r1)*0.5/dz ;
856  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
857  pRMax = tanRMax*p.z() + (r1+r2)*0.5 ;
858  safeR = (pRMax - rho)/secRMax ;
859 
860  if (safeZ < safeR) { safe = safeZ; }
861  else { safe = safeR; }
862  if ( safe < 0.5 * kCarTolerance ) { safe = 0; }
863  return safe ;
864 }
static constexpr double mm
Definition: G4SIunits.hh:115
double x() const
int G4int
Definition: G4Types.hh:78
double z() const
void DumpInfo() const
G4GLOB_DLL std::ostream G4cout
EInside Inside(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
double G4double
Definition: G4Types.hh:76
double perp() const

Here is the call graph for this function:

void G4Paraboloid::Extent ( G4ThreeVector pMin,
G4ThreeVector pMax 
) const
virtual

Reimplemented from G4VSolid.

Definition at line 170 of file G4Paraboloid.cc.

171 {
172  pMin.set(-r2,-r2,-dz);
173  pMax.set( r2, r2, dz);
174 
175  // Check correctness of the bounding box
176  //
177  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
178  {
179  std::ostringstream message;
180  message << "Bad bounding box (min >= max) for solid: "
181  << GetName() << " !"
182  << "\npMin = " << pMin
183  << "\npMax = " << pMax;
184  G4Exception("G4Paraboloid::Extent()", "GeomMgt0001", JustWarning, message);
185  DumpInfo();
186  }
187 }
void set(double x, double y, double z)
G4String GetName() const
double x() const
double z() const
void DumpInfo() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4Paraboloid::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

G4GeometryType G4Paraboloid::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 870 of file G4Paraboloid.cc.

871 {
872  return G4String("G4Paraboloid");
873 }
G4ThreeVector G4Paraboloid::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 909 of file G4Paraboloid.cc.

910 {
911  G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea;
912  G4double z = G4RandFlat::shoot(0.,1.);
913  G4double phi = G4RandFlat::shoot(0., twopi);
914  if(pi*(sqr(r1) + sqr(r2))/A >= z)
915  {
916  G4double rho;
917  if(pi * sqr(r1) / A > z)
918  {
919  rho = r1 * std::sqrt(G4RandFlat::shoot(0., 1.));
920  return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz);
921  }
922  else
923  {
924  rho = r2 * std::sqrt(G4RandFlat::shoot(0., 1));
925  return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz);
926  }
927  }
928  else
929  {
930  z = G4RandFlat::shoot(0., 1.)*2*dz - dz;
931  return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi),
932  std::sqrt(z*k1 + k2)*std::sin(phi), z);
933  }
934 }
ThreeVector shoot(const G4int Ap, const G4int Af)
CLHEP::Hep3Vector G4ThreeVector
G4double CalculateSurfaceArea() const
static constexpr double twopi
Definition: G4SIunits.hh:76
double A(double temperature)
tuple z
Definition: test.py:28
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4Polyhedron * G4Paraboloid::GetPolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 951 of file G4Paraboloid.cc.

952 {
953  if (!fpPolyhedron ||
957  {
958  G4AutoLock l(&polyhedronMutex);
959  delete fpPolyhedron;
961  fRebuildPolyhedron = false;
962  l.unlock();
963  }
964  return fpPolyhedron;
965 }
G4Polyhedron * fpPolyhedron
static G4int GetNumberOfRotationSteps()
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4Polyhedron * CreatePolyhedron() const
G4bool fRebuildPolyhedron

Here is the call graph for this function:

G4double G4Paraboloid::GetRadiusMinusZ ( ) const
inline

Here is the caller graph for this function:

G4double G4Paraboloid::GetRadiusPlusZ ( ) const
inline

Here is the caller graph for this function:

G4double G4Paraboloid::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

G4double G4Paraboloid::GetZHalfLength ( ) const
inline

Here is the caller graph for this function:

EInside G4Paraboloid::Inside ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 213 of file G4Paraboloid.cc.

214 {
215  // First check is the point is above or below the solid.
216  //
217  if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; }
218 
219  G4double rho2 = p.perp2(),
220  rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance),
221  A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
222 
223  if(A < 0 && sqr(A) > rhoSurfTimesTol2)
224  {
225  // Actually checking rho < radius of paraboloid at z = p.z().
226  // We're either inside or in lower/upper cutoff area.
227 
228  if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
229  {
230  // We're in the upper/lower cutoff area, sides have a paraboloid shape
231  // maybe further checks should be made to make these nicer
232 
233  return kSurface;
234  }
235  else
236  {
237  return kInside;
238  }
239  }
240  else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
241  {
242  // We're in the parabolic surface.
243 
244  return kSurface;
245  }
246  else
247  {
248  return kOutside;
249  }
250 }
double perp2() const
double z() const
double A(double temperature)
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4Paraboloid & G4Paraboloid::operator= ( const G4Paraboloid rhs)

Definition at line 134 of file G4Paraboloid.cc.

135 {
136  // Check assignment to self
137  //
138  if (this == &rhs) { return *this; }
139 
140  // Copy base class data
141  //
142  G4VSolid::operator=(rhs);
143 
144  // Copy data
145  //
146  fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
147  dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2;
148  fRebuildPolyhedron = false;
149  delete fpPolyhedron; fpPolyhedron = 0;
150 
151  return *this;
152 }
G4Polyhedron * fpPolyhedron
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111
G4bool fRebuildPolyhedron

Here is the call graph for this function:

void G4Paraboloid::SetRadiusMinusZ ( G4double  R1)
inline
void G4Paraboloid::SetRadiusPlusZ ( G4double  R2)
inline
void G4Paraboloid::SetZHalfLength ( G4double  dz)
inline
std::ostream & G4Paraboloid::StreamInfo ( std::ostream &  os) const
virtual

Implements G4VSolid.

Definition at line 888 of file G4Paraboloid.cc.

889 {
890  G4int oldprc = os.precision(16);
891  os << "-----------------------------------------------------------\n"
892  << " *** Dump for solid - " << GetName() << " ***\n"
893  << " ===================================================\n"
894  << " Solid type: G4Paraboloid\n"
895  << " Parameters: \n"
896  << " z half-axis: " << dz/mm << " mm \n"
897  << " radius at -dz: " << r1/mm << " mm \n"
898  << " radius at dz: " << r2/mm << " mm \n"
899  << "-----------------------------------------------------------\n";
900  os.precision(oldprc);
901 
902  return os;
903 }
G4String GetName() const
static constexpr double mm
Definition: G4SIunits.hh:115
int G4int
Definition: G4Types.hh:78

Here is the call graph for this function:

G4ThreeVector G4Paraboloid::SurfaceNormal ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 255 of file G4Paraboloid.cc.

256 {
257  G4ThreeVector n(0, 0, 0);
258  if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance)
259  {
260  // If above or below just return normal vector for the cutoff plane.
261 
262  n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z()));
263  }
264  else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
265  {
266  // This means we're somewhere in the plane z = dz or z = -dz.
267  // (As far as the program is concerned anyway.
268 
269  if(p.z() < 0) // Are we in upper or lower plane?
270  {
271  if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance))
272  {
273  n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit();
274  }
275  else if(r1 < 0.5 * kCarTolerance
276  || p.perp2() > sqr(r1 - 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  else
288  {
289  if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance))
290  {
291  n = G4ThreeVector(p.x(), p.y(), 0.).unit();
292  }
293  else if(r2 < 0.5 * kCarTolerance
294  || p.perp2() > sqr(r2 - 0.5 * kCarTolerance))
295  {
296  n = G4ThreeVector(p.x(), p.y(), 0.).unit()
297  + G4ThreeVector(0., 0., 1.).unit();
298  n = n.unit();
299  }
300  else
301  {
302  n = G4ThreeVector(0., 0., 1.);
303  }
304  }
305  }
306  else
307  {
308  G4double rho2 = p.perp2();
309  G4double rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance);
310  G4double A = rho2 - ((k1 *p.z() + k2)
311  + 0.25 * kCarTolerance * kCarTolerance);
312 
313  if(A < 0 && sqr(A) > rhoSurfTimesTol2)
314  {
315  // Actually checking rho < radius of paraboloid at z = p.z().
316  // We're inside.
317 
318  if(p.mag2() != 0) { n = p.unit(); }
319  }
320  else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
321  {
322  // We're in the parabolic surface.
323 
324  n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
325  }
326  else
327  {
328  n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
329  }
330  }
331 
332  if(n.mag2() == 0)
333  {
334  std::ostringstream message;
335  message << "No normal defined for this point p." << G4endl
336  << " p = " << 1 / mm * p << " mm";
337  G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002",
338  JustWarning, message);
339  }
340  return n;
341 }
static constexpr double mm
Definition: G4SIunits.hh:115
double perp2() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
double z() const
double A(double temperature)
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
Hep3Vector unit() const
double y() const
double mag2() const
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Member Data Documentation

G4Polyhedron* G4Paraboloid::fpPolyhedron
mutableprotected

Definition at line 152 of file G4Paraboloid.hh.

G4bool G4Paraboloid::fRebuildPolyhedron
mutableprotected

Definition at line 151 of file G4Paraboloid.hh.


The documentation for this class was generated from the following files: