G4Cons Class Reference

#include <G4Cons.hh>

Inheritance diagram for G4Cons:

G4CSGSolid G4VSolid

Public Member Functions

 G4Cons (const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
 ~G4Cons ()
G4double GetInnerRadiusMinusZ () const
G4double GetOuterRadiusMinusZ () const
G4double GetInnerRadiusPlusZ () const
G4double GetOuterRadiusPlusZ () const
G4double GetZHalfLength () const
G4double GetStartPhiAngle () const
G4double GetDeltaPhiAngle () const
void SetInnerRadiusMinusZ (G4double Rmin1)
void SetOuterRadiusMinusZ (G4double Rmax1)
void SetInnerRadiusPlusZ (G4double Rmin2)
void SetOuterRadiusPlusZ (G4double Rmax2)
void SetZHalfLength (G4double newDz)
void SetStartPhiAngle (G4double newSPhi, G4bool trig=true)
void SetDeltaPhiAngle (G4double newDPhi)
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
void DescribeYourselfTo (G4VGraphicsScene &scene) const
G4PolyhedronCreatePolyhedron () const
G4NURBSCreateNURBS () const
 G4Cons (__void__ &)
 G4Cons (const G4Cons &rhs)
G4Consoperator= (const G4Cons &rhs)
G4double GetRmin1 () const
G4double GetRmax1 () const
G4double GetRmin2 () const
G4double GetRmax2 () const
G4double GetDz () const
G4double GetSPhi () const
G4double GetDPhi () const

Detailed Description

Definition at line 74 of file G4Cons.hh.


Constructor & Destructor Documentation

G4Cons::G4Cons ( const G4String pName,
G4double  pRmin1,
G4double  pRmax1,
G4double  pRmin2,
G4double  pRmax2,
G4double  pDz,
G4double  pSPhi,
G4double  pDPhi 
)

Definition at line 79 of file G4Cons.cc.

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

Referenced by Clone().

00084   : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
00085     fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
00086 {
00087   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00088   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00089 
00090   // Check z-len
00091   //
00092   if ( pDz < 0 )
00093   {
00094     std::ostringstream message;
00095     message << "Invalid Z half-length for Solid: " << GetName() << G4endl
00096             << "        hZ = " << pDz;
00097     G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
00098                 FatalException, message);
00099   }
00100 
00101   // Check radii
00102   //
00103   if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
00104   {
00105     std::ostringstream message;
00106     message << "Invalid values of radii for Solid: " << GetName() << G4endl
00107             << "        pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
00108             << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
00109     G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
00110                 FatalException, message) ;
00111   }
00112   if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
00113   if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
00114 
00115   // Check angles
00116   //
00117   CheckPhiAngles(pSPhi, pDPhi);
00118 }

G4Cons::~G4Cons (  ) 

Definition at line 138 of file G4Cons.cc.

00139 {
00140 }

G4Cons::G4Cons ( __void__ &   ) 

Definition at line 125 of file G4Cons.cc.

00126   : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
00127     fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
00128     fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
00129     cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
00130     fPhiFullCone(false)
00131 {
00132 }

G4Cons::G4Cons ( const G4Cons rhs  ) 

Definition at line 146 of file G4Cons.cc.

00147   : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
00148     kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
00149     fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
00150     fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
00151     cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
00152     sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
00153     cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone)
00154 {
00155 }


Member Function Documentation

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

Implements G4VSolid.

Definition at line 261 of file G4Cons.cc.

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

