G4Sphere Class Reference

#include <G4Sphere.hh>

Inheritance diagram for G4Sphere:

G4CSGSolid G4VSolid

Public Member Functions

 G4Sphere (const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
 ~G4Sphere ()
G4double GetInnerRadius () const
G4double GetOuterRadius () const
G4double GetStartPhiAngle () const
G4double GetDeltaPhiAngle () const
G4double GetStartThetaAngle () const
G4double GetDeltaThetaAngle () const
void SetInnerRadius (G4double newRMin)
void SetOuterRadius (G4double newRmax)
void SetStartPhiAngle (G4double newSphi, G4bool trig=true)
void SetDeltaPhiAngle (G4double newDphi)
void SetStartThetaAngle (G4double newSTheta)
void SetDeltaThetaAngle (G4double newDTheta)
G4double GetCubicVolume ()
G4double GetSurfaceArea ()
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
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
G4ThreeVector GetPointOnSurface () const
G4VSolidClone () const
std::ostream & StreamInfo (std::ostream &os) const
G4VisExtent GetExtent () const
void DescribeYourselfTo (G4VGraphicsScene &scene) const
G4PolyhedronCreatePolyhedron () const
G4NURBSCreateNURBS () const
 G4Sphere (__void__ &)
 G4Sphere (const G4Sphere &rhs)
G4Sphereoperator= (const G4Sphere &rhs)
G4double GetRmin () const
G4double GetRmax () const
G4double GetSPhi () const
G4double GetDPhi () const
G4double GetSTheta () const
G4double GetDTheta () const
G4double GetInsideRadius () const
void SetInsideRadius (G4double newRmin)

Detailed Description

Definition at line 77 of file G4Sphere.hh.


Constructor & Destructor Documentation

G4Sphere::G4Sphere ( const G4String pName,
G4double  pRmin,
G4double  pRmax,
G4double  pSPhi,
G4double  pDPhi,
G4double  pSTheta,
G4double  pDTheta 
)

Definition at line 90 of file G4Sphere.cc.

References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), and G4GeometryTolerance::GetRadialTolerance().

Referenced by Clone().

00094   : G4CSGSolid(pName), fEpsilon(2.e-11),
00095     fFullPhiSphere(true), fFullThetaSphere(true)
00096 {
00097   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00098 
00099   // Check radii and set radial tolerances
00100 
00101   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00102   if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
00103   {
00104     std::ostringstream message;
00105     message << "Invalid radii for Solid: " << GetName() << G4endl
00106             << "        pRmin = " << pRmin << ", pRmax = " << pRmax;
00107     G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
00108                 FatalException, message);
00109   }
00110   fRmin=pRmin; fRmax=pRmax;
00111   fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
00112   fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
00113 
00114   // Check angles
00115 
00116   CheckPhiAngles(pSPhi, pDPhi);
00117   CheckThetaAngles(pSTheta, pDTheta);
00118 }

G4Sphere::~G4Sphere (  ) 

Definition at line 141 of file G4Sphere.cc.

00142 {
00143 }

G4Sphere::G4Sphere ( __void__ &   ) 

Definition at line 125 of file G4Sphere.cc.

00126   : G4CSGSolid(a), fRminTolerance(0.), fRmaxTolerance(0.),
00127     kAngTolerance(0.), kRadTolerance(0.), fEpsilon(0.),
00128     fRmin(0.), fRmax(0.), fSPhi(0.), fDPhi(0.), fSTheta(0.),
00129     fDTheta(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
00130     sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), hDPhi(0.), cPhi(0.),
00131     ePhi(0.), sinSTheta(0.), cosSTheta(0.), sinETheta(0.), cosETheta(0.),
00132     tanSTheta(0.), tanSTheta2(0.), tanETheta(0.), tanETheta2(0.), eTheta(0.),
00133     fFullPhiSphere(false), fFullThetaSphere(false), fFullSphere(true)
00134 {
00135 }

G4Sphere::G4Sphere ( const G4Sphere rhs  ) 

Definition at line 149 of file G4Sphere.cc.

00150   : G4CSGSolid(rhs), fRminTolerance(rhs.fRminTolerance),
00151     fRmaxTolerance(rhs.fRmaxTolerance), kAngTolerance(rhs.kAngTolerance),
00152     kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
00153     fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
00154     fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
00155     sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
00156     cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
00157     sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
00158     sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
00159     hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
00160     sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
00161     sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
00162     tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
00163     tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
00164     fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
00165     fFullSphere(rhs.fFullSphere)
00166 {
00167 }


Member Function Documentation

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

Implements G4VSolid.

Definition at line 220 of file G4Sphere.cc.

References G4VSolid::CalculateClippedPolygonExtent(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), and G4AffineTransform::TransformPoint().

