G4Torus Class Reference

#include <G4Torus.hh>

Inheritance diagram for G4Torus:

G4CSGSolid G4VSolid

Public Member Functions

 G4Torus (const G4String &pName, G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
 ~G4Torus ()
G4double GetRmin () const
G4double GetRmax () const
G4double GetRtor () const
G4double GetSPhi () const
G4double GetDPhi () const
G4double GetCubicVolume ()
G4double GetSurfaceArea ()
EInside Inside (const G4ThreeVector &p) const
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
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
void SetAllParameters (G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
 G4Torus (__void__ &)
 G4Torus (const G4Torus &rhs)
G4Torusoperator= (const G4Torus &rhs)

Detailed Description

Definition at line 102 of file G4Torus.hh.


Constructor & Destructor Documentation

G4Torus::G4Torus ( const G4String pName,
G4double  pRmin,
G4double  pRmax,
G4double  pRtor,
G4double  pSPhi,
G4double  pDPhi 
)

Definition at line 84 of file G4Torus.cc.

References SetAllParameters().

Referenced by Clone().

00090   : G4CSGSolid(pName)
00091 {
00092   SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
00093 }

G4Torus::~G4Torus (  ) 

Definition at line 193 of file G4Torus.cc.

00194 {}

G4Torus::G4Torus ( __void__ &   ) 

Definition at line 182 of file G4Torus.cc.

00183   : G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
00184     fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
00185     kRadTolerance(0.), kAngTolerance(0.)
00186 {
00187 }

G4Torus::G4Torus ( const G4Torus rhs  ) 

Definition at line 200 of file G4Torus.cc.

00201   : G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
00202     fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
00203     fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
00204     kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance)
00205 {
00206 }


Member Function Documentation

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

Implements G4VSolid.

Definition at line 412 of file G4Torus.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().

00416 {
00417   if ((!pTransform.IsRotated()) && (fDPhi==twopi) && (fRmin==0))
00418   {
00419     // Special case handling for unrotated solid torus
00420     // Compute x/y/z mins and maxs for bounding box respecting limits,
00421     // with early returns if outside limits. Then switch() on pAxis,
00422     // and compute exact x and y limit for x/y case
00423       
00424     G4double xoffset,xMin,xMax;
00425     G4double yoffset,yMin,yMax;
00426     G4double zoffset,zMin,zMax;
00427 
00428     G4double RTorus,delta,diff1,diff2,maxDiff,newMin,newMax;
00429     G4double xoff1,xoff2,yoff1,yoff2;
00430 
00431     xoffset = pTransform.NetTranslation().x();
00432     xMin    = xoffset - fRmax - fRtor ;
00433     xMax    = xoffset + fRmax + fRtor ;
00434 
00435     if (pVoxelLimit.IsXLimited())
00436     {
00437       if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00438         || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00439         return false ;
00440       else
00441       {
00442         if (xMin < pVoxelLimit.GetMinXExtent())
00443         {
00444           xMin = pVoxelLimit.GetMinXExtent() ;
00445         }
00446         if (xMax > pVoxelLimit.GetMaxXExtent())
00447         {
00448           xMax = pVoxelLimit.GetMaxXExtent() ;
00449         }
00450       }
00451     }
00452     yoffset = pTransform.NetTranslation().y();
00453     yMin    = yoffset - fRmax - fRtor ;
00454     yMax    = yoffset + fRmax + fRtor ;
00455 
00456     if (pVoxelLimit.IsYLimited())
00457     {
00458       if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00459         || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00460       {
00461         return false ;
00462       }
00463       else
00464       {
00465         if (yMin < pVoxelLimit.GetMinYExtent() )
00466         { 
00467           yMin = pVoxelLimit.GetMinYExtent() ;
00468         }
00469         if (yMax > pVoxelLimit.GetMaxYExtent() )
00470         {
00471           yMax = pVoxelLimit.GetMaxYExtent() ;
00472         }
00473       }
00474     }
00475     zoffset = pTransform.NetTranslation().z() ;
00476     zMin    = zoffset - fRmax ;
00477     zMax    = zoffset + fRmax ;
00478 
00479     if (pVoxelLimit.IsZLimited())
00480     {
00481       if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00482         || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) )
00483       {
00484         return false ;
00485       }
00486       else
00487       {
00488         if (zMin < pVoxelLimit.GetMinZExtent() )
00489         {
00490           zMin = pVoxelLimit.GetMinZExtent() ;
00491         }
00492         if (zMax > pVoxelLimit.GetMaxZExtent() )
00493         {
00494           zMax = pVoxelLimit.GetMaxZExtent() ;
00495         }
00496       }
00497     }
00498 
00499     // Known to cut cylinder
00500     
00501     switch (pAxis)
00502     {
00503       case kXAxis:
00504         yoff1=yoffset-yMin;
00505         yoff2=yMax-yoffset;
00506         if ( yoff1 >= 0 && yoff2 >= 0 )
00507         {
00508           // Y limits cross max/min x => no change
00509           //
00510           pMin = xMin ;
00511           pMax = xMax ;
00512         }
00513         else
00514         {
00515           // Y limits don't cross max/min x => compute max delta x,
00516           // hence new mins/maxs
00517           //
00518        
00519           RTorus=fRmax+fRtor;
00520           delta   = RTorus*RTorus - yoff1*yoff1;
00521           diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
00522           delta   = RTorus*RTorus - yoff2*yoff2;
00523           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
00524           maxDiff = (diff1 > diff2) ? diff1:diff2 ;
00525           newMin  = xoffset - maxDiff ;
00526           newMax  = xoffset + maxDiff ;
00527           pMin    = (newMin < xMin) ? xMin : newMin ;
00528           pMax    = (newMax > xMax) ? xMax : newMax ;
00529         }
00530         break;
00531 
00532       case kYAxis:
00533         xoff1 = xoffset - xMin ;
00534         xoff2 = xMax - xoffset ;
00535         if (xoff1 >= 0 && xoff2 >= 0 )
00536         {
00537           // X limits cross max/min y => no change
00538           //
00539           pMin = yMin ;
00540           pMax = yMax ;
00541         } 
00542         else
00543         {
00544           // X limits don't cross max/min y => compute max delta y,
00545           // hence new mins/maxs
00546           //
00547           RTorus=fRmax+fRtor;
00548           delta   = RTorus*RTorus - xoff1*xoff1;
00549           diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
00550           delta   = RTorus*RTorus - xoff2*xoff2;
00551           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
00552           maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
00553           newMin  = yoffset - maxDiff ;
00554           newMax  = yoffset + maxDiff ;
00555           pMin    = (newMin < yMin) ? yMin : newMin ;
00556           pMax    = (newMax > yMax) ? yMax : newMax ;
00557         }
00558         break;
00559 
00560       case kZAxis:
00561         pMin=zMin;
00562         pMax=zMax;
00563         break;
00564 
00565       default:
00566         break;
00567     }
00568     pMin -= kCarTolerance ;
00569     pMax += kCarTolerance ;
00570 
00571     return true;
00572   }
00573   else
00574   {
00575     G4int i, noEntries, noBetweenSections4 ;
00576     G4bool existsAfterClip = false ;
00577 
00578     // Calculate rotated vertex coordinates
00579 
00580     G4ThreeVectorList *vertices ;
00581     G4int noPolygonVertices ;  // will be 4 
00582     vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ;
00583 
00584     pMin = +kInfinity ;
00585     pMax = -kInfinity ;
00586 
00587     noEntries          = vertices->size() ;
00588     noBetweenSections4 = noEntries - noPolygonVertices ;
00589       
00590     for (i=0;i<noEntries;i+=noPolygonVertices)
00591     {
00592       ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
00593     }    
00594     for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
00595     {
00596       ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
00597     }
00598     if (pMin!=kInfinity||pMax!=-kInfinity)
00599     {
00600       existsAfterClip = true ; // Add 2*tolerance to avoid precision troubles
00601       pMin           -= kCarTolerance ;
00602       pMax           += kCarTolerance ;
00603     }
00604     else
00605     {
00606       // Check for case where completely enveloping clipping volume
00607       // If point inside then we are confident that the solid completely
00608       // envelopes the clipping volume. Hence set min/max extents according
00609       // to clipping volume extents along the specified axis.
00610 
00611       G4ThreeVector clipCentre(
00612           (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00613           (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00614           (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5  ) ;
00615         
00616       if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
00617       {
00618         existsAfterClip = true ;
00619         pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
00620         pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
00621       }
00622     }
00623     delete vertices;
00624     return existsAfterClip;
00625   }
00626 }

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

Reimplemented from G4VSolid.

Definition at line 1703 of file G4Torus.cc.

References G4Torus().

01704 {
01705   return new G4Torus(*this);
01706 }

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

Reimplemented from G4VSolid.

Definition at line 237 of file G4Torus.cc.

References G4VPVParameterisation::ComputeDimensions().

00240 {
00241   p->ComputeDimensions(*this,n,pRep);
00242 }

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

Reimplemented from G4VSolid.

Definition at line 1793 of file G4Torus.cc.

01794 {
01795   G4NURBS* pNURBS;
01796   if (fRmin != 0) 
01797   {
01798     if (fDPhi >= twopi) 
01799     {
01800       pNURBS = new G4NURBStube(fRmin, fRmax, fRtor);
01801     }
01802     else 
01803     {
01804       pNURBS = new G4NURBStubesector(fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi);
01805     }
01806   }
01807   else 
01808   {
01809     if (fDPhi >= twopi) 
01810     {
01811       pNURBS = new G4NURBScylinder (fRmax, fRtor);
01812     }
01813     else 
01814     {
01815       const G4double epsilon = 1.e-4; // Cylinder sector not yet available!
01816       pNURBS = new G4NURBStubesector (epsilon, fRmax, fRtor,
01817                                       fSPhi, fSPhi + fDPhi);
01818     }
01819   }
01820   return pNURBS;
01821 }

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

Reimplemented from G4VSolid.

Definition at line 1788 of file G4Torus.cc.

01789 {
01790   return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
01791 }

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

Implements G4VSolid.

Definition at line 1783 of file G4Torus.cc.

References G4VGraphicsScene::AddSolid().

01784 {
01785   scene.AddSolid (*this);
01786 }

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

Implements G4VSolid.

Definition at line 1139 of file G4Torus.cc.

01140 {
01141   G4double safe=0.0, safe1, safe2 ;
01142   G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
01143   G4double rho2, rho, pt2, pt ;
01144     
01145   rho2 = p.x()*p.x() + p.y()*p.y() ;
01146   rho  = std::sqrt(rho2) ;
01147   pt2  = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
01148   pt   = std::sqrt(pt2) ;
01149 
01150   safe1 = fRmin - pt ;
01151   safe2 = pt - fRmax ;
01152 
01153   if (safe1 > safe2)  { safe = safe1; }
01154   else                { safe = safe2; }
01155 
01156   if ( fDPhi < twopi && rho )
01157   {
01158     phiC    = fSPhi + fDPhi*0.5 ;
01159     cosPhiC = std::cos(phiC) ;
01160     sinPhiC = std::sin(phiC) ;
01161     cosPsi  = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
01162 
01163     if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point
01164     {                                  // Point lies outside phi range
01165       if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
01166       {
01167         safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
01168       }
01169       else
01170       {
01171         ePhi    = fSPhi + fDPhi ;
01172         safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
01173       }
01174       if (safePhi > safe)  { safe = safePhi ; }
01175     }
01176   }
01177   if (safe < 0 )  { safe = 0 ; }
01178   return safe;
01179 }

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

Implements G4VSolid.

Definition at line 987 of file G4Torus.cc.

References G4VSolid::kCarTolerance.

Referenced by SurfaceNormal().

00989 {
00990 
00991   G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value
00992 
00993   G4double  sd[4] ;
00994 
00995   // Precalculated trig for phi intersections - used by r,z intersections to
00996   //                                            check validity
00997 
00998   G4bool seg;        // true if segmented
00999   G4double hDPhi;    // half dphi
01000   G4double cPhi,sinCPhi=0.,cosCPhi=0.;  // central phi
01001 
01002   G4double tolORMin2;  // `generous' radii squared
01003   G4double tolORMax2;
01004 
01005   static const G4double halfCarTolerance = 0.5*kCarTolerance;
01006  
01007   G4double Dist,xi,yi,zi,rhoi2,it2; // Intersection point variables
01008 
01009   G4double Comp;
01010   G4double cosSPhi,sinSPhi;       // Trig for phi start intersect
01011   G4double ePhi,cosEPhi,sinEPhi;  // for phi end intersect
01012 
01013   // Set phi divided flag and precalcs
01014   //
01015   if ( fDPhi < twopi )
01016   {
01017     seg        = true ;
01018     hDPhi      = 0.5*fDPhi ;    // half delta phi
01019     cPhi       = fSPhi + hDPhi ;
01020     sinCPhi    = std::sin(cPhi) ;
01021     cosCPhi    = std::cos(cPhi) ;
01022   }
01023   else
01024   {
01025     seg = false ;
01026   }
01027 
01028   if (fRmin > fRminTolerance) // Calculate tolerant rmin and rmax
01029   {
01030     tolORMin2 = (fRmin - fRminTolerance)*(fRmin - fRminTolerance) ;
01031   }
01032   else
01033   {
01034     tolORMin2 = 0 ;
01035   }
01036   tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax + fRmaxTolerance) ;
01037 
01038   // Intersection with Rmax (possible return) and Rmin (must also check phi)
01039 
01040   G4double Rtor2 = fRtor*fRtor ;
01041 
01042   snxt = SolveNumericJT(p,v,fRmax,true);
01043 
01044   if (fRmin)  // Possible Rmin intersection
01045   {
01046     sd[0] = SolveNumericJT(p,v,fRmin,true);
01047     if ( sd[0] < snxt )  { snxt = sd[0] ; }
01048   }
01049 
01050   //
01051   // Phi segment intersection
01052   //
01053   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
01054   //
01055   // o NOTE: Large duplication of code between sphi & ephi checks
01056   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
01057   //            intersection check <=0 -> >=0
01058   //         -> use some form of loop Construct ?
01059 
01060   if (seg)
01061   {
01062     sinSPhi = std::sin(fSPhi) ; // First phi surface ('S'tarting phi)
01063     cosSPhi = std::cos(fSPhi) ;
01064     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;  // Component in outwards
01065                                                // normal direction
01066     if (Comp < 0 )
01067     {
01068       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
01069 
01070       if (Dist < halfCarTolerance)
01071       {
01072         sphi = Dist/Comp ;
01073         if (sphi < snxt)
01074         {
01075           if ( sphi < 0 )  { sphi = 0 ; }
01076 
01077           xi    = p.x() + sphi*v.x() ;
01078           yi    = p.y() + sphi*v.y() ;
01079           zi    = p.z() + sphi*v.z() ;
01080           rhoi2 = xi*xi + yi*yi ;
01081           it2   = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
01082 
01083           if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
01084           {
01085             // r intersection is good - check intersecting
01086             // with correct half-plane
01087             //
01088             if ((yi*cosCPhi-xi*sinCPhi)<=0)  { snxt=sphi; }
01089           }
01090         }
01091       }
01092     }
01093     ePhi=fSPhi+fDPhi;    // Second phi surface ('E'nding phi)
01094     sinEPhi=std::sin(ePhi);
01095     cosEPhi=std::cos(ePhi);
01096     Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
01097         
01098     if ( Comp < 0 )   // Component in outwards normal dirn
01099     {
01100       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
01101 
01102       if (Dist < halfCarTolerance )
01103       {
01104         sphi = Dist/Comp ;
01105 
01106         if (sphi < snxt )
01107         {
01108           if (sphi < 0 )  { sphi = 0 ; }
01109        
01110           xi    = p.x() + sphi*v.x() ;
01111           yi    = p.y() + sphi*v.y() ;
01112           zi    = p.z() + sphi*v.z() ;
01113           rhoi2 = xi*xi + yi*yi ;
01114           it2   = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
01115 
01116           if (it2 >= tolORMin2 && it2 <= tolORMax2)
01117           {
01118             // z and r intersections good - check intersecting
01119             // with correct half-plane
01120             //
01121             if ((yi*cosCPhi-xi*sinCPhi)>=0)  { snxt=sphi; }
01122           }    
01123         }
01124       }
01125     }
01126   }
01127   if(snxt < halfCarTolerance)  { snxt = 0.0 ; }
01128 
01129   return snxt ;
01130 }

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

Implements G4VSolid.

Definition at line 1542 of file G4Torus.cc.

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

01543 {
01544   G4double safe=0.0,safeR1,safeR2;
01545   G4double rho2,rho,pt2,pt ;
01546   G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
01547   rho2 = p.x()*p.x() + p.y()*p.y() ;
01548   rho  = std::sqrt(rho2) ;
01549   pt2  = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
01550   pt   = std::sqrt(pt2) ;
01551 
01552 #ifdef G4CSGDEBUG
01553   if( Inside(p) == kOutside )
01554   {
01555      G4int oldprc = G4cout.precision(16) ;
01556      G4cout << G4endl ;
01557      DumpInfo();
01558      G4cout << "Position:"  << G4endl << G4endl ;
01559      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
01560      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
01561      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
01562      G4cout.precision(oldprc);
01563      G4Exception("G4Torus::DistanceToOut(p)", "GeomSolids1002",
01564                  JustWarning, "Point p is outside !?" );
01565   }
01566 #endif
01567 
01568   if (fRmin)
01569   {
01570     safeR1 = pt - fRmin ;
01571     safeR2 = fRmax - pt ;
01572 
01573     if (safeR1 < safeR2)  { safe = safeR1 ; }
01574     else                  { safe = safeR2 ; }
01575   }
01576   else
01577   {
01578     safe = fRmax - pt ;
01579   }  
01580 
01581   // Check if phi divided, Calc distances closest phi plane
01582   //
01583   if (fDPhi<twopi) // Above/below central phi of Torus?
01584   {
01585     phiC    = fSPhi + fDPhi*0.5 ;
01586     cosPhiC = std::cos(phiC) ;
01587     sinPhiC = std::sin(phiC) ;
01588 
01589     if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
01590     {
01591       safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
01592     }
01593     else
01594     {
01595       ePhi    = fSPhi + fDPhi ;
01596       safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
01597     }
01598     if (safePhi < safe)  { safe = safePhi ; }
01599   }
01600   if (safe < 0)  { safe = 0 ; }
01601   return safe ;  
01602 }

G4double G4Torus::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 1187 of file G4Torus.cc.

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

Referenced by SurfaceNormal().

01192 {
01193   ESide    side = kNull, sidephi = kNull ;
01194   G4double snxt = kInfinity, sphi, sd[4] ;
01195 
01196   static const G4double halfCarTolerance = 0.5*kCarTolerance;
01197   static const G4double halfAngTolerance = 0.5*kAngTolerance;
01198 
01199   // Vars for phi intersection
01200   //
01201   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
01202   G4double cPhi, sinCPhi, cosCPhi ;
01203   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
01204 
01205   // Radial Intersections Defenitions & General Precals
01206 
01208 
01209 #if 1
01210 
01211   // This is the version with the calculation of CalcNorm = true 
01212   // To be done: Check the precision of this calculation.
01213   // If you want return always validNorm = false, then take the version below
01214   
01215   G4double rho2  = p.x()*p.x()+p.y()*p.y();
01216   G4double rho   = std::sqrt(rho2) ;
01217 
01218   G4double pt2   = rho2 + p.z()*p.z() + fRtor * (fRtor - 2.0*rho);
01219     // Regroup for slightly better FP accuracy
01220 
01221   if( pt2 < 0.0)
01222   {
01223      pt2= std::fabs( pt2 );
01224   }
01225      
01226   G4double pt    = std::sqrt(pt2) ;
01227 
01228   G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
01229 
01230   G4double tolRMax = fRmax - fRmaxTolerance ;
01231    
01232   G4double vDotNmax   = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
01233   G4double pDotxyNmax = (1 - fRtor/rho) ;
01234 
01235   if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) )
01236   {
01237     // On tolerant boundary & heading outwards (or perpendicular to) outer
01238     // radial surface -> leaving immediately with *n for really convex part
01239     // only
01240       
01241     if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) ) 
01242     {
01243       *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
01244                           p.y()*(1 - fRtor/rho)/pt,
01245                           p.z()/pt                  ) ;
01246       *validNorm = true ;
01247     }
01248      
01249     return snxt = 0 ; // Leaving by Rmax immediately
01250   }
01251   
01252   snxt = SolveNumericJT(p,v,fRmax,false);  
01253   side = kRMax ;
01254 
01255   // rmin
01256 
01257   if ( fRmin )
01258   {
01259     G4double tolRMin = fRmin + fRminTolerance ;
01260 
01261     if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) )
01262     {
01263       if (calcNorm)  { *validNorm = false ; } // Concave surface of the torus
01264       return  snxt = 0 ;                      // Leaving by Rmin immediately
01265     }
01266     
01267     sd[0] = SolveNumericJT(p,v,fRmin,false);
01268     if ( sd[0] < snxt )
01269     {
01270       snxt = sd[0] ;
01271       side = kRMin ;
01272     }
01273   }
01274 
01275 #else
01276 
01277   // this is the "conservative" version which return always validnorm = false
01278   // NOTE: using this version the unit test testG4Torus will break
01279 
01280   snxt = SolveNumericJT(p,v,fRmax,false);  
01281   side = kRMax ;
01282 
01283   if ( fRmin )
01284   {
01285     sd[0] = SolveNumericJT(p,v,fRmin,false);
01286     if ( sd[0] < snxt )
01287     {
01288       snxt = sd[0] ;
01289       side = kRMin ;
01290     }
01291   }
01292 
01293   if ( calcNorm && (snxt == 0.0) )
01294   {
01295     *validNorm = false ;    // Leaving solid, but possible re-intersection
01296     return snxt  ;
01297   }
01298 
01299 #endif
01300   
01301   if (fDPhi < twopi)  // Phi Intersections
01302   {
01303     sinSPhi = std::sin(fSPhi) ;
01304     cosSPhi = std::cos(fSPhi) ;
01305     ePhi    = fSPhi + fDPhi ;
01306     sinEPhi = std::sin(ePhi) ;
01307     cosEPhi = std::cos(ePhi) ;
01308     cPhi    = fSPhi + fDPhi*0.5 ;
01309     sinCPhi = std::sin(cPhi) ;
01310     cosCPhi = std::cos(cPhi) ;
01311     
01312     // angle calculation with correction 
01313     // of difference in domain of atan2 and Sphi
01314     //
01315     vphi = std::atan2(v.y(),v.x()) ;
01316      
01317     if ( vphi < fSPhi - halfAngTolerance  )    { vphi += twopi; }
01318     else if ( vphi > ePhi + halfAngTolerance ) { vphi -= twopi; }
01319 
01320     if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
01321     {
01322       pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
01323       pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
01324 
01325       // Comp -ve when in direction of outwards normal
01326       //
01327       compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
01328       compE   = sinEPhi*v.x() - cosEPhi*v.y() ;
01329       sidephi = kNull ;
01330      
01331       if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
01332                             && (pDistE <= halfCarTolerance) ) )
01333        || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
01334                             && (pDistE >  halfCarTolerance) ) )  )
01335       {
01336         // Inside both phi *full* planes
01337 
01338         if ( compS < 0 )
01339         {
01340           sphi = pDistS/compS ;
01341             
01342           if (sphi >= -halfCarTolerance)
01343           {
01344             xi = p.x() + sphi*v.x() ;
01345             yi = p.y() + sphi*v.y() ;
01346               
01347             // Check intersecting with correct half-plane
01348             // (if not -> no intersect)
01349             //
01350             if ( (std::fabs(xi)<=kCarTolerance)
01351               && (std::fabs(yi)<=kCarTolerance) )
01352             {
01353               sidephi = kSPhi;
01354               if ( ((fSPhi-halfAngTolerance)<=vphi)
01355                 && ((ePhi+halfAngTolerance)>=vphi) )
01356               {
01357                 sphi = kInfinity;
01358               }
01359             }
01360             else if ( yi*cosCPhi-xi*sinCPhi >=0 )
01361             {
01362               sphi = kInfinity ;
01363             }
01364             else
01365             {
01366               sidephi = kSPhi ;
01367             }       
01368           }
01369           else
01370           {
01371             sphi = kInfinity ;
01372           }
01373         }
01374         else
01375         {
01376           sphi = kInfinity ;
01377         }
01378 
01379         if ( compE < 0 )
01380         {
01381           sphi2 = pDistE/compE ;
01382             
01383           // Only check further if < starting phi intersection
01384           //
01385           if ( (sphi2 > -kCarTolerance) && (sphi2 < sphi) )
01386           {
01387             xi = p.x() + sphi2*v.x() ;
01388             yi = p.y() + sphi2*v.y() ;
01389               
01390             if ( (std::fabs(xi)<=kCarTolerance)
01391               && (std::fabs(yi)<=kCarTolerance) )
01392             {
01393               // Leaving via ending phi
01394               //
01395               if( !( (fSPhi-halfAngTolerance <= vphi)
01396                   && (ePhi+halfAngTolerance >= vphi) ) )
01397               {
01398                 sidephi = kEPhi ;
01399                 sphi = sphi2;
01400               }
01401             } 
01402             else    // Check intersecting with correct half-plane 
01403             {
01404               if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
01405               {
01406                 // Leaving via ending phi
01407                 //
01408                 sidephi = kEPhi ;
01409                 sphi = sphi2;
01410                
01411               }
01412             }
01413           }
01414         }
01415       }
01416       else
01417       {
01418         sphi = kInfinity ;
01419       }
01420     } 
01421     else
01422     {
01423       // On z axis + travel not || to z axis -> if phi of vector direction
01424       // within phi of shape, Step limited by rmax, else Step =0
01425 
01426       vphi = std::atan2(v.y(),v.x());
01427  
01428       if ( ( fSPhi-halfAngTolerance <= vphi ) && 
01429            ( vphi <= ( ePhi+halfAngTolerance ) ) )
01430       {
01431         sphi = kInfinity;
01432       }
01433       else
01434       {
01435         sidephi = kSPhi ; // arbitrary 
01436         sphi=0;
01437       }
01438     }
01439 
01440     // Order intersections
01441 
01442     if (sphi<snxt)
01443     {
01444       snxt=sphi;
01445       side=sidephi;
01446     }     
01447   }
01448 
01449   G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
01450   // Note: by numerical computation we know where the ray hits the torus
01451   // So I propose to return the side where the ray hits
01452 
01453   if (calcNorm)
01454   {
01455     switch(side)
01456     {
01457       case kRMax:                     // n is unit vector 
01458         xi    = p.x() + snxt*v.x() ;
01459         yi    =p.y() + snxt*v.y() ;
01460         zi    = p.z() + snxt*v.z() ;
01461         rhoi2 = xi*xi + yi*yi ;
01462         rhoi  = std::sqrt(rhoi2) ;
01463         it2   = std::fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ;
01464         it    = std::sqrt(it2) ;
01465         iDotxyNmax = (1-fRtor/rhoi) ;
01466         if(iDotxyNmax >= -2.*fRmaxTolerance) // really convex part of Rmax
01467         {                       
01468           *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
01469                               yi*(1-fRtor/rhoi)/it,
01470                               zi/it                 ) ;
01471           *validNorm = true ;
01472         }
01473         else
01474         {
01475           *validNorm = false ; // concave-convex part of Rmax
01476         }
01477         break ;
01478 
01479       case kRMin:
01480         *validNorm = false ;  // Rmin is concave or concave-convex
01481         break;
01482 
01483       case kSPhi:
01484         if (fDPhi <= pi )
01485         {
01486           *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
01487           *validNorm=true;
01488         }
01489         else
01490         {
01491           *validNorm = false ;
01492         }
01493         break ;
01494 
01495       case kEPhi:
01496         if (fDPhi <= pi)
01497         {
01498           *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
01499           *validNorm=true;
01500         }
01501         else
01502         {
01503           *validNorm = false ;
01504         }
01505         break;
01506 
01507       default:
01508 
01509         // It seems we go here from time to time ...
01510 
01511         G4cout << G4endl;
01512         DumpInfo();
01513         std::ostringstream message;
01514         G4int oldprc = message.precision(16);
01515         message << "Undefined side for valid surface normal to solid."
01516                 << G4endl
01517                 << "Position:"  << G4endl << G4endl
01518                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
01519                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
01520                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
01521                 << "Direction:" << G4endl << G4endl
01522                 << "v.x() = "   << v.x() << G4endl
01523                 << "v.y() = "   << v.y() << G4endl
01524                 << "v.z() = "   << v.z() << G4endl << G4endl
01525                 << "Proposed distance :" << G4endl << G4endl
01526                 << "snxt = " << snxt/mm << " mm" << G4endl;
01527         message.precision(oldprc);
01528         G4Exception("G4Torus::DistanceToOut(p,v,..)",
01529                     "GeomSolids1002",JustWarning, message);
01530         break;
01531     }
01532   }
01533   if ( snxt<halfCarTolerance )  { snxt=0 ; }
01534 
01535   return snxt;
01536 }

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

