G4Paraboloid Class Reference

#include <G4Paraboloid.hh>

Inheritance diagram for G4Paraboloid:

G4VSolid

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)
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
G4NURBSCreateNURBS () const
 G4Paraboloid (__void__ &)
 G4Paraboloid (const G4Paraboloid &rhs)
G4Paraboloidoperator= (const G4Paraboloid &rhs)

Protected Attributes

G4PolyhedronfpPolyhedron

Detailed Description

Definition at line 64 of file G4Paraboloid.hh.


Constructor & Destructor Documentation

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

Definition at line 59 of file G4Paraboloid.cc.

References FatalErrorInArgument, G4Exception(), and G4VSolid::GetName().

Referenced by Clone().

00063  : G4VSolid(pName), fpPolyhedron(0), fSurfaceArea(0.), fCubicVolume(0.) 
00064 
00065 {
00066   if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) )
00067   {
00068     std::ostringstream message;
00069     message << "Invalid dimensions. Negative Input Values or R1>=R2 - "
00070             << GetName();
00071     G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002", 
00072                 FatalErrorInArgument, message,
00073                 "Z half-length must be larger than zero or R1>=R2.");
00074   }
00075 
00076   r1 = pR1;
00077   r2 = pR2;
00078   dz = pDz;
00079 
00080   // r1^2 = k1 * (-dz) + k2
00081   // r2^2 = k1 * ( dz) + k2
00082   // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
00083   // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
00084 
00085   k1 = (r2 * r2 - r1 * r1) / 2 / dz;
00086   k2 = (r2 * r2 + r1 * r1) / 2;
00087 }

G4Paraboloid::~G4Paraboloid (  )  [virtual]

Definition at line 104 of file G4Paraboloid.cc.

00105 {
00106 }

G4Paraboloid::G4Paraboloid ( __void__ &   ) 

Definition at line 94 of file G4Paraboloid.cc.

00095   : G4VSolid(a), fpPolyhedron(0), fSurfaceArea(0.), fCubicVolume(0.),
00096     dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
00097 {
00098 }

G4Paraboloid::G4Paraboloid ( const G4Paraboloid rhs  ) 

Definition at line 112 of file G4Paraboloid.cc.

00113   : G4VSolid(rhs), fpPolyhedron(0),
00114     fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume),
00115     dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2)
00116 {
00117 }


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 161 of file G4Paraboloid.cc.