00266 {
00267   if ( !pTransform.IsRotated() && (fDPhi == twopi)
00268     && (fRmin1 == 0) && (fRmin2 == 0) )
00269   {
00270     // Special case handling for unrotated solid cones
00271     // Compute z/x/y mins and maxs for bounding box respecting limits,
00272     // with early returns if outside limits. Then switch() on pAxis,
00273     // and compute exact x and y limit for x/y case
00274       
00275     G4double xoffset, xMin, xMax ;
00276     G4double yoffset, yMin, yMax ;
00277     G4double zoffset, zMin, zMax ;
00278 
00279     G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
00280     G4double xoff1, xoff2, yoff1, yoff2 ;
00281       
00282     zoffset = pTransform.NetTranslation().z();
00283     zMin    = zoffset - fDz ;
00284     zMax    = zoffset + fDz ;
00285 
00286     if (pVoxelLimit.IsZLimited())
00287     {
00288       if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) || 
00289           (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance)  )
00290       {
00291         return false ;
00292       }
00293       else
00294       {
00295         if ( zMin < pVoxelLimit.GetMinZExtent() )
00296         {
00297           zMin = pVoxelLimit.GetMinZExtent() ;
00298         }
00299         if ( zMax > pVoxelLimit.GetMaxZExtent() )
00300         {
00301           zMax = pVoxelLimit.GetMaxZExtent() ;
00302         }
00303       }
00304     }
00305     xoffset = pTransform.NetTranslation().x() ;
00306     RMax    = (fRmax2 >= fRmax1) ?  zMax : zMin  ;                  
00307     xMax    = xoffset + (fRmax1 + fRmax2)*0.5 + 
00308               (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
00309     xMin    = 2*xoffset-xMax ;
00310 
00311     if (pVoxelLimit.IsXLimited())
00312     {
00313       if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) || 
00314            (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance)  )
00315       {
00316         return false ;
00317       }
00318       else
00319       {
00320         if ( xMin < pVoxelLimit.GetMinXExtent() )
00321         {
00322           xMin = pVoxelLimit.GetMinXExtent() ;
00323         }
00324         if ( xMax > pVoxelLimit.GetMaxXExtent() )
00325         {
00326           xMax=pVoxelLimit.GetMaxXExtent() ;
00327         }
00328       }
00329     }
00330     yoffset = pTransform.NetTranslation().y() ;
00331     yMax    = yoffset + (fRmax1 + fRmax2)*0.5 + 
00332               (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
00333     yMin    = 2*yoffset-yMax ;
00334     RMax    = yMax - yoffset ;  // = max radius due to Zmax/Zmin cuttings
00335 
00336     if (pVoxelLimit.IsYLimited())
00337     {
00338       if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) || 
00339            (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance)  )
00340       {
00341         return false ;
00342       }
00343       else
00344       {
00345         if ( yMin < pVoxelLimit.GetMinYExtent() )
00346         {
00347           yMin = pVoxelLimit.GetMinYExtent() ;
00348         }
00349         if ( yMax > pVoxelLimit.GetMaxYExtent() )
00350         {
00351           yMax = pVoxelLimit.GetMaxYExtent() ;
00352         }
00353       }
00354     }    
00355     switch (pAxis) // Known to cut cones
00356     {
00357       case kXAxis:
00358         yoff1 = yoffset - yMin ;
00359         yoff2 = yMax - yoffset ;
00360 
00361         if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
00362         {                                 // => no change
00363           pMin = xMin ;
00364           pMax = xMax ;
00365         }
00366         else
00367         {
00368           // Y limits don't cross max/min x => compute max delta x,
00369           // hence new mins/maxs
00370           delta=RMax*RMax-yoff1*yoff1;
00371           diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00372           delta=RMax*RMax-yoff2*yoff2;
00373           diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00374           maxDiff = (diff1>diff2) ? diff1:diff2 ;
00375           newMin  = xoffset - maxDiff ;
00376           newMax  = xoffset + maxDiff ;
00377           pMin    = ( newMin < xMin ) ? xMin : newMin  ;
00378           pMax    = ( newMax > xMax) ? xMax : newMax ;
00379         } 
00380       break ;
00381 
00382       case kYAxis:
00383         xoff1 = xoffset - xMin ;
00384         xoff2 = xMax - xoffset ;
00385 
00386         if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
00387         {                                  // => no change
00388           pMin = yMin ;
00389           pMax = yMax ;
00390         }
00391         else
00392         {
00393           // X limits don't cross max/min y => compute max delta y,
00394           // hence new mins/maxs
00395           delta=RMax*RMax-xoff1*xoff1;
00396           diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00397           delta=RMax*RMax-xoff2*xoff2;
00398           diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00399           maxDiff = (diff1 > diff2) ? diff1:diff2 ;
00400           newMin  = yoffset - maxDiff ;
00401           newMax  = yoffset + maxDiff ;
00402           pMin    = (newMin < yMin) ? yMin : newMin ;
00403           pMax    = (newMax > yMax) ? yMax : newMax ;
00404         }
00405       break ;
00406 
00407       case kZAxis:
00408         pMin = zMin ;
00409         pMax = zMax ;
00410       break ;
00411       
00412       default:
00413       break ;
00414     }
00415     pMin -= kCarTolerance ;
00416     pMax += kCarTolerance ;
00417 
00418     return true ;
00419   }
00420   else   // Calculate rotated vertex coordinates
00421   {
00422     G4int i, noEntries, noBetweenSections4 ;
00423     G4bool existsAfterClip = false ;
00424     G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
00425 
00426     pMin = +kInfinity ;
00427     pMax = -kInfinity ;
00428 
00429     noEntries          = vertices->size() ;
00430     noBetweenSections4 = noEntries-4 ;
00431       
00432     for ( i = 0 ; i < noEntries ; i += 4 )
00433     {
00434       ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
00435     }
00436     for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
00437     {
00438       ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
00439     }    
00440     if ( (pMin != kInfinity) || (pMax != -kInfinity) )
00441     {
00442       existsAfterClip = true ;
00443         
00444       // Add 2*tolerance to avoid precision troubles
00445 
00446       pMin -= kCarTolerance ;
00447       pMax += kCarTolerance ;
00448     }
00449     else
00450     {
00451       // Check for case where completely enveloping clipping volume
00452       // If point inside then we are confident that the solid completely
00453       // envelopes the clipping volume. Hence set min/max extents according
00454       // to clipping volume extents along the specified axis.
00455        
00456       G4ThreeVector clipCentre(
00457       (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
00458       (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
00459       (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5  ) ;
00460         
00461       if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
00462       {
00463         existsAfterClip = true ;
00464         pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
00465         pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
00466       }
00467     }
00468     delete vertices ;
00469     return existsAfterClip ;
00470   }
00471 }

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

Reimplemented from G4VSolid.

Definition at line 2254 of file G4Cons.cc.

References G4Cons().

02255 {
02256   return new G4Cons(*this);
02257 }

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

Reimplemented from G4VSolid.

Definition at line 249 of file G4Cons.cc.

References G4VPVParameterisation::ComputeDimensions().

00252 {
00253   p->ComputeDimensions(*this,n,pRep) ;
00254 }

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

Reimplemented from G4VSolid.

Definition at line 2387 of file G4Cons.cc.

02388 {
02389   G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;
02390   return new G4NURBSbox (RMax, RMax, fDz);       // Box for now!!!
02391 }

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

Reimplemented from G4VSolid.

Definition at line 2382 of file G4Cons.cc.

02383 {
02384   return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
02385 }

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

Implements G4VSolid.

Definition at line 2377 of file G4Cons.cc.

References G4VGraphicsScene::AddSolid().

02378 {
02379   scene.AddSolid (*this);
02380 }

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

Implements G4VSolid.

Definition at line 1384 of file G4Cons.cc.

01385 {
01386   G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
01387   G4double tanRMin, secRMin, pRMin ;
01388   G4double tanRMax, secRMax, pRMax ;
01389 
01390   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
01391   safeZ = std::fabs(p.z()) - fDz ;
01392 
01393   if ( fRmin1 || fRmin2 )
01394   {
01395     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
01396     secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
01397     pRMin   = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
01398     safeR1  = (pRMin - rho)/secRMin ;
01399 
01400     tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
01401     secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
01402     pRMax   = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
01403     safeR2  = (rho - pRMax)/secRMax ;
01404 
01405     if ( safeR1 > safeR2) { safe = safeR1; }
01406     else                  { safe = safeR2; }
01407   }
01408   else
01409   {
01410     tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
01411     secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
01412     pRMax   = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
01413     safe    = (rho - pRMax)/secRMax ;
01414   }
01415   if ( safeZ > safe )  { safe = safeZ; }
01416 
01417   if ( !fPhiFullCone && rho )
01418   {
01419     // Psi=angle from central phi to point
01420 
01421     cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
01422 
01423     if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
01424     {
01425       if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
01426       {
01427         safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
01428       }
01429       else
01430       {
01431         safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
01432       }
01433       if ( safePhi > safe )  { safe = safePhi; }
01434     }
01435   }
01436   if ( safe < 0.0 )  { safe = 0.0; }
01437 
01438   return safe ;
01439 }

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

Implements G4VSolid.

Definition at line 712 of file G4Cons.cc.

References G4VSolid::kCarTolerance.

00714 {
00715   G4double snxt = kInfinity ;      // snxt = default return value
00716   const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2.
00717   static const G4double halfCarTolerance=kCarTolerance*0.5;
00718   static const G4double halfRadTolerance=kRadTolerance*0.5;
00719 
00720   G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;  // Data for cones
00721   G4double tanRMin,secRMin,rMinAv,rMinOAv ;
00722   G4double rout,rin ;
00723 
00724   G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared
00725   G4double tolORMax2,tolIRMax,tolIRMax2 ;
00726   G4double tolODz,tolIDz ;
00727 
00728   G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars
00729 
00730   G4double t1,t2,t3,b,c,d ;    // Quadratic solver variables 
00731   G4double nt1,nt2,nt3 ;
00732   G4double Comp ;
00733 
00734   G4ThreeVector Normal;
00735 
00736   // Cone Precalcs
00737 
00738   tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
00739   secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
00740   rMinAv  = (fRmin1 + fRmin2)*0.5 ;
00741 
00742   if (rMinAv > halfRadTolerance)
00743   {
00744     rMinOAv = rMinAv - halfRadTolerance ;
00745   }
00746   else
00747   {
00748     rMinOAv = 0.0 ;
00749   }  
00750   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
00751   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
00752   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
00753   rMaxOAv = rMaxAv + halfRadTolerance ;
00754    
00755   // Intersection with z-surfaces
00756 
00757   tolIDz = fDz - halfCarTolerance ;
00758   tolODz = fDz + halfCarTolerance ;
00759 
00760   if (std::fabs(p.z()) >= tolIDz)
00761   {
00762     if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
00763     {
00764       sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
00765 
00766       if( sd < 0.0 )  { sd = 0.0; }                    // negative dist -> zero
00767 
00768       xi   = p.x() + sd*v.x() ;  // Intersection coords
00769       yi   = p.y() + sd*v.y() ;
00770       rhoi2 = xi*xi + yi*yi  ;
00771 
00772       // Check validity of intersection
00773       // Calculate (outer) tolerant radi^2 at intersecion
00774 
00775       if (v.z() > 0)
00776       {
00777         tolORMin  = fRmin1 - halfRadTolerance*secRMin ;
00778         tolIRMin  = fRmin1 + halfRadTolerance*secRMin ;
00779         tolIRMax  = fRmax1 - halfRadTolerance*secRMin ;
00780         tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
00781                     (fRmax1 + halfRadTolerance*secRMax) ;
00782       }
00783       else
00784       {
00785         tolORMin  = fRmin2 - halfRadTolerance*secRMin ;
00786         tolIRMin  = fRmin2 + halfRadTolerance*secRMin ;
00787         tolIRMax  = fRmax2 - halfRadTolerance*secRMin ;
00788         tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
00789                     (fRmax2 + halfRadTolerance*secRMax) ;
00790       }
00791       if ( tolORMin > 0 ) 
00792       {
00793         tolORMin2 = tolORMin*tolORMin ;
00794         tolIRMin2 = tolIRMin*tolIRMin ;
00795       }
00796       else                
00797       {
00798         tolORMin2 = 0.0 ;
00799         tolIRMin2 = 0.0 ;
00800       }
00801       if ( tolIRMax > 0 )  { tolIRMax2 = tolIRMax*tolIRMax; }     
00802       else                 { tolIRMax2 = 0.0; }
00803       
00804       if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
00805       {
00806         if ( !fPhiFullCone && rhoi2 )
00807         {
00808           // Psi = angle made with central (average) phi of shape
00809 
00810           cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
00811 
00812           if (cosPsi >= cosHDPhiIT)  { return sd; }
00813         }
00814         else
00815         {
00816           return sd;
00817         }
00818       }
00819     }
00820     else  // On/outside extent, and heading away  -> cannot intersect
00821     {
00822       return snxt ;  
00823     }
00824   }
00825     
00826 // ----> Can not intersect z surfaces
00827 
00828 
00829 // Intersection with outer cone (possible return) and
00830 //                   inner cone (must also check phi)
00831 //
00832 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
00833 //
00834 // Intersects with x^2+y^2=(a*z+b)^2
00835 //
00836 // where a=tanRMax or tanRMin
00837 //       b=rMaxAv  or rMinAv
00838 //
00839 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
00840 //     t1                        t2                      t3  
00841 //
00842 //  \--------u-------/       \-----------v----------/ \---------w--------/
00843 //
00844 
00845   t1   = 1.0 - v.z()*v.z() ;
00846   t2   = p.x()*v.x() + p.y()*v.y() ;
00847   t3   = p.x()*p.x() + p.y()*p.y() ;
00848   rin  = tanRMin*p.z() + rMinAv ;
00849   rout = tanRMax*p.z() + rMaxAv ;
00850 
00851   // Outer Cone Intersection
00852   // Must be outside/on outer cone for valid intersection
00853 
00854   nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
00855   nt2 = t2 - tanRMax*v.z()*rout ;
00856   nt3 = t3 - rout*rout ;
00857 
00858   if (std::fabs(nt1) > kRadTolerance)  // Equation quadratic => 2 roots
00859   {
00860     b = nt2/nt1;
00861     c = nt3/nt1;
00862     d = b*b-c  ;
00863     if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
00864       || (rout < 0) )
00865     {
00866       // If outside real cone (should be rho-rout>kRadTolerance*0.5
00867       // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
00868 
00869       if (d >= 0)
00870       {
00871           
00872         if ((rout < 0) && (nt3 <= 0))
00873         {
00874           // Inside `shadow cone' with -ve radius
00875           // -> 2nd root could be on real cone
00876 
00877           if (b>0) { sd = c/(-b-std::sqrt(d)); }
00878           else     { sd = -b + std::sqrt(d);   }
00879         }
00880         else
00881         {
00882           if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
00883           {
00884             sd=c/(-b+std::sqrt(d));
00885           }
00886           else
00887           {
00888             if ( c <= 0 ) // second >=0
00889             {
00890               sd = -b + std::sqrt(d) ;
00891               if((sd<0) & (sd>-halfRadTolerance)) sd=0;
00892             }
00893             else  // both negative, travel away
00894             {
00895               return kInfinity ;
00896             }
00897           }
00898         }
00899         if ( sd > 0 )  // If 'forwards'. Check z intersection
00900         {
00901           if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
00902           {               // 64 bits systems. Split long distances and recompute
00903             G4double fTerm = sd-std::fmod(sd,dRmax);
00904             sd = fTerm + DistanceToIn(p+fTerm*v,v);
00905           } 
00906           zi = p.z() + sd*v.z() ;
00907 
00908           if (std::fabs(zi) <= tolODz)
00909           {
00910             // Z ok. Check phi intersection if reqd
00911 
00912             if ( fPhiFullCone )  { return sd; }
00913             else
00914             {
00915               xi     = p.x() + sd*v.x() ;
00916               yi     = p.y() + sd*v.y() ;
00917               ri     = rMaxAv + zi*tanRMax ;
00918               cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
00919 
00920               if ( cosPsi >= cosHDPhiIT )  { return sd; }
00921             }
00922           }
00923         }                // end if (sd>0)
00924       }
00925     }
00926     else
00927     {
00928       // Inside outer cone
00929       // check not inside, and heading through G4Cons (-> 0 to in)
00930 
00931       if ( ( t3  > (rin + halfRadTolerance*secRMin)*
00932                    (rin + halfRadTolerance*secRMin) )
00933         && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
00934       {
00935         // Inside cones, delta r -ve, inside z extent
00936         // Point is on the Surface => check Direction using  Normal.dot(v)
00937 
00938         xi     = p.x() ;
00939         yi     = p.y()  ;
00940         risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
00941         Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
00942         if ( !fPhiFullCone )
00943         {
00944           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
00945           if ( cosPsi >= cosHDPhiIT )
00946           {
00947             if ( Normal.dot(v) <= 0 )  { return 0.0; }
00948           }
00949         }
00950         else
00951         {             
00952           if ( Normal.dot(v) <= 0 )  { return 0.0; }
00953         }
00954       }
00955     }
00956   }
00957   else  //  Single root case 
00958   {
00959     if ( std::fabs(nt2) > kRadTolerance )
00960     {
00961       sd = -0.5*nt3/nt2 ;
00962 
00963       if ( sd < 0 )  { return kInfinity; }   // travel away
00964       else  // sd >= 0,  If 'forwards'. Check z intersection
00965       {
00966         zi = p.z() + sd*v.z() ;
00967 
00968         if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
00969         {
00970           // Z ok. Check phi intersection if reqd
00971 
00972           if ( fPhiFullCone )  { return sd; }
00973           else
00974           {
00975             xi     = p.x() + sd*v.x() ;
00976             yi     = p.y() + sd*v.y() ;
00977             ri     = rMaxAv + zi*tanRMax ;
00978             cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
00979 
00980             if (cosPsi >= cosHDPhiIT)  { return sd; }
00981           }
00982         }
00983       }
00984     }
00985     else  //    travel || cone surface from its origin
00986     {
00987       sd = kInfinity ;
00988     }
00989   }
00990 
00991   // Inner Cone Intersection
00992   // o Space is divided into 3 areas:
00993   //   1) Radius greater than real inner cone & imaginary cone & outside
00994   //      tolerance
00995   //   2) Radius less than inner or imaginary cone & outside tolarance
00996   //   3) Within tolerance of real or imaginary cones
00997   //      - Extra checks needed for 3's intersections
00998   //        => lots of duplicated code
00999 
01000   if (rMinAv)
01001   { 
01002     nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
01003     nt2 = t2 - tanRMin*v.z()*rin ;
01004     nt3 = t3 - rin*rin ;
01005  
01006     if ( nt1 )
01007     {
01008       if ( nt3 > rin*kRadTolerance*secRMin )
01009       {
01010         // At radius greater than real & imaginary cones
01011         // -> 2nd root, with zi check
01012 
01013         b = nt2/nt1 ;
01014         c = nt3/nt1 ;
01015         d = b*b-c ;
01016         if (d >= 0)   // > 0
01017         {
01018            if(b>0){sd = c/( -b-std::sqrt(d));}
01019            else   {sd = -b + std::sqrt(d) ;}
01020 
01021           if ( sd >= 0 )   // > 0
01022           {
01023             if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
01024             {               // 64 bits systems. Split long distance and recompute
01025               G4double fTerm = sd-std::fmod(sd,dRmax);
01026               sd = fTerm + DistanceToIn(p+fTerm*v,v);
01027             } 
01028             zi = p.z() + sd*v.z() ;
01029 
01030             if ( std::fabs(zi) <= tolODz )
01031             {
01032               if ( !fPhiFullCone )
01033               {
01034                 xi     = p.x() + sd*v.x() ;
01035                 yi     = p.y() + sd*v.y() ;
01036                 ri     = rMinAv + zi*tanRMin ;
01037                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01038 
01039                 if (cosPsi >= cosHDPhiIT)
01040                 { 
01041                   if ( sd > halfRadTolerance )  { snxt=sd; }
01042                   else
01043                   {
01044                     // Calculate a normal vector in order to check Direction
01045 
01046                     risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01047                     Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
01048                     if ( Normal.dot(v) <= 0 )  { snxt = sd; }
01049                   } 
01050                 }
01051               }
01052               else
01053               {
01054                 if ( sd > halfRadTolerance )  { return sd; }
01055                 else
01056                 {
01057                   // Calculate a normal vector in order to check Direction
01058 
01059                   xi     = p.x() + sd*v.x() ;
01060                   yi     = p.y() + sd*v.y() ;
01061                   risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01062                   Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01063                   if ( Normal.dot(v) <= 0 )  { return sd; }
01064                 }
01065               }
01066             }
01067           }
01068         }
01069       }
01070       else  if ( nt3 < -rin*kRadTolerance*secRMin )
01071       {
01072         // Within radius of inner cone (real or imaginary)
01073         // -> Try 2nd root, with checking intersection is with real cone
01074         // -> If check fails, try 1st root, also checking intersection is
01075         //    on real cone
01076 
01077         b = nt2/nt1 ;
01078         c = nt3/nt1 ;
01079         d = b*b - c ;
01080 
01081         if ( d >= 0 )  // > 0
01082         {
01083           if (b>0) { sd = c/(-b-std::sqrt(d)); }
01084           else     { sd = -b + std::sqrt(d);   }
01085           zi = p.z() + sd*v.z() ;
01086           ri = rMinAv + zi*tanRMin ;
01087 
01088           if ( ri > 0 )
01089           {
01090             if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )  // sd > 0
01091             {
01092               if ( sd>dRmax ) // Avoid rounding errors due to precision issues
01093               {               // seen on 64 bits systems. Split and recompute
01094                 G4double fTerm = sd-std::fmod(sd,dRmax);
01095                 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01096               } 
01097               if ( !fPhiFullCone )
01098               {
01099                 xi     = p.x() + sd*v.x() ;
01100                 yi     = p.y() + sd*v.y() ;
01101                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01102 
01103                 if (cosPsi >= cosHDPhiOT)
01104                 {
01105                   if ( sd > halfRadTolerance )  { snxt=sd; }
01106                   else
01107                   {
01108                     // Calculate a normal vector in order to check Direction
01109 
01110                     risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01111                     Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
01112                     if ( Normal.dot(v) <= 0 )  { snxt = sd; } 
01113                   }
01114                 }
01115               }
01116               else
01117               {
01118                 if( sd > halfRadTolerance )  { return sd; }
01119                 else
01120                 {
01121                   // Calculate a normal vector in order to check Direction
01122 
01123                   xi     = p.x() + sd*v.x() ;
01124                   yi     = p.y() + sd*v.y() ;
01125                   risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01126                   Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01127                   if ( Normal.dot(v) <= 0 )  { return sd; }
01128                 } 
01129               }
01130             }
01131           }
01132           else
01133           {
01134             if (b>0) { sd = -b - std::sqrt(d);   }
01135             else     { sd = c/(-b+std::sqrt(d)); }
01136             zi = p.z() + sd*v.z() ;
01137             ri = rMinAv + zi*tanRMin ;
01138 
01139             if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
01140             {
01141               if ( sd>dRmax ) // Avoid rounding errors due to precision issues
01142               {               // seen on 64 bits systems. Split and recompute
01143                 G4double fTerm = sd-std::fmod(sd,dRmax);
01144                 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01145               } 
01146               if ( !fPhiFullCone )
01147               {
01148                 xi     = p.x() + sd*v.x() ;
01149                 yi     = p.y() + sd*v.y() ;
01150                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01151 
01152                 if (cosPsi >= cosHDPhiIT)
01153                 {
01154                   if ( sd > halfRadTolerance )  { snxt=sd; }
01155                   else
01156                   {
01157                     // Calculate a normal vector in order to check Direction
01158 
01159                     risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01160                     Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
01161                     if ( Normal.dot(v) <= 0 )  { snxt = sd; } 
01162                   }
01163                 }
01164               }
01165               else
01166               {
01167                 if ( sd > halfRadTolerance )  { return sd; }
01168                 else
01169                 {
01170                   // Calculate a normal vector in order to check Direction
01171 
01172                   xi     = p.x() + sd*v.x() ;
01173                   yi     = p.y() + sd*v.y() ;
01174                   risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01175                   Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01176                   if ( Normal.dot(v) <= 0 )  { return sd; }
01177                 } 
01178               }
01179             }
01180           }
01181         }
01182       }
01183       else
01184       {
01185         // Within kRadTol*0.5 of inner cone (real OR imaginary)
01186         // ----> Check not travelling through (=>0 to in)
01187         // ----> if not:
01188         //    -2nd root with validity check
01189 
01190         if ( std::fabs(p.z()) <= tolODz )
01191         {
01192           if ( nt2 > 0 )
01193           {
01194             // Inside inner real cone, heading outwards, inside z range
01195 
01196             if ( !fPhiFullCone )
01197             {
01198               cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
01199 
01200               if (cosPsi >= cosHDPhiIT)  { return 0.0; }
01201             }
01202             else  { return 0.0; }
01203           }
01204           else
01205           {
01206             // Within z extent, but not travelling through
01207             // -> 2nd root or kInfinity if 1st root on imaginary cone
01208 
01209             b = nt2/nt1 ;
01210             c = nt3/nt1 ;
01211             d = b*b - c ;
01212 
01213             if ( d >= 0 )   // > 0
01214             {
01215               if (b>0) { sd = -b - std::sqrt(d);   }
01216               else     { sd = c/(-b+std::sqrt(d)); }
01217               zi = p.z() + sd*v.z() ;
01218               ri = rMinAv + zi*tanRMin ;
01219               
01220               if ( ri > 0 )   // 2nd root
01221               {
01222                 if (b>0) { sd = c/(-b-std::sqrt(d)); }
01223                 else     { sd = -b + std::sqrt(d);   }
01224                 
01225                 zi = p.z() + sd*v.z() ;
01226 
01227                 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )  // sd>0
01228                 {
01229                   if ( sd>dRmax ) // Avoid rounding errors due to precision issue
01230                   {               // seen on 64 bits systems. Split and recompute
01231                     G4double fTerm = sd-std::fmod(sd,dRmax);
01232                     sd = fTerm + DistanceToIn(p+fTerm*v,v);
01233                   } 
01234                   if ( !fPhiFullCone )
01235                   {
01236                     xi     = p.x() + sd*v.x() ;
01237                     yi     = p.y() + sd*v.y() ;
01238                     ri     = rMinAv + zi*tanRMin ;
01239                     cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01240 
01241                     if ( cosPsi >= cosHDPhiIT )  { snxt = sd; }
01242                   }
01243                   else  { return sd; }
01244                 }
01245               }
01246               else  { return kInfinity; }
01247             }
01248           }
01249         }
01250         else   // 2nd root
01251         {
01252           b = nt2/nt1 ;
01253           c = nt3/nt1 ;
01254           d = b*b - c ;
01255 
01256           if ( d > 0 )
01257           {  
01258             if (b>0) { sd = c/(-b-std::sqrt(d)); }
01259             else     { sd = -b + std::sqrt(d) ;  }
01260             zi = p.z() + sd*v.z() ;
01261 
01262             if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )  // sd>0
01263             {
01264               if ( sd>dRmax ) // Avoid rounding errors due to precision issues
01265               {               // seen on 64 bits systems. Split and recompute
01266                 G4double fTerm = sd-std::fmod(sd,dRmax);
01267                 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01268               } 
01269               if ( !fPhiFullCone )
01270               {
01271                 xi     = p.x() + sd*v.x();
01272                 yi     = p.y() + sd*v.y();
01273                 ri     = rMinAv + zi*tanRMin ;
01274                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
01275 
01276                 if (cosPsi >= cosHDPhiIT)  { snxt = sd; }
01277               }
01278               else  { return sd; }
01279             }
01280           }
01281         }
01282       }
01283     }
01284   }
01285 
01286   // Phi segment intersection
01287   //
01288   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
01289   //
01290   // o NOTE: Large duplication of code between sphi & ephi checks
01291   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
01292   //            intersection check <=0 -> >=0
01293   //         -> Should use some form of loop Construct
01294 
01295   if ( !fPhiFullCone )
01296   {
01297     // First phi surface (starting phi)
01298 
01299     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
01300                     
01301     if ( Comp < 0 )    // Component in outwards normal dirn
01302     {
01303       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
01304 
01305       if (Dist < halfCarTolerance)
01306       {
01307         sd = Dist/Comp ;
01308 
01309         if ( sd < snxt )
01310         {
01311           if ( sd < 0 )  { sd = 0.0; }
01312 
01313           zi = p.z() + sd*v.z() ;
01314 
01315           if ( std::fabs(zi) <= tolODz )
01316           {
01317             xi        = p.x() + sd*v.x() ;
01318             yi        = p.y() + sd*v.y() ;
01319             rhoi2     = xi*xi + yi*yi ;
01320             tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
01321             tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
01322 
01323             if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
01324             {
01325               // z and r intersections good - check intersecting with
01326               // correct half-plane
01327 
01328               if ((yi*cosCPhi - xi*sinCPhi) <= 0 )  { snxt = sd; }
01329             }
01330           }
01331         }
01332       }
01333     }
01334 
01335     // Second phi surface (Ending phi)
01336 
01337     Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
01338         
01339     if ( Comp < 0 )   // Component in outwards normal dirn
01340     {
01341       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
01342       if (Dist < halfCarTolerance)
01343       {
01344         sd = Dist/Comp ;
01345 
01346         if ( sd < snxt )
01347         {
01348           if ( sd < 0 )  { sd = 0.0; }
01349 
01350           zi = p.z() + sd*v.z() ;
01351 
01352           if (std::fabs(zi) <= tolODz)
01353           {
01354             xi        = p.x() + sd*v.x() ;
01355             yi        = p.y() + sd*v.y() ;
01356             rhoi2     = xi*xi + yi*yi ;
01357             tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
01358             tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
01359 
01360             if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
01361             {
01362               // z and r intersections good - check intersecting with
01363               // correct half-plane
01364 
01365               if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  { snxt = sd; }
01366             }
01367           }
01368         }
01369       }
01370     }
01371   }
01372   if (snxt < halfCarTolerance)  { snxt = 0.; }
01373 
01374   return snxt ;
01375 }

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