00224 {
00225   if ( fFullSphere )
00226   {
00227     // Special case handling for solid spheres-shells
00228     // (rotation doesn't influence).
00229     // Compute x/y/z mins and maxs for bounding box respecting limits,
00230     // with early returns if outside limits. Then switch() on pAxis,
00231     // and compute exact x and y limit for x/y case
00232       
00233     G4double xoffset,xMin,xMax;
00234     G4double yoffset,yMin,yMax;
00235     G4double zoffset,zMin,zMax;
00236 
00237     G4double diff1,diff2,delta,maxDiff,newMin,newMax;
00238     G4double xoff1,xoff2,yoff1,yoff2;
00239 
00240     xoffset=pTransform.NetTranslation().x();
00241     xMin=xoffset-fRmax;
00242     xMax=xoffset+fRmax;
00243     if (pVoxelLimit.IsXLimited())
00244     {
00245       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00246         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00247       {
00248         return false;
00249       }
00250       else
00251       {
00252         if (xMin<pVoxelLimit.GetMinXExtent())
00253         {
00254           xMin=pVoxelLimit.GetMinXExtent();
00255         }
00256         if (xMax>pVoxelLimit.GetMaxXExtent())
00257         {
00258           xMax=pVoxelLimit.GetMaxXExtent();
00259         }
00260       }
00261     }
00262 
00263     yoffset=pTransform.NetTranslation().y();
00264     yMin=yoffset-fRmax;
00265     yMax=yoffset+fRmax;
00266     if (pVoxelLimit.IsYLimited())
00267     {
00268       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00269         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00270       {
00271         return false;
00272       }
00273       else
00274       {
00275         if (yMin<pVoxelLimit.GetMinYExtent())
00276         {
00277           yMin=pVoxelLimit.GetMinYExtent();
00278         }
00279         if (yMax>pVoxelLimit.GetMaxYExtent())
00280         {
00281           yMax=pVoxelLimit.GetMaxYExtent();
00282         }
00283       }
00284     }
00285 
00286     zoffset=pTransform.NetTranslation().z();
00287     zMin=zoffset-fRmax;
00288     zMax=zoffset+fRmax;
00289     if (pVoxelLimit.IsZLimited())
00290     {
00291       if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00292         || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
00293       {
00294         return false;
00295       }
00296       else
00297       {
00298         if (zMin<pVoxelLimit.GetMinZExtent())
00299         {
00300           zMin=pVoxelLimit.GetMinZExtent();
00301         }
00302         if (zMax>pVoxelLimit.GetMaxZExtent())
00303         {
00304           zMax=pVoxelLimit.GetMaxZExtent();
00305         }
00306       }
00307     }
00308 
00309     // Known to cut sphere
00310 
00311     switch (pAxis)
00312     {
00313       case kXAxis:
00314         yoff1=yoffset-yMin;
00315         yoff2=yMax-yoffset;
00316         if ((yoff1>=0) && (yoff2>=0))
00317         {
00318           // Y limits cross max/min x => no change
00319           //
00320           pMin=xMin;
00321           pMax=xMax;
00322         }
00323         else
00324         {
00325           // Y limits don't cross max/min x => compute max delta x,
00326           // hence new mins/maxs
00327           //
00328           delta=fRmax*fRmax-yoff1*yoff1;
00329           diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00330           delta=fRmax*fRmax-yoff2*yoff2;
00331           diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00332           maxDiff=(diff1>diff2) ? diff1:diff2;
00333           newMin=xoffset-maxDiff;
00334           newMax=xoffset+maxDiff;
00335           pMin=(newMin<xMin) ? xMin : newMin;
00336           pMax=(newMax>xMax) ? xMax : newMax;
00337         }
00338         break;
00339       case kYAxis:
00340         xoff1=xoffset-xMin;
00341         xoff2=xMax-xoffset;
00342         if ((xoff1>=0) && (xoff2>=0))
00343         {
00344           // X limits cross max/min y => no change
00345           //
00346           pMin=yMin;
00347           pMax=yMax;
00348         }
00349         else
00350         {
00351           // X limits don't cross max/min y => compute max delta y,
00352           // hence new mins/maxs
00353           //
00354           delta=fRmax*fRmax-xoff1*xoff1;
00355           diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00356           delta=fRmax*fRmax-xoff2*xoff2;
00357           diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00358           maxDiff=(diff1>diff2) ? diff1:diff2;
00359           newMin=yoffset-maxDiff;
00360           newMax=yoffset+maxDiff;
00361           pMin=(newMin<yMin) ? yMin : newMin;
00362           pMax=(newMax>yMax) ? yMax : newMax;
00363         }
00364         break;
00365       case kZAxis:
00366         pMin=zMin;
00367         pMax=zMax;
00368         break;
00369       default:
00370         break;
00371     }
00372     pMin-=kCarTolerance;
00373     pMax+=kCarTolerance;
00374 
00375     return true;  
00376   }
00377   else       // Transformed cutted sphere
00378   {
00379     G4int i,j,noEntries,noBetweenSections;
00380     G4bool existsAfterClip=false;
00381 
00382     // Calculate rotated vertex coordinates
00383 
00384     G4ThreeVectorList* vertices;
00385     G4int  noPolygonVertices ;
00386     vertices=CreateRotatedVertices(pTransform,noPolygonVertices);
00387 
00388     pMin=+kInfinity;
00389     pMax=-kInfinity;
00390 
00391     noEntries=vertices->size();  // noPolygonVertices*noPhiCrossSections
00392     noBetweenSections=noEntries-noPolygonVertices;
00393 
00394     G4ThreeVectorList ThetaPolygon ;
00395     for (i=0;i<noEntries;i+=noPolygonVertices)
00396     {
00397       for(j=0;j<(noPolygonVertices/2)-1;j++)
00398       {
00399         ThetaPolygon.push_back((*vertices)[i+j]) ;      
00400         ThetaPolygon.push_back((*vertices)[i+j+1]) ;      
00401         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;      
00402         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;      
00403         CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
00404         ThetaPolygon.clear() ;
00405       }
00406     }
00407     for (i=0;i<noBetweenSections;i+=noPolygonVertices)
00408     {
00409       for(j=0;j<noPolygonVertices-1;j++)
00410       {
00411         ThetaPolygon.push_back((*vertices)[i+j]) ;      
00412         ThetaPolygon.push_back((*vertices)[i+j+1]) ;      
00413         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;      
00414         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;      
00415         CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
00416         ThetaPolygon.clear() ;
00417       }
00418       ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;      
00419       ThetaPolygon.push_back((*vertices)[i]) ;  
00420       ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;      
00421       ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;      
00422       CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
00423       ThetaPolygon.clear() ;
00424     }
00425       
00426     if ((pMin!=kInfinity) || (pMax!=-kInfinity))
00427     {
00428       existsAfterClip=true;
00429 
00430       // Add 2*tolerance to avoid precision troubles
00431       //
00432       pMin-=kCarTolerance;
00433       pMax+=kCarTolerance;
00434     }
00435     else
00436     {
00437       // Check for case where completely enveloping clipping volume
00438       // If point inside then we are confident that the solid completely
00439       // envelopes the clipping volume. Hence set min/max extents according
00440       // to clipping volume extents along the specified axis.
00441 
00442       G4ThreeVector clipCentre(
00443           (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00444           (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00445           (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
00446         
00447       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
00448       {
00449         existsAfterClip=true;
00450         pMin=pVoxelLimit.GetMinExtent(pAxis);
00451         pMax=pVoxelLimit.GetMaxExtent(pAxis);
00452       }
00453     }
00454     delete vertices;
00455     return existsAfterClip;
00456   }
00457 }

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

Reimplemented from G4VSolid.

Definition at line 3009 of file G4Sphere.cc.

References G4Sphere().

03010 {
03011   return new G4Sphere(*this);
03012 }

void G4Sphere::ComputeDimensions ( G4VPVParameterisation p,
const G4int  n,
const G4VPhysicalVolume pRep 
) [virtual]

Reimplemented from G4VSolid.

Definition at line 209 of file G4Sphere.cc.

References G4VPVParameterisation::ComputeDimensions().

00212 {
00213   p->ComputeDimensions(*this,n,pRep);
00214 }

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

Reimplemented from G4VSolid.

Definition at line 3187 of file G4Sphere.cc.

03188 {
03189   return new G4NURBSbox (fRmax, fRmax, fRmax);       // Box for now!!!
03190 }

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

Reimplemented from G4VSolid.

Definition at line 3182 of file G4Sphere.cc.

03183 {
03184   return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
03185 }

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

Implements G4VSolid.

Definition at line 3177 of file G4Sphere.cc.

References G4VGraphicsScene::AddSolid().

03178 {
03179   scene.AddSolid (*this);
03180 }

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

Implements G4VSolid.

Definition at line 1801 of file G4Sphere.cc.

References G4INCL::Math::pi.

01802 {
01803   G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
01804   G4double rho2,rds,rho;
01805   G4double cosPsi;
01806   G4double pTheta,dTheta1,dTheta2;
01807   rho2=p.x()*p.x()+p.y()*p.y();
01808   rds=std::sqrt(rho2+p.z()*p.z());
01809   rho=std::sqrt(rho2);
01810 
01811   //
01812   // Distance to r shells
01813   //    
01814   if (fRmin)
01815   {
01816     safeRMin=fRmin-rds;
01817     safeRMax=rds-fRmax;
01818     if (safeRMin>safeRMax)
01819     {
01820       safe=safeRMin;
01821     }
01822     else
01823     {
01824       safe=safeRMax;
01825     }
01826   }
01827   else
01828   {
01829     safe=rds-fRmax;
01830   }
01831 
01832   //
01833   // Distance to phi extent
01834   //
01835   if (!fFullPhiSphere && rho)
01836   {
01837     // Psi=angle from central phi to point
01838     //
01839     cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
01840     if (cosPsi<std::cos(hDPhi))
01841     {
01842       // Point lies outside phi range
01843       //
01844       if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
01845       {
01846         safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
01847       }
01848       else
01849       {
01850         safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
01851       }
01852       if (safePhi>safe)  { safe=safePhi; }
01853     }
01854   }
01855   //
01856   // Distance to Theta extent
01857   //    
01858   if ((rds!=0.0) && (!fFullThetaSphere))
01859   {
01860     pTheta=std::acos(p.z()/rds);
01861     if (pTheta<0)  { pTheta+=pi; }
01862     dTheta1=fSTheta-pTheta;
01863     dTheta2=pTheta-eTheta;
01864     if (dTheta1>dTheta2)
01865     {
01866       if (dTheta1>=0)             // WHY ???????????
01867       {
01868         safeTheta=rds*std::sin(dTheta1);
01869         if (safe<=safeTheta)
01870         {
01871           safe=safeTheta;
01872         }
01873       }
01874     }
01875     else
01876     {
01877       if (dTheta2>=0)
01878       {
01879         safeTheta=rds*std::sin(dTheta2);
01880         if (safe<=safeTheta)
01881         {
01882           safe=safeTheta;
01883         }
01884       }
01885     }
01886   }
01887 
01888   if (safe<0)  { safe=0; }
01889   return safe;
01890 }

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

Implements G4VSolid.

Definition at line 867 of file G4Sphere.cc.

References G4VSolid::kCarTolerance, and G4INCL::Math::pi.

00869 {
00870   G4double snxt = kInfinity ;      // snxt = default return value
00871   G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
00872   G4double tolSTheta=0., tolETheta=0. ;
00873   const G4double dRmax = 100.*fRmax;
00874 
00875   static const G4double halfCarTolerance = kCarTolerance*0.5;
00876   static const G4double halfAngTolerance = kAngTolerance*0.5;
00877   const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
00878   const G4double halfRminTolerance = fRminTolerance*0.5;
00879   const G4double tolORMin2 = (fRmin>halfRminTolerance)
00880                ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
00881   const G4double tolIRMin2 =
00882                (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
00883   const G4double tolORMax2 =
00884                (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
00885   const G4double tolIRMax2 =
00886                (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
00887 
00888   // Intersection point
00889   //
00890   G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
00891 
00892   // Phi intersection
00893   //
00894   G4double Comp ; 
00895 
00896   // Phi precalcs
00897   //
00898   G4double Dist, cosPsi ;
00899 
00900   // Theta precalcs
00901   //
00902   G4double dist2STheta, dist2ETheta ;
00903   G4double t1, t2, b, c, d2, d, sd = kInfinity ;
00904 
00905   // General Precalcs
00906   //
00907   rho2 = p.x()*p.x() + p.y()*p.y() ;
00908   rad2 = rho2 + p.z()*p.z() ;
00909   pTheta = std::atan2(std::sqrt(rho2),p.z()) ;
00910 
00911   pDotV2d = p.x()*v.x() + p.y()*v.y() ;
00912   pDotV3d = pDotV2d + p.z()*v.z() ;
00913 
00914   // Theta precalcs
00915   //
00916   if (!fFullThetaSphere)
00917   {
00918     tolSTheta = fSTheta - halfAngTolerance ;
00919     tolETheta = eTheta + halfAngTolerance ;
00920   }
00921 
00922   // Outer spherical shell intersection
00923   // - Only if outside tolerant fRmax
00924   // - Check for if inside and outer G4Sphere heading through solid (-> 0)
00925   // - No intersect -> no intersection with G4Sphere
00926   //
00927   // Shell eqn: x^2+y^2+z^2=RSPH^2
00928   //
00929   // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
00930   //
00931   // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
00932   // =>      rad2        +2sd(pDotV3d)       +sd^2                =R^2
00933   //
00934   // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
00935 
00936   c = rad2 - fRmax*fRmax ;
00937 
00938   if (c > fRmaxTolerance*fRmax)
00939   {
00940     // If outside tolerant boundary of outer G4Sphere
00941     // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]
00942 
00943     d2 = pDotV3d*pDotV3d - c ;
00944 
00945     if ( d2 >= 0 )
00946     {
00947       sd = -pDotV3d - std::sqrt(d2) ;
00948 
00949       if (sd >= 0 )
00950       {
00951         if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen on
00952         {               // 64 bits systems. Split long distances and recompute
00953           G4double fTerm = sd-std::fmod(sd,dRmax);
00954           sd = fTerm + DistanceToIn(p+fTerm*v,v);
00955         } 
00956         xi   = p.x() + sd*v.x() ;
00957         yi   = p.y() + sd*v.y() ;
00958         rhoi = std::sqrt(xi*xi + yi*yi) ;
00959 
00960         if (!fFullPhiSphere && rhoi)    // Check phi intersection
00961         {
00962           cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
00963 
00964           if (cosPsi >= cosHDPhiOT)
00965           {
00966             if (!fFullThetaSphere)   // Check theta intersection
00967             {
00968               zi = p.z() + sd*v.z() ;
00969 
00970               // rhoi & zi can never both be 0
00971               // (=>intersect at origin =>fRmax=0)
00972               //
00973               iTheta = std::atan2(rhoi,zi) ;
00974               if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
00975               {
00976                 return snxt = sd ;
00977               }
00978             }
00979             else
00980             {
00981               return snxt=sd;
00982             }
00983           }
00984         }
00985         else
00986         {
00987           if (!fFullThetaSphere)    // Check theta intersection
00988           {
00989             zi = p.z() + sd*v.z() ;
00990 
00991             // rhoi & zi can never both be 0
00992             // (=>intersect at origin => fRmax=0 !)
00993             //
00994             iTheta = std::atan2(rhoi,zi) ;
00995             if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
00996             {
00997               return snxt=sd;
00998             }
00999           }
01000           else
01001           {
01002             return snxt = sd;
01003           }
01004         }          
01005       }
01006     }
01007     else    // No intersection with G4Sphere
01008     {
01009       return snxt=kInfinity;
01010     }
01011   }
01012   else
01013   {
01014     // Inside outer radius
01015     // check not inside, and heading through G4Sphere (-> 0 to in)
01016 
01017     d2 = pDotV3d*pDotV3d - c ;
01018 
01019     if ( (rad2 > tolIRMax2)
01020       && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
01021     {
01022       if (!fFullPhiSphere)
01023       {
01024         // Use inner phi tolerant boundary -> if on tolerant
01025         // phi boundaries, phi intersect code handles leaving/entering checks
01026 
01027         cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
01028 
01029         if (cosPsi>=cosHDPhiIT)
01030         { 
01031           // inside radii, delta r -ve, inside phi
01032 
01033           if ( !fFullThetaSphere )
01034           {
01035             if ( (pTheta >= tolSTheta + kAngTolerance)
01036               && (pTheta <= tolETheta - kAngTolerance) )
01037             {
01038               return snxt=0;
01039             }
01040           }
01041           else    // strictly inside Theta in both cases
01042           {
01043             return snxt=0;
01044           }
01045         }
01046       }
01047       else
01048       {
01049         if ( !fFullThetaSphere )
01050         {
01051           if ( (pTheta >= tolSTheta + kAngTolerance)
01052             && (pTheta <= tolETheta - kAngTolerance) )
01053           {
01054             return snxt=0;
01055           }
01056         }
01057         else   // strictly inside Theta in both cases
01058         {
01059           return snxt=0;
01060         }
01061       }
01062     }
01063   }
01064 
01065   // Inner spherical shell intersection
01066   // - Always farthest root, because would have passed through outer
01067   //   surface first.
01068   // - Tolerant check if travelling through solid
01069 
01070   if (fRmin)
01071   {
01072     c  = rad2 - fRmin*fRmin ;
01073     d2 = pDotV3d*pDotV3d - c ;
01074 
01075     // Within tolerance inner radius of inner G4Sphere
01076     // Check for immediate entry/already inside and travelling outwards
01077 
01078     if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
01079       && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
01080     {
01081       if ( !fFullPhiSphere )
01082       {
01083         // Use inner phi tolerant boundary -> if on tolerant
01084         // phi boundaries, phi intersect code handles leaving/entering checks
01085 
01086         cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
01087         if (cosPsi >= cosHDPhiIT)
01088         { 
01089           // inside radii, delta r -ve, inside phi
01090           //
01091           if ( !fFullThetaSphere )
01092           {
01093             if ( (pTheta >= tolSTheta + kAngTolerance)
01094               && (pTheta <= tolETheta - kAngTolerance) )
01095             {
01096               return snxt=0;
01097             }
01098           }
01099           else
01100           {
01101             return snxt = 0 ;
01102           }
01103         }
01104       }
01105       else
01106       {
01107         if ( !fFullThetaSphere )
01108         {
01109           if ( (pTheta >= tolSTheta + kAngTolerance)
01110             && (pTheta <= tolETheta - kAngTolerance) )
01111           {
01112             return snxt = 0 ;
01113           }
01114         }
01115         else
01116         {
01117           return snxt=0;
01118         }
01119       }
01120     }
01121     else   // Not special tolerant case
01122     {
01123       if (d2 >= 0)
01124       {
01125         sd = -pDotV3d + std::sqrt(d2) ;
01126         if ( sd >= halfRminTolerance )  // It was >= 0 ??
01127         {
01128           xi   = p.x() + sd*v.x() ;
01129           yi   = p.y() + sd*v.y() ;
01130           rhoi = std::sqrt(xi*xi+yi*yi) ;
01131 
01132           if ( !fFullPhiSphere && rhoi )   // Check phi intersection
01133           {
01134             cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
01135 
01136             if (cosPsi >= cosHDPhiOT)
01137             {
01138               if ( !fFullThetaSphere )  // Check theta intersection
01139               {
01140                 zi = p.z() + sd*v.z() ;
01141 
01142                 // rhoi & zi can never both be 0
01143                 // (=>intersect at origin =>fRmax=0)
01144                 //
01145                 iTheta = std::atan2(rhoi,zi) ;
01146                 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
01147                 {
01148                   snxt = sd;
01149                 }
01150               }
01151               else
01152               {
01153                 snxt=sd;
01154               }
01155             }
01156           }
01157           else
01158           {
01159             if ( !fFullThetaSphere )   // Check theta intersection
01160             {
01161               zi = p.z() + sd*v.z() ;
01162 
01163               // rhoi & zi can never both be 0
01164               // (=>intersect at origin => fRmax=0 !)
01165               //
01166               iTheta = std::atan2(rhoi,zi) ;
01167               if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
01168               {
01169                 snxt = sd;
01170               }
01171             }
01172             else
01173             {
01174               snxt = sd;
01175             }
01176           }
01177         }
01178       }
01179     }
01180   }
01181 
01182   // Phi segment intersection
01183   //
01184   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
01185   //
01186   // o NOTE: Large duplication of code between sphi & ephi checks
01187   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
01188   //            intersection check <=0 -> >=0
01189   //         -> Should use some form of loop Construct
01190   //
01191   if ( !fFullPhiSphere )
01192   {
01193     // First phi surface ('S'tarting phi)
01194     // Comp = Component in outwards normal dirn
01195     //
01196     Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
01197                     
01198     if ( Comp < 0 )
01199     {
01200       Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
01201 
01202       if (Dist < halfCarTolerance)
01203       {
01204         sd = Dist/Comp ;
01205 
01206         if (sd < snxt)
01207         {
01208           if ( sd > 0 )
01209           {
01210             xi    = p.x() + sd*v.x() ;
01211             yi    = p.y() + sd*v.y() ;
01212             zi    = p.z() + sd*v.z() ;
01213             rhoi2 = xi*xi + yi*yi   ;
01214             radi2 = rhoi2 + zi*zi   ;
01215           }
01216           else
01217           {
01218             sd    = 0     ;
01219             xi    = p.x() ;
01220             yi    = p.y() ;
01221             zi    = p.z() ;
01222             rhoi2 = rho2  ;
01223             radi2 = rad2  ;
01224           }
01225           if ( (radi2 <= tolORMax2)
01226             && (radi2 >= tolORMin2)
01227             && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
01228           {
01229             // Check theta intersection
01230             // rhoi & zi can never both be 0
01231             // (=>intersect at origin =>fRmax=0)
01232             //
01233             if ( !fFullThetaSphere )
01234             {
01235               iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
01236               if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
01237               {
01238                 // r and theta intersections good
01239                 // - check intersecting with correct half-plane
01240 
01241                 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
01242                 {
01243                   snxt = sd;
01244                 }
01245               }
01246             }
01247             else
01248             {
01249               snxt = sd;
01250             }
01251           }
01252         }
01253       }
01254     }
01255 
01256     // Second phi surface ('E'nding phi)
01257     // Component in outwards normal dirn
01258 
01259     Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
01260         
01261     if (Comp < 0)
01262     {
01263       Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
01264       if ( Dist < halfCarTolerance )
01265       {
01266         sd = Dist/Comp ;
01267 
01268         if ( sd < snxt )
01269         {
01270           if (sd > 0)
01271           {
01272             xi    = p.x() + sd*v.x() ;
01273             yi    = p.y() + sd*v.y() ;
01274             zi    = p.z() + sd*v.z() ;
01275             rhoi2 = xi*xi + yi*yi   ;
01276             radi2 = rhoi2 + zi*zi   ;
01277           }
01278           else
01279           {
01280             sd    = 0     ;
01281             xi    = p.x() ;
01282             yi    = p.y() ;
01283             zi    = p.z() ;
01284             rhoi2 = rho2  ;
01285             radi2 = rad2  ;
01286           }
01287           if ( (radi2 <= tolORMax2)
01288             && (radi2 >= tolORMin2)
01289             && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
01290           {
01291             // Check theta intersection
01292             // rhoi & zi can never both be 0
01293             // (=>intersect at origin =>fRmax=0)
01294             //
01295             if ( !fFullThetaSphere )
01296             {
01297               iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
01298               if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
01299               {
01300                 // r and theta intersections good
01301                 // - check intersecting with correct half-plane
01302 
01303                 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
01304                 {
01305                   snxt = sd;
01306                 }
01307               }
01308             }
01309             else
01310             {
01311               snxt = sd;
01312             }
01313           }
01314         }
01315       }
01316     }
01317   }
01318 
01319   // Theta segment intersection
01320 
01321   if ( !fFullThetaSphere )
01322   {
01323 
01324     // Intersection with theta surfaces
01325     // Known failure cases:
01326     // o  Inside tolerance of stheta surface, skim
01327     //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
01328     //
01329     //    To solve: Check 2nd root of etheta surface in addition to stheta
01330     //
01331     // o  start/end theta is exactly pi/2 
01332     // Intersections with cones
01333     //
01334     // Cone equation: x^2+y^2=z^2tan^2(t)
01335     //
01336     // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
01337     //
01338     // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
01339     //       + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
01340     //
01341     // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
01342 
01343     if (fSTheta)
01344     {
01345       dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
01346     }
01347     else
01348     {
01349       dist2STheta = kInfinity ;
01350     }
01351     if ( eTheta < pi )
01352     {
01353       dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
01354     }
01355     else
01356     {
01357       dist2ETheta=kInfinity;
01358     }      
01359     if ( pTheta < tolSTheta )
01360     {
01361       // Inside (theta<stheta-tol) stheta cone
01362       // First root of stheta cone, second if first root -ve
01363 
01364       t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01365       t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01366       if (t1)
01367       {   
01368         b  = t2/t1 ;
01369         c  = dist2STheta/t1 ;
01370         d2 = b*b - c ;
01371 
01372         if ( d2 >= 0 )
01373         {
01374           d  = std::sqrt(d2) ;
01375           sd = -b - d ;    // First root
01376           zi = p.z() + sd*v.z();
01377 
01378           if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) )
01379           {
01380             sd = -b+d;    // Second root
01381           }
01382           if ((sd >= 0) && (sd < snxt))
01383           {
01384             xi    = p.x() + sd*v.x();
01385             yi    = p.y() + sd*v.y();
01386             zi    = p.z() + sd*v.z();
01387             rhoi2 = xi*xi + yi*yi;
01388             radi2 = rhoi2 + zi*zi;
01389             if ( (radi2 <= tolORMax2)
01390               && (radi2 >= tolORMin2)
01391               && (zi*(fSTheta - halfpi) <= 0) )
01392             {
01393               if ( !fFullPhiSphere && rhoi2 )  // Check phi intersection
01394               {
01395                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01396                 if (cosPsi >= cosHDPhiOT)
01397                 {
01398                   snxt = sd;
01399                 }
01400               }
01401               else
01402               {
01403                 snxt = sd;
01404               }
01405             }
01406           }
01407         }
01408       }
01409 
01410       // Possible intersection with ETheta cone. 
01411       // Second >= 0 root should be considered
01412         
01413       if ( eTheta < pi )
01414       {
01415         t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01416         t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01417         if (t1)
01418         { 
01419           b  = t2/t1 ;
01420           c  = dist2ETheta/t1 ;
01421           d2 = b*b - c ;
01422 
01423           if (d2 >= 0)
01424           {
01425             d  = std::sqrt(d2) ;
01426             sd = -b + d ;    // Second root
01427 
01428             if ( (sd >= 0) && (sd < snxt) )
01429             {
01430               xi    = p.x() + sd*v.x() ;
01431               yi    = p.y() + sd*v.y() ;
01432               zi    = p.z() + sd*v.z() ;
01433               rhoi2 = xi*xi + yi*yi   ;
01434               radi2 = rhoi2 + zi*zi   ;
01435 
01436               if ( (radi2 <= tolORMax2)
01437                 && (radi2 >= tolORMin2)
01438                 && (zi*(eTheta - halfpi) <= 0) )
01439               {
01440                 if (!fFullPhiSphere && rhoi2)   // Check phi intersection
01441                 {
01442                   cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01443                   if (cosPsi >= cosHDPhiOT)
01444                   {
01445                     snxt = sd;
01446                   }
01447                 }
01448                 else
01449                 {
01450                   snxt = sd;
01451                 }
01452               }
01453             }
01454           }
01455         }
01456       }
01457     }  
01458     else if ( pTheta > tolETheta ) 
01459     { 
01460       // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
01461       // Inside (theta > etheta+tol) e-theta cone
01462       // First root of etheta cone, second if first root 'imaginary'
01463 
01464       t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01465       t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01466       if (t1)
01467       {  
01468         b  = t2/t1 ;
01469         c  = dist2ETheta/t1 ;
01470         d2 = b*b - c ;
01471 
01472         if (d2 >= 0)
01473         {
01474           d  = std::sqrt(d2) ;
01475           sd = -b - d ;    // First root
01476           zi = p.z() + sd*v.z();
01477 
01478           if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) )
01479           {
01480             sd = -b + d ;           // second root
01481           }
01482           if ( (sd >= 0) && (sd < snxt) )
01483           {
01484             xi    = p.x() + sd*v.x() ;
01485             yi    = p.y() + sd*v.y() ;
01486             zi    = p.z() + sd*v.z() ;
01487             rhoi2 = xi*xi + yi*yi   ;
01488             radi2 = rhoi2 + zi*zi   ;
01489 
01490             if ( (radi2 <= tolORMax2)
01491               && (radi2 >= tolORMin2) 
01492               && (zi*(eTheta - halfpi) <= 0) )
01493             {
01494               if (!fFullPhiSphere && rhoi2)  // Check phi intersection
01495               {
01496                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01497                 if (cosPsi >= cosHDPhiOT)
01498                 {
01499                   snxt = sd;
01500                 }
01501               }
01502               else
01503               {
01504                 snxt = sd;
01505               }
01506             }
01507           }
01508         }
01509       }
01510 
01511       // Possible intersection with STheta cone. 
01512       // Second >= 0 root should be considered
01513         
01514       if ( fSTheta )
01515       {
01516         t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01517         t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01518         if (t1)
01519         { 
01520           b  = t2/t1 ;
01521           c  = dist2STheta/t1 ;
01522           d2 = b*b - c ;
01523 
01524           if (d2 >= 0)
01525           {
01526             d  = std::sqrt(d2) ;
01527             sd = -b + d ;    // Second root
01528 
01529             if ( (sd >= 0) && (sd < snxt) )
01530             {
01531               xi    = p.x() + sd*v.x() ;
01532               yi    = p.y() + sd*v.y() ;
01533               zi    = p.z() + sd*v.z() ;
01534               rhoi2 = xi*xi + yi*yi   ;
01535               radi2 = rhoi2 + zi*zi   ;
01536 
01537               if ( (radi2 <= tolORMax2)
01538                 && (radi2 >= tolORMin2)
01539                 && (zi*(fSTheta - halfpi) <= 0) )
01540               {
01541                 if (!fFullPhiSphere && rhoi2)   // Check phi intersection
01542                 {
01543                   cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01544                   if (cosPsi >= cosHDPhiOT)
01545                   {
01546                     snxt = sd;
01547                   }
01548                 }
01549                 else
01550                 {
01551                   snxt = sd;
01552                 }
01553               }
01554             }
01555           }
01556         }
01557       }  
01558     }     
01559     else if ( (pTheta < tolSTheta + kAngTolerance)
01560            && (fSTheta > halfAngTolerance) )
01561     {
01562       // In tolerance of stheta
01563       // If entering through solid [r,phi] => 0 to in
01564       // else try 2nd root
01565 
01566       t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01567       if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
01568         || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
01569         || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
01570       {
01571         if (!fFullPhiSphere && rho2)  // Check phi intersection
01572         {
01573           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
01574           if (cosPsi >= cosHDPhiIT)
01575           {
01576             return 0 ;
01577           }
01578         }
01579         else
01580         {
01581           return 0 ;
01582         }
01583       }
01584 
01585       // Not entering immediately/travelling through
01586 
01587       t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01588       if (t1)
01589       { 
01590         b  = t2/t1 ;
01591         c  = dist2STheta/t1 ;
01592         d2 = b*b - c ;
01593 
01594         if (d2 >= 0)
01595         {
01596           d  = std::sqrt(d2) ;
01597           sd = -b + d ;
01598           if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
01599           {   // ^^^^^^^^^^^^^^^^^^^^^  shouldn't it be >=0 instead ?
01600             xi    = p.x() + sd*v.x() ;
01601             yi    = p.y() + sd*v.y() ;
01602             zi    = p.z() + sd*v.z() ;
01603             rhoi2 = xi*xi + yi*yi   ;
01604             radi2 = rhoi2 + zi*zi   ;
01605 
01606             if ( (radi2 <= tolORMax2)
01607               && (radi2 >= tolORMin2)
01608               && (zi*(fSTheta - halfpi) <= 0) )
01609             {
01610               if ( !fFullPhiSphere && rhoi2 )    // Check phi intersection
01611               {
01612                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01613                 if ( cosPsi >= cosHDPhiOT )
01614                 {
01615                   snxt = sd;
01616                 }
01617               }
01618               else
01619               {
01620                 snxt = sd;
01621               }
01622             }
01623           }
01624         }
01625       }
01626     }   
01627     else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
01628     {
01629 
01630       // In tolerance of etheta
01631       // If entering through solid [r,phi] => 0 to in
01632       // else try 2nd root
01633 
01634       t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01635 
01636       if (   ((t2<0) && (eTheta < halfpi)
01637           && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
01638         ||   ((t2>=0) && (eTheta > halfpi)
01639           && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
01640         ||   ((v.z()>0) && (eTheta == halfpi)
01641           && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))  )
01642       {
01643         if (!fFullPhiSphere && rho2)   // Check phi intersection
01644         {
01645           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
01646           if (cosPsi >= cosHDPhiIT)
01647           {
01648             return 0 ;
01649           }
01650         }
01651         else
01652         {
01653           return 0 ;
01654         }
01655       }
01656 
01657       // Not entering immediately/travelling through
01658 
01659       t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01660       if (t1)
01661       { 
01662         b  = t2/t1 ;
01663         c  = dist2ETheta/t1 ;
01664         d2 = b*b - c ;
01665 
01666         if (d2 >= 0)
01667         {
01668           d  = std::sqrt(d2) ;
01669           sd = -b + d ;
01670         
01671           if ( (sd >= halfCarTolerance)
01672             && (sd < snxt) && (eTheta > halfpi) )
01673           {
01674             xi    = p.x() + sd*v.x() ;
01675             yi    = p.y() + sd*v.y() ;
01676             zi    = p.z() + sd*v.z() ;
01677             rhoi2 = xi*xi + yi*yi   ;
01678             radi2 = rhoi2 + zi*zi   ;
01679 
01680             if ( (radi2 <= tolORMax2)
01681               && (radi2 >= tolORMin2)
01682               && (zi*(eTheta - halfpi) <= 0) )
01683             {
01684               if (!fFullPhiSphere && rhoi2)   // Check phi intersection
01685               {
01686                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01687                 if (cosPsi >= cosHDPhiOT)
01688                 {
01689                   snxt = sd;
01690                 }
01691               }
01692               else
01693               {
01694                 snxt = sd;
01695               }
01696             }
01697           }
01698         } 
01699       }   
01700     }  
01701     else
01702     {
01703       // stheta+tol<theta<etheta-tol
01704       // For BOTH stheta & etheta check 2nd root for validity [r,phi]
01705 
01706       t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01707       t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01708       if (t1)
01709       { 
01710         b  = t2/t1;
01711         c  = dist2STheta/t1 ;
01712         d2 = b*b - c ;
01713 
01714         if (d2 >= 0)
01715         {
01716           d  = std::sqrt(d2) ;
01717           sd = -b + d ;    // second root
01718 
01719           if ((sd >= 0) && (sd < snxt))
01720           {
01721             xi    = p.x() + sd*v.x() ;
01722             yi    = p.y() + sd*v.y() ;
01723             zi    = p.z() + sd*v.z() ;
01724             rhoi2 = xi*xi + yi*yi   ;
01725             radi2 = rhoi2 + zi*zi   ;
01726 
01727             if ( (radi2 <= tolORMax2)
01728               && (radi2 >= tolORMin2)
01729               && (zi*(fSTheta - halfpi) <= 0) )
01730             {
01731               if (!fFullPhiSphere && rhoi2)   // Check phi intersection
01732               {
01733                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01734                 if (cosPsi >= cosHDPhiOT)
01735                 {
01736                   snxt = sd;
01737                 }
01738               }
01739               else
01740               {
01741                 snxt = sd;
01742               }
01743             }
01744           }
01745         }
01746       }        
01747       t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01748       t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01749       if (t1)
01750       {   
01751         b  = t2/t1 ;
01752         c  = dist2ETheta/t1 ;
01753         d2 = b*b - c ;
01754 
01755         if (d2 >= 0)
01756         {
01757           d  = std::sqrt(d2) ;
01758           sd = -b + d;    // second root
01759 
01760           if ((sd >= 0) && (sd < snxt))
01761           {
01762             xi    = p.x() + sd*v.x() ;
01763             yi    = p.y() + sd*v.y() ;
01764             zi    = p.z() + sd*v.z() ;
01765             rhoi2 = xi*xi + yi*yi   ;
01766             radi2 = rhoi2 + zi*zi   ;
01767 
01768             if ( (radi2 <= tolORMax2)
01769               && (radi2 >= tolORMin2)
01770               && (zi*(eTheta - halfpi) <= 0) )
01771             {
01772               if (!fFullPhiSphere && rhoi2)   // Check phi intersection
01773               {
01774                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01775                 if ( cosPsi >= cosHDPhiOT )
01776                 {
01777                   snxt = sd;
01778                 }
01779               }
01780               else
01781               {
01782                 snxt = sd;
01783               }
01784             }
01785           }
01786         }
01787       }
01788     }  
01789   }
01790   return snxt;
01791 }

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

Implements G4VSolid.

Definition at line 2781 of file G4Sphere.cc.

References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, kOutside, and G4INCL::Math::pi.

02782 {
02783   G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
02784   G4double rho2,rds,rho;
02785   G4double pTheta,dTheta1,dTheta2;
02786   rho2=p.x()*p.x()+p.y()*p.y();
02787   rds=std::sqrt(rho2+p.z()*p.z());
02788   rho=std::sqrt(rho2);
02789 
02790 #ifdef G4CSGDEBUG
02791   if( Inside(p) == kOutside )
02792   {
02793      G4int old_prc = G4cout.precision(16);
02794      G4cout << G4endl;
02795      DumpInfo();
02796      G4cout << "Position:"  << G4endl << G4endl ;
02797      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
02798      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
02799      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
02800      G4cout.precision(old_prc) ;
02801      G4Exception("G4Sphere::DistanceToOut(p)",
02802                  "GeomSolids1002", JustWarning, "Point p is outside !?" );
02803   }
02804 #endif
02805 
02806   //
02807   // Distance to r shells
02808   //    
02809   if (fRmin)
02810   {
02811     safeRMin=rds-fRmin;
02812     safeRMax=fRmax-rds;
02813     if (safeRMin<safeRMax)
02814     {
02815       safe=safeRMin;
02816     }
02817     else
02818     {
02819       safe=safeRMax;
02820     }
02821   }
02822   else
02823   {
02824     safe=fRmax-rds;
02825   }
02826 
02827   //
02828   // Distance to phi extent
02829   //
02830   if (!fFullPhiSphere && rho)
02831   {
02832     if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
02833     {
02834       safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
02835     }
02836     else
02837     {
02838       safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
02839     }
02840     if (safePhi<safe)  { safe=safePhi; }
02841   }
02842 
02843   //
02844   // Distance to Theta extent
02845   //    
02846   if (rds)
02847   {
02848     pTheta=std::acos(p.z()/rds);
02849     if (pTheta<0)  { pTheta+=pi; }
02850     dTheta1=pTheta-fSTheta;
02851     dTheta2=eTheta-pTheta;
02852     if (dTheta1<dTheta2)  { safeTheta=rds*std::sin(dTheta1); }
02853     else                  { safeTheta=rds*std::sin(dTheta2); }
02854     if (safe>safeTheta)   { safe=safeTheta; }
02855   }
02856 
02857   if (safe<0)  { safe=0; }
02858   return safe;
02859 }

G4double G4Sphere::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 1897 of file G4Sphere.cc.

References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, and G4INCL::Math::pi.

01902 {
01903   G4double snxt = kInfinity;     // snxt is default return value
01904   G4double sphi= kInfinity,stheta= kInfinity;
01905   ESide side=kNull,sidephi=kNull,sidetheta=kNull;  
01906 
01907   static const G4double halfCarTolerance = kCarTolerance*0.5;
01908   static const G4double halfAngTolerance = kAngTolerance*0.5;
01909   const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
01910   const G4double halfRminTolerance = fRminTolerance*0.5;
01911   const G4double Rmax_plus  = fRmax + halfRmaxTolerance;
01912   const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
01913   G4double t1,t2;
01914   G4double b,c,d;
01915 
01916   // Variables for phi intersection:
01917 
01918   G4double pDistS,compS,pDistE,compE,sphi2,vphi;
01919     
01920   G4double rho2,rad2,pDotV2d,pDotV3d;
01921 
01922   G4double xi,yi,zi;      // Intersection point
01923 
01924   // Theta precals
01925   //
01926   G4double rhoSecTheta;
01927   G4double dist2STheta, dist2ETheta, distTheta;
01928   G4double d2,sd;
01929 
01930   // General Precalcs
01931   //
01932   rho2 = p.x()*p.x()+p.y()*p.y();
01933   rad2 = rho2+p.z()*p.z();
01934 
01935   pDotV2d = p.x()*v.x()+p.y()*v.y();
01936   pDotV3d = pDotV2d+p.z()*v.z();
01937 
01938   // Radial Intersections from G4Sphere::DistanceToIn
01939   //
01940   // Outer spherical shell intersection
01941   // - Only if outside tolerant fRmax
01942   // - Check for if inside and outer G4Sphere heading through solid (-> 0)
01943   // - No intersect -> no intersection with G4Sphere
01944   //
01945   // Shell eqn: x^2+y^2+z^2=RSPH^2
01946   //
01947   // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
01948   //
01949   // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
01950   // =>      rad2        +2sd(pDotV3d)       +sd^2                =R^2
01951   //
01952   // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
01953 
01954   if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
01955   {
01956     c = rad2 - fRmax*fRmax;
01957 
01958     if (c < fRmaxTolerance*fRmax) 
01959     {
01960       // Within tolerant Outer radius 
01961       // 
01962       // The test is
01963       //     rad  - fRmax < 0.5*kRadTolerance
01964       // =>  rad  < fRmax + 0.5*kRadTol
01965       // =>  rad2 < (fRmax + 0.5*kRadTol)^2
01966       // =>  rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
01967       // =>  rad2 - fRmax^2    <~    fRmax*kRadTol 
01968 
01969       d2 = pDotV3d*pDotV3d - c;
01970 
01971       if( (c >- fRmaxTolerance*fRmax)       // on tolerant surface
01972        && ((pDotV3d >=0) || (d2 < 0)) )     // leaving outside from Rmax 
01973                                             // not re-entering
01974       {
01975         if(calcNorm)
01976         {
01977           *validNorm = true ;
01978           *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
01979         }
01980         return snxt = 0;
01981       }
01982       else 
01983       {
01984         snxt = -pDotV3d+std::sqrt(d2);    // second root since inside Rmax
01985         side =  kRMax ; 
01986       }
01987     }
01988 
01989     // Inner spherical shell intersection:
01990     // Always first >=0 root, because would have passed
01991     // from outside of Rmin surface .
01992 
01993     if (fRmin)
01994     {
01995       c  = rad2 - fRmin*fRmin;
01996       d2 = pDotV3d*pDotV3d - c;
01997 
01998       if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
01999       {
02000         if ( (c < fRminTolerance*fRmin)              // leaving from Rmin
02001           && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
02002         {
02003           if(calcNorm)  { *validNorm = false; }  // Rmin surface is concave         
02004           return snxt = 0 ;
02005         }
02006         else
02007         {  
02008           if ( d2 >= 0. )
02009           {
02010             sd = -pDotV3d-std::sqrt(d2);
02011 
02012             if ( sd >= 0. )     // Always intersect Rmin first
02013             {
02014               snxt = sd ;
02015               side = kRMin ;
02016             }
02017           }
02018         }
02019       }
02020     }
02021   }
02022 
02023   // Theta segment intersection
02024 
02025   if ( !fFullThetaSphere )
02026   {
02027     // Intersection with theta surfaces
02028     //
02029     // Known failure cases:
02030     // o  Inside tolerance of stheta surface, skim
02031     //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
02032     //
02033     //    To solve: Check 2nd root of etheta surface in addition to stheta
02034     //
02035     // o  start/end theta is exactly pi/2 
02036     //
02037     // Intersections with cones
02038     //
02039     // Cone equation: x^2+y^2=z^2tan^2(t)
02040     //
02041     // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
02042     //
02043     // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
02044     //       + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
02045     //
02046     // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
02047     //
02048   
02049     if(fSTheta) // intersection with first cons
02050     {
02051       if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
02052       {
02053         if( v.z() > 0. ) 
02054         {
02055           if ( std::fabs( p.z() ) <= halfRmaxTolerance )
02056           {
02057             if(calcNorm)
02058             {
02059               *validNorm = true;
02060               *n = G4ThreeVector(0.,0.,1.);
02061             }
02062             return snxt = 0 ;
02063           }  
02064           stheta    = -p.z()/v.z();
02065           sidetheta = kSTheta;
02066         }
02067       }
02068       else // kons is not plane 
02069       {
02070         t1          = 1-v.z()*v.z()*(1+tanSTheta2);
02071         t2          = pDotV2d-p.z()*v.z()*tanSTheta2;  // ~vDotN if p on cons
02072         dist2STheta = rho2-p.z()*p.z()*tanSTheta2;     // t3
02073 
02074         distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
02075 
02076         if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
02077         {                                      // v parallel to kons
02078           if( v.z() > 0. )
02079           {
02080             if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
02081             {
02082               if( (fSTheta < halfpi) && (p.z() > 0.) )
02083               {
02084                 if( calcNorm )  { *validNorm = false; }
02085                 return snxt = 0.;
02086               }
02087               else if( (fSTheta > halfpi) && (p.z() <= 0) )
02088               {
02089                 if( calcNorm ) 
02090                 {
02091                   *validNorm = true;
02092                   if (rho2)
02093                   {
02094                     rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
02095                    
02096                     *n = G4ThreeVector( p.x()/rhoSecTheta,   
02097                                         p.y()/rhoSecTheta,
02098                                         std::sin(fSTheta)  );
02099                   }
02100                   else *n = G4ThreeVector(0.,0.,1.);
02101                 }
02102                 return snxt = 0.;               
02103               }
02104             }
02105             stheta    = -0.5*dist2STheta/t2;
02106             sidetheta = kSTheta;
02107           }  
02108         }      // 2nd order equation, 1st root of fSTheta cone,
02109         else   // 2nd if 1st root -ve
02110         {
02111           if( std::fabs(distTheta) < halfRmaxTolerance )
02112           {
02113             if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
02114             {
02115               if( calcNorm ) 
02116               {
02117                 *validNorm = true;
02118                 if (rho2)
02119                 {
02120                   rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
02121                    
02122                   *n = G4ThreeVector( p.x()/rhoSecTheta,   
02123                                       p.y()/rhoSecTheta,
02124                                       std::sin(fSTheta)  );
02125                 }
02126                 else  { *n = G4ThreeVector(0.,0.,1.); }
02127               }
02128               return snxt = 0.;
02129             }
02130             else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
02131             {
02132               if( calcNorm )  { *validNorm = false; }
02133               return snxt = 0.;
02134             }                               
02135           }
02136           b  = t2/t1;
02137           c  = dist2STheta/t1;
02138           d2 = b*b - c ;
02139 
02140           if ( d2 >= 0. )
02141           {
02142             d = std::sqrt(d2);
02143 
02144             if( fSTheta > halfpi )
02145             {
02146               sd = -b - d;         // First root
02147 
02148               if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
02149                ||  (sd < 0.)  || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) )     ) 
02150               {
02151                 sd = -b + d ; // 2nd root
02152               }
02153               if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )  
02154               {
02155                 stheta    = sd;
02156                 sidetheta = kSTheta;
02157               }
02158             }
02159             else // sTheta < pi/2, concave surface, no normal
02160             {
02161               sd = -b - d;         // First root
02162 
02163               if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
02164                 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) )   )
02165               {
02166                 sd = -b + d ; // 2nd root
02167               }
02168               if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )  
02169               {
02170                 stheta    = sd;
02171                 sidetheta = kSTheta;
02172               }            
02173             }
02174           }
02175         }
02176       }
02177     }
02178     if (eTheta < pi) // intersection with second cons
02179     {
02180       if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
02181       {
02182         if( v.z() < 0. ) 
02183         {
02184           if ( std::fabs( p.z() ) <= halfRmaxTolerance )
02185           {
02186             if(calcNorm)
02187             {
02188               *validNorm = true;
02189               *n = G4ThreeVector(0.,0.,-1.);
02190             }
02191             return snxt = 0 ;
02192           }  
02193           sd = -p.z()/v.z();
02194 
02195           if( sd < stheta )
02196           {
02197             stheta    = sd;
02198             sidetheta = kETheta;
02199           }
02200         }
02201       }
02202       else // kons is not plane 
02203       {
02204         t1          = 1-v.z()*v.z()*(1+tanETheta2);
02205         t2          = pDotV2d-p.z()*v.z()*tanETheta2;  // ~vDotN if p on cons
02206         dist2ETheta = rho2-p.z()*p.z()*tanETheta2;     // t3
02207 
02208         distTheta = std::sqrt(rho2)-p.z()*tanETheta;
02209 
02210         if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
02211         {                                      // v parallel to kons
02212           if( v.z() < 0. )
02213           {
02214             if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
02215             {
02216               if( (eTheta > halfpi) && (p.z() < 0.) )
02217               {
02218                 if( calcNorm )  { *validNorm = false; }
02219                 return snxt = 0.;
02220               }
02221               else if ( (eTheta < halfpi) && (p.z() >= 0) )
02222               {
02223                 if( calcNorm ) 
02224                 {
02225                   *validNorm = true;
02226                   if (rho2)
02227                   {
02228                     rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
02229                     *n = G4ThreeVector( p.x()/rhoSecTheta,   
02230                                         p.y()/rhoSecTheta,
02231                                         -sinETheta  );
02232                   }
02233                   else  { *n = G4ThreeVector(0.,0.,-1.); }
02234                 }
02235                 return snxt = 0.;               
02236               }
02237             }
02238             sd = -0.5*dist2ETheta/t2;
02239 
02240             if( sd < stheta )
02241             {
02242               stheta    = sd;
02243               sidetheta = kETheta;
02244             }
02245           }  
02246         }      // 2nd order equation, 1st root of fSTheta cone
02247         else   // 2nd if 1st root -ve
02248         {
02249           if ( std::fabs(distTheta) < halfRmaxTolerance )
02250           {
02251             if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
02252             {
02253               if( calcNorm ) 
02254               {
02255                 *validNorm = true;
02256                 if (rho2)
02257                 {
02258                     rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
02259                     *n = G4ThreeVector( p.x()/rhoSecTheta,   
02260                                         p.y()/rhoSecTheta,
02261                                         -sinETheta  );
02262                 }
02263                 else *n = G4ThreeVector(0.,0.,-1.);
02264               }                           
02265               return snxt = 0.;
02266             }
02267             else if ( (eTheta > halfpi)
02268                    && (t2 < 0.) && (p.z() <=0.) ) // leave
02269             {
02270               if( calcNorm )  { *validNorm = false; }
02271               return snxt = 0.;
02272             }                               
02273           }
02274           b  = t2/t1;
02275           c  = dist2ETheta/t1;
02276           d2 = b*b - c ;
02277 
02278           if ( d2 >= 0. )
02279           {
02280             d = std::sqrt(d2);
02281 
02282             if( eTheta < halfpi )
02283             {
02284               sd = -b - d;         // First root
02285 
02286               if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
02287                || (sd < 0.) ) 
02288               {
02289                 sd = -b + d ; // 2nd root
02290               }
02291               if( sd > halfRmaxTolerance )  
02292               {
02293                 if( sd < stheta )
02294                 {
02295                   stheta    = sd;
02296                   sidetheta = kETheta;
02297                 }
02298               }
02299             }
02300             else // sTheta+fDTheta > pi/2, concave surface, no normal
02301             {
02302               sd = -b - d;         // First root
02303 
02304               if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
02305                 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
02306               {
02307                 sd = -b + d ; // 2nd root
02308               }
02309               if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )  
02310               {
02311                 if( sd < stheta )
02312                 {
02313                   stheta    = sd;
02314                   sidetheta = kETheta;
02315                 }
02316               }            
02317             }
02318           }
02319         }
02320       }
02321     }
02322 
02323   } // end theta intersections
02324 
02325   // Phi Intersection
02326     
02327   if ( !fFullPhiSphere )
02328   {
02329     if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
02330     {
02331       // pDist -ve when inside
02332 
02333       pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
02334       pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
02335 
02336       // Comp -ve when in direction of outwards normal
02337 
02338       compS   = -sinSPhi*v.x()+cosSPhi*v.y() ;
02339       compE   =  sinEPhi*v.x()-cosEPhi*v.y() ;
02340       sidephi = kNull ;
02341 
02342       if ( (pDistS <= 0) && (pDistE <= 0) )
02343       {
02344         // Inside both phi *full* planes
02345 
02346         if ( compS < 0 )
02347         {
02348           sphi = pDistS/compS ;
02349           xi   = p.x()+sphi*v.x() ;
02350           yi   = p.y()+sphi*v.y() ;
02351 
02352           // Check intersection with correct half-plane (if not -> no intersect)
02353           //
02354           if( (std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance) )
02355           {
02356             vphi = std::atan2(v.y(),v.x());
02357             sidephi = kSPhi;
02358             if ( ( (fSPhi-halfAngTolerance) <= vphi)
02359               && ( (ePhi+halfAngTolerance)  >= vphi) )
02360             {
02361               sphi = kInfinity;
02362             }
02363           }
02364           else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
02365           {
02366             sphi=kInfinity;
02367           }
02368           else
02369           {
02370             sidephi = kSPhi ;
02371             if ( pDistS > -halfCarTolerance)  { sphi = 0; } // Leave by sphi 
02372           }
02373         }
02374         else  { sphi = kInfinity; }
02375 
02376         if ( compE < 0 )
02377         {
02378           sphi2=pDistE/compE ;
02379           if (sphi2 < sphi) // Only check further if < starting phi intersection
02380           {
02381             xi = p.x()+sphi2*v.x() ;
02382             yi = p.y()+sphi2*v.y() ;
02383 
02384             // Check intersection with correct half-plane
02385             //
02386             if ((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
02387             {
02388               // Leaving via ending phi
02389               //
02390               vphi = std::atan2(v.y(),v.x()) ;
02391                
02392               if( !((fSPhi-halfAngTolerance <= vphi)
02393                   &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
02394               { 
02395                 sidephi = kEPhi;
02396                 if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
02397                 else                                { sphi = 0.0;   }
02398               }
02399             }
02400             else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
02401             {
02402               sidephi = kEPhi ;
02403               if ( pDistE <= -halfCarTolerance )
02404               {
02405                 sphi=sphi2;
02406               }
02407               else 
02408               {
02409                 sphi = 0 ;
02410               }
02411             }
02412           }
02413         }        
02414       }
02415       else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
02416       {
02417         if ( pDistS <= pDistE )
02418         {
02419           sidephi = kSPhi ;
02420         }
02421         else
02422         {
02423           sidephi = kEPhi ;
02424         }
02425         if ( fDPhi > pi )
02426         {
02427           if ( (compS < 0) && (compE < 0) )  { sphi = 0; }
02428           else                               { sphi = kInfinity; }
02429         }
02430         else
02431         {
02432           // if towards both >=0 then once inside (after error)
02433           // will remain inside
02434 
02435           if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
02436           else                                { sphi = 0; }
02437         }    
02438       }
02439       else if ( (pDistS > 0) && (pDistE < 0) )
02440       {
02441         // Outside full starting plane, inside full ending plane
02442 
02443         if ( fDPhi > pi )
02444         {
02445           if ( compE < 0 )
02446           {
02447             sphi = pDistE/compE ;
02448             xi   = p.x() + sphi*v.x() ;
02449             yi   = p.y() + sphi*v.y() ;
02450 
02451             // Check intersection in correct half-plane
02452             // (if not -> not leaving phi extent)
02453             //
02454             if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
02455             {
02456               vphi = std::atan2(v.y(),v.x());
02457               sidephi = kSPhi;
02458               if ( ( (fSPhi-halfAngTolerance) <= vphi)
02459                 && ( (ePhi+halfAngTolerance)  >= vphi) )
02460               {
02461                 sphi = kInfinity;
02462               }
02463             }
02464             else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
02465             {
02466               sphi = kInfinity ;
02467             }
02468             else // Leaving via Ending phi
02469             {
02470               sidephi = kEPhi ;
02471               if ( pDistE > -halfCarTolerance )  { sphi = 0.; }
02472             }
02473           }
02474           else
02475           {
02476             sphi = kInfinity ;
02477           }
02478         }
02479         else
02480         {
02481           if ( compS >= 0 )
02482           {
02483             if ( compE < 0 )
02484             {            
02485               sphi = pDistE/compE ;
02486               xi   = p.x() + sphi*v.x() ;
02487               yi   = p.y() + sphi*v.y() ;
02488 
02489               // Check intersection in correct half-plane
02490               // (if not -> remain in extent)
02491               //
02492               if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
02493               {
02494                 vphi = std::atan2(v.y(),v.x());
02495                 sidephi = kSPhi;
02496                 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02497                   && ( (ePhi+halfAngTolerance)  >= vphi) )
02498                 {
02499                   sphi = kInfinity;
02500                 }
02501               }
02502               else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
02503               {
02504                 sphi=kInfinity;
02505               }
02506               else // otherwise leaving via Ending phi
02507               {
02508                 sidephi = kEPhi ;
02509               }
02510             }
02511             else sphi=kInfinity;
02512           }
02513           else // leaving immediately by starting phi
02514           {
02515             sidephi = kSPhi ;
02516             sphi    = 0 ;
02517           }
02518         }
02519       }
02520       else
02521       {
02522         // Must be pDistS < 0 && pDistE > 0
02523         // Inside full starting plane, outside full ending plane
02524 
02525         if ( fDPhi > pi )
02526         {
02527           if ( compS < 0 )
02528           {
02529             sphi=pDistS/compS;
02530             xi=p.x()+sphi*v.x();
02531             yi=p.y()+sphi*v.y();
02532             
02533             // Check intersection in correct half-plane
02534             // (if not -> not leaving phi extent)
02535             //
02536             if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
02537             {
02538               vphi = std::atan2(v.y(),v.x()) ;
02539               sidephi = kSPhi;
02540               if ( ( (fSPhi-halfAngTolerance) <= vphi)
02541                 && ( (ePhi+halfAngTolerance)  >= vphi) )
02542               {
02543               sphi = kInfinity;
02544               }
02545             }
02546             else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
02547             {
02548               sphi = kInfinity ;
02549             }
02550             else  // Leaving via Starting phi
02551             {
02552               sidephi = kSPhi ;
02553               if ( pDistS > -halfCarTolerance )  { sphi = 0; }
02554             }
02555           }
02556           else
02557           {
02558             sphi = kInfinity ;
02559           }
02560         }
02561         else
02562         {
02563           if ( compE >= 0 )
02564           {
02565             if ( compS < 0 )
02566             {
02567               sphi = pDistS/compS ;
02568               xi   = p.x()+sphi*v.x() ;
02569               yi   = p.y()+sphi*v.y() ;
02570               
02571               // Check intersection in correct half-plane
02572               // (if not -> remain in extent)
02573               //
02574               if((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
02575               {
02576                 vphi = std::atan2(v.y(),v.x()) ;
02577                 sidephi = kSPhi;
02578                 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02579                   && ( (ePhi+halfAngTolerance)  >= vphi) )
02580                 {
02581                   sphi = kInfinity;
02582                 }
02583               }
02584               else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
02585               {
02586                 sphi = kInfinity ;
02587               }
02588               else // otherwise leaving via Starting phi
02589               {
02590                 sidephi = kSPhi ;
02591               }
02592             }
02593             else
02594             {
02595               sphi = kInfinity ;
02596             }
02597           }
02598           else // leaving immediately by ending
02599           {
02600             sidephi = kEPhi ;
02601             sphi    = 0     ;
02602           }
02603         }
02604       }      
02605     }
02606     else
02607     {
02608       // On z axis + travel not || to z axis -> if phi of vector direction
02609       // within phi of shape, Step limited by rmax, else Step =0
02610 
02611       if ( v.x() || v.y() )
02612       {
02613         vphi = std::atan2(v.y(),v.x()) ;
02614         if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
02615         {
02616           sphi = kInfinity;
02617         }
02618         else
02619         {
02620           sidephi = kSPhi ; // arbitrary 
02621           sphi    = 0     ;
02622         }
02623       }
02624       else  // travel along z - no phi intersection
02625       {
02626         sphi = kInfinity ;
02627       }
02628     }
02629     if ( sphi < snxt )  // Order intersecttions
02630     {
02631       snxt = sphi ;
02632       side = sidephi ;
02633     }
02634   }
02635   if (stheta < snxt ) // Order intersections
02636   {
02637     snxt = stheta ;
02638     side = sidetheta ;
02639   }
02640 
02641   if (calcNorm)    // Output switch operator
02642   {
02643     switch( side )
02644     {
02645       case kRMax:
02646         xi=p.x()+snxt*v.x();
02647         yi=p.y()+snxt*v.y();
02648         zi=p.z()+snxt*v.z();
02649         *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
02650         *validNorm=true;
02651         break;
02652 
02653       case kRMin:
02654         *validNorm=false;  // Rmin is concave
02655         break;
02656 
02657       case kSPhi:
02658         if ( fDPhi <= pi )     // Normal to Phi-
02659         {
02660           *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
02661           *validNorm=true;
02662         }
02663         else  { *validNorm=false; }
02664         break ;
02665 
02666       case kEPhi:
02667         if ( fDPhi <= pi )      // Normal to Phi+
02668         {
02669           *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
02670           *validNorm=true;
02671         }
02672         else  { *validNorm=false; }
02673         break;
02674 
02675       case kSTheta:
02676         if( fSTheta == halfpi )
02677         {
02678           *n=G4ThreeVector(0.,0.,1.);
02679           *validNorm=true;
02680         }
02681         else if ( fSTheta > halfpi )
02682         {
02683           xi = p.x() + snxt*v.x();
02684           yi = p.y() + snxt*v.y();
02685           rho2=xi*xi+yi*yi;
02686           if (rho2)
02687           { 
02688             rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
02689             *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
02690                                -tanSTheta/std::sqrt(1+tanSTheta2));
02691           }
02692           else
02693           {
02694             *n = G4ThreeVector(0.,0.,1.);
02695           }
02696           *validNorm=true;
02697         }
02698         else  { *validNorm=false; }  // Concave STheta cone
02699         break;
02700 
02701       case kETheta:
02702         if( eTheta == halfpi )
02703         {
02704           *n         = G4ThreeVector(0.,0.,-1.);
02705           *validNorm = true;
02706         }
02707         else if ( eTheta < halfpi )
02708         {
02709           xi=p.x()+snxt*v.x();
02710           yi=p.y()+snxt*v.y();
02711           rho2=xi*xi+yi*yi;
02712           if (rho2)
02713           { 
02714             rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
02715             *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
02716                                -tanETheta/std::sqrt(1+tanETheta2) );
02717           }
02718           else
02719           {
02720             *n = G4ThreeVector(0.,0.,-1.);
02721           }
02722           *validNorm=true;
02723         }
02724         else  { *validNorm=false; }   // Concave ETheta cone
02725         break;
02726 
02727       default:
02728         G4cout << G4endl;
02729         DumpInfo();
02730         std::ostringstream message;
02731         G4int oldprc = message.precision(16);
02732         message << "Undefined side for valid surface normal to solid."
02733                 << G4endl
02734                 << "Position:"  << G4endl << G4endl
02735                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
02736                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
02737                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
02738                 << "Direction:" << G4endl << G4endl
02739                 << "v.x() = "   << v.x() << G4endl
02740                 << "v.y() = "   << v.y() << G4endl
02741                 << "v.z() = "   << v.z() << G4endl << G4endl
02742                 << "Proposed distance :" << G4endl << G4endl
02743                 << "snxt = "    << snxt/mm << " mm" << G4endl;
02744         message.precision(oldprc);
02745         G4Exception("G4Sphere::DistanceToOut(p,v,..)",
02746                     "GeomSolids1002", JustWarning, message);
02747         break;
02748     }
02749   }
02750   if (snxt == kInfinity)
02751   {
02752     G4cout << G4endl;
02753     DumpInfo();
02754     std::ostringstream message;
02755     G4int oldprc = message.precision(16);
02756     message << "Logic error: snxt = kInfinity  ???" << G4endl
02757             << "Position:"  << G4endl << G4endl
02758             << "p.x() = "   << p.x()/mm << " mm" << G4endl
02759             << "p.y() = "   << p.y()/mm << " mm" << G4endl
02760             << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
02761             << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
02762             << " mm" << G4endl << G4endl
02763             << "Direction:" << G4endl << G4endl
02764             << "v.x() = "   << v.x() << G4endl
02765             << "v.y() = "   << v.y() << G4endl
02766             << "v.z() = "   << v.z() << G4endl << G4endl
02767             << "Proposed distance :" << G4endl << G4endl
02768             << "snxt = "    << snxt/mm << " mm" << G4endl;
02769     message.precision(oldprc);
02770     G4Exception("G4Sphere::DistanceToOut(p,v,..)",
02771                 "GeomSolids1002", JustWarning, message);
02772   }
02773 
02774   return snxt;
02775 }

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

Reimplemented from G4VSolid.

Definition at line 309 of file G4Sphere.icc.

References G4CSGSolid::fCubicVolume.

00310 {
00311   if(fCubicVolume != 0.) {;}
00312   else { fCubicVolume = fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta))*
00313                               (fRmax*fRmax*fRmax-fRmin*fRmin*fRmin)/3.; }
00314   return fCubicVolume;
00315 }