References G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), G4AffineTransform::IsRotated(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetRotation(), and G4AffineTransform::NetTranslation().

00165 {
00166   G4double xMin = -r2 + pTransform.NetTranslation().x(),
00167            xMax = r2 + pTransform.NetTranslation().x(),
00168            yMin = -r2 + pTransform.NetTranslation().y(),
00169            yMax = r2 + pTransform.NetTranslation().y(),
00170            zMin = -dz + pTransform.NetTranslation().z(),
00171            zMax = dz + pTransform.NetTranslation().z();
00172 
00173   if(!pTransform.IsRotated()
00174   || pTransform.NetRotation()(G4ThreeVector(0, 0, 1)) == G4ThreeVector(0, 0, 1))
00175   {
00176     if(pVoxelLimit.IsXLimited())
00177     {
00178       if(pVoxelLimit.GetMaxXExtent() < xMin - 0.5 * kCarTolerance
00179       || pVoxelLimit.GetMinXExtent() > xMax + 0.5 * kCarTolerance)
00180       {
00181         return false;
00182       }
00183       else
00184       {
00185         if(pVoxelLimit.GetMinXExtent() > xMin)
00186         {
00187           xMin = pVoxelLimit.GetMinXExtent();
00188         }
00189         if(pVoxelLimit.GetMaxXExtent() < xMax)
00190         {
00191           xMax = pVoxelLimit.GetMaxXExtent();
00192         }
00193       }
00194     }
00195     if(pVoxelLimit.IsYLimited())
00196     {
00197       if(pVoxelLimit.GetMaxYExtent() < yMin - 0.5 * kCarTolerance
00198       || pVoxelLimit.GetMinYExtent() > yMax + 0.5 * kCarTolerance)
00199       {
00200         return false;
00201       }
00202       else
00203       {
00204         if(pVoxelLimit.GetMinYExtent() > yMin)
00205         {
00206           yMin = pVoxelLimit.GetMinYExtent();
00207         }
00208         if(pVoxelLimit.GetMaxYExtent() < yMax)
00209         {
00210           yMax = pVoxelLimit.GetMaxYExtent();
00211         }
00212       }
00213     }
00214     if(pVoxelLimit.IsZLimited())
00215     {
00216       if(pVoxelLimit.GetMaxZExtent() < zMin - 0.5 * kCarTolerance
00217       || pVoxelLimit.GetMinZExtent() > zMax + 0.5 * kCarTolerance)
00218       {
00219         return false;
00220       }
00221       else
00222       {
00223         if(pVoxelLimit.GetMinZExtent() > zMin)
00224         {
00225           zMin = pVoxelLimit.GetMinZExtent();
00226         }
00227         if(pVoxelLimit.GetMaxZExtent() < zMax)
00228         {
00229           zMax = pVoxelLimit.GetMaxZExtent();
00230         }
00231       }
00232     }
00233     switch(pAxis)
00234     {
00235       case kXAxis:
00236         pMin = xMin;
00237         pMax = xMax;
00238         break;
00239       case kYAxis:
00240         pMin = yMin;
00241         pMax = yMax;
00242         break;
00243       case kZAxis:
00244         pMin = zMin;
00245         pMax = zMax;
00246         break;
00247       default:
00248         pMin = 0;
00249         pMax = 0;
00250         return false;
00251     }
00252   }
00253   else
00254   {
00255     G4bool existsAfterClip=true;
00256 
00257     // Calculate rotated vertex coordinates
00258 
00259     G4int noPolygonVertices=0;
00260     G4ThreeVectorList* vertices
00261       = CreateRotatedVertices(pTransform,noPolygonVertices);
00262 
00263     if(pAxis == kXAxis || pAxis == kYAxis || pAxis == kZAxis)
00264     {
00265 
00266       pMin =  kInfinity;
00267       pMax = -kInfinity;
00268 
00269       for(G4ThreeVectorList::iterator it = vertices->begin();
00270           it < vertices->end(); it++)
00271       {
00272         if(pMin > (*it)[pAxis]) pMin = (*it)[pAxis];
00273         if((*it)[pAxis] < pVoxelLimit.GetMinExtent(pAxis))
00274         {
00275           pMin = pVoxelLimit.GetMinExtent(pAxis);
00276         }
00277         if(pMax < (*it)[pAxis])
00278         {
00279           pMax = (*it)[pAxis];
00280         }
00281         if((*it)[pAxis] > pVoxelLimit.GetMaxExtent(pAxis))
00282         {
00283           pMax = pVoxelLimit.GetMaxExtent(pAxis);
00284         }
00285       }
00286 
00287       if(pMin > pVoxelLimit.GetMaxExtent(pAxis)
00288       || pMax < pVoxelLimit.GetMinExtent(pAxis)) { existsAfterClip = false; }
00289     }
00290     else
00291     {
00292       pMin = 0;
00293       pMax = 0;
00294       existsAfterClip = false;
00295     }
00296     delete vertices;
00297     return existsAfterClip;
00298   }
00299   return true;
00300 }

G4double G4Paraboloid::CalculateSurfaceArea (  )  const [inline]

Definition at line 135 of file G4Paraboloid.icc.

References G4INCL::Math::pi, and sqr().

Referenced by GetPointOnSurface(), and GetSurfaceArea().

00136 {
00137   G4double h1, h2, A1, A2;
00138 
00139   h1 = k2/k1 + dz;
00140   h2 = k2/k1 - dz;
00141 
00142   // Calculate surface area for the paraboloid full paraboloid
00143   // cutoff at z = dz (not the cutoff area though).
00144 
00145   A1 = sqr(r2) + 4 * sqr(h1);
00146   A1 *= sqr(A1); // Sets A1 = A1^3
00147   A1 = CLHEP::pi * r2 /6 / sqr(h1) * ( std::sqrt(A1) - r2 * r2 * r2);
00148 
00149   // Calculate surface area for the paraboloid full paraboloid
00150   // cutoff at z = -dz (not the cutoff area though).
00151 
00152   A2 = sqr(r1) + 4 * sqr(h2);
00153   A2 *= sqr(A2);// Sets A2 = A2^3
00154 
00155   if(h2 != 0)
00156     { A2 = CLHEP::pi * r1 /6 / sqr(h2) * ( std::sqrt(A2) - r1 * r1 * r1); }
00157   else
00158     { A2 = 0.; }
00159 
00160   return fSurfaceArea = A1 - A2 + (sqr(r1) + sqr(r2))*CLHEP::pi;
00161 }

G4VSolid * G4Paraboloid::Clone (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 972 of file G4Paraboloid.cc.

References G4Paraboloid().

00973 {
00974   return new G4Paraboloid(*this);
00975 }

G4NURBS * G4Paraboloid::CreateNURBS (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 1149 of file G4Paraboloid.cc.

01150 {
01151   // Box for now!!!
01152   //
01153   return new G4NURBSbox(r1, r1, dz);
01154 }

G4Polyhedron * G4Paraboloid::CreatePolyhedron (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 1156 of file G4Paraboloid.cc.

Referenced by GetPolyhedron().

01157 {
01158   return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
01159 }

void G4Paraboloid::DescribeYourselfTo ( G4VGraphicsScene scene  )  const [virtual]

Implements G4VSolid.

Definition at line 1144 of file G4Paraboloid.cc.

References G4VGraphicsScene::AddSolid().

01145 {
01146   scene.AddSolid(*this);
01147 }

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

Implements G4VSolid.

Definition at line 575 of file G4Paraboloid.cc.

00576 {
00577   G4double safz = -dz+std::fabs(p.z());
00578   if(safz<0) { safz=0; }
00579   G4double safr = kInfinity;
00580 
00581   G4double rho = p.x()*p.x()+p.y()*p.y();
00582   G4double paraRho = (p.z()-k2)/k1;
00583   G4double sqrho = std::sqrt(rho);
00584 
00585   if(paraRho<0)
00586   {
00587     safr=sqrho-r2;
00588     if(safr>safz) { safz=safr; }
00589     return safz;
00590   }
00591 
00592   G4double sqprho = std::sqrt(paraRho);
00593   G4double dRho = sqrho-sqprho;
00594   if(dRho<0) { return safz; }
00595 
00596   G4double talf = -2.*k1*sqprho;
00597   G4double tmp  = 1+talf*talf;
00598   if(tmp<0.) { return safz; }
00599 
00600   G4double salf = talf/std::sqrt(tmp);
00601   safr = std::fabs(dRho*salf);
00602   if(safr>safz) { safz=safr; }
00603 
00604   return safz;
00605 }

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

Implements G4VSolid.

Definition at line 442 of file G4Paraboloid.cc.

References G4endl, G4Exception(), G4VSolid::GetName(), Inside(), JustWarning, G4VSolid::kCarTolerance, kInside, and sqr().

00444 {
00445   G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
00446   G4double tol2 = kCarTolerance*kCarTolerance;
00447   G4double tolh = 0.5*kCarTolerance;
00448 
00449   if(r2 && p.z() > - tolh + dz) 
00450   {
00451     // If the points is above check for intersection with upper edge.
00452 
00453     if(v.z() < 0)
00454     {
00455       G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz.
00456       if(sqr(p.x() + v.x()*intersection)
00457        + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance))
00458       {
00459         if(p.z() < tolh + dz)
00460           { return 0; }
00461         else
00462           { return intersection; }
00463       }
00464     }
00465     else  // Direction away, no possibility of intersection
00466     {
00467       return kInfinity;
00468     }
00469   }
00470   else if(r1 && p.z() < tolh - dz)
00471   {
00472     // If the points is belove check for intersection with lower edge.
00473 
00474     if(v.z() > 0)
00475     {
00476       G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz.
00477       if(sqr(p.x() + v.x()*intersection)
00478        + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance))
00479       {
00480         if(p.z() > -tolh - dz)
00481         {
00482           return 0;
00483         }
00484         else
00485         {
00486           return intersection;
00487         }
00488       }
00489     }
00490     else  // Direction away, no possibility of intersection
00491     {
00492       return kInfinity;
00493     }
00494   }
00495 
00496   G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(),
00497            vRho2 = v.perp2(), intersection,
00498            B = (k1 * p.z() + k2 - rho2) * vRho2;
00499 
00500   if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) )
00501     || (p.z() < - dz+kCarTolerance)
00502     || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside.
00503   {
00504     // Is there a problem with squaring rho twice?
00505 
00506     if(vRho2<tol2) // Needs to be treated seperately.
00507     {
00508       intersection = ((rho2 - k2)/k1 - p.z())/v.z();
00509       if(intersection < 0) { return kInfinity; }
00510       else if(std::fabs(p.z() + v.z() * intersection) <= dz)
00511       {
00512         return intersection;
00513       }
00514       else
00515       {
00516         return kInfinity;
00517       }
00518     }
00519     else if(A*A + B < 0) // No real intersections.
00520     {
00521       return kInfinity;
00522     }
00523     else
00524     {
00525       intersection = (A - std::sqrt(B + sqr(A))) / vRho2;
00526       if(intersection < 0)
00527       {
00528         return kInfinity;
00529       }
00530       else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh)
00531       {
00532         return intersection;
00533       }
00534       else
00535       {
00536         return kInfinity;
00537       }
00538     }
00539   }
00540   else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2)
00541   {
00542     // If this is true we're somewhere in the border.
00543 
00544     G4ThreeVector normal(p.x(), p.y(), -k1/2);
00545     if(normal.dot(v) <= 0)
00546       { return 0; }
00547     else
00548       { return kInfinity; }
00549   }
00550   else
00551   {
00552     std::ostringstream message;
00553     if(Inside(p) == kInside)
00554     {
00555       message << "Point p is inside! - " << GetName() << G4endl;
00556     }
00557     else
00558     {
00559       message << "Likely a problem in this function, for solid: " << GetName()
00560               << G4endl;
00561     }
00562     message << "          p = " << p * (1/mm) << " mm" << G4endl
00563             << "          v = " << v * (1/mm) << " mm";
00564     G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002",
00565                 JustWarning, message);
00566     return 0;
00567   }
00568 }

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