Implements G4VSolid.

Definition at line 2073 of file G4Cons.cc.

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

02074 {
02075   G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
02076   G4double tanRMin, secRMin, pRMin;
02077   G4double tanRMax, secRMax, pRMax;
02078 
02079 #ifdef G4CSGDEBUG
02080   if( Inside(p) == kOutside )
02081   {
02082     G4int oldprc=G4cout.precision(16) ;
02083     G4cout << G4endl ;
02084     DumpInfo();
02085     G4cout << "Position:"  << G4endl << G4endl ;
02086     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
02087     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
02088     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
02089     G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
02090            << " mm" << G4endl << G4endl ;
02091     if( (p.x() != 0.) || (p.x() != 0.) )
02092     {
02093       G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree
02094              << " degree" << G4endl << G4endl ; 
02095     }
02096     G4cout.precision(oldprc) ;
02097     G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
02098                 JustWarning, "Point p is outside !?" );
02099   }
02100 #endif
02101 
02102   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
02103   safeZ = fDz - std::fabs(p.z()) ;
02104 
02105   if (fRmin1 || fRmin2)
02106   {
02107     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
02108     secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
02109     pRMin   = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
02110     safeR1  = (rho - pRMin)/secRMin ;
02111   }
02112   else
02113   {
02114     safeR1 = kInfinity ;
02115   }
02116 
02117   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
02118   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
02119   pRMax   = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
02120   safeR2  = (pRMax - rho)/secRMax ;
02121 
02122   if (safeR1 < safeR2)  { safe = safeR1; }
02123   else                  { safe = safeR2; }
02124   if (safeZ < safe)     { safe = safeZ ; }
02125 
02126   // Check if phi divided, Calc distances closest phi plane
02127 
02128   if (!fPhiFullCone)
02129   {
02130     // Above/below central phi of G4Cons?
02131 
02132     if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
02133     {
02134       safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
02135     }
02136     else
02137     {
02138       safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
02139     }
02140     if (safePhi < safe)  { safe = safePhi; }
02141   }
02142   if ( safe < 0 )  { safe = 0; }
02143 
02144   return safe ;
02145 }