G4double G4Sphere::GetDeltaPhiAngle (  )  const [inline]

Definition at line 62 of file G4Sphere.icc.

Referenced by GetDPhi(), G4tgbGeometryDumper::GetSolidParams(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

00063 {
00064   return fDPhi;
00065 }

G4double G4Sphere::GetDeltaThetaAngle (  )  const [inline]

Definition at line 73 of file G4Sphere.icc.

Referenced by GetDTheta(), G4tgbGeometryDumper::GetSolidParams(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

00074 {
00075   return fDTheta;
00076 }

G4double G4Sphere::GetDPhi (  )  const [inline]

Definition at line 291 of file G4Sphere.icc.

References GetDeltaPhiAngle().

00292 {
00293   return GetDeltaPhiAngle();
00294 }

G4double G4Sphere::GetDTheta (  )  const [inline]

Definition at line 303 of file G4Sphere.icc.

References GetDeltaThetaAngle().

00304 {
00305   return GetDeltaThetaAngle();
00306 }

G4GeometryType G4Sphere::GetEntityType (  )  const [virtual]

Implements G4VSolid.

Definition at line 3000 of file G4Sphere.cc.

03001 {
03002   return G4String("G4Sphere");
03003 }

G4VisExtent G4Sphere::GetExtent (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 3171 of file G4Sphere.cc.

03172 {
03173   return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
03174 }

G4double G4Sphere::GetInnerRadius (  )  const [inline]

Definition at line 44 of file G4Sphere.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams().

00045 {
00046   return fRmin;
00047 }

G4double G4Sphere::GetInsideRadius (  )  const [inline]

Definition at line 38 of file G4Sphere.icc.

Referenced by GetRmin(), G4PSSphereSurfaceFlux::IsSelectedSurface(), G4PSSphereSurfaceCurrent::IsSelectedSurface(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

00039 {
00040   return fRmin;
00041 }

G4double G4Sphere::GetOuterRadius (  )  const [inline]

Definition at line 50 of file G4Sphere.icc.

Referenced by GetRmax(), G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

00051 {
00052   return fRmax;
00053 }

G4ThreeVector G4Sphere::GetPointOnSurface (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 3042 of file G4Sphere.cc.

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

03043 {
03044   G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
03045   G4double height1, height2, slant1, slant2, costheta, sintheta, rRand;
03046 
03047   height1 = (fRmax-fRmin)*cosSTheta;
03048   height2 = (fRmax-fRmin)*cosETheta;
03049   slant1  = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
03050   slant2  = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
03051   rRand   = GetRadiusInRing(fRmin,fRmax);
03052   
03053   aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
03054   aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
03055   aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
03056   aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
03057   aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
03058   
03059   phi = RandFlat::shoot(fSPhi, ePhi); 
03060   cosphi = std::cos(phi); 
03061   sinphi = std::sin(phi);
03062   costheta = RandFlat::shoot(cosETheta,cosSTheta);
03063   sintheta = std::sqrt(1.-sqr(costheta));
03064 
03065   if(fFullPhiSphere) { aFiv = 0; }
03066   if(fSTheta == 0)   { aThr=0; }
03067   if(eTheta == pi) { aFou = 0; }
03068   if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
03069   if(eTheta == halfpi)  { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }
03070 
03071   chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
03072   if( (chose>=0.) && (chose<aOne) )
03073   {
03074     return G4ThreeVector(fRmax*sintheta*cosphi,
03075                          fRmax*sintheta*sinphi, fRmax*costheta);
03076   }
03077   else if( (chose>=aOne) && (chose<aOne+aTwo) )
03078   {
03079     return G4ThreeVector(fRmin*sintheta*cosphi,
03080                          fRmin*sintheta*sinphi, fRmin*costheta);
03081   }
03082   else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
03083   {
03084     if (fSTheta != halfpi)
03085     {
03086       zRand = RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
03087       return G4ThreeVector(tanSTheta*zRand*cosphi,
03088                            tanSTheta*zRand*sinphi,zRand);
03089     }
03090     else
03091     {
03092       return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
03093     }    
03094   }
03095   else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
03096   {
03097     if(eTheta != halfpi)
03098     {
03099       zRand = RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
03100       return G4ThreeVector  (tanETheta*zRand*cosphi,
03101                              tanETheta*zRand*sinphi,zRand);
03102     }
03103     else
03104     {
03105       return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
03106     }
03107   }
03108   else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
03109   {
03110     return G4ThreeVector(rRand*sintheta*cosSPhi,
03111                          rRand*sintheta*sinSPhi,rRand*costheta);
03112   }
03113   else
03114   {
03115     return G4ThreeVector(rRand*sintheta*cosEPhi,
03116                          rRand*sintheta*sinEPhi,rRand*costheta);
03117   }
03118 }

G4double G4Sphere::GetRmax (  )  const [inline]

Definition at line 279 of file G4Sphere.icc.

References GetOuterRadius().

00280 {
00281   return GetOuterRadius();
00282 }

G4double G4Sphere::GetRmin (  )  const [inline]

Definition at line 273 of file G4Sphere.icc.

References GetInsideRadius().

00274 {
00275   return GetInsideRadius();
00276 }

G4double G4Sphere::GetSPhi (  )  const [inline]

Definition at line 285 of file G4Sphere.icc.

References GetStartPhiAngle().

00286 {
00287   return GetStartPhiAngle();
00288 }

G4double G4Sphere::GetStartPhiAngle (  )  const [inline]

Definition at line 56 of file G4Sphere.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), GetSPhi(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

00057 {
00058   return fSPhi;
00059 }

G4double G4Sphere::GetStartThetaAngle (  )  const [inline]

Definition at line 68 of file G4Sphere.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), GetSTheta(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

00069 {
00070   return fSTheta;
00071 }

G4double G4Sphere::GetSTheta (  )  const [inline]

Definition at line 297 of file G4Sphere.icc.

References GetStartThetaAngle().

00298 {
00299   return GetStartThetaAngle();
00300 }

G4double G4Sphere::GetSurfaceArea (  )  [virtual]

Reimplemented from G4VSolid.

Definition at line 3124 of file G4Sphere.cc.

References G4CSGSolid::fSurfaceArea, and G4INCL::Math::pi.

03125 {
03126   if(fSurfaceArea != 0.) {;}
03127   else
03128   {   
03129     G4double Rsq=fRmax*fRmax;
03130     G4double rsq=fRmin*fRmin;
03131          
03132     fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
03133     if(!fFullPhiSphere)
03134     {
03135       fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
03136     }
03137     if(fSTheta >0)
03138     {
03139       G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
03140                               + std::pow(cosSTheta,2));
03141       if(fDPhi>pi)
03142       { 
03143         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
03144       }
03145       else
03146       {
03147         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
03148       }
03149     }
03150     if(eTheta < pi)
03151     {
03152       G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
03153                               + std::pow(cosETheta,2));
03154       if(fDPhi>pi)
03155       { 
03156         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
03157       }
03158       else
03159       {
03160         fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
03161       }
03162     }
03163   }
03164   return fSurfaceArea;
03165 }

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

Implements G4VSolid.

Definition at line 465 of file G4Sphere.cc.

References kInside, kOutside, and kSurface.

Referenced by CalculateExtent(), and DistanceToOut().

00466 {
00467   G4double rho,rho2,rad2,tolRMin,tolRMax;
00468   G4double pPhi,pTheta;
00469   EInside in = kOutside;
00470   static const G4double halfAngTolerance = kAngTolerance*0.5;
00471   const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
00472   const G4double halfRminTolerance = fRminTolerance*0.5;
00473   const G4double Rmax_minus = fRmax - halfRmaxTolerance;
00474   const G4double Rmin_plus  = (fRmin > 0) ? fRmin+halfRminTolerance : 0;
00475 
00476   rho2 = p.x()*p.x() + p.y()*p.y() ;
00477   rad2 = rho2 + p.z()*p.z() ;
00478 
00479   // Check radial surfaces. Sets 'in'
00480 
00481   tolRMin = Rmin_plus;
00482   tolRMax = Rmax_minus;
00483 
00484   if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
00485   {
00486     in = kInside;
00487   }
00488   else
00489   {
00490     tolRMax = fRmax + halfRmaxTolerance;                  // outside case
00491     tolRMin = std::max(fRmin-halfRminTolerance, 0.);      // outside case
00492     if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
00493     {
00494       in = kSurface;
00495     }
00496     else
00497     {
00498       return in = kOutside;
00499     }
00500   }
00501 
00502   // Phi boundaries   : Do not check if it has no phi boundary!
00503 
00504   if ( !fFullPhiSphere && rho2 )  // [fDPhi < twopi] and [p.x or p.y]
00505   {
00506     pPhi = std::atan2(p.y(),p.x()) ;
00507 
00508     if      ( pPhi < fSPhi - halfAngTolerance  ) { pPhi += twopi; }
00509     else if ( pPhi > ePhi + halfAngTolerance )   { pPhi -= twopi; }
00510     
00511     if ( (pPhi < fSPhi - halfAngTolerance)
00512       || (pPhi > ePhi + halfAngTolerance) )      { return in = kOutside; }
00513     
00514     else if (in == kInside)  // else it's kSurface anyway already
00515     {
00516       if ( (pPhi < fSPhi + halfAngTolerance)
00517         || (pPhi > ePhi - halfAngTolerance) )    { in = kSurface; }     
00518     }
00519   }
00520 
00521   // Theta bondaries
00522   
00523   if ( (rho2 || p.z()) && (!fFullThetaSphere) )
00524   {
00525     rho    = std::sqrt(rho2);
00526     pTheta = std::atan2(rho,p.z());
00527 
00528     if ( in == kInside )
00529     {
00530       if ( (pTheta < fSTheta + halfAngTolerance)
00531         || (pTheta > eTheta - halfAngTolerance) )
00532       {
00533         if ( (pTheta >= fSTheta - halfAngTolerance)
00534           && (pTheta <= eTheta + halfAngTolerance) )
00535         {
00536           in = kSurface;
00537         }
00538         else
00539         {
00540           in = kOutside;
00541         }
00542       }
00543     }
00544     else
00545     {
00546       if ( (pTheta < fSTheta - halfAngTolerance)
00547         || (pTheta > eTheta + halfAngTolerance) )
00548       {
00549         in = kOutside;
00550       }
00551     }
00552   }
00553   return in;
00554 }

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

Definition at line 173 of file G4Sphere.cc.

References cosCPhi, cosEPhi, cosETheta, cosHDPhiIT, cosHDPhiOT, cosSPhi, cosSTheta, cPhi, ePhi, eTheta, fDPhi, fDTheta, fEpsilon, fFullPhiSphere, fFullSphere, fFullThetaSphere, fRmax, fRmaxTolerance, fRmin, fRminTolerance, fSPhi, fSTheta, hDPhi, kAngTolerance, kRadTolerance, G4CSGSolid::operator=(), sinCPhi, sinEPhi, sinETheta, sinSPhi, sinSTheta, tanETheta, tanETheta2, tanSTheta, and tanSTheta2.

00174 {
00175    // Check assignment to self
00176    //
00177    if (this == &rhs)  { return *this; }
00178 
00179    // Copy base class data
00180    //
00181    G4CSGSolid::operator=(rhs);
00182 
00183    // Copy data
00184    //
00185    fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
00186    kAngTolerance = rhs.kAngTolerance; kRadTolerance = rhs.kRadTolerance;
00187    fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; fRmax = rhs.fRmax;
00188    fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSTheta = rhs.fSTheta;
00189    fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
00190    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
00191    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
00192    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
00193    hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = rhs.ePhi;
00194    sinSTheta = rhs.sinSTheta; cosSTheta = rhs.cosSTheta;
00195    sinETheta = rhs.sinETheta; cosETheta = rhs.cosETheta;
00196    tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs.tanSTheta2;
00197    tanETheta = rhs.tanETheta; tanETheta2 = rhs.tanETheta2;
00198    eTheta = rhs.eTheta; fFullPhiSphere = rhs.fFullPhiSphere;
00199    fFullThetaSphere = rhs.fFullThetaSphere; fFullSphere = rhs.fFullSphere;
00200 
00201    return *this;
00202 }

void G4Sphere::SetDeltaPhiAngle ( G4double  newDphi  )  [inline]

Definition at line 250 of file G4Sphere.icc.

00251 {
00252   CheckPhiAngles(fSPhi, newDPhi);
00253   Initialize();
00254 }

void G4Sphere::SetDeltaThetaAngle ( G4double  newDTheta  )  [inline]

Definition at line 264 of file G4Sphere.icc.

00265 {
00266   CheckThetaAngles(fSTheta, newDTheta);
00267   Initialize();
00268 }

void G4Sphere::SetInnerRadius ( G4double  newRMin  )  [inline]

Definition at line 224 of file G4Sphere.icc.

References SetInsideRadius().

00225 {
00226   SetInsideRadius(newRmin);
00227 }

void G4Sphere::SetInsideRadius ( G4double  newRmin  )  [inline]

Definition at line 216 of file G4Sphere.icc.

Referenced by SetInnerRadius().

00217 {
00218   fRmin= newRmin;
00219   fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
00220   Initialize();
00221 }

void G4Sphere::SetOuterRadius ( G4double  newRmax  )  [inline]

Definition at line 230 of file G4Sphere.icc.

00231 {
00232   fRmax= newRmax;
00233   fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
00234   Initialize();
00235 }

void G4Sphere::SetStartPhiAngle ( G4double  newSphi,
G4bool  trig = true 
) [inline]

Definition at line 238 of file G4Sphere.icc.

00239 {
00240   // Flag 'compute' can be used to explicitely avoid recomputation of
00241   // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
00242 
00243   CheckSPhiAngle(newSPhi);
00244   fFullPhiSphere = false;
00245   if (compute)  { InitializePhiTrigonometry(); }
00246   Initialize();
00247 }

void G4Sphere::SetStartThetaAngle ( G4double  newSTheta  )  [inline]

Definition at line 257 of file G4Sphere.icc.

00258 {
00259   CheckThetaAngles(newSTheta, fDTheta);
00260   Initialize();
00261 }

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

Reimplemented from G4CSGSolid.

Definition at line 3018 of file G4Sphere.cc.

References G4VSolid::GetName().

03019 {
03020   G4int oldprc = os.precision(16);
03021   os << "-----------------------------------------------------------\n"
03022      << "    *** Dump for solid - " << GetName() << " ***\n"
03023      << "    ===================================================\n"
03024      << " Solid type: G4Sphere\n"
03025      << " Parameters: \n"
03026      << "    inner radius: " << fRmin/mm << " mm \n"
03027      << "    outer radius: " << fRmax/mm << " mm \n"
03028      << "    starting phi of segment  : " << fSPhi/degree << " degrees \n"
03029      << "    delta phi of segment     : " << fDPhi/degree << " degrees \n"
03030      << "    starting theta of segment: " << fSTheta/degree << " degrees \n"
03031      << "    delta theta of segment   : " << fDTheta/degree << " degrees \n"
03032      << "-----------------------------------------------------------\n";
03033   os.precision(oldprc);
03034 
03035   return os;
03036 }

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

Implements G4VSolid.

Definition at line 562 of file G4Sphere.cc.

References G4Exception(), JustWarning, G4VSolid::kCarTolerance, and G4INCL::Math::pi.

00563 {
00564   G4int noSurfaces = 0;  
00565   G4double rho, rho2, radius, pTheta, pPhi=0.;
00566   G4double distRMin = kInfinity;
00567   G4double distSPhi = kInfinity, distEPhi = kInfinity;
00568   G4double distSTheta = kInfinity, distETheta = kInfinity;
00569   G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
00570   G4ThreeVector norm, sumnorm(0.,0.,0.);
00571 
00572   static const G4double halfCarTolerance = 0.5*kCarTolerance;
00573   static const G4double halfAngTolerance = 0.5*kAngTolerance;
00574 
00575   rho2 = p.x()*p.x()+p.y()*p.y();
00576   radius = std::sqrt(rho2+p.z()*p.z());
00577   rho  = std::sqrt(rho2);
00578 
00579   G4double    distRMax = std::fabs(radius-fRmax);
00580   if (fRmin)  distRMin = std::fabs(radius-fRmin);
00581     
00582   if ( rho && !fFullSphere )
00583   {
00584     pPhi = std::atan2(p.y(),p.x());
00585 
00586     if (pPhi < fSPhi-halfAngTolerance)     { pPhi += twopi; }
00587     else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
00588   }
00589   if ( !fFullPhiSphere )
00590   {
00591     if ( rho )
00592     {
00593       distSPhi = std::fabs( pPhi-fSPhi ); 
00594       distEPhi = std::fabs( pPhi-ePhi ); 
00595     }
00596     else if( !fRmin )
00597     {
00598       distSPhi = 0.; 
00599       distEPhi = 0.; 
00600     }
00601     nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
00602     nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
00603   }        
00604   if ( !fFullThetaSphere )
00605   {
00606     if ( rho )
00607     {
00608       pTheta     = std::atan2(rho,p.z());
00609       distSTheta = std::fabs(pTheta-fSTheta); 
00610       distETheta = std::fabs(pTheta-eTheta);
00611  
00612       nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
00613                           -cosSTheta*p.y()/rho,
00614                            sinSTheta          );
00615 
00616       nTe = G4ThreeVector( cosETheta*p.x()/rho,
00617                            cosETheta*p.y()/rho,
00618                           -sinETheta          );    
00619     }
00620     else if( !fRmin )
00621     {
00622       if ( fSTheta )  
00623       {              
00624         distSTheta = 0.;
00625         nTs = G4ThreeVector(0.,0.,-1.);
00626       }
00627       if ( eTheta < pi )
00628       {              
00629         distETheta = 0.;
00630         nTe = G4ThreeVector(0.,0.,1.);
00631       }
00632     }    
00633   }
00634   if( radius )  { nR = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius); }
00635 
00636   if( distRMax <= halfCarTolerance )
00637   {
00638     noSurfaces ++;
00639     sumnorm += nR;
00640   }
00641   if( fRmin && (distRMin <= halfCarTolerance) )
00642   {
00643     noSurfaces ++;
00644     sumnorm -= nR;
00645   }
00646   if( !fFullPhiSphere )   
00647   {
00648     if (distSPhi <= halfAngTolerance)
00649     {
00650       noSurfaces ++;
00651       sumnorm += nPs;
00652     }
00653     if (distEPhi <= halfAngTolerance) 
00654     {
00655       noSurfaces ++;
00656       sumnorm += nPe;
00657     }
00658   }
00659   if ( !fFullThetaSphere )
00660   {
00661     if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
00662     {
00663       noSurfaces ++;
00664       if ((radius <= halfCarTolerance) && fFullPhiSphere)  { sumnorm += nZ;  }
00665       else                                                 { sumnorm += nTs; }
00666     }
00667     if ((distETheta <= halfAngTolerance) && (eTheta < pi)) 
00668     {
00669       noSurfaces ++;
00670       if ((radius <= halfCarTolerance) && fFullPhiSphere)  { sumnorm -= nZ;  }
00671       else                                                 { sumnorm += nTe; }
00672       if(sumnorm.z() == 0.)  { sumnorm += nZ; }
00673     }
00674   }
00675   if ( noSurfaces == 0 )
00676   {
00677 #ifdef G4CSGDEBUG
00678     G4Exception("G4Sphere::SurfaceNormal(p)", "GeomSolids1002",
00679                 JustWarning, "Point p is not on surface !?" ); 
00680 #endif
00681      norm = ApproxSurfaceNormal(p);
00682   }
00683   else if ( noSurfaces == 1 )  { norm = sumnorm; }
00684   else                         { norm = sumnorm.unit(); }
00685   return norm;
00686 }


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