Implements G4VSolid.

Definition at line 922 of file G4Paraboloid.cc.

References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, G4VSolid::kCarTolerance, and kOutside.

00923 {
00924   G4double safe=0.0,rho,safeR,safeZ ;
00925   G4double tanRMax,secRMax,pRMax ;
00926 
00927 #ifdef G4SPECSDEBUG
00928   if( Inside(p) == kOutside )
00929   {
00930      G4cout << G4endl ;
00931      DumpInfo();
00932      std::ostringstream message;
00933      G4int oldprc = message.precision(16);
00934      message << "Point p is outside !?" << G4endl
00935              << "Position:" << G4endl
00936              << "   p.x() = "   << p.x()/mm << " mm" << G4endl
00937              << "   p.y() = "   << p.y()/mm << " mm" << G4endl
00938              << "   p.z() = "   << p.z()/mm << " mm";
00939      message.precision(oldprc) ;
00940      G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002",
00941                  JustWarning, message);
00942   }
00943 #endif
00944 
00945   rho = p.perp();
00946   safeZ = dz - std::fabs(p.z()) ;
00947 
00948   tanRMax = (r2 - r1)*0.5/dz ;
00949   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
00950   pRMax   = tanRMax*p.z() + (r1+r2)*0.5 ;
00951   safeR  = (pRMax - rho)/secRMax ;
00952 
00953   if (safeZ < safeR) { safe = safeZ; }
00954   else { safe = safeR; }
00955   if ( safe < 0.5 * kCarTolerance ) { safe = 0; }
00956   return safe ;
00957 }

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 611 of file G4Paraboloid.cc.