G4double G4Cons::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 1446 of file G4Cons.cc.

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

01451 {
01452   ESide side = kNull, sider = kNull, sidephi = kNull;
01453 
01454   static const G4double halfCarTolerance=kCarTolerance*0.5;
01455   static const G4double halfRadTolerance=kRadTolerance*0.5;
01456   static const G4double halfAngTolerance=kAngTolerance*0.5;
01457 
01458   G4double snxt,srd,sphi,pdist ;
01459 
01460   G4double tanRMax, secRMax, rMaxAv ;  // Data for outer cone
01461   G4double tanRMin, secRMin, rMinAv ;  // Data for inner cone
01462 
01463   G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
01464   G4double b, c, d, sr2, sr3 ;
01465 
01466   // Vars for intersection within tolerance
01467 
01468   ESide    sidetol = kNull ;
01469   G4double slentol = kInfinity ;
01470 
01471   // Vars for phi intersection:
01472 
01473   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
01474   G4double zi, ri, deltaRoi2 ;
01475 
01476   // Z plane intersection
01477 
01478   if ( v.z() > 0.0 )
01479   {
01480     pdist = fDz - p.z() ;
01481 
01482     if (pdist > halfCarTolerance)
01483     {
01484       snxt = pdist/v.z() ;
01485       side = kPZ ;
01486     }
01487     else
01488     {
01489       if (calcNorm)
01490       {
01491         *n         = G4ThreeVector(0,0,1) ;
01492         *validNorm = true ;
01493       }
01494       return  snxt = 0.0;
01495     }
01496   }
01497   else if ( v.z() < 0.0 )
01498   {
01499     pdist = fDz + p.z() ;
01500 
01501     if ( pdist > halfCarTolerance)
01502     {
01503       snxt = -pdist/v.z() ;
01504       side = kMZ ;
01505     }
01506     else
01507     {
01508       if ( calcNorm )
01509       {
01510         *n         = G4ThreeVector(0,0,-1) ;
01511         *validNorm = true ;
01512       }
01513       return snxt = 0.0 ;
01514     }
01515   }
01516   else     // Travel perpendicular to z axis
01517   {
01518     snxt = kInfinity ;    
01519     side = kNull ;
01520   }
01521 
01522   // Radial Intersections
01523   //
01524   // Intersection with outer cone (possible return) and
01525   //                   inner cone (must also check phi)
01526   //
01527   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
01528   //
01529   // Intersects with x^2+y^2=(a*z+b)^2
01530   //
01531   // where a=tanRMax or tanRMin
01532   //       b=rMaxAv  or rMinAv
01533   //
01534   // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
01535   //     t1                        t2                      t3  
01536   //
01537   //  \--------u-------/       \-----------v----------/ \---------w--------/
01538 
01539   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
01540   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
01541   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
01542 
01543 
01544   t1   = 1.0 - v.z()*v.z() ;      // since v normalised
01545   t2   = p.x()*v.x() + p.y()*v.y() ;
01546   t3   = p.x()*p.x() + p.y()*p.y() ;
01547   rout = tanRMax*p.z() + rMaxAv ;
01548 
01549   nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
01550   nt2 = t2 - tanRMax*v.z()*rout ;
01551   nt3 = t3 - rout*rout ;
01552 
01553   if (v.z() > 0.0)
01554   {
01555     deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
01556                 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
01557   }
01558   else if ( v.z() < 0.0 )
01559   {
01560     deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
01561                 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
01562   }
01563   else
01564   {
01565     deltaRoi2 = 1.0;
01566   }
01567 
01568   if ( nt1 && (deltaRoi2 > 0.0) )  
01569   {
01570     // Equation quadratic => 2 roots : second root must be leaving
01571 
01572     b = nt2/nt1 ;
01573     c = nt3/nt1 ;
01574     d = b*b - c ;
01575 
01576     if ( d >= 0 )
01577     {
01578       // Check if on outer cone & heading outwards
01579       // NOTE: Should use rho-rout>-kRadTolerance*0.5
01580         
01581       if (nt3 > -halfRadTolerance && nt2 >= 0 )
01582       {
01583         if (calcNorm)
01584         {
01585           risec      = std::sqrt(t3)*secRMax ;
01586           *validNorm = true ;
01587           *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
01588         }
01589         return snxt=0 ;
01590       }
01591       else
01592       {
01593         sider = kRMax  ;
01594         if (b>0) { srd = -b - std::sqrt(d);    }
01595         else     { srd = c/(-b+std::sqrt(d)) ; }
01596 
01597         zi    = p.z() + srd*v.z() ;
01598         ri    = tanRMax*zi + rMaxAv ;
01599           
01600         if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
01601         {
01602           // An intersection within the tolerance
01603           //   we will Store it in case it is good -
01604           // 
01605           slentol = srd ;
01606           sidetol = kRMax ;
01607         }            
01608         if ( (ri < 0) || (srd < halfRadTolerance) )
01609         {
01610           // Safety: if both roots -ve ensure that srd cannot `win'
01611           //         distance to out
01612 
01613           if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
01614           else     { sr2 = -b + std::sqrt(d);   }
01615           zi  = p.z() + sr2*v.z() ;
01616           ri  = tanRMax*zi + rMaxAv ;
01617 
01618           if ((ri >= 0) && (sr2 > halfRadTolerance))
01619           {
01620             srd = sr2;
01621           }
01622           else
01623           {
01624             srd = kInfinity ;
01625 
01626             if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
01627             {
01628               // An intersection within the tolerance.
01629               // Storing it in case it is good.
01630 
01631               slentol = sr2 ;
01632               sidetol = kRMax ;
01633             }
01634           }
01635         }
01636       }
01637     }
01638     else
01639     {
01640       // No intersection with outer cone & not parallel
01641       // -> already outside, no intersection
01642 
01643       if ( calcNorm )
01644       {
01645         risec      = std::sqrt(t3)*secRMax;
01646         *validNorm = true;
01647         *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
01648       }
01649       return snxt = 0.0 ;
01650     }
01651   }
01652   else if ( nt2 && (deltaRoi2 > 0.0) )
01653   {
01654     // Linear case (only one intersection) => point outside outer cone
01655 
01656     if ( calcNorm )
01657     {
01658       risec      = std::sqrt(t3)*secRMax;
01659       *validNorm = true;
01660       *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
01661     }
01662     return snxt = 0.0 ;
01663   }
01664   else
01665   {
01666     // No intersection -> parallel to outer cone
01667     // => Z or inner cone intersection
01668 
01669     srd = kInfinity ;
01670   }
01671 
01672   // Check possible intersection within tolerance
01673 
01674   if ( slentol <= halfCarTolerance )
01675   {
01676     // An intersection within the tolerance was found.  
01677     // We must accept it only if the momentum points outwards.  
01678     //
01679     // G4ThreeVector ptTol ;  // The point of the intersection  
01680     // ptTol= p + slentol*v ;
01681     // ri=tanRMax*zi+rMaxAv ;
01682     //
01683     // Calculate a normal vector,  as below
01684 
01685     xi    = p.x() + slentol*v.x();
01686     yi    = p.y() + slentol*v.y();
01687     risec = std::sqrt(xi*xi + yi*yi)*secRMax;
01688     G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
01689 
01690     if ( Normal.dot(v) > 0 )    // We will leave the Cone immediatelly
01691     {
01692       if ( calcNorm ) 
01693       {
01694         *n         = Normal.unit() ;
01695         *validNorm = true ;
01696       }
01697       return snxt = 0.0 ;
01698     }
01699     else // On the surface, but not heading out so we ignore this intersection
01700     {    //                                        (as it is within tolerance).
01701       slentol = kInfinity ;
01702     }
01703   }
01704 
01705   // Inner Cone intersection
01706 
01707   if ( fRmin1 || fRmin2 )
01708   {
01709     tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
01710     nt1     = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
01711 
01712     if ( nt1 )
01713     {
01714       secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
01715       rMinAv  = (fRmin1 + fRmin2)*0.5 ;    
01716       rin     = tanRMin*p.z() + rMinAv ;
01717       nt2     = t2 - tanRMin*v.z()*rin ;
01718       nt3     = t3 - rin*rin ;
01719       
01720       // Equation quadratic => 2 roots : first root must be leaving
01721 
01722       b = nt2/nt1 ;
01723       c = nt3/nt1 ;
01724       d = b*b - c ;
01725 
01726       if ( d >= 0.0 )
01727       {
01728         // NOTE: should be rho-rin<kRadTolerance*0.5,
01729         //       but using squared versions for efficiency
01730 
01731         if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25)) 
01732         {
01733           if ( nt2 < 0.0 )
01734           {
01735             if (calcNorm)  { *validNorm = false; }
01736             return          snxt      = 0.0;
01737           }
01738         }
01739         else
01740         {
01741           if (b>0) { sr2 = -b - std::sqrt(d);   }
01742           else     { sr2 = c/(-b+std::sqrt(d)); }
01743           zi  = p.z() + sr2*v.z() ;
01744           ri  = tanRMin*zi + rMinAv ;
01745 
01746           if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
01747           {
01748             // An intersection within the tolerance
01749             // storing it in case it is good.
01750 
01751             slentol = sr2 ;
01752             sidetol = kRMax ;
01753           }
01754           if( (ri<0) || (sr2 < halfRadTolerance) )
01755           {
01756             if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
01757             else     { sr3 = -b + std::sqrt(d) ;  }
01758 
01759             // Safety: if both roots -ve ensure that srd cannot `win'
01760             //         distancetoout
01761 
01762             if  ( sr3 > halfRadTolerance )
01763             {
01764               if( sr3 < srd )
01765               {
01766                 zi = p.z() + sr3*v.z() ;
01767                 ri = tanRMin*zi + rMinAv ;
01768 
01769                 if ( ri >= 0.0 )
01770                 {
01771                   srd=sr3 ;
01772                   sider=kRMin ;
01773                 }
01774               } 
01775             }
01776             else if ( sr3 > -halfRadTolerance )
01777             {
01778               // Intersection in tolerance. Store to check if it's good
01779 
01780               slentol = sr3 ;
01781               sidetol = kRMin ;
01782             }
01783           }
01784           else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
01785           {
01786             srd   = sr2 ;
01787             sider = kRMin ;
01788           }
01789           else if (sr2 > -halfCarTolerance)
01790           {
01791             // Intersection in tolerance. Store to check if it's good
01792 
01793             slentol = sr2 ;
01794             sidetol = kRMin ;
01795           }    
01796           if( slentol <= halfCarTolerance  )
01797           {
01798             // An intersection within the tolerance was found. 
01799             // We must accept it only if  the momentum points outwards. 
01800 
01801             G4ThreeVector Normal ; 
01802             
01803             // Calculate a normal vector,  as below
01804 
01805             xi     = p.x() + slentol*v.x() ;
01806             yi     = p.y() + slentol*v.y() ;
01807             if( sidetol==kRMax )
01808             {
01809               risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
01810               Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
01811             }
01812             else
01813             {
01814               risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
01815               Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01816             }
01817             if( Normal.dot(v) > 0 )
01818             {
01819               // We will leave the cone immediately
01820 
01821               if( calcNorm ) 
01822               {
01823                 *n         = Normal.unit() ;
01824                 *validNorm = true ;
01825               }
01826               return snxt = 0.0 ;
01827             }
01828             else 
01829             { 
01830               // On the surface, but not heading out so we ignore this
01831               // intersection (as it is within tolerance). 
01832 
01833               slentol = kInfinity ;
01834             }        
01835           }
01836         }
01837       }
01838     }
01839   }
01840 
01841   // Linear case => point outside inner cone ---> outer cone intersect
01842   //
01843   // Phi Intersection
01844   
01845   if ( !fPhiFullCone )
01846   {
01847     // add angle calculation with correction 
01848     // of the difference in domain of atan2 and Sphi
01849 
01850     vphi = std::atan2(v.y(),v.x()) ;
01851 
01852     if ( vphi < fSPhi - halfAngTolerance  )              { vphi += twopi; }
01853     else if ( vphi > fSPhi + fDPhi + halfAngTolerance )  { vphi -= twopi; }
01854 
01855     if ( p.x() || p.y() )   // Check if on z axis (rho not needed later)
01856     {
01857       // pDist -ve when inside
01858 
01859       pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
01860       pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
01861 
01862       // Comp -ve when in direction of outwards normal
01863 
01864       compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
01865       compE = sinEPhi*v.x() - cosEPhi*v.y() ;
01866 
01867       sidephi = kNull ;
01868 
01869       if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
01870                             && (pDistE <= halfCarTolerance) ) )
01871          || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
01872                               && (pDistE >  halfCarTolerance) ) )  )
01873       {
01874         // Inside both phi *full* planes
01875         if ( compS < 0 )
01876         {
01877           sphi = pDistS/compS ;
01878           if (sphi >= -halfCarTolerance)
01879           {
01880             xi = p.x() + sphi*v.x() ;
01881             yi = p.y() + sphi*v.y() ;
01882 
01883             // Check intersecting with correct half-plane
01884             // (if not -> no intersect)
01885             //
01886             if ( (std::fabs(xi)<=kCarTolerance)
01887               && (std::fabs(yi)<=kCarTolerance) )
01888             {
01889               sidephi= kSPhi;
01890               if ( ( fSPhi-halfAngTolerance <= vphi )
01891                 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
01892               {
01893                 sphi = kInfinity;
01894               }
01895             }
01896             else
01897             if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
01898             {
01899               sphi = kInfinity ;
01900             }
01901             else
01902             {
01903               sidephi = kSPhi ;
01904               if ( pDistS > -halfCarTolerance )
01905               {
01906                 sphi = 0.0 ; // Leave by sphi immediately
01907               }    
01908             }       
01909           }
01910           else
01911           {
01912             sphi = kInfinity ;
01913           }
01914         }
01915         else
01916         {
01917           sphi = kInfinity ;
01918         }
01919 
01920         if ( compE < 0 )
01921         {
01922           sphi2 = pDistE/compE ;
01923 
01924           // Only check further if < starting phi intersection
01925           //
01926           if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
01927           {
01928             xi = p.x() + sphi2*v.x() ;
01929             yi = p.y() + sphi2*v.y() ;
01930 
01931             // Check intersecting with correct half-plane
01932 
01933             if ( (std::fabs(xi)<=kCarTolerance)
01934               && (std::fabs(yi)<=kCarTolerance) )
01935             {
01936               // Leaving via ending phi
01937 
01938               if(!( (fSPhi-halfAngTolerance <= vphi)
01939                  && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
01940               {
01941                 sidephi = kEPhi ;
01942                 if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
01943                 else                                { sphi = 0.0; }
01944               }
01945             }
01946             else // Check intersecting with correct half-plane
01947             if ( yi*cosCPhi-xi*sinCPhi >= 0 )
01948             {
01949               // Leaving via ending phi
01950 
01951               sidephi = kEPhi ;
01952               if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
01953               else                                { sphi = 0.0; }
01954             }
01955           }
01956         }
01957       }
01958       else
01959       {
01960         sphi = kInfinity ;
01961       }
01962     }
01963     else
01964     {
01965       // On z axis + travel not || to z axis -> if phi of vector direction
01966       // within phi of shape, Step limited by rmax, else Step =0
01967 
01968       if ( (fSPhi-halfAngTolerance <= vphi)
01969         && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
01970       {
01971         sphi = kInfinity ;
01972       }
01973       else
01974       {
01975         sidephi = kSPhi  ;   // arbitrary 
01976         sphi    = 0.0 ;
01977       }
01978     }      
01979     if ( sphi < snxt )  // Order intersecttions
01980     {
01981       snxt=sphi ;
01982       side=sidephi ;
01983     }
01984   }
01985   if ( srd < snxt )  // Order intersections
01986   {
01987     snxt = srd   ;
01988     side = sider ;
01989   }
01990   if (calcNorm)
01991   {
01992     switch(side)
01993     {                     // Note: returned vector not normalised
01994       case kRMax:         // (divide by frmax for unit vector)
01995         xi         = p.x() + snxt*v.x() ;
01996         yi         = p.y() + snxt*v.y() ;
01997         risec      = std::sqrt(xi*xi + yi*yi)*secRMax ;
01998         *n         = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
01999         *validNorm = true ;
02000         break ;
02001       case kRMin:
02002         *validNorm = false ;  // Rmin is inconvex
02003         break ;
02004       case kSPhi:
02005         if ( fDPhi <= pi )
02006         {
02007           *n         = G4ThreeVector(sinSPhi, -cosSPhi, 0);
02008           *validNorm = true ;
02009         }
02010         else
02011         {
02012           *validNorm = false ;
02013         }
02014         break ;
02015       case kEPhi:
02016         if ( fDPhi <= pi )
02017         {
02018           *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
02019           *validNorm = true ;
02020         }
02021         else
02022         {
02023           *validNorm = false ;
02024         }
02025         break ;
02026       case kPZ:
02027         *n         = G4ThreeVector(0,0,1) ;
02028         *validNorm = true ;
02029         break ;
02030       case kMZ:
02031         *n         = G4ThreeVector(0,0,-1) ;
02032         *validNorm = true ;
02033         break ;
02034       default:
02035         G4cout << G4endl ;
02036         DumpInfo();
02037         std::ostringstream message;
02038         G4int oldprc = message.precision(16) ;
02039         message << "Undefined side for valid surface normal to solid."
02040                 << G4endl
02041                 << "Position:"  << G4endl << G4endl
02042                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
02043                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
02044                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
02045                 << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
02046                 << " mm" << G4endl << G4endl ;
02047         if( p.x() != 0. || p.x() != 0.)
02048         {
02049            message << "point phi = "   << std::atan2(p.y(),p.x())/degree
02050                    << " degree" << G4endl << G4endl ; 
02051         }
02052         message << "Direction:" << G4endl << G4endl
02053                 << "v.x() = "   << v.x() << G4endl
02054                 << "v.y() = "   << v.y() << G4endl
02055                 << "v.z() = "   << v.z() << G4endl<< G4endl
02056                 << "Proposed distance :" << G4endl<< G4endl
02057                 << "snxt = "    << snxt/mm << " mm" << G4endl ;
02058         message.precision(oldprc) ;
02059         G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
02060                     JustWarning, message) ;
02061         break ;
02062     }
02063   }
02064   if (snxt < halfCarTolerance)  { snxt = 0.; }
02065 
02066   return snxt ;
02067 }

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