Reimplemented from G4VSolid.

Definition at line 68 of file G4Torus.icc.

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

00069 {
00070   if(fCubicVolume != 0.) {;}
00071   else  { fCubicVolume = fDPhi*CLHEP::pi*fRtor*(fRmax*fRmax-fRmin*fRmin); }
00072   return fCubicVolume;
00073 }

G4double G4Torus::GetDPhi (  )  const [inline]

Definition at line 62 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

00063 {
00064   return fDPhi ;
00065 }

G4GeometryType G4Torus::GetEntityType (  )  const [virtual]

Implements G4VSolid.

Definition at line 1694 of file G4Torus.cc.

01695 {
01696   return G4String("G4Torus");
01697 }

G4ThreeVector G4Torus::GetPointOnSurface (  )  const [virtual]

Reimplemented from G4VSolid.

Definition at line 1735 of file G4Torus.cc.

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

01736 {
01737   G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
01738    
01739   phi   = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
01740   theta = RandFlat::shoot(0.,twopi);
01741   
01742   cosu   = std::cos(phi);    sinu = std::sin(phi);
01743   cosv   = std::cos(theta);  sinv = std::sin(theta); 
01744 
01745   // compute the areas
01746 
01747   aOut   = (fDPhi)*twopi*fRtor*fRmax;
01748   aIn    = (fDPhi)*twopi*fRtor*fRmin;
01749   aSide  = pi*(fRmax*fRmax-fRmin*fRmin);
01750   
01751   if ((fSPhi == 0) && (fDPhi == twopi)){ aSide = 0; }
01752   chose = RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
01753 
01754   if(chose < aOut)
01755   {
01756     return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
01757                           (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
01758   }
01759   else if( (chose >= aOut) && (chose < aOut + aIn) )
01760   {
01761     return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
01762                           (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
01763   }
01764   else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
01765   {
01766     rRand = GetRadiusInRing(fRmin,fRmax);
01767     return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
01768                           (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
01769   }
01770   else
01771   {   
01772     rRand = GetRadiusInRing(fRmin,fRmax);
01773     return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
01774                           (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi), 
01775                           rRand*sinv);
01776    }
01777 }

G4double G4Torus::GetRmax (  )  const [inline]

Definition at line 44 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

00045 {
00046   return fRmax ;
00047 }

G4double G4Torus::GetRmin (  )  const [inline]

Definition at line 38 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

00039 {
00040   return fRmin ;
00041 }

G4double G4Torus::GetRtor (  )  const [inline]

Definition at line 50 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

00051 {
00052   return fRtor ;
00053 }

G4double G4Torus::GetSPhi (  )  const [inline]

Definition at line 56 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

00057 {
00058   return fSPhi ;
00059 }

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

Reimplemented from G4VSolid.

Definition at line 76 of file G4Torus.icc.

References G4CSGSolid::fSurfaceArea.

00077 {
00078   if(fSurfaceArea != 0.) {;}
00079   else
00080   { 
00081     fSurfaceArea = fDPhi*CLHEP::twopi*fRtor*(fRmax+fRmin);
00082     if(fDPhi < CLHEP::twopi)
00083     {
00084       fSurfaceArea = fSurfaceArea + CLHEP::twopi*(fRmax*fRmax-fRmin*fRmin);
00085     } 
00086   }
00087   return fSurfaceArea;
00088 }

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

Implements G4VSolid.

Definition at line 632 of file G4Torus.cc.

References kInside, kOutside, and kSurface.

Referenced by CalculateExtent(), DistanceToOut(), and SurfaceNormal().

00633 {
00634   G4double r2, pt2, pPhi, tolRMin, tolRMax ;
00635 
00636   EInside in = kOutside ;
00637   static const G4double halfAngTolerance = 0.5*kAngTolerance;
00638 
00639                                                // General precals
00640   r2  = p.x()*p.x() + p.y()*p.y() ;
00641   pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*std::sqrt(r2) ;
00642 
00643   if (fRmin) tolRMin = fRmin + fRminTolerance ;
00644   else       tolRMin = 0 ;
00645 
00646   tolRMax = fRmax - fRmaxTolerance;
00647       
00648   if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
00649   {
00650     if ( fDPhi == twopi || pt2 == 0 )  // on torus swept axis
00651     {
00652       in = kInside ;
00653     }
00654     else
00655     {
00656       // Try inner tolerant phi boundaries (=>inside)
00657       // if not inside, try outer tolerant phi boundaries
00658 
00659       pPhi = std::atan2(p.y(),p.x()) ;
00660 
00661       if ( pPhi < -halfAngTolerance )  { pPhi += twopi ; }  // 0<=pPhi<2pi
00662       if ( fSPhi >= 0 )
00663       {
00664         if ( (std::fabs(pPhi) < halfAngTolerance)
00665             && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
00666         { 
00667             pPhi += twopi ; // 0 <= pPhi < 2pi
00668         }
00669         if ( (pPhi >= fSPhi + halfAngTolerance)
00670             && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
00671         {
00672           in = kInside ;
00673         }
00674           else if ( (pPhi >= fSPhi - halfAngTolerance)
00675                  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
00676         {
00677           in = kSurface ;
00678         }
00679       }
00680       else  // fSPhi < 0
00681       {
00682           if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
00683             && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) )  {;}
00684           else
00685           {
00686             in = kSurface ;
00687           }
00688       }
00689     }
00690   }
00691   else   // Try generous boundaries
00692   {
00693     tolRMin = fRmin - fRminTolerance ;
00694     tolRMax = fRmax + fRmaxTolerance ;
00695 
00696     if (tolRMin < 0 )  { tolRMin = 0 ; }
00697 
00698     if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
00699     {
00700       if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis
00701       {
00702         in = kSurface ;
00703       }
00704       else // Try outer tolerant phi boundaries only
00705       {
00706         pPhi = std::atan2(p.y(),p.x()) ;
00707 
00708         if ( pPhi < -halfAngTolerance )  { pPhi += twopi ; }  // 0<=pPhi<2pi
00709         if ( fSPhi >= 0 )
00710         {
00711           if ( (std::fabs(pPhi) < halfAngTolerance)
00712             && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
00713           { 
00714             pPhi += twopi ; // 0 <= pPhi < 2pi
00715           }
00716           if ( (pPhi >= fSPhi - halfAngTolerance)
00717             && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
00718           {
00719             in = kSurface;
00720           }
00721         }
00722         else  // fSPhi < 0
00723         {
00724           if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
00725             && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) )  {;}
00726           else
00727           {
00728             in = kSurface ;
00729           }
00730         }
00731       }
00732     }
00733   }
00734   return in ;
00735 }

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

Definition at line 212 of file G4Torus.cc.

References fDPhi, fRmax, fRmaxTolerance, fRmin, fRminTolerance, fRtor, fSPhi, kAngTolerance, kRadTolerance, and G4CSGSolid::operator=().

00213 {
00214    // Check assignment to self
00215    //
00216    if (this == &rhs)  { return *this; }
00217 
00218    // Copy base class data
00219    //
00220    G4CSGSolid::operator=(rhs);
00221 
00222    // Copy data
00223    //
00224    fRmin = rhs.fRmin; fRmax = rhs.fRmax;
00225    fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
00226    fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
00227    kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance;
00228 
00229    return *this;
00230 }

void G4Torus::SetAllParameters ( G4double  pRmin,
G4double  pRmax,
G4double  pRtor,
G4double  pSPhi,
G4double  pDPhi 
)

Definition at line 100 of file G4Torus.cc.

References FatalException, G4CSGSolid::fCubicVolume, G4CSGSolid::fpPolyhedron, G4CSGSolid::fSurfaceArea, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), and G4VSolid::kCarTolerance.

Referenced by G4Torus().

00105 {
00106   const G4double fEpsilon = 4.e-11;  // relative tolerance of radii
00107 
00108   fCubicVolume = 0.;
00109   fSurfaceArea = 0.;
00110   fpPolyhedron = 0;
00111 
00112   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00113   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00114   
00115   if ( pRtor >= pRmax+1.e3*kCarTolerance )  // Check swept radius, as in G4Cons
00116   {
00117     fRtor = pRtor ;
00118   }
00119   else
00120   {
00121     std::ostringstream message;
00122     message << "Invalid swept radius for Solid: " << GetName() << G4endl
00123             << "        pRtor = " << pRtor << ", pRmax = " << pRmax;
00124     G4Exception("G4Torus::SetAllParameters()",
00125                 "GeomSolids0002", FatalException, message);
00126   }
00127 
00128   // Check radii, as in G4Cons
00129   //
00130   if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
00131   {
00132     if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; }
00133     else                             { fRmin = 0.0   ; }
00134     fRmax = pRmax ;
00135   }
00136   else
00137   {
00138     std::ostringstream message;
00139     message << "Invalid values of radii for Solid: " << GetName() << G4endl
00140             << "        pRmin = " << pRmin << ", pRmax = " << pRmax;
00141     G4Exception("G4Torus::SetAllParameters()",
00142                 "GeomSolids0002", FatalException, message);
00143   }
00144 
00145   // Relative tolerances
00146   //
00147   fRminTolerance = (fRmin)
00148                  ? 0.5*std::max( kRadTolerance, fEpsilon*(fRtor-fRmin )) : 0;
00149   fRmaxTolerance = 0.5*std::max( kRadTolerance, fEpsilon*(fRtor+fRmax) );
00150 
00151   // Check angles
00152   //
00153   if ( pDPhi >= twopi )  { fDPhi = twopi ; }
00154   else
00155   {
00156     if (pDPhi > 0)       { fDPhi = pDPhi ; }
00157     else
00158     {
00159       std::ostringstream message;
00160       message << "Invalid Z delta-Phi for Solid: " << GetName() << G4endl
00161               << "        pDPhi = " << pDPhi;
00162       G4Exception("G4Torus::SetAllParameters()",
00163                   "GeomSolids0002", FatalException, message);
00164     }
00165   }
00166   
00167   // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
00168   //
00169   fSPhi = pSPhi;
00170 
00171   if (fSPhi < 0)  { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; }
00172   else            { fSPhi = std::fmod(fSPhi,twopi) ; }
00173 
00174   if (fSPhi+fDPhi > twopi)  { fSPhi-=twopi ; }
00175 }

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

Reimplemented from G4CSGSolid.

Definition at line 1712 of file G4Torus.cc.

References G4VSolid::GetName().

01713 {
01714   G4int oldprc = os.precision(16);
01715   os << "-----------------------------------------------------------\n"
01716      << "    *** Dump for solid - " << GetName() << " ***\n"
01717      << "    ===================================================\n"
01718      << " Solid type: G4Torus\n"
01719      << " Parameters: \n"
01720      << "    inner radius: " << fRmin/mm << " mm \n"
01721      << "    outer radius: " << fRmax/mm << " mm \n"
01722      << "    swept radius: " << fRtor/mm << " mm \n"
01723      << "    starting phi: " << fSPhi/degree << " degrees \n"
01724      << "    delta phi   : " << fDPhi/degree << " degrees \n"
01725      << "-----------------------------------------------------------\n";
01726   os.precision(oldprc);
01727 
01728   return os;
01729 }

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

Implements G4VSolid.

Definition at line 743 of file G4Torus.cc.

References DistanceToIn(), DistanceToOut(), G4endl, G4Exception(), Inside(), JustWarning, G4VSolid::kCarTolerance, kInside, kOutside, and kSurface.

00744 {
00745   G4int noSurfaces = 0;  
00746   G4double rho2, rho, pt2, pt, pPhi;
00747   G4double distRMin = kInfinity;
00748   G4double distSPhi = kInfinity, distEPhi = kInfinity;
00749 
00750   // To cope with precision loss
00751   //
00752   const G4double delta = std::max(10.0*kCarTolerance,
00753                                   1.0e-8*(fRtor+fRmax));
00754   const G4double dAngle = 10.0*kAngTolerance;
00755 
00756   G4ThreeVector nR, nPs, nPe;
00757   G4ThreeVector norm, sumnorm(0.,0.,0.);
00758 
00759   rho2 = p.x()*p.x() + p.y()*p.y();
00760   rho = std::sqrt(rho2);
00761   pt2 = rho2+p.z()*p.z() +fRtor * (fRtor-2*rho);
00762   pt2 = std::max(pt2, 0.0); // std::fabs(pt2);
00763   pt = std::sqrt(pt2) ;
00764 
00765   G4double  distRMax = std::fabs(pt - fRmax);
00766   if(fRmin) distRMin = std::fabs(pt - fRmin);
00767 
00768   if( rho > delta && pt != 0.0 )
00769   {
00770     G4double redFactor= (rho-fRtor)/rho;
00771     nR = G4ThreeVector( p.x()*redFactor,  // p.x()*(1.-fRtor/rho),
00772                         p.y()*redFactor,  // p.y()*(1.-fRtor/rho),
00773                         p.z()          );
00774     nR *= 1.0/pt;
00775   }
00776 
00777   if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
00778   {
00779     if ( rho )
00780     {
00781       pPhi = std::atan2(p.y(),p.x());
00782 
00783       if(pPhi < fSPhi-delta)            { pPhi += twopi; }
00784       else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
00785 
00786       distSPhi = std::fabs( pPhi - fSPhi );
00787       distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
00788     }
00789     nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
00790     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
00791   } 
00792   if( distRMax <= delta )
00793   {
00794     noSurfaces ++;
00795     sumnorm += nR;
00796   }
00797   else if( fRmin && (distRMin <= delta) ) // Must not be on both Outer and Inner
00798   {
00799     noSurfaces ++;
00800     sumnorm -= nR;
00801   }
00802 
00803   //  To be on one of the 'phi' surfaces,
00804   //  it must be within the 'tube' - with tolerance
00805 
00806   if( (fDPhi < twopi) && (fRmin-delta <= pt) && (pt <= (fRmax+delta)) )
00807   {
00808     if (distSPhi <= dAngle)
00809     {
00810       noSurfaces ++;
00811       sumnorm += nPs;
00812     }
00813     if (distEPhi <= dAngle) 
00814     {
00815       noSurfaces ++;
00816       sumnorm += nPe;
00817     }
00818   }
00819   if ( noSurfaces == 0 )
00820   {
00821 #ifdef G4CSGDEBUG
00822      G4ExceptionDescription ed;
00823      ed.precision(16);
00824 
00825      EInside  inIt= Inside( p );
00826      
00827      if( inIt != kSurface )
00828      {
00829         ed << " ERROR>  Surface Normal was called for Torus,"
00830            << " with point not on surface." << G4endl;
00831      }
00832      else
00833      {
00834         ed << " ERROR>  Surface Normal has not found a surface, "
00835            << " despite the point being on the surface. " <<G4endl;
00836      }
00837 
00838      if( inIt != kInside)
00839      {
00840          ed << " Safety (Dist To In)  = " << DistanceToIn(p) << G4endl;
00841      }
00842      if( inIt != kOutside)
00843      {
00844          ed << " Safety (Dist to Out) = " << DistanceToOut(p) << G4endl;
00845      }
00846      ed << " Coordinates of point : " << p << G4endl;
00847      ed << " Parameters  of solid : " << G4endl << *this << G4endl;
00848 
00849      if( inIt == kSurface )
00850      {
00851         G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
00852                     JustWarning, ed,
00853                     "Failing to find normal, even though point is on surface!" );
00854      }
00855      else
00856      {
00857         static const char* NameInside[3]= { "Inside", "Surface", "Outside" };
00858         ed << "  The point is " << NameInside[inIt] << " the solid. "<< G4endl;
00859         G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
00860                     JustWarning, ed, "Point p is not on surface !?" );
00861      }
00862 #endif
00863      norm = ApproxSurfaceNormal(p);
00864   }
00865   else if ( noSurfaces == 1 )  { norm = sumnorm; }
00866   else                         { norm = sumnorm.unit(); }
00867 
00868   // G4cout << "G4Torus::SurfaceNormal p= " << p << " returns norm= " << norm << G4endl;
00869 
00870   return norm ;
00871 }


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