References G4endl, G4Exception(), Inside(), JustWarning, G4VSolid::kCarTolerance, kOutside, and sqr().

00616 {
00617   G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
00618   G4double vRho2 = v.perp2(), intersection;
00619   G4double tol2 = kCarTolerance*kCarTolerance;
00620   G4double tolh = 0.5*kCarTolerance;
00621 
00622   if(calcNorm) { *validNorm = false; }
00623 
00624   // We have that the particle p follows the line x = p + s * v
00625   // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and
00626   // z = p.z() + s * v.z()
00627   // The equation for all points on the surface (surface expanded for
00628   // to include all z) x^2 + y^2 = k1 * z + k2 => .. =>
00629   // => s = (A +- std::sqrt(A^2 + B)) / vRho2
00630   // where:
00631   //
00632   G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
00633   //
00634   // and:
00635   //
00636   G4double B = (-rho2 + paraRho2) * vRho2;
00637 
00638   if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2
00639     && std::fabs(p.z()) < dz - kCarTolerance)
00640   {
00641     // Make sure it's safely inside.
00642 
00643     if(v.z() > 0)
00644     {
00645       // It's heading upwards, check where it colides with the plane z = dz.
00646       // When it does, is that in the surface of the paraboloid.
00647       // z = p.z() + variable * v.z() for all points the particle can go.
00648       // => variable = (z - p.z()) / v.z() so intersection must be:
00649 
00650       intersection = (dz - p.z()) / v.z();
00651       G4ThreeVector ip = p + intersection * v; // Point of intersection.
00652 
00653       if(ip.perp2() < sqr(r2 + kCarTolerance))
00654       {
00655         if(calcNorm)
00656         {
00657           *n = G4ThreeVector(0, 0, 1);
00658           if(r2 < tolh || ip.perp2() > sqr(r2 - tolh))
00659           {
00660             *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00661             *n = n->unit();
00662           }
00663           *validNorm = true;
00664         }
00665         return intersection;
00666       }
00667     }
00668     else if(v.z() < 0)
00669     {
00670       // It's heading downwards, check were it colides with the plane z = -dz.
00671       // When it does, is that in the surface of the paraboloid.
00672       // z = p.z() + variable * v.z() for all points the particle can go.
00673       // => variable = (z - p.z()) / v.z() so intersection must be:
00674 
00675       intersection = (-dz - p.z()) / v.z();
00676       G4ThreeVector ip = p + intersection * v; // Point of intersection.
00677 
00678       if(ip.perp2() < sqr(r1 + tolh))
00679       {
00680         if(calcNorm)
00681         {
00682           *n = G4ThreeVector(0, 0, -1);
00683           if(r1 < tolh || ip.perp2() > sqr(r1 - tolh))
00684           {
00685             *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00686             *n = n->unit();
00687           }
00688           *validNorm = true;
00689         }
00690         return intersection;
00691       }
00692     }
00693 
00694     // Now check for collisions with paraboloid surface.
00695 
00696     if(vRho2 == 0) // Needs to be treated seperately.
00697     {
00698       intersection = ((rho2 - k2)/k1 - p.z())/v.z();
00699       if(calcNorm)
00700       {
00701         G4ThreeVector intersectionP = p + v * intersection;
00702         *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
00703         *n = n->unit();
00704 
00705         *validNorm = true;
00706       }
00707       return intersection;
00708     }
00709     else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
00710     {
00711       // intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
00712       // The above calculation has a precision problem:
00713       // known problem of solving quadratic equation with small A  
00714 
00715       A = A/vRho2;
00716       B = (k1 * p.z() + k2 - rho2)/vRho2;
00717       intersection = B/(-A + std::sqrt(B + sqr(A)));
00718       if(calcNorm)
00719       {
00720         G4ThreeVector intersectionP = p + v * intersection;
00721         *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
00722         *n = n->unit();
00723         *validNorm = true;
00724       }
00725       return intersection;
00726     }
00727     std::ostringstream message;
00728     message << "There is no intersection between given line and solid!"
00729             << G4endl
00730             << "          p = " << p << G4endl
00731             << "          v = " << v;
00732     G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
00733                 JustWarning, message);
00734 
00735     return kInfinity;
00736   }
00737   else if ( (rho2 < paraRho2 + kCarTolerance
00738          || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 )
00739          && std::fabs(p.z()) < dz + tolh) 
00740   {
00741     // If this is true we're somewhere in the border.
00742     
00743     G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
00744 
00745     if(std::fabs(p.z()) > dz - tolh)
00746     {
00747       // We're in the lower or upper edge
00748       //
00749       if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
00750       {             // If we're heading out of the object that is treated here
00751         if(calcNorm)
00752         {
00753           *validNorm = true;
00754           if(p.z() > 0)
00755             { *n = G4ThreeVector(0, 0, 1); }
00756           else
00757             { *n = G4ThreeVector(0, 0, -1); }
00758         }
00759         return 0;
00760       }
00761 
00762       if(v.z() == 0)
00763       {
00764         // Case where we're moving inside the surface needs to be
00765         // treated separately.
00766         // Distance until it goes out through a side is returned.
00767 
00768         G4double r = (p.z() > 0)? r2 : r1;
00769         G4double pDotV = p.dot(v);
00770         A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y()));
00771         intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2;
00772 
00773         if(calcNorm)
00774         {
00775           *validNorm = true;
00776 
00777           *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z()))
00778               + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y()
00779               * intersection, -k1/2).unit()).unit();
00780         }
00781         return intersection;
00782       }
00783     }
00784     //
00785     // Problem in the Logic :: Following condition for point on upper surface 
00786     //                         and Vz<0  will return 0 (Problem #1015), but
00787     //                         it has to return intersection with parabolic
00788     //                         surface or with lower plane surface (z = -dz)
00789     // The logic has to be  :: If not found intersection until now,
00790     // do not exit but continue to search for possible intersection.
00791     // Only for point situated on both borders (Z and parabolic)
00792     // this condition has to be taken into account and done later
00793     //
00794     //
00795     // else if(normal.dot(v) >= 0)
00796     // {
00797     //   if(calcNorm)
00798     //   {
00799     //     *validNorm = true;
00800     //     *n = normal.unit();
00801     //   }
00802     //   return 0;
00803     // }
00804 
00805     if(v.z() > 0)
00806     {
00807       // Check for collision with upper edge.
00808 
00809       intersection = (dz - p.z()) / v.z();
00810       G4ThreeVector ip = p + intersection * v;
00811 
00812       if(ip.perp2() < sqr(r2 - tolh))
00813       {
00814         if(calcNorm)
00815         {
00816           *validNorm = true;
00817           *n = G4ThreeVector(0, 0, 1);
00818         }
00819         return intersection;
00820       }
00821       else if(ip.perp2() < sqr(r2 + tolh))
00822       {
00823         if(calcNorm)
00824         {
00825           *validNorm = true;
00826           *n = G4ThreeVector(0, 0, 1)
00827              + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00828           *n = n->unit();
00829         }
00830         return intersection;
00831       }
00832     }
00833     if( v.z() < 0)
00834     {
00835       // Check for collision with lower edge.
00836 
00837       intersection = (-dz - p.z()) / v.z();
00838       G4ThreeVector ip = p + intersection * v;
00839 
00840       if(ip.perp2() < sqr(r1 - tolh))
00841       {
00842         if(calcNorm)
00843         {
00844           *validNorm = true;
00845           *n = G4ThreeVector(0, 0, -1);
00846         }
00847         return intersection;
00848       }
00849       else if(ip.perp2() < sqr(r1 + tolh))
00850       {
00851         if(calcNorm)
00852         {
00853           *validNorm = true;
00854           *n = G4ThreeVector(0, 0, -1)
00855              + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00856           *n = n->unit();
00857         }
00858         return intersection;
00859       }
00860     }
00861 
00862     // Note: comparison with zero below would not be correct !
00863     //
00864     if(std::fabs(vRho2) > tol2) // precision error in the calculation of
00865     {                           // intersection = (A+std::sqrt(B+sqr(A)))/vRho2
00866       A = A/vRho2;
00867       B = (k1 * p.z() + k2 - rho2);
00868       if(std::fabs(B)>kCarTolerance)
00869       {
00870         B = (B)/vRho2;
00871         intersection = B/(-A + std::sqrt(B + sqr(A)));
00872       }
00873       else                      // Point is On both borders: Z and parabolic
00874       {                         // solution depends on normal.dot(v) sign
00875         if(normal.dot(v) >= 0)
00876         {
00877           if(calcNorm)
00878           {
00879             *validNorm = true;
00880             *n = normal.unit();
00881           }
00882           return 0;
00883         }
00884         intersection = 2.*A;
00885       }
00886     }
00887     else
00888     {
00889       intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
00890     }
00891 
00892     if(calcNorm)
00893     {
00894       *validNorm = true;
00895       *n = G4ThreeVector(p.x() + intersection * v.x(), p.y()
00896          + intersection * v.y(), - k1 / 2);
00897       *n = n->unit();
00898     }
00899     return intersection;
00900   }
00901   else
00902   {
00903 #ifdef G4SPECSDEBUG
00904     if(kOutside == Inside(p))
00905     {
00906       G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
00907                   JustWarning, "Point p is outside!");
00908     }
00909     else
00910       G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
00911                   JustWarning, "There's an error in this functions code.");
00912 #endif
00913     return kInfinity;
00914   }
00915   return 0;
00916 } 