Reimplemented from G4VSolid.

Definition at line 254 of file G4Cons.icc.

References G4CSGSolid::fCubicVolume.

00255 {
00256   if(fCubicVolume != 0.) {;}
00257   else   
00258   {
00259     G4double Rmean, rMean, deltaR, deltar;
00260 
00261     Rmean  = 0.5*(fRmax1+fRmax2);
00262     deltaR = fRmax1-fRmax2;
00263 
00264     rMean  = 0.5*(fRmin1+fRmin2);
00265     deltar = fRmin1-fRmin2;
00266     fCubicVolume = fDPhi*fDz*(Rmean*Rmean-rMean*rMean
00267                             +(deltaR*deltaR-deltar*deltar)/12); 
00268   }
00269   return fCubicVolume;
00270 }

G4double G4Cons::GetDeltaPhiAngle (  )  const [inline]

Definition at line 74 of file G4Cons.icc.

Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsPhi::G4ParameterisationConsPhi(), G4VParameterisationCons::G4VParameterisationCons(), GetDPhi(), G4ParameterisationConsPhi::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

00075 {
00076   return fDPhi;
00077 }

G4double G4Cons::GetDPhi (  )  const [inline]

Definition at line 248 of file G4Cons.icc.

References GetDeltaPhiAngle().

