G4Cons.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4Cons.cc 69788 2013-05-15 12:06:57Z gcosmo $
00028 // GEANT4 tag $Name: $
00029 //
00030 //
00031 // class G4Cons
00032 //
00033 // Implementation for G4Cons class
00034 //
00035 // History:
00036 //
00037 // 05.04.12 M.Kelsey:   GetPointOnSurface() throw flat in sqrt(r)
00038 // 12.10.09 T.Nikitina: Added to DistanceToIn(p,v) check on the direction in
00039 //                      case of point on surface
00040 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
00041 // 13.09.96 V.Grichine: Review and final modifications
00042 // ~1994    P.Kent: Created, as main part of the geometry prototype
00043 // --------------------------------------------------------------------
00044 
00045 #include "G4Cons.hh"
00046 
00047 #include "G4VoxelLimits.hh"
00048 #include "G4AffineTransform.hh"
00049 #include "G4GeometryTolerance.hh"
00050 
00051 #include "G4VPVParameterisation.hh"
00052 
00053 #include "meshdefs.hh"
00054 
00055 #include "Randomize.hh"
00056 
00057 #include "G4VGraphicsScene.hh"
00058 #include "G4Polyhedron.hh"
00059 #include "G4NURBS.hh"
00060 #include "G4NURBSbox.hh"
00061 
00062 using namespace CLHEP;
00063  
00065 //
00066 // Private enum: Not for external use - used by distanceToOut
00067 
00068 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
00069 
00070 // used by normal
00071 
00072 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
00073 
00075 //
00076 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
00077 //               - note if pDPhi>2PI then reset to 2PI
00078 
00079 G4Cons::G4Cons( const G4String& pName,
00080                       G4double  pRmin1, G4double pRmax1,
00081                       G4double  pRmin2, G4double pRmax2,
00082                       G4double pDz,
00083                       G4double pSPhi, G4double pDPhi)
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 }
00119 
00121 //
00122 // Fake default constructor - sets only member data and allocates memory
00123 //                            for usage restricted to object persistency.
00124 //
00125 G4Cons::G4Cons( __void__& a )
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 }
00133 
00135 //
00136 // Destructor
00137 
00138 G4Cons::~G4Cons()
00139 {
00140 }
00141 
00143 //
00144 // Copy constructor
00145 
00146 G4Cons::G4Cons(const G4Cons& rhs)
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 }
00156 
00158 //
00159 // Assignment operator
00160 
00161 G4Cons& G4Cons::operator = (const G4Cons& rhs) 
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 }
00186 
00188 //
00189 // Return whether point inside/outside/on surface
00190 
00191 EInside G4Cons::Inside(const G4ThreeVector& p) const
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 }
00243 
00245 //
00246 // Dispatch to parameterisation for replication mechanism dimension
00247 // computation & modification.
00248 
00249 void G4Cons::ComputeDimensions(      G4VPVParameterisation* p,
00250                                const G4int                  n,
00251                                const G4VPhysicalVolume*     pRep    )
00252 {
00253   p->ComputeDimensions(*this,n,pRep) ;
00254 }
00255 
00256 
00258 //
00259 // Calculate extent under transform and specified limit
00260 
00261 G4bool G4Cons::CalculateExtent( const EAxis              pAxis,
00262               const G4VoxelLimits&     pVoxelLimit,
00263               const G4AffineTransform& pTransform,
00264                     G4double&          pMin,
00265                     G4double&          pMax  ) const
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 }
00472 
00474 //
00475 // Return unit normal of surface closest to p
00476 // - note if point on z axis, ignore phi divided sides
00477 // - unsafe if point close to z axis a rmin=0 - no explicit checks
00478 
00479 G4ThreeVector G4Cons::SurfaceNormal( const G4ThreeVector& p) const
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 }
00580 
00582 //
00583 // Algorithm for SurfaceNormal() following the original specification
00584 // for points not on the surface
00585 
00586 G4ThreeVector G4Cons::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00587 {
00588   ENorm side ;
00589   G4ThreeVector norm ;
00590   G4double rho, phi ;
00591   G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
00592   G4double tanRMin, secRMin, pRMin, widRMin ;
00593   G4double tanRMax, secRMax, pRMax, widRMax ;
00594 
00595   distZ = std::fabs(std::fabs(p.z()) - fDz) ;
00596   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
00597 
00598   tanRMin  = (fRmin2 - fRmin1)*0.5/fDz ;
00599   secRMin  = std::sqrt(1 + tanRMin*tanRMin) ;
00600   pRMin    = rho - p.z()*tanRMin ;
00601   widRMin  = fRmin2 - fDz*tanRMin ;
00602   distRMin = std::fabs(pRMin - widRMin)/secRMin ;
00603 
00604   tanRMax  = (fRmax2 - fRmax1)*0.5/fDz ;
00605   secRMax  = std::sqrt(1+tanRMax*tanRMax) ;
00606   pRMax    = rho - p.z()*tanRMax ;
00607   widRMax  = fRmax2 - fDz*tanRMax ;
00608   distRMax = std::fabs(pRMax - widRMax)/secRMax ;
00609   
00610   if (distRMin < distRMax)  // First minimum
00611   {
00612     if (distZ < distRMin)
00613     {
00614       distMin = distZ ;
00615       side    = kNZ ;
00616     }
00617     else
00618     {
00619       distMin = distRMin ;
00620       side    = kNRMin ;
00621     }
00622   }
00623   else
00624   {
00625     if (distZ < distRMax)
00626     {
00627       distMin = distZ ;
00628       side    = kNZ ;
00629     }
00630     else
00631     {
00632       distMin = distRMax ;
00633       side    = kNRMax ;
00634     }
00635   }
00636   if ( !fPhiFullCone && rho )  // Protected against (0,0,z) 
00637   {
00638     phi = std::atan2(p.y(),p.x()) ;
00639 
00640     if (phi < 0)  { phi += twopi; }
00641 
00642     if (fSPhi < 0)  { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
00643     else            { distSPhi = std::fabs(phi - fSPhi)*rho; }
00644 
00645     distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
00646 
00647     // Find new minimum
00648 
00649     if (distSPhi < distEPhi)
00650     {
00651       if (distSPhi < distMin)  { side = kNSPhi; }
00652     }
00653     else 
00654     {
00655       if (distEPhi < distMin)  { side = kNEPhi; }
00656     }
00657   }    
00658   switch (side)
00659   {
00660     case kNRMin:      // Inner radius
00661       rho *= secRMin ;
00662       norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
00663       break ;
00664     case kNRMax:      // Outer radius
00665       rho *= secRMax ;
00666       norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
00667       break ;
00668     case kNZ:         // +/- dz
00669       if (p.z() > 0)  { norm = G4ThreeVector(0,0,1);  }
00670       else            { norm = G4ThreeVector(0,0,-1); }
00671       break ;
00672     case kNSPhi:
00673       norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
00674       break ;
00675     case kNEPhi:
00676       norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
00677       break ;
00678     default:          // Should never reach this case...
00679       DumpInfo();
00680       G4Exception("G4Cons::ApproxSurfaceNormal()",
00681                   "GeomSolids1002", JustWarning,
00682                   "Undefined side for valid surface normal to solid.");
00683       break ;    
00684   }
00685   return norm ;
00686 }
00687 
00689 //
00690 // Calculate distance to shape from outside, along normalised vector
00691 // - return kInfinity if no intersection, or intersection distance <= tolerance
00692 //
00693 // - Compute the intersection with the z planes 
00694 //        - if at valid r, phi, return
00695 //
00696 // -> If point is outside cone, compute intersection with rmax1*0.5
00697 //        - if at valid phi,z return
00698 //        - if inside outer cone, handle case when on tolerant outer cone
00699 //          boundary and heading inwards(->0 to in)
00700 //
00701 // -> Compute intersection with inner cone, taking largest +ve root
00702 //        - if valid (in z,phi), save intersction
00703 //
00704 //    -> If phi segmented, compute intersections with phi half planes
00705 //        - return smallest of valid phi intersections and
00706 //          inner radius intersection
00707 //
00708 // NOTE:
00709 // - `if valid' implies tolerant checking of intersection points
00710 // - z, phi intersection from Tubs
00711 
00712 G4double G4Cons::DistanceToIn( const G4ThreeVector& p,
00713                                const G4ThreeVector& v   ) const
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 }
01376 
01378 // 
01379 // Calculate distance (<= actual) to closest surface of shape from outside
01380 // - Calculate distance to z, radial planes
01381 // - Only to phi planes if outside phi extent
01382 // - Return 0 if point inside
01383 
01384 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
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 }
01440 
01442 //
01443 // Calculate distance to surface of shape from 'inside', allowing for tolerance
01444 // - Only Calc rmax intersection if no valid rmin intersection
01445 
01446 G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
01447                                 const G4ThreeVector& v,
01448                                 const G4bool calcNorm,
01449                                       G4bool *validNorm,
01450                                       G4ThreeVector *n) const
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 }
02068 
02070 //
02071 // Calculate distance (<=actual) to closest surface of shape from inside
02072 
02073 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
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 }
02146 
02148 //
02149 // Create a List containing the transformed vertices
02150 // Ordering [0-3] -fDz cross section
02151 //          [4-7] +fDz cross section such that [0] is below [4],
02152 //                                             [1] below [5] etc.
02153 // Note:
02154 //  Caller has deletion resposibility
02155 //  Potential improvement: For last slice, use actual ending angle
02156 //                         to avoid rounding error problems.
02157 
02158 G4ThreeVectorList*
02159 G4Cons::CreateRotatedVertices(const G4AffineTransform& pTransform) const
02160 {
02161   G4ThreeVectorList* vertices ;
02162   G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
02163   G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
02164   G4double cosCrossAngle, sinCrossAngle, sAngle ;
02165   G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
02166   G4int crossSection, noCrossSections ;
02167 
02168   // Compute no of cross-sections necessary to mesh cone
02169     
02170   noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
02171 
02172   if (noCrossSections < kMinMeshSections)
02173   {
02174     noCrossSections = kMinMeshSections ;
02175   }
02176   else if (noCrossSections > kMaxMeshSections)
02177   {
02178     noCrossSections = kMaxMeshSections ;
02179   }
02180   meshAngle = fDPhi/(noCrossSections - 1) ;
02181 
02182   meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
02183   meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;
02184 
02185   // If complete in phi, set start angle such that mesh will be at RMax
02186   // on the x axis. Will give better extent calculations when not rotated.
02187 
02188   if ( fPhiFullCone && (fSPhi == 0.0) )
02189   {
02190     sAngle = -meshAngle*0.5 ;
02191   }
02192   else
02193   {
02194     sAngle = fSPhi ;
02195   } 
02196   vertices = new G4ThreeVectorList();
02197 
02198   if (vertices)
02199   {
02200     vertices->reserve(noCrossSections*4) ;
02201     for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
02202     {
02203       // Compute coordinates of cross section at section crossSection
02204 
02205       crossAngle    = sAngle + crossSection*meshAngle ;
02206       cosCrossAngle = std::cos(crossAngle) ;
02207       sinCrossAngle = std::sin(crossAngle) ;
02208 
02209       rMaxX1 = meshRMax1*cosCrossAngle ;
02210       rMaxY1 = meshRMax1*sinCrossAngle ;
02211       rMaxX2 = meshRMax2*cosCrossAngle ;
02212       rMaxY2 = meshRMax2*sinCrossAngle ;
02213         
02214       rMinX1 = fRmin1*cosCrossAngle ;
02215       rMinY1 = fRmin1*sinCrossAngle ;
02216       rMinX2 = fRmin2*cosCrossAngle ;
02217       rMinY2 = fRmin2*sinCrossAngle ;
02218         
02219       vertex0 = G4ThreeVector(rMinX1,rMinY1,-fDz) ;
02220       vertex1 = G4ThreeVector(rMaxX1,rMaxY1,-fDz) ;
02221       vertex2 = G4ThreeVector(rMaxX2,rMaxY2,+fDz) ;
02222       vertex3 = G4ThreeVector(rMinX2,rMinY2,+fDz) ;
02223 
02224       vertices->push_back(pTransform.TransformPoint(vertex0)) ;
02225       vertices->push_back(pTransform.TransformPoint(vertex1)) ;
02226       vertices->push_back(pTransform.TransformPoint(vertex2)) ;
02227       vertices->push_back(pTransform.TransformPoint(vertex3)) ;
02228     }
02229   }
02230   else
02231   {
02232     DumpInfo();
02233     G4Exception("G4Cons::CreateRotatedVertices()",
02234                 "GeomSolids0003", FatalException,
02235                 "Error in allocation of vertices. Out of memory !");
02236   }
02237 
02238   return vertices ;
02239 }
02240 
02242 //
02243 // GetEntityType
02244 
02245 G4GeometryType G4Cons::GetEntityType() const
02246 {
02247   return G4String("G4Cons");
02248 }
02249 
02251 //
02252 // Make a clone of the object
02253 //
02254 G4VSolid* G4Cons::Clone() const
02255 {
02256   return new G4Cons(*this);
02257 }
02258 
02260 //
02261 // Stream object contents to an output stream
02262 
02263 std::ostream& G4Cons::StreamInfo(std::ostream& os) const
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 }
02283 
02284 
02285 
02287 //
02288 // GetPointOnSurface
02289 
02290 G4ThreeVector G4Cons::GetPointOnSurface() const
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 }
02372 
02374 //
02375 // Methods for visualisation
02376 
02377 void G4Cons::DescribeYourselfTo (G4VGraphicsScene& scene) const
02378 {
02379   scene.AddSolid (*this);
02380 }
02381 
02382 G4Polyhedron* G4Cons::CreatePolyhedron () const
02383 {
02384   return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
02385 }
02386 
02387 G4NURBS* G4Cons::CreateNURBS () const
02388 {
02389   G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;
02390   return new G4NURBSbox (RMax, RMax, fDz);       // Box for now!!!
02391 }

Generated on Mon May 27 17:47:57 2013 for Geant4 by  doxygen 1.4.7