G4double G4Paraboloid::GetCubicVolume (  )  [inline, virtual]

Reimplemented from G4VSolid.

Definition at line 123 of file G4Paraboloid.icc.

00124 {
00125   if(fCubicVolume != 0 ) {;}
00126   else
00127   {
00128     fCubicVolume = CLHEP::twopi * k2 * dz;
00129   }
00130   return fCubicVolume;
00131 }

G4GeometryType G4Paraboloid::GetEntityType (  )  const [virtual]

Implements G4VSolid.

Definition at line 963 of file G4Paraboloid.cc.

00964 {
00965   return G4String("G4Paraboloid");
00966 }

G4ThreeVector G4Paraboloid::GetPointOnSurface (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 1002 of file G4Paraboloid.cc.

References CalculateSurfaceArea(), G4INCL::Math::pi, and sqr().

01003 {
01004   G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea;
01005   G4double z = RandFlat::shoot(0.,1.);
01006   G4double phi = RandFlat::shoot(0., twopi);
01007   if(pi*(sqr(r1) + sqr(r2))/A >= z)
01008   {
01009     G4double rho;
01010     if(pi * sqr(r1) / A > z)
01011     {
01012       rho = r1 * std::sqrt(RandFlat::shoot(0., 1.));
01013       return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz);
01014     }
01015     else
01016     {
01017       rho = r2 * std::sqrt(RandFlat::shoot(0., 1));
01018       return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz);
01019     }
01020   }
01021   else
01022   {
01023     z = RandFlat::shoot(0., 1.)*2*dz - dz;
01024     return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi),
01025                          std::sqrt(z*k1 + k2)*std::sin(phi), z);
01026   }
01027 }