00249 {
00250   return GetDeltaPhiAngle();
00251 }

G4double G4Cons::GetDz (  )  const [inline]

Definition at line 236 of file G4Cons.icc.

References GetZHalfLength().

00237 {
00238   return GetZHalfLength();
00239 }

G4GeometryType G4Cons::GetEntityType (  )  const [virtual]

Implements G4VSolid.

Definition at line 2245 of file G4Cons.cc.

02246 {
02247   return G4String("G4Cons");
02248 }

G4double G4Cons::GetInnerRadiusMinusZ (  )  const [inline]

Definition at line 38 of file G4Cons.icc.

Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), GetRmin1(), and G4tgbGeometryDumper::GetSolidParams().

00039 {
00040   return fRmin1 ;
00041 }

G4double G4Cons::GetInnerRadiusPlusZ (  )  const [inline]

Definition at line 50 of file G4Cons.icc.

Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), GetRmin2(), and G4tgbGeometryDumper::GetSolidParams().

00051 {
00052   return fRmin2 ;
00053 }

G4double G4Cons::GetOuterRadiusMinusZ (  )  const [inline]

Definition at line 44 of file G4Cons.icc.

Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), GetRmax1(), and G4tgbGeometryDumper::GetSolidParams().

00045 {
00046   return fRmax1 ;
00047 }