G4Polyhedron * G4Paraboloid::GetPolyhedron (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 1162 of file G4Paraboloid.cc.

References CreatePolyhedron(), fpPolyhedron, and G4Polyhedron::GetNumberOfRotationStepsAtTimeOfCreation().

01163 {
01164   if (!fpPolyhedron ||
01165       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01166       fpPolyhedron->GetNumberOfRotationSteps())
01167   {
01168     delete fpPolyhedron;
01169     fpPolyhedron = CreatePolyhedron();
01170   }
01171   return fpPolyhedron;
01172 }

G4double G4Paraboloid::GetRadiusMinusZ (  )  const [inline]

Definition at line 51 of file G4Paraboloid.icc.

Referenced by G4GDMLWriteSolids::ParaboloidWrite().

00052 {
00053   return r1;
00054 }

G4double G4Paraboloid::GetRadiusPlusZ (  )  const [inline]

Definition at line 45 of file G4Paraboloid.icc.

Referenced by G4GDMLWriteSolids::ParaboloidWrite().

00046 {
00047   return r2;
00048 }

G4double G4Paraboloid::GetSurfaceArea (  )  [inline, virtual]

Reimplemented from G4VSolid.

Definition at line 164 of file G4Paraboloid.icc.

References CalculateSurfaceArea().

00165 {
00166   if(fSurfaceArea == 0.) CalculateSurfaceArea();
00167 
00168   return fSurfaceArea;
00169 }

G4double G4Paraboloid::GetZHalfLength (  )  const [inline]

Definition at line 39 of file G4Paraboloid.icc.

Referenced by G4GDMLWriteSolids::ParaboloidWrite().

00040 {
00041   return dz;
00042 }

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

Implements G4VSolid.

Definition at line 306 of file G4Paraboloid.cc.

References G4VSolid::kCarTolerance, kInside, kOutside, kSurface, and sqr().

Referenced by DistanceToIn(), and DistanceToOut().

00307 {
00308   // First check is  the point is above or below the solid.
00309   //
00310   if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; }
00311 
00312   G4double rho2 = p.perp2(),
00313            rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance),
00314            A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
00315  
00316   if(A < 0 && sqr(A) > rhoSurfTimesTol2)
00317   {
00318     // Actually checking rho < radius of paraboloid at z = p.z().
00319     // We're either inside or in lower/upper cutoff area.
00320    
00321     if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
00322     {
00323       // We're in the upper/lower cutoff area, sides have a paraboloid shape
00324       // maybe further checks should be made to make these nicer
00325 
00326       return kSurface;
00327     }
00328     else
00329     {
00330       return kInside;
00331     }
00332   }
00333   else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
00334   {
00335     // We're in the parabolic surface.
00336 
00337     return kSurface;
00338   }
00339   else
00340   {
00341     return kOutside;
00342   }
00343 }

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

Definition at line 124 of file G4Paraboloid.cc.

References dz, fCubicVolume, fpPolyhedron, fSurfaceArea, k1, k2, G4VSolid::operator=(), r1, and r2.

00125 {
00126    // Check assignment to self
00127    //
00128    if (this == &rhs)  { return *this; }
00129 
00130    // Copy base class data
00131    //
00132    G4VSolid::operator=(rhs);
00133 
00134    // Copy data
00135    //
00136    fpPolyhedron = 0; 
00137    fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
00138    dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2;
00139 
00140    return *this;
00141 }

void G4Paraboloid::SetRadiusMinusZ ( G4double  R1  )  [inline]

Definition at line 101 of file G4Paraboloid.icc.

References FatalException, G4Exception(), and sqr().

00102 {
00103   if(pR1 < 0 || pR1 >= r2)
00104   {
00105     G4Exception("G4Paraboloid::SetRadiusMinusZ()", "GeomSolids0002", 
00106                 FatalException, "Invalid dimensions.");
00107   }
00108   else
00109   {
00110     r1 = pR1;
00111     k1 = (sqr(r2) - sqr(r1)) / (2 * dz);
00112     k2 = (sqr(r2) + sqr(r1)) / 2;
00113 
00114     // This informs GetSurfaceArea() and GetCubicVolume() that it needs
00115     // to recalculate buffered value.
00116     //
00117     fSurfaceArea = 0.; 
00118     fCubicVolume = 0.;
00119   }
00120 }

void G4Paraboloid::SetRadiusPlusZ ( G4double  R2  )  [inline]

Definition at line 79 of file G4Paraboloid.icc.

References FatalException, G4Exception(), and sqr().

00080 {
00081   if(pR2 <= 0 || pR2 <= r1)
00082   {
00083     G4Exception("G4Paraboloid::SetRadiusPlusZ()", "GeomSolids0002", 
00084                 FatalException, "Invalid dimensions.");
00085   }
00086   else
00087   {
00088     r2 = pR2;
00089     k1 = (sqr(r2) - sqr(r1)) / (2 * dz);
00090     k2 = (sqr(r2) + sqr(r1)) / 2;
00091 
00092     // This informs GetSurfaceArea() and GetCubicVolume() that it needs
00093     // to recalculate buffered value.
00094     //
00095     fSurfaceArea = 0.; 
00096     fCubicVolume = 0.;
00097   }
00098 }

void G4Paraboloid::SetZHalfLength ( G4double  dz  )  [inline]