G4double G4Cons::GetOuterRadiusPlusZ (  )  const [inline]

Definition at line 56 of file G4Cons.icc.

Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), GetRmax2(), and G4tgbGeometryDumper::GetSolidParams().

00057 {
00058   return fRmax2 ;
00059 }

G4ThreeVector G4Cons::GetPointOnSurface (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 2290 of file G4Cons.cc.

References G4CSGSolid::GetRadiusInRing(), and sqr().

02291 {   
02292   // declare working variables
02293   //
02294   G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
02295   G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
02296   rone = (fRmax1-fRmax2)/(2.*fDz);
02297   rtwo = (fRmin1-fRmin2)/(2.*fDz);
02298   qone=0.; qtwo=0.;
02299   if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
02300   if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
02301   slin   = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz));
02302   slout  = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz));
02303   Aone   = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;       
02304   Atwo   = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
02305   Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); 
02306   Afour  = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
02307   Afive  = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
02308   
02309   phi    = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
02310   cosu   = std::cos(phi);  sinu = std::sin(phi);
02311   rRand1 = GetRadiusInRing(fRmin1, fRmin2);
02312   rRand2 = GetRadiusInRing(fRmax1, fRmax2);
02313   
02314   if ( (fSPhi == 0.) && fPhiFullCone )  { Afive = 0.; }
02315   chose  = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
02316  
02317   if( (chose >= 0.) && (chose < Aone) )
02318   {
02319     if(fRmin1 != fRmin2)
02320     {
02321       zRand = RandFlat::shoot(-1.*fDz,fDz); 
02322       return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
02323                             rtwo*sinu*(qtwo-zRand), zRand);
02324     }
02325     else
02326     {
02327       return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
02328                            RandFlat::shoot(-1.*fDz,fDz));
02329     }
02330   }
02331   else if( (chose >= Aone) && (chose <= Aone + Atwo) )
02332   {
02333     if(fRmax1 != fRmax2)
02334     {
02335       zRand = RandFlat::shoot(-1.*fDz,fDz); 
02336       return G4ThreeVector (rone*cosu*(qone-zRand),
02337                             rone*sinu*(qone-zRand), zRand);
02338     }    
02339     else
02340     {
02341       return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
02342                            RandFlat::shoot(-1.*fDz,fDz));
02343     }
02344   }
02345   else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
02346   {
02347     return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
02348   }
02349   else if( (chose >= Aone + Atwo + Athree)
02350         && (chose < Aone + Atwo + Athree + Afour) )
02351   {
02352     return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz);
02353   }
02354   else if( (chose >= Aone + Atwo + Athree + Afour)
02355         && (chose < Aone + Atwo + Athree + Afour + Afive) )
02356   {
02357     zRand  = RandFlat::shoot(-1.*fDz,fDz);
02358     rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
02359                              fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 
02360     return G4ThreeVector (rRand1*std::cos(fSPhi),
02361                           rRand1*std::sin(fSPhi), zRand);
02362   }
02363   else
02364   { 
02365     zRand  = RandFlat::shoot(-1.*fDz,fDz);
02366     rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
02367                              fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 
02368     return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi),
02369                           rRand1*std::sin(fSPhi+fDPhi), zRand);
02370   }
02371 }

G4double G4Cons::GetRmax1 (  )  const [inline]

Definition at line 218 of file G4Cons.icc.

References GetOuterRadiusMinusZ().

00219 {
00220   return GetOuterRadiusMinusZ();
00221 }

G4double G4Cons::GetRmax2 (  )  const [inline]

Definition at line 230 of file G4Cons.icc.

References GetOuterRadiusPlusZ().

00231 {
00232   return GetOuterRadiusPlusZ();
00233 }

G4double G4Cons::GetRmin1 (  )  const [inline]

Definition at line 212 of file G4Cons.icc.

References GetInnerRadiusMinusZ().

00213 {
00214   return GetInnerRadiusMinusZ();
00215 }

G4double G4Cons::GetRmin2 (  )  const [inline]

Definition at line 224 of file G4Cons.icc.

References GetInnerRadiusPlusZ().

00225 {
00226   return GetInnerRadiusPlusZ();
00227 }

G4double G4Cons::GetSPhi (  )  const [inline]

Definition at line 242 of file G4Cons.icc.

References GetStartPhiAngle().

00243 {
00244   return GetStartPhiAngle();
00245 }

G4double G4Cons::GetStartPhiAngle (  )  const [inline]

Definition at line 68 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), G4tgbGeometryDumper::GetSolidParams(), and GetSPhi().

00069 {
00070   return fSPhi ;
00071 }

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

Reimplemented from G4VSolid.

Definition at line 273 of file G4Cons.icc.

References G4CSGSolid::fSurfaceArea.

00274 {
00275   if(fSurfaceArea != 0.) {;}
00276   else   
00277   {
00278     G4double mmin, mmax, dmin, dmax;
00279 
00280     mmin= (fRmin1+fRmin2)*0.5;
00281     mmax= (fRmax1+fRmax2)*0.5;
00282     dmin= (fRmin2-fRmin1);
00283     dmax= (fRmax2-fRmax1);
00284     
00285     fSurfaceArea = fDPhi*( mmin * std::sqrt(dmin*dmin+4*fDz*fDz)
00286                          + mmax * std::sqrt(dmax*dmax+4*fDz*fDz)
00287                          + 0.5*(fRmax1*fRmax1-fRmin1*fRmin1
00288                                +fRmax2*fRmax2-fRmin2*fRmin2 ));
00289     if(!fPhiFullCone)
00290     {
00291       fSurfaceArea = fSurfaceArea+4*fDz*(mmax-mmin);
00292     }
00293   }
00294   return fSurfaceArea;
00295 }

G4double G4Cons::GetZHalfLength (  )  const [inline]