Definition at line 57 of file G4Paraboloid.icc.

References FatalException, G4Exception(), and sqr().

00058 {
00059   if(pDz <= 0)
00060   {
00061     G4Exception("G4Paraboloid::SetZHalfLength()", "GeomSolids0002", 
00062                 FatalException, "Invalid dimensions.");
00063   }
00064   else
00065   {
00066     dz = pDz;
00067     k1 = (sqr(r2) - sqr(r1)) / (2 * dz);
00068     k2 = (sqr(r2) + sqr(r1)) / 2;
00069 
00070     // This informs GetSurfaceArea() and GetCubicVolume() that it needs
00071     // to recalculate buffered value.
00072     //
00073     fSurfaceArea = 0.; 
00074     fCubicVolume = 0.;
00075   }
00076 }

std::ostream & G4Paraboloid::StreamInfo ( std::ostream &  os  )  const [virtual]

Implements G4VSolid.

Definition at line 981 of file G4Paraboloid.cc.

References G4VSolid::GetName().

00982 {
00983   G4int oldprc = os.precision(16);
00984   os << "-----------------------------------------------------------\n"
00985      << "    *** Dump for solid - " << GetName() << " ***\n"
00986      << "    ===================================================\n"
00987      << " Solid type: G4Paraboloid\n"
00988      << " Parameters: \n"
00989      << "    z half-axis:   " << dz/mm << " mm \n"
00990      << "    radius at -dz: " << r1/mm << " mm \n"
00991      << "    radius at dz:  " << r2/mm << " mm \n"
00992      << "-----------------------------------------------------------\n";
00993   os.precision(oldprc);
00994 
00995   return os;
00996 }

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

Implements G4VSolid.

Definition at line 348 of file G4Paraboloid.cc.

References G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, CLHEP::detail::n, and sqr().

00349 {
00350   G4ThreeVector n(0, 0, 0);
00351   if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance)
00352   {
00353     // If above or below just return normal vector for the cutoff plane.
00354 
00355     n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z()));
00356   }
00357   else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
00358   {
00359     // This means we're somewhere in the plane z = dz or z = -dz.
00360     // (As far as the program is concerned anyway.
00361 
00362     if(p.z() < 0) // Are we in upper or lower plane?
00363     {
00364       if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance))
00365       {
00366         n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit();
00367       }
00368       else if(r1 < 0.5 * kCarTolerance
00369            || p.perp2() > sqr(r1 - 0.5 * kCarTolerance))
00370       {
00371         n = G4ThreeVector(p.x(), p.y(), 0.).unit()
00372           + G4ThreeVector(0., 0., -1.).unit();
00373         n = n.unit();
00374       }
00375       else
00376       {
00377         n = G4ThreeVector(0., 0., -1.);
00378       }
00379     }
00380     else
00381     {
00382       if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance))
00383       {
00384         n = G4ThreeVector(p.x(), p.y(), 0.).unit();
00385       }
00386       else if(r2 < 0.5 * kCarTolerance
00387            || p.perp2() > sqr(r2 - 0.5 * kCarTolerance))
00388       {
00389         n = G4ThreeVector(p.x(), p.y(), 0.).unit()
00390           + G4ThreeVector(0., 0., 1.).unit();
00391         n = n.unit();
00392       }
00393       else
00394       {
00395         n = G4ThreeVector(0., 0., 1.);
00396       }
00397     }
00398   }
00399   else
00400   {
00401     G4double rho2 = p.perp2();
00402     G4double rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance);
00403     G4double A = rho2 - ((k1 *p.z() + k2)
00404                + 0.25 * kCarTolerance * kCarTolerance);
00405 
00406     if(A < 0 && sqr(A) > rhoSurfTimesTol2)
00407     {
00408       // Actually checking rho < radius of paraboloid at z = p.z().
00409       // We're inside.
00410 
00411       if(p.mag2() != 0) { n = p.unit(); }
00412     }
00413     else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
00414     {
00415       // We're in the parabolic surface.
00416 
00417       n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
00418     }
00419     else
00420     {
00421       n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
00422     }
00423   }
00424 
00425   if(n.mag2() == 0)
00426   {
00427     std::ostringstream message;
00428     message << "No normal defined for this point p." << G4endl
00429             << "          p = " << 1 / mm * p << " mm";
00430     G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002",
00431                 JustWarning, message);
00432   }
00433   return n;
00434 }


Field Documentation

G4Polyhedron* G4Paraboloid::fpPolyhedron [mutable, protected]

Definition at line 140 of file G4Paraboloid.hh.

Referenced by GetPolyhedron(), and operator=().


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:49 2013 for Geant4 by  doxygen 1.4.7