Definition at line 62 of file G4Cons.icc.

Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsZ::ComputeTransformation(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsZ::G4ParameterisationConsZ(), G4VParameterisationCons::G4VParameterisationCons(), GetDz(), G4ParameterisationConsZ::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

00063 {
00064   return fDz ;
00065 }

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

Implements G4VSolid.

Definition at line 191 of file G4Cons.cc.

References G4VSolid::kCarTolerance, kInside, kOutside, and kSurface.

Referenced by CalculateExtent(), and DistanceToOut().

00192 {
00193   G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
00194   EInside in;
00195   static const G4double halfCarTolerance=kCarTolerance*0.5;
00196   static const G4double halfRadTolerance=kRadTolerance*0.5;
00197   static const G4double halfAngTolerance=kAngTolerance*0.5;
00198 
00199   if (std::fabs(p.z()) > fDz + halfCarTolerance )  { return in = kOutside; }
00200   else if(std::fabs(p.z()) >= fDz - halfCarTolerance )    { in = kSurface; }
00201   else                                                    { in = kInside;  }
00202 
00203   r2 = p.x()*p.x() + p.y()*p.y() ;
00204   rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
00205   rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
00206 
00207   // rh2 = rh*rh;
00208 
00209   tolRMin = rl - halfRadTolerance;
00210   if ( tolRMin < 0 )  { tolRMin = 0; }
00211   tolRMax = rh + halfRadTolerance;
00212 
00213   if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
00214 
00215   if (rl) { tolRMin = rl + halfRadTolerance; }
00216   else    { tolRMin = 0.0; }
00217   tolRMax = rh - halfRadTolerance;
00218       
00219   if (in == kInside) // else it's kSurface already
00220   {
00221      if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
00222   }
00223   if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
00224   {
00225     pPhi = std::atan2(p.y(),p.x()) ;
00226 
00227     if ( pPhi < fSPhi - halfAngTolerance  )             { pPhi += twopi; }
00228     else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
00229     
00230     if ( (pPhi < fSPhi - halfAngTolerance) ||          
00231          (pPhi > fSPhi + fDPhi + halfAngTolerance) )  { return in = kOutside; }
00232       
00233     else if (in == kInside)  // else it's kSurface anyway already
00234     {
00235        if ( (pPhi < fSPhi + halfAngTolerance) || 
00236             (pPhi > fSPhi + fDPhi - halfAngTolerance) )  { in = kSurface; }
00237     }
00238   }
00239   else if ( !fPhiFullCone )  { in = kSurface; }
00240 
00241   return in ;
00242 }

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

Definition at line 161 of file G4Cons.cc.

References cosCPhi, cosEPhi, cosHDPhiIT, cosHDPhiOT, cosSPhi, fDPhi, fDz, fPhiFullCone, fRmax1, fRmax2, fRmin1, fRmin2, fSPhi, kAngTolerance, kRadTolerance, G4CSGSolid::operator=(), sinCPhi, sinEPhi, and sinSPhi.

00162 {
00163    // Check assignment to self
00164    //
00165    if (this == &rhs)  { return *this; }
00166 
00167    // Copy base class data
00168    //
00169    G4CSGSolid::operator=(rhs);
00170 
00171    // Copy data
00172    //
00173    kRadTolerance = rhs.kRadTolerance;
00174    kAngTolerance = rhs.kAngTolerance;
00175    fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
00176    fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
00177    fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
00178    sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
00179    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
00180    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
00181    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
00182    fPhiFullCone = rhs.fPhiFullCone;
00183 
00184    return *this;
00185 }

void G4Cons::SetDeltaPhiAngle ( G4double  newDPhi  )  [inline]

Definition at line 203 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00204 {
00205   CheckPhiAngles(fSPhi, newDPhi);
00206   Initialize();
00207 }

void G4Cons::SetInnerRadiusMinusZ ( G4double  Rmin1  )  [inline]

Definition at line 157 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00158 {
00159   fRmin1= Rmin1 ;
00160   Initialize();
00161 }

void G4Cons::SetInnerRadiusPlusZ ( G4double  Rmin2  )  [inline]

Definition at line 171 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00172 {
00173   fRmin2= Rmin2 ;
00174   Initialize();
00175 }

void G4Cons::SetOuterRadiusMinusZ ( G4double  Rmax1  )  [inline]

Definition at line 164 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00165 {
00166   fRmax1= Rmax1 ;
00167   Initialize();
00168 }

void G4Cons::SetOuterRadiusPlusZ ( G4double  Rmax2  )  [inline]

Definition at line 178 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00179 {
00180   fRmax2= Rmax2 ;
00181   Initialize();
00182 }

void G4Cons::SetStartPhiAngle ( G4double  newSPhi,
G4bool  trig = true 
) [inline]

Definition at line 192 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00193 {
00194   // Flag 'compute' can be used to explicitely avoid recomputation of
00195   // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
00196 
00197   CheckSPhiAngle(newSPhi);
00198   fPhiFullCone = false;
00199   if (compute)  { InitializeTrigonometry(); }
00200   Initialize();
00201 }

void G4Cons::SetZHalfLength ( G4double  newDz  )  [inline]

Definition at line 185 of file G4Cons.icc.

Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().

00186 {
00187   fDz= newDz ;
00188   Initialize();
00189 }

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

Reimplemented from G4CSGSolid.

Definition at line 2263 of file G4Cons.cc.

References G4VSolid::GetName().

02264 {
02265   G4int oldprc = os.precision(16);
02266   os << "-----------------------------------------------------------\n"
02267      << "    *** Dump for solid - " << GetName() << " ***\n"
02268      << "    ===================================================\n"
02269      << " Solid type: G4Cons\n"
02270      << " Parameters: \n"
02271      << "   inside  -fDz radius: "  << fRmin1/mm << " mm \n"
02272      << "   outside -fDz radius: "  << fRmax1/mm << " mm \n"
02273      << "   inside  +fDz radius: "  << fRmin2/mm << " mm \n"
02274      << "   outside +fDz radius: "  << fRmax2/mm << " mm \n"
02275      << "   half length in Z   : "  << fDz/mm << " mm \n"
02276      << "   starting angle of segment: " << fSPhi/degree << " degrees \n"
02277      << "   delta angle of segment   : " << fDPhi/degree << " degrees \n"
02278      << "-----------------------------------------------------------\n";
02279   os.precision(oldprc);
02280 
02281   return os;
02282 }

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

Implements G4VSolid.

Definition at line 479 of file G4Cons.cc.

References G4Exception(), JustWarning, and G4VSolid::kCarTolerance.

00480 {
00481   G4int noSurfaces = 0;
00482   G4double rho, pPhi;
00483   G4double distZ, distRMin, distRMax;
00484   G4double distSPhi = kInfinity, distEPhi = kInfinity;
00485   G4double tanRMin, secRMin, pRMin, widRMin;
00486   G4double tanRMax, secRMax, pRMax, widRMax;
00487 
00488   static const G4double delta  = 0.5*kCarTolerance;
00489   static const G4double dAngle = 0.5*kAngTolerance;
00490   
00491   G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
00492   G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
00493 
00494   distZ = std::fabs(std::fabs(p.z()) - fDz);
00495   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y());
00496 
00497   tanRMin  = (fRmin2 - fRmin1)*0.5/fDz;
00498   secRMin  = std::sqrt(1 + tanRMin*tanRMin);
00499   pRMin    = rho - p.z()*tanRMin;
00500   widRMin  = fRmin2 - fDz*tanRMin;
00501   distRMin = std::fabs(pRMin - widRMin)/secRMin;
00502 
00503   tanRMax  = (fRmax2 - fRmax1)*0.5/fDz;
00504   secRMax  = std::sqrt(1+tanRMax*tanRMax);
00505   pRMax    = rho - p.z()*tanRMax;
00506   widRMax  = fRmax2 - fDz*tanRMax;
00507   distRMax = std::fabs(pRMax - widRMax)/secRMax;
00508 
00509   if (!fPhiFullCone)   // Protected against (0,0,z) 
00510   {
00511     if ( rho )
00512     {
00513       pPhi = std::atan2(p.y(),p.x());
00514 
00515       if (pPhi  < fSPhi-delta)            { pPhi += twopi; }
00516       else if (pPhi > fSPhi+fDPhi+delta)  { pPhi -= twopi; }
00517 
00518       distSPhi = std::fabs( pPhi - fSPhi ); 
00519       distEPhi = std::fabs( pPhi - fSPhi - fDPhi ); 
00520     }
00521     else if( !(fRmin1) || !(fRmin2) )
00522     {
00523       distSPhi = 0.; 
00524       distEPhi = 0.; 
00525     }
00526     nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
00527     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
00528   }
00529   if ( rho > delta )   
00530   {
00531     nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
00532     if (fRmin1 || fRmin2)
00533     {
00534       nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
00535     }
00536   }
00537 
00538   if( distRMax <= delta )
00539   {
00540     noSurfaces ++;
00541     sumnorm += nR;
00542   }
00543   if( (fRmin1 || fRmin2) && (distRMin <= delta) )
00544   {
00545     noSurfaces ++;
00546     sumnorm += nr;
00547   }
00548   if( !fPhiFullCone )   
00549   {
00550     if (distSPhi <= dAngle)
00551     {
00552       noSurfaces ++;
00553       sumnorm += nPs;
00554     }
00555     if (distEPhi <= dAngle) 
00556     {
00557       noSurfaces ++;
00558       sumnorm += nPe;
00559     }
00560   }
00561   if (distZ <= delta)  
00562   {
00563     noSurfaces ++;
00564     if ( p.z() >= 0.)  { sumnorm += nZ; }
00565     else               { sumnorm -= nZ; }
00566   }
00567   if ( noSurfaces == 0 )
00568   {
00569 #ifdef G4CSGDEBUG
00570     G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
00571                 JustWarning, "Point p is not on surface !?" );
00572 #endif 
00573      norm = ApproxSurfaceNormal(p);
00574   }
00575   else if ( noSurfaces == 1 )  { norm = sumnorm; }
00576   else                         { norm = sumnorm.unit(); }
00577 
00578   return norm ;
00579 }


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