Changeset 1228 for trunk/source/geometry/solids/CSG/src

Timestamp:

Jan 8, 2010, 11:56:51 AM (14 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

Location:

trunk/source/geometry/solids/CSG/src

Files:

• G4Box.cc (modified) (1 diff)
• G4CSGSolid.cc (modified) (1 diff)
• G4Cons.cc (modified) (21 diffs)
• G4Orb.cc (modified) (10 diffs)
• G4Para.cc (modified) (1 diff)
• G4Sphere.cc (modified) (126 diffs)
• G4Torus.cc (modified) (2 diffs)
• G4Trap.cc (modified) (1 diff)
• G4Trd.cc (modified) (1 diff)
• G4Tubs.cc (modified) (18 diffs)

trunk/source/geometry/solids/CSG/src/G4Box.cc

@@ -26 +26 @@
 //
 // $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // 

trunk/source/geometry/solids/CSG/src/G4CSGSolid.cc

@@ -26 +26 @@
 //
 // $Id: G4CSGSolid.cc,v 1.13 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // --------------------------------------------------------------------

trunk/source/geometry/solids/CSG/src/G4Cons.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4Cons.cc,v 1.60 2008/11/06 15:26:53 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// $Id: G4Cons.cc,v 1.67 2009/11/12 11:53:11 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 //
@@ -35 +35 @@
 // History:
 //
+// 12.10.09 T.Nikitina: Added to DistanceToIn(p,v) check on the direction in
+//                      case of point on surface
 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
 // 13.09.96 V.Grichine: Review and final modifications
@@ -79 +81 @@
                       G4double pDz,
                       G4double pSPhi, G4double pDPhi)
-  : G4CSGSolid(pName)
+  : G4CSGSolid(pName), fSPhi(0), fDPhi(0)
 {
-  // Check z-len
-
   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 
+  // Check z-len
+  //
   if ( pDz > 0 )
   {
@@ -100 +102 @@
 
   // Check radii
-
+  //
   if ( (pRmin1<pRmax1) && (pRmin2<pRmax2) && (pRmin1>=0) && (pRmin2>=0) )
   {
@@ -121 +123 @@
   }
 
-  fPhiFullCone = true;
-  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
-  {
-    fDPhi=twopi;
-    fSPhi=0;
-  }
-  else
-  {
-    fPhiFullCone = false;
-    if ( pDPhi > 0 )
-    {
-      fDPhi = pDPhi;
-    }
-    else
-    {
-      G4cerr << "ERROR - G4Cons()::G4Cons(): " << GetName() << G4endl
-             << "        Negative delta-Phi ! - "
-             << pDPhi << G4endl;
-      G4Exception("G4Cons::G4Cons()", "InvalidSetup",
-                  FatalException, "Invalid dphi.");
-    }
-
-    // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
-
-    if ( pSPhi < 0 )
-    {
-      fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi);
-    }
-    else
-    {
-      fSPhi = std::fmod(pSPhi,twopi) ;
-    }
-    if ( fSPhi+fDPhi > twopi )
-    {
-      fSPhi -= twopi ;
-    }
-  }
-  InitializeTrigonometry();
+  // Check angles
+  //
+  CheckPhiAngles(pSPhi, pDPhi);
 }
 
@@ -705 +672 @@
 {
   G4double snxt = kInfinity ;      // snxt = default return value
-
+  const G4double dRmax = 100*std::min(fRmax1,fRmax2);
   static const G4double halfCarTolerance=kCarTolerance*0.5;
   static const G4double halfRadTolerance=kRadTolerance*0.5;
@@ -717 +684 @@
   G4double tolODz,tolIDz ;
 
-  G4double Dist,s,xi,yi,zi,ri=0.,rhoi2,cosPsi ; // Intersection point variables
+  G4double Dist,s,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars
 
   G4double t1,t2,t3,b,c,d ;    // Quadratic solver variables 
   G4double nt1,nt2,nt3 ;
   G4double Comp ;
+
+  G4ThreeVector Normal;
 
   // Cone Precalcs
@@ -887 +856 @@
         if ( s > 0 )  // If 'forwards'. Check z intersection
         {
+          if ( s>dRmax ) // Avoid rounding errors due to precision issues on
+          {              // 64 bits systems. Split long distances and recompute
+            G4double fTerm = s-std::fmod(s,dRmax);
+            s = fTerm + DistanceToIn(p+fTerm*v,v);
+          } 
           zi = p.z() + s*v.z() ;
 
@@ -917 +891 @@
       {
         // Inside cones, delta r -ve, inside z extent
-
+        // Point is on the Surface => check Direction using  Normal.dot(v)
+
+        xi     = p.x() ;
+        yi     = p.y()  ;
+        risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
+        Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
         if ( !fPhiFullCone )
         {
           cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
-
-          if (cosPsi >= cosHDPhiIT)  { return 0.0; }
-        }
-        else  { return 0.0; }
+          if ( cosPsi >= cosHDPhiIT )
+          {
+            if ( Normal.dot(v) <= 0 )  { return 0.0; }
+          }
+        }
+        else
+        {             
+          if ( Normal.dot(v) <= 0 )  { return 0.0; }
+        }
       }
     }
@@ -972 +956 @@
 
   if (rMinAv)
-  {
+  { 
     nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
     nt2 = t2 - tanRMin*v.z()*rin ;
@@ -993 +977 @@
           if ( s >= 0 )   // > 0
           {
+            if ( s>dRmax ) // Avoid rounding errors due to precision issues on
+            {              // 64 bits systems. Split long distance and recompute
+              G4double fTerm = s-std::fmod(s,dRmax);
+              s = fTerm + DistanceToIn(p+fTerm*v,v);
+            } 
             zi = p.z() + s*v.z() ;
 
@@ -1004 +993 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
+                if (cosPsi >= cosHDPhiIT)
+                { 
+                  if ( s > halfRadTolerance )  { snxt=s; }
+                  else
+                  {
+                    // Calculate a normal vector in order to check Direction
+
+                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
+                    if ( Normal.dot(v) <= 0 )  { snxt = s; }
+                  } 
+                }
               }
-              else  { return s; }
+              else
+              {
+                if ( s > halfRadTolerance )  { return s; }
+                else
+                {
+                  // Calculate a normal vector in order to check Direction
+
+                  xi     = p.x() + s*v.x() ;
+                  yi     = p.y() + s*v.y() ;
+                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
+                  if ( Normal.dot(v) <= 0 )  { return s; }
+                }
+              }
             }
           }
@@ -1032 +1045 @@
             if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s > 0
             {
+              if ( s>dRmax ) // Avoid rounding errors due to precision issues
+              {              // seen on 64 bits systems. Split and recompute
+                G4double fTerm = s-std::fmod(s,dRmax);
+                s = fTerm + DistanceToIn(p+fTerm*v,v);
+              } 
               if ( !fPhiFullCone )
               {
@@ -1038 +1056 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiOT)  { snxt = s; }
+                if (cosPsi >= cosHDPhiOT)
+                {
+                  if ( s > halfRadTolerance )  { snxt=s; }
+                  else
+                  {
+                    // Calculate a normal vector in order to check Direction
+
+                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
+                    if ( Normal.dot(v) <= 0 )  { snxt = s; } 
+                  }
+                }
               }
-              else  { return s; }
+              else
+              {
+                if( s > halfRadTolerance )  { return s; }
+                else
+                {
+                  // Calculate a normal vector in order to check Direction
+
+                  xi     = p.x() + s*v.x() ;
+                  yi     = p.y() + s*v.y() ;
+                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
+                  if ( Normal.dot(v) <= 0 )  { return s; }
+                } 
+              }
             }
           }
@@ -1051 +1093 @@
             if ( (s >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // s>0
             {
+              if ( s>dRmax ) // Avoid rounding errors due to precision issues
+              {              // seen on 64 bits systems. Split and recompute
+                G4double fTerm = s-std::fmod(s,dRmax);
+                s = fTerm + DistanceToIn(p+fTerm*v,v);
+              } 
               if ( !fPhiFullCone )
               {
@@ -1057 +1104 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
+                if (cosPsi >= cosHDPhiIT)
+                {
+                  if ( s > halfRadTolerance )  { snxt=s; }
+                  else
+                  {
+                    // Calculate a normal vector in order to check Direction
+
+                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
+                    if ( Normal.dot(v) <= 0 )  { snxt = s; } 
+                  }
+                }
               }
-              else  { return s; }
+              else
+              {
+                if ( s > halfRadTolerance )  { return s; }
+                else
+                {
+                  // Calculate a normal vector in order to check Direction
+
+                  xi     = p.x() + s*v.x() ;
+                  yi     = p.y() + s*v.y() ;
+                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
+                  if ( Normal.dot(v) <= 0 )  { return s; }
+                } 
+              }
             }
           }
@@ -1099 +1170 @@
               zi = p.z() + s*v.z() ;
               ri = rMinAv + zi*tanRMin ;
-
+              
               if ( ri > 0 )   // 2nd root
               {
@@ -1107 +1178 @@
                 if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
                 {
+                  if ( s>dRmax ) // Avoid rounding errors due to precision issue
+                  {              // seen on 64 bits systems. Split and recompute
+                    G4double fTerm = s-std::fmod(s,dRmax);
+                    s = fTerm + DistanceToIn(p+fTerm*v,v);
+                  } 
                   if ( !fPhiFullCone )
                   {
@@ -1136 +1212 @@
             if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
             {
+              if ( s>dRmax ) // Avoid rounding errors due to precision issues
+              {              // seen on 64 bits systems. Split and recompute
+                G4double fTerm = s-std::fmod(s,dRmax);
+                s = fTerm + DistanceToIn(p+fTerm*v,v);
+              } 
               if ( !fPhiFullCone )
               {
@@ -1335 +1416 @@
   // Vars for intersection within tolerance
 
-  ESide    sidetol ;
+  ESide    sidetol = kNull ;
   G4double slentol = kInfinity ;
 
@@ -1669 +1750 @@
             xi     = p.x() + slentol*v.x() ;
             yi     = p.y() + slentol*v.y() ;
-            risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
-            Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMin/secRMin) ;
-            
+            if( sidetol==kRMax )
+            {
+              risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
+              Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
+            }
+            else
+            {
+              risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
+              Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
+            }
             if( Normal.dot(v) > 0 )
             {

trunk/source/geometry/solids/CSG/src/G4Orb.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Orb.cc,v 1.24 2007/05/18 07:38:01 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// $Id: G4Orb.cc,v 1.30 2009/11/30 10:20:38 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // class G4Orb
@@ -298 +298 @@
   rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z() ;
 
+  G4double rad = std::sqrt(rad2);
+
   // G4double rad = std::sqrt(rad2);
   // Check radial surface
@@ -304 +306 @@
   tolRMax = fRmax - fRmaxTolerance*0.5 ;
     
-  if ( rad2 <= tolRMax*tolRMax )  in = kInside ;
+  if ( rad <= tolRMax )  { in = kInside ; }
   else
   {
     tolRMax = fRmax + fRmaxTolerance*0.5 ;       
-    if ( rad2 <= tolRMax*tolRMax ) in = kSurface ;
-    else                           in = kOutside ;
+    if ( rad <= tolRMax )  { in = kSurface ; }
+    else                   { in = kOutside ; }
   }
   return in;
@@ -360 +362 @@
   G4double c, d2, s = kInfinity ;
 
+  const G4double dRmax = 100.*fRmax;
+
   // General Precalcs
 
@@ -384 +388 @@
   // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
 
-  c = rad2 - fRmax*fRmax ;
+
+  G4double rad = std::sqrt(rad2);
+  c = (rad - fRmax)*(rad + fRmax);
 
   if ( c > fRmaxTolerance*fRmax )
@@ -396 +402 @@
     {
       s = -pDotV3d - std::sqrt(d2) ;
-
-      if (s >= 0 ) return snxt = s;
-           
+      if ( s >= 0 )
+      {
+        if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
+        {              // 64 bits systems. Split long distances and recompute
+          G4double fTerm = s-std::fmod(s,dRmax);
+          s = fTerm + DistanceToIn(p+fTerm*v,v);
+        } 
+        return snxt = s;
+      }
     }
     else    // No intersection with G4Orb
@@ -410 +422 @@
     {
       d2 = pDotV3d*pDotV3d - c ;             
-      //  if ( pDotV3d >= 0 ) return snxt = kInfinity;
-      if ( d2 < fRmaxTolerance*fRmax || pDotV3d >= 0 ) return snxt = kInfinity;
-      else                return snxt = 0.;
-    }
+      if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
+      {
+        return snxt = kInfinity;
+      }
+      else
+      {
+        return snxt = 0.;
+      }
+    }
+#ifdef G4CSGDEBUG
     else // inside ???
     {
@@ -419 +437 @@
                   JustWarning, "Point p is inside !?");
     }
+#endif
   }
   return snxt;
@@ -431 +450 @@
 G4double G4Orb::DistanceToIn( const G4ThreeVector& p ) const
 {
-  G4double safe=0.0, rad  = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
-                 safe = rad - fRmax;
-  if( safe < 0 ) safe = 0. ;
+  G4double safe = 0.0,
+           rad  = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
+  safe = rad - fRmax;
+  if( safe < 0 ) { safe = 0.; }
   return safe;
 }
@@ -475 +495 @@
   
   const G4double  Rmax_plus = fRmax + fRmaxTolerance*0.5;
-
-  if( rad2 <= Rmax_plus*Rmax_plus )
-  {
-    c = rad2-fRmax*fRmax ;
-
-    if ( c < fRmaxTolerance*fRmax) 
+  G4double rad = std::sqrt(rad2);
+
+  if ( rad <= Rmax_plus )
+  {
+    c = (rad - fRmax)*(rad + fRmax);
+
+    if ( c < fRmaxTolerance*fRmax ) 
     {
       // Within tolerant Outer radius 

trunk/source/geometry/solids/CSG/src/G4Para.cc

@@ -26 +26 @@
 //
 // $Id: G4Para.cc,v 1.39 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // class G4Para

trunk/source/geometry/solids/CSG/src/G4Sphere.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4Sphere.cc,v 1.68 2008/07/07 09:35:16 grichine Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// $Id: G4Sphere.cc,v 1.84 2009/08/07 15:56:23 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // class G4Sphere
@@ -34 +34 @@
 // History:
 //
+// 14.09.09 T.Nikitina: fix for phi section in DistanceToOut(p,v,..),as for G4Tubs,G4Cons 
+// 26.03.09 G.Cosmo   : optimisations and uniform use of local radial tolerance
 // 12.06.08 V.Grichine: fix for theta intersections in DistanceToOut(p,v,...)
 // 22.07.05 O.Link    : Added check for intersection with double cone
@@ -41 +43 @@
 // 02.06.04 V.Grichine: bug fixed in DistanceToIn(p,v), on Rmax,Rmin go inside
 // 30.10.03 J.Apostolakis: new algorithm in Inside for SPhi-sections
-// 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi < SPhi, SPhi<0
+// 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi<SPhi, SPhi<0
 // 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance
 // 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451
@@ -53 +55 @@
 // --------------------------------------------------------------------
 
-#include <assert.h>
-
 #include "G4Sphere.hh"
 
@@ -92 +92 @@
                           G4double pSPhi, G4double pDPhi,
                           G4double pSTheta, G4double pDTheta )
-  : G4CSGSolid(pName)
+  : G4CSGSolid(pName), fFullPhiSphere(true), fFullThetaSphere(true)
 {
-  fEpsilon = 1.0e-14;
-
-  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
+  fEpsilon = 2.0e-11;  // relative radial tolerance constant
+
   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 
-  // Check radii
-
-  if (pRmin<pRmax&&pRmin>=0)
+  // Check radii and set radial tolerances
+
+  G4double kRadTolerance = G4GeometryTolerance::GetInstance()
+                         ->GetRadialTolerance();
+  if ( (pRmin < pRmax) && (pRmax >= 10*kRadTolerance) && (pRmin >= 0) )
   {
     fRmin=pRmin; fRmax=pRmax;
+    fRminTolerance = (pRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
+    fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
   }
   else
@@ -116 +119 @@
   // Check angles
 
-  if (pDPhi>=twopi)
-  {
-    fDPhi=twopi;
-  }
-  else if (pDPhi>0)
-  {
-    fDPhi=pDPhi;
-  }
-  else
-  {
-    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
-           << "        Negative Z delta-Phi ! - "
-           << pDPhi << G4endl;
-    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
-                "Invalid DPhi.");
-  }
-
-  // Convert fSPhi to 0-2PI
-
-  if (pSPhi<0)
-  {
-    fSPhi=twopi-std::fmod(std::fabs(pSPhi),twopi);
-  }
-  else
-  {
-    fSPhi=std::fmod(pSPhi,twopi);
-  }
-
-  // Sphere is placed such that fSPhi+fDPhi>twopi !
-  // fSPhi could be < 0 !!?
-  //
-  if (fSPhi+fDPhi>twopi) fSPhi-=twopi;
-
-  // Check theta angles
-
-  if (pSTheta<0 || pSTheta>pi)
-  {
-    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl;
-    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
-                "stheta outside 0-PI range.");
-  }
-  else
-  {
-    fSTheta=pSTheta;
-  }
-
-  if (pDTheta+pSTheta>=pi)
-  {
-    fDTheta=pi-pSTheta;
-  }
-  else if (pDTheta>0)
-  {
-    fDTheta=pDTheta;
-  }
-  else
-  {
-    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
-           << "        Negative delta-Theta ! - "
-           << pDTheta << G4endl;
-    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
-                "Invalid pDTheta.");
-  }
+  CheckPhiAngles(pSPhi, pDPhi);
+  CheckThetaAngles(pSTheta, pDTheta);
 }
 
@@ -219 +162 @@
                                         G4double& pMin, G4double& pMax ) const
 {
-  if ( fDPhi==twopi && fDTheta==pi)  // !pTransform.IsRotated() &&
+  if ( fFullSphere )
   {
     // Special case handling for solid spheres-shells
@@ -310 +253 @@
         yoff1=yoffset-yMin;
         yoff2=yMax-yoffset;
-        if (yoff1>=0&&yoff2>=0)
+        if ((yoff1>=0) && (yoff2>=0))
         {
           // Y limits cross max/min x => no change
@@ -334 +277 @@
         xoff1=xoffset-xMin;
         xoff2=xMax-xoffset;
-        if (xoff1>=0&&xoff2>=0)
+        if ((xoff1>=0) && (xoff2>=0))
         {
           // X limits cross max/min y => no change
@@ -416 +359 @@
     }
       
-    if (pMin!=kInfinity || pMax!=-kInfinity)
+    if ((pMin!=kInfinity) || (pMax!=-kInfinity))
     {
       existsAfterClip=true;
@@ -453 +396 @@
 // Return whether point inside/outside/on surface
 // Split into radius, phi, theta checks
-// Each check modifies `in', or returns as approprate
+// Each check modifies 'in', or returns as approprate
 
 EInside G4Sphere::Inside( const G4ThreeVector& p ) const
@@ -459 +402 @@
   G4double rho,rho2,rad2,tolRMin,tolRMax;
   G4double pPhi,pTheta;
-  EInside in=kOutside;
+  EInside in = kOutside;
+  static const G4double halfAngTolerance = kAngTolerance*0.5;
+  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
+  const G4double halfRminTolerance = fRminTolerance*0.5;
+  const G4double Rmax_minus = fRmax - halfRmaxTolerance;
+  const G4double Rmin_plus  = (fRmin > 0) ? fRmin+halfRminTolerance : 0;
 
   rho2 = p.x()*p.x() + p.y()*p.y() ;
   rad2 = rho2 + p.z()*p.z() ;
 
-  //  if(rad2 >= 1.369e+19) DBG();
-  //  G4double rad = std::sqrt(rad2);
-  // Check radial surfaces
-  // sets `in'
-
-  if ( fRmin ) tolRMin = fRmin + kRadTolerance*0.5;
-  else         tolRMin = 0 ;
-  
-  tolRMax = fRmax - kRadTolerance*0.5 ;
-  //  const G4double  fractionTolerance = 1.0e-12;
-  const G4double  flexRadMaxTolerance = // kRadTolerance;
-    std::max(kRadTolerance, fEpsilon * fRmax);
-
-  const G4double  Rmax_minus = fRmax - flexRadMaxTolerance*0.5;
-  const G4double  flexRadMinTolerance = std::max(kRadTolerance, 
-                     fEpsilon * fRmin);
-  const G4double  Rmin_plus = (fRmin > 0) ? fRmin + flexRadMinTolerance*0.5 : 0 ;
+  // Check radial surfaces. Sets 'in'
+
+  tolRMin = Rmin_plus;
+  tolRMax = Rmax_minus;
+
+  if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
+  {
+    in = kInside;
+  }
+  else
+  {
+    tolRMax = fRmax + halfRmaxTolerance;                  // outside case
+    tolRMin = std::max(fRmin-halfRminTolerance, 0.);      // outside case
+    if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
+    {
+      in = kSurface;
+    }
+    else
+    {
+      return in = kOutside;
+    }
+  }
+
+  // Phi boundaries   : Do not check if it has no phi boundary!
+
+  if ( !fFullPhiSphere && rho2 )  // [fDPhi < twopi] and [p.x or p.y]
+  {
+    pPhi = std::atan2(p.y(),p.x()) ;
+
+    if      ( pPhi < fSPhi - halfAngTolerance  ) { pPhi += twopi; }
+    else if ( pPhi > ePhi + halfAngTolerance )   { pPhi -= twopi; }
     
-if(rad2 <= Rmax_minus*Rmax_minus && rad2 >= Rmin_plus*Rmin_plus) in = kInside ; 
-
-// if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kInside ;
-  // if ( rad <= tolRMax && rad >= tolRMin )  in = kInside ;
-  else
-  {
-    tolRMax = fRmax + kRadTolerance*0.5 ;
-    tolRMin = fRmin - kRadTolerance*0.5 ;
-
-    if ( tolRMin < 0.0 ) tolRMin = 0.0 ;
-    
-     if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kSurface ;
-    //  if ( rad <= tolRMax && rad >= tolRMin )  in = kSurface ;
-    else                                                return in = kOutside ;
-  }
-
-  // Phi boundaries   : Do not check if it has no phi boundary!
-  // (in != kOutside). It is new J.Apostolakis proposal of 30.10.03
-
-  if ( ( fDPhi < twopi - kAngTolerance ) &&
-       ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
-  {
-    pPhi = std::atan2(p.y(),p.x()) ;
-
-    if      ( pPhi < fSPhi - kAngTolerance*0.5  )         pPhi += twopi ; 
-    else if ( pPhi > fSPhi + fDPhi + kAngTolerance*0.5 )  pPhi -= twopi; 
-    
-    if ((pPhi < fSPhi - kAngTolerance*0.5) ||  
-        (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) )  return in = kOutside ;
+    if ( (pPhi < fSPhi - halfAngTolerance)
+      || (pPhi > ePhi + halfAngTolerance) )      { return in = kOutside; }
     
     else if (in == kInside)  // else it's kSurface anyway already
     {
-      if ( (pPhi < fSPhi + kAngTolerance*0.5) || 
-           (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) )      in = kSurface ;       
+      if ( (pPhi < fSPhi + halfAngTolerance)
+        || (pPhi > ePhi - halfAngTolerance) )    { in = kSurface; }     
     }
   }
 
   // Theta bondaries
-  // (in!=kOutside)
   
-  if ( (rho2 || p.z()) && fDTheta < pi - kAngTolerance*0.5 )
+  if ( (rho2 || p.z()) && (!fFullThetaSphere) )
   {
     rho    = std::sqrt(rho2);
@@ -529 +463 @@
     if ( in == kInside )
     {
-      if ( (pTheta < fSTheta + kAngTolerance*0.5)
-        || (pTheta > fSTheta + fDTheta - kAngTolerance*0.5) )
-      {
-        if ( (pTheta >= fSTheta - kAngTolerance*0.5)
-          && (pTheta <= fSTheta + fDTheta + kAngTolerance*0.5) )
-        {
-          in = kSurface ;
+      if ( (pTheta < fSTheta + halfAngTolerance)
+        || (pTheta > eTheta - halfAngTolerance) )
+      {
+        if ( (pTheta >= fSTheta - halfAngTolerance)
+          && (pTheta <= eTheta + halfAngTolerance) )
+        {
+          in = kSurface;
         }
         else
         {
-          in = kOutside ;
+          in = kOutside;
         }
       }
@@ -545 +479 @@
     else
     {
-      if ( (pTheta < fSTheta - kAngTolerance*0.5)
-        || (pTheta > fSTheta + fDTheta + kAngTolerance*0.5) )
-      {
-        in = kOutside ;
+      if ( (pTheta < fSTheta - halfAngTolerance)
+        || (pTheta > eTheta + halfAngTolerance) )
+      {
+        in = kOutside;
       }
     }
@@ -568 +502 @@
   G4double distSPhi = kInfinity, distEPhi = kInfinity;
   G4double distSTheta = kInfinity, distETheta = kInfinity;
-  G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
   G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
   G4ThreeVector norm, sumnorm(0.,0.,0.);
+
+  static const G4double halfCarTolerance = 0.5*kCarTolerance;
+  static const G4double halfAngTolerance = 0.5*kAngTolerance;
 
   rho2 = p.x()*p.x()+p.y()*p.y();
@@ -579 +515 @@
   if (fRmin)  distRMin = std::fabs(rad-fRmin);
     
-  if ( rho && (fDPhi < twopi || fDTheta < pi) )
+  if ( rho && !fFullSphere )
   {
     pPhi = std::atan2(p.y(),p.x());
 
-    if(pPhi  < fSPhi-dAngle)           pPhi     += twopi;
-    else if(pPhi > fSPhi+fDPhi+dAngle) pPhi     -= twopi;
-  }
-  if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
+    if (pPhi < fSPhi-halfAngTolerance)     { pPhi += twopi; }
+    else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
+  }
+  if ( !fFullPhiSphere )
   {
     if ( rho )
     {
-      distSPhi = std::fabs( pPhi - fSPhi ); 
-      distEPhi = std::fabs(pPhi-fSPhi-fDPhi); 
+      distSPhi = std::fabs( pPhi-fSPhi ); 
+      distEPhi = std::fabs( pPhi-ePhi ); 
     }
     else if( !fRmin )
@@ -598 +534 @@
       distEPhi = 0.; 
     }
-    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
-    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
+    nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
+    nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
   }        
-  if ( fDTheta < pi ) // && rad ) // old limitation against (0,0,0)
+  if ( !fFullThetaSphere )
   {
     if ( rho )
@@ -607 +543 @@
       pTheta     = std::atan2(rho,p.z());
       distSTheta = std::fabs(pTheta-fSTheta); 
-      distETheta = std::fabs(pTheta-fSTheta-fDTheta);
+      distETheta = std::fabs(pTheta-eTheta);
  
-      nTs = G4ThreeVector(-std::cos(fSTheta)*p.x()/rho, //  *std::cos(pPhi),
-                          -std::cos(fSTheta)*p.y()/rho, //  *std::sin(pPhi),
-                           std::sin(fSTheta)                   );
-
-      nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*p.x()/rho, // *std::cos(pPhi),
-                           std::cos(fSTheta+fDTheta)*p.y()/rho, // *std::sin(pPhi),
-                          -std::sin(fSTheta+fDTheta)                  );    
+      nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
+                          -cosSTheta*p.y()/rho,
+                           sinSTheta          );
+
+      nTe = G4ThreeVector( cosETheta*p.x()/rho,
+                           cosETheta*p.y()/rho,
+                          -sinETheta          );    
     }
     else if( !fRmin )
@@ -622 +558 @@
       {              
         distSTheta = 0.;
-        nTs = G4ThreeVector(0.,0.,-1.);
-      }
-      if ( fSTheta + fDTheta < pi ) // distETheta = 0.;
+        nTs = G4ThreeVector(0.,0.,-1.);
+      }
+      if ( eTheta < pi )
       {              
         distETheta = 0.;
-        nTe = G4ThreeVector(0.,0.,1.);
+        nTe = G4ThreeVector(0.,0.,1.);
       }
     }    
   }
-  if( rad )  nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);
-
-  if( distRMax <= delta )
+  if( rad )  { nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad); }
+
+  if( distRMax <= halfCarTolerance )
   {
     noSurfaces ++;
     sumnorm += nR;
   }
-  if( fRmin && distRMin <= delta )
+  if( fRmin && (distRMin <= halfCarTolerance) )
   {
     noSurfaces ++;
     sumnorm -= nR;
   }
-  if( fDPhi < twopi )   
-  {
-    if (distSPhi <= dAngle)
+  if( !fFullPhiSphere )   
+  {
+    if (distSPhi <= halfAngTolerance)
     {
       noSurfaces ++;
       sumnorm += nPs;
     }
-    if (distEPhi <= dAngle) 
+    if (distEPhi <= halfAngTolerance) 
     {
       noSurfaces ++;
@@ -656 +592 @@
     }
   }
-  if ( fDTheta < pi )
-  {
-    if (distSTheta <= dAngle && fSTheta > 0.)
+  if ( !fFullThetaSphere )
+  {
+    if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
     {
       noSurfaces ++;
-      if( rad <= delta && fDPhi >= twopi) sumnorm += nZ;
-      else                                sumnorm += nTs;
-    }
-    if (distETheta <= dAngle && fSTheta+fDTheta < pi) 
+      if ((rad <= halfCarTolerance) && fFullPhiSphere)  { sumnorm += nZ;  }
+      else                                              { sumnorm += nTs; }
+    }
+    if ((distETheta <= halfAngTolerance) && (eTheta < pi)) 
     {
       noSurfaces ++;
-      if( rad <= delta && fDPhi >= twopi) sumnorm -= nZ;
-      else                                sumnorm += nTe;
-      if(sumnorm.z() == 0.)               sumnorm += nZ;
+      if ((rad <= halfCarTolerance) && fFullPhiSphere)  { sumnorm -= nZ;  }
+      else                                              { sumnorm += nTe; }
+      if(sumnorm.z() == 0.)  { sumnorm += nZ; }
     }
   }
@@ -680 +616 @@
      norm = ApproxSurfaceNormal(p);
   }
-  else if ( noSurfaces == 1 ) norm = sumnorm;
-  else                        norm = sumnorm.unit();
+  else if ( noSurfaces == 1 )  { norm = sumnorm; }
+  else                         { norm = sumnorm.unit(); }
   return norm;
 }
@@ -735 +671 @@
     
   pPhi = std::atan2(p.y(),p.x());
-  if (pPhi<0) pPhi += twopi;
-
-  if (fDPhi<twopi&&rho)
+  if (pPhi<0) { pPhi += twopi; }
+
+  if (!fFullPhiSphere && rho)
   {
     if (fSPhi<0)
@@ -774 +710 @@
   //
 
-  if (fDTheta<pi&&rad)
+  if (!fFullThetaSphere && rad)
   {
     pTheta=std::atan2(rho,p.z());
@@ -809 +745 @@
       break;
     case kNSPhi:
-      norm=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
+      norm=G4ThreeVector(sinSPhi,-cosSPhi,0);
       break;
     case kNEPhi:
-      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
+      norm=G4ThreeVector(-sinEPhi,cosEPhi,0);
       break;
     case kNSTheta:
-      norm=G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi),
-                         -std::cos(fSTheta)*std::sin(pPhi),
-                          std::sin(fSTheta)            );
-      //  G4cout<<G4endl<<" case kNSTheta:"<<G4endl;
-      //  G4cout<<"pPhi = "<<pPhi<<G4endl;
-      //  G4cout<<"rad  = "<<rad<<G4endl;
-      //  G4cout<<"pho  = "<<rho<<G4endl;
-      //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
-      //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
+      norm=G4ThreeVector(-cosSTheta*std::cos(pPhi),
+                         -cosSTheta*std::sin(pPhi),
+                          sinSTheta            );
       break;
     case kNETheta:
-      norm=G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi),
-                          std::cos(fSTheta+fDTheta)*std::sin(pPhi),
-                         -std::sin(fSTheta+fDTheta)              );
-
-      //  G4cout<<G4endl<<" case kNETheta:"<<G4endl;
-      //  G4cout<<"pPhi = "<<pPhi<<G4endl;
-      //  G4cout<<"rad  = "<<rad<<G4endl;
-      //  G4cout<<"pho  = "<<rho<<G4endl;
-      //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
-      //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
+      norm=G4ThreeVector( cosETheta*std::cos(pPhi),
+                          cosETheta*std::sin(pPhi),
+                         -sinETheta              );
       break;
     default:
       DumpInfo();
-      G4Exception("G4Sphere::ApproxSurfaceNormal()", "Notification", JustWarning,
+      G4Exception("G4Sphere::ApproxSurfaceNormal()","Notification",JustWarning,
                   "Undefined side for valid surface normal to solid.");
       break;    
-  } // end case
+  }
 
   return norm;
@@ -880 +803 @@
 {
   G4double snxt = kInfinity ;      // snxt = default return value
-
   G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
-
-  G4double tolIRMin2, tolORMin2, tolORMax2, tolIRMax2 ;
   G4double tolSTheta=0., tolETheta=0. ;
+  const G4double dRmax = 100.*fRmax;
+
+  static const G4double halfCarTolerance = kCarTolerance*0.5;
+  static const G4double halfAngTolerance = kAngTolerance*0.5;
+  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
+  const G4double halfRminTolerance = fRminTolerance*0.5;
+  const G4double tolORMin2 = (fRmin>halfRminTolerance)
+               ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
+  const G4double tolIRMin2 =
+               (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
+  const G4double tolORMax2 =
+               (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
+  const G4double tolIRMax2 =
+               (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
 
   // Intersection point
-
+  //
   G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
 
   // Phi intersection
-
-  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi , Comp ; 
-
-  // Phi flag and precalcs
-
-  G4bool segPhi ;        
-  G4double hDPhi, hDPhiOT, hDPhiIT, cPhi, sinCPhi=0., cosCPhi=0. ; 
-  G4double cosHDPhiOT=0., cosHDPhiIT=0. ;
+  //
+  G4double Comp ; 
+
+  // Phi precalcs
+  //
   G4double Dist, cosPsi ;
 
-  G4bool segTheta ;                             // Theta flag and precals
-  G4double tanSTheta, tanETheta ;
-  G4double tanSTheta2, tanETheta2 ;
+  // Theta precalcs
+  //
   G4double dist2STheta, dist2ETheta ;
   G4double t1, t2, b, c, d2, d, s = kInfinity ;
 
   // General Precalcs
-
+  //
   rho2 = p.x()*p.x() + p.y()*p.y() ;
   rad2 = rho2 + p.z()*p.z() ;
@@ -916 +846 @@
   pDotV3d = pDotV2d + p.z()*v.z() ;
 
-  // Radial Precalcs
-
-  if (fRmin > kRadTolerance*0.5)
-  {
-    tolORMin2=(fRmin-kRadTolerance*0.5)*(fRmin-kRadTolerance*0.5);
-  }
-  else
-  {
-    tolORMin2 = 0 ;
-  }
-  tolIRMin2 = (fRmin+kRadTolerance*0.5)*(fRmin+kRadTolerance*0.5) ;
-  tolORMax2 = (fRmax+kRadTolerance*0.5)*(fRmax+kRadTolerance*0.5) ;
-  tolIRMax2 = (fRmax-kRadTolerance*0.5)*(fRmax-kRadTolerance*0.5) ;
-
-  // Set phi divided flag and precalcs
-
-  if (fDPhi < twopi)
-  {
-    segPhi = true ;
-    hDPhi = 0.5*fDPhi ;    // half delta phi
-    cPhi = fSPhi + hDPhi ;
-
-    hDPhiOT = hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi 
-    hDPhiIT = hDPhi-0.5*kAngTolerance;
-
-    sinCPhi    = std::sin(cPhi) ;
-    cosCPhi    = std::cos(cPhi) ;
-    cosHDPhiOT = std::cos(hDPhiOT) ;
-    cosHDPhiIT = std::cos(hDPhiIT) ;
-  }
-  else
-  {
-    segPhi = false ;
-  }
-
   // Theta precalcs
-    
-  if (fDTheta < pi )
-  {
-    segTheta  = true ;
-    tolSTheta = fSTheta - kAngTolerance*0.5 ;
-    tolETheta = fSTheta + fDTheta + kAngTolerance*0.5 ;
-  }
-  else
-  {
-    segTheta = false ;
+  //
+  if (!fFullThetaSphere)
+  {
+    tolSTheta = fSTheta - halfAngTolerance ;
+    tolETheta = eTheta + halfAngTolerance ;
   }
 
@@ -979 +869 @@
 
   c = rad2 - fRmax*fRmax ;
-  const G4double  flexRadMaxTolerance = // kRadTolerance;
-    std::max(kRadTolerance, fEpsilon * fRmax);
-
-  //  if (c > kRadTolerance*fRmax)
-  if (c > flexRadMaxTolerance*fRmax)
-  {
-    // If outside toleranct boundary of outer G4Sphere
-    // [should be std::sqrt(rad2)-fRmax > kRadTolerance*0.5]
+
+  if (c > fRmaxTolerance*fRmax)
+  {
+    // If outside tolerant boundary of outer G4Sphere
+    // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]
 
     d2 = pDotV3d*pDotV3d - c ;
@@ -996 +883 @@
       if (s >= 0 )
       {
+        if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
+        {              // 64 bits systems. Split long distances and recompute
+          G4double fTerm = s-std::fmod(s,dRmax);
+          s = fTerm + DistanceToIn(p+fTerm*v,v);
+        } 
         xi   = p.x() + s*v.x() ;
         yi   = p.y() + s*v.y() ;
         rhoi = std::sqrt(xi*xi + yi*yi) ;
 
-        if (segPhi && rhoi)    // Check phi intersection
+        if (!fFullPhiSphere && rhoi)    // Check phi intersection
         {
           cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
@@ -1006 +898 @@
           if (cosPsi >= cosHDPhiOT)
           {
-            if (segTheta)   // Check theta intersection
+            if (!fFullThetaSphere)   // Check theta intersection
             {
               zi = p.z() + s*v.z() ;
@@ -1027 +919 @@
         else
         {
-          if (segTheta)    // Check theta intersection
+          if (!fFullThetaSphere)    // Check theta intersection
           {
             zi = p.z() + s*v.z() ;
@@ -1059 +951 @@
     d2 = pDotV3d*pDotV3d - c ;
 
-    // if (rad2 > tolIRMin2 && pDotV3d < 0 )
-
-    if (rad2 > tolIRMax2 && ( d2 >= flexRadMaxTolerance*fRmax && pDotV3d < 0 ) )
-    {
-      if (segPhi)
+    if ( (rad2 > tolIRMax2)
+      && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
+    {
+      if (!fFullPhiSphere)
       {
         // Use inner phi tolerant boundary -> if on tolerant
@@ -1074 +965 @@
           // inside radii, delta r -ve, inside phi
 
-          if (segTheta)
+          if ( !fFullThetaSphere )
           {
             if ( (pTheta >= tolSTheta + kAngTolerance)
@@ -1090 +981 @@
       else
       {
-        if ( segTheta )
+        if ( !fFullThetaSphere )
         {
           if ( (pTheta >= tolSTheta + kAngTolerance)
@@ -1109 +1000 @@
   // - Always farthest root, because would have passed through outer
   //   surface first.
-  // - Tolerant check for if travelling through solid
+  // - Tolerant check if travelling through solid
 
   if (fRmin)
@@ -1119 +1010 @@
     // Check for immediate entry/already inside and travelling outwards
 
-    // if (c >- kRadTolerance*0.5 && pDotV3d >= 0 && rad2 < tolIRMin2 )
-
-    if ( c > -kRadTolerance*0.5 && rad2 < tolIRMin2 && 
-         ( d2 < fRmin*kCarTolerance || pDotV3d >= 0 ) )
-    {
-      if (segPhi)
+    if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
+      && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
+    {
+      if ( !fFullPhiSphere )
       {
         // Use inner phi tolerant boundary -> if on tolerant
@@ -1134 +1023 @@
           // inside radii, delta r -ve, inside phi
           //
-          if (segTheta)
+          if ( !fFullThetaSphere )
           {
             if ( (pTheta >= tolSTheta + kAngTolerance)
@@ -1150 +1039 @@
       else
       {
-        if (segTheta)
+        if ( !fFullThetaSphere )
         {
           if ( (pTheta >= tolSTheta + kAngTolerance)
@@ -1166 +1055 @@
     else   // Not special tolerant case
     {
-      //  d2 = pDotV3d*pDotV3d - c ;
-
       if (d2 >= 0)
       {
         s = -pDotV3d + std::sqrt(d2) ;
-        if ( s >= kRadTolerance*0.5 )  // It was >= 0 ??
+        if ( s >= halfRminTolerance )  // It was >= 0 ??
         {
           xi   = p.x() + s*v.x() ;
@@ -1177 +1064 @@
           rhoi = std::sqrt(xi*xi+yi*yi) ;
 
-          if ( segPhi && rhoi )   // Check phi intersection
+          if ( !fFullPhiSphere && rhoi )   // Check phi intersection
           {
             cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
@@ -1183 +1070 @@
             if (cosPsi >= cosHDPhiOT)
             {
-              if (segTheta)  // Check theta intersection
+              if ( !fFullThetaSphere )  // Check theta intersection
               {
                 zi = p.z() + s*v.z() ;
@@ -1204 +1091 @@
           else
           {
-            if (segTheta)   // Check theta intersection
+            if ( !fFullThetaSphere )   // Check theta intersection
             {
               zi = p.z() + s*v.z() ;
@@ -1214 +1101 @@
               if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
               {
-                snxt = s ;
+                snxt = s;
               }
             }
             else
             {
-              snxt=s;
+              snxt = s;
             }
           }
@@ -1236 +1123 @@
   //         -> Should use some form of loop Construct
   //
-  if ( segPhi )
-  {
-    // First phi surface (`S'tarting phi)
-
-    sinSPhi = std::sin(fSPhi) ;
-    cosSPhi = std::cos(fSPhi) ;
-
+  if ( !fFullPhiSphere )
+  {
+    // First phi surface ('S'tarting phi)
     // Comp = Component in outwards normal dirn
     //
-    Comp    = v.x()*sinSPhi - v.y()*cosSPhi  ;
+    Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
                     
     if ( Comp < 0 )
@@ -1251 +1134 @@
       Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
 
-      if (Dist < kCarTolerance*0.5)
+      if (Dist < halfCarTolerance)
       {
         s = Dist/Comp ;
@@ -1282 +1165 @@
             // (=>intersect at origin =>fRmax=0)
             //
-            if ( segTheta )
+            if ( !fFullThetaSphere )
             {
               iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
@@ -1305 +1188 @@
     }
 
-    // Second phi surface (`E'nding phi)
-
-    ePhi    = fSPhi + fDPhi ;
-    sinEPhi = std::sin(ePhi)     ;
-    cosEPhi = std::cos(ePhi)     ;
-
-    // Compnent in outwards normal dirn
-
-    Comp    = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
+    // Second phi surface ('E'nding phi)
+    // Component in outwards normal dirn
+
+    Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
         
     if (Comp < 0)
     {
       Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
-      if ( Dist < kCarTolerance*0.5 )
+      if ( Dist < halfCarTolerance )
       {
         s = Dist/Comp ;
@@ -1340 +1218 @@
             rhoi2 = rho2  ;
             radi2 = rad2  ;
-          } if ( (radi2 <= tolORMax2)
+          }
+          if ( (radi2 <= tolORMax2)
             && (radi2 >= tolORMin2)
             && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
@@ -1348 +1227 @@
             // (=>intersect at origin =>fRmax=0)
             //
-            if ( segTheta )
+            if ( !fFullThetaSphere )
             {
               iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
@@ -1374 +1253 @@
   // Theta segment intersection
 
-  if ( segTheta )
+  if ( !fFullThetaSphere )
   {
 
@@ -1396 +1275 @@
     // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
 
-    tanSTheta  = std::tan(fSTheta)         ;
-    tanSTheta2 = tanSTheta*tanSTheta  ;
-    tanETheta  = std::tan(fSTheta+fDTheta) ;
-    tanETheta2 = tanETheta*tanETheta  ;
-      
     if (fSTheta)
     {
@@ -1409 +1283 @@
       dist2STheta = kInfinity ;
     }
-    if ( fSTheta + fDTheta < pi )
+    if ( eTheta < pi )
     {
       dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
     }
-      else
+    else
     {
       dist2ETheta=kInfinity;
     }      
-    if ( pTheta < tolSTheta) // dist2STheta<-kRadTolerance*0.5 && dist2ETheta>0)
+    if ( pTheta < tolSTheta )
     {
       // Inside (theta<stheta-tol) s theta cone
@@ -1424 +1298 @@
       t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
       t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
-        
-      b  = t2/t1 ;
-      c  = dist2STheta/t1 ;
-      d2 = b*b - c ;
-
-      if ( d2 >= 0 )
-      {
-        d = std::sqrt(d2) ;
-        s = -b - d ;    // First root
-        zi    = p.z() + s*v.z();
-
-        if ( s < 0 || zi*(fSTheta - halfpi) > 0 )
-        {
-          s = -b+d;    // Second root
-        }
-        if (s >= 0 && s < snxt)
-        {
-          xi    = p.x() + s*v.x();
-          yi    = p.y() + s*v.y();
+      if (t1)
+      {   
+        b  = t2/t1 ;
+        c  = dist2STheta/t1 ;
+        d2 = b*b - c ;
+
+        if ( d2 >= 0 )
+        {
+          d = std::sqrt(d2) ;
+          s = -b - d ;    // First root
           zi    = p.z() + s*v.z();
-          rhoi2 = xi*xi + yi*yi;
-          radi2 = rhoi2 + zi*zi;
-          if ( (radi2 <= tolORMax2)
-            && (radi2 >= tolORMin2)
-            && (zi*(fSTheta - halfpi) <= 0) )
-          {
-            if ( segPhi && rhoi2 )  // Check phi intersection
-            {
-              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-              if (cosPsi >= cosHDPhiOT)
+
+          if ( (s < 0) || (zi*(fSTheta - halfpi) > 0) )
+          {
+            s = -b+d;    // Second root
+          }
+          if ((s >= 0) && (s < snxt))
+          {
+            xi    = p.x() + s*v.x();
+            yi    = p.y() + s*v.y();
+            zi    = p.z() + s*v.z();
+            rhoi2 = xi*xi + yi*yi;
+            radi2 = rhoi2 + zi*zi;
+            if ( (radi2 <= tolORMax2)
+              && (radi2 >= tolORMin2)
+              && (zi*(fSTheta - halfpi) <= 0) )
+            {
+              if ( !fFullPhiSphere && rhoi2 )  // Check phi intersection
+              {
+                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+                if (cosPsi >= cosHDPhiOT)
+                {
+                  snxt = s ;
+                }
+              }
+              else
               {
                 snxt = s ;
               }
-            }
-            else
-            {
-              snxt = s ;
             }
           }
@@ -1469 +1345 @@
       // Second >= 0 root should be considered
         
-      if ( fSTheta + fDTheta < pi )
+      if ( eTheta < pi )
       {
         t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
         t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
-        
+        if (t1)
+        { 
+          b  = t2/t1 ;
+          c  = dist2ETheta/t1 ;
+          d2 = b*b - c ;
+
+          if (d2 >= 0)
+          {
+            d = std::sqrt(d2) ;
+            s = -b + d ;    // Second root
+
+            if ( (s >= 0) && (s < snxt) )
+            {
+              xi    = p.x() + s*v.x() ;
+              yi    = p.y() + s*v.y() ;
+              zi    = p.z() + s*v.z() ;
+              rhoi2 = xi*xi + yi*yi   ;
+              radi2 = rhoi2 + zi*zi   ;
+
+              if ( (radi2 <= tolORMax2)
+                && (radi2 >= tolORMin2)
+                && (zi*(eTheta - halfpi) <= 0) )
+              {
+                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
+                {
+                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+                  if (cosPsi >= cosHDPhiOT)
+                  {
+                    snxt = s ;
+                  }
+                }
+                else
+                {
+                  snxt = s ;
+                }
+              }
+            }
+          }
+        }
+      }
+    }  
+    else if ( pTheta > tolETheta ) 
+    { 
+      // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
+      // Inside (theta > etheta+tol) e-theta cone
+      // First root of etheta cone, second if first root 'imaginary'
+
+      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
+      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
+      if (t1)
+      {  
         b  = t2/t1 ;
         c  = dist2ETheta/t1 ;
@@ -1481 +1407 @@
         {
           d = std::sqrt(d2) ;
-          s = -b + d ;    // Second root
-
-          if (s >= 0 && s < snxt)
+          s = -b - d ;    // First root
+          zi    = p.z() + s*v.z();
+
+          if ( (s < 0) || (zi*(eTheta - halfpi) > 0) )
+          {
+            s = -b + d ;           // second root
+          }
+          if ( (s >= 0) && (s < snxt) )
           {
             xi    = p.x() + s*v.x() ;
@@ -1492 +1423 @@
 
             if ( (radi2 <= tolORMax2)
-              && (radi2 >= tolORMin2)
-              && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
-            {
-              if (segPhi && rhoi2)   // Check phi intersection
+              && (radi2 >= tolORMin2) 
+              && (zi*(eTheta - halfpi) <= 0) )
+            {
+              if (!fFullPhiSphere && rhoi2)  // Check phi intersection
               {
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
@@ -1511 +1442 @@
         }
       }
-    }  
-    else if ( pTheta > tolETheta ) 
-    { 
-      // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
-      // Inside (theta > etheta+tol) e-theta cone
-      // First root of etheta cone, second if first root `imaginary'
-
-      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
-      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
-        
-      b  = t2/t1 ;
-      c  = dist2ETheta/t1 ;
-      d2 = b*b - c ;
-
-      if (d2 >= 0)
-      {
-        d = std::sqrt(d2) ;
-        s = -b - d ;    // First root
-        zi    = p.z() + s*v.z();
-
-        if (s < 0 || zi*(fSTheta + fDTheta - halfpi) > 0)
-        {
-          s = -b + d ;           // second root
-        }
-        if (s >= 0 && s < snxt)
-        {
-          xi    = p.x() + s*v.x() ;
-          yi    = p.y() + s*v.y() ;
-          zi    = p.z() + s*v.z() ;
-          rhoi2 = xi*xi + yi*yi   ;
-          radi2 = rhoi2 + zi*zi   ;
-
-          if ( (radi2 <= tolORMax2)
-            && (radi2 >= tolORMin2) 
-            && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
-          {
-            if (segPhi && rhoi2)  // Check phi intersection
-            {
-              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-              if (cosPsi >= cosHDPhiOT)
-              {
-                snxt = s ;
-              }
-            }
-            else
-            {
-              snxt = s ;
-            }
-          }
-        }
-      }
 
       // Possible intersection with STheta cone. 
@@ -1570 +1450 @@
         t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
         t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
-
+        if (t1)
+        { 
+          b  = t2/t1 ;
+          c  = dist2STheta/t1 ;
+          d2 = b*b - c ;
+
+          if (d2 >= 0)
+          {
+            d = std::sqrt(d2) ;
+            s = -b + d ;    // Second root
+
+            if ( (s >= 0) && (s < snxt) )
+            {
+              xi    = p.x() + s*v.x() ;
+              yi    = p.y() + s*v.y() ;
+              zi    = p.z() + s*v.z() ;
+              rhoi2 = xi*xi + yi*yi   ;
+              radi2 = rhoi2 + zi*zi   ;
+
+              if ( (radi2 <= tolORMax2)
+                && (radi2 >= tolORMin2)
+                && (zi*(fSTheta - halfpi) <= 0) )
+              {
+                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
+                {
+                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+                  if (cosPsi >= cosHDPhiOT)
+                  {
+                    snxt = s ;
+                  }
+                }
+                else
+                {
+                  snxt = s ;
+                }
+              }
+            }
+          }
+        }
+      }  
+    }     
+    else if ( (pTheta < tolSTheta + kAngTolerance)
+           && (fSTheta > halfAngTolerance) )
+    {
+      // In tolerance of stheta
+      // If entering through solid [r,phi] => 0 to in
+      // else try 2nd root
+
+      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
+      if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
+        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
+        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
+      {
+        if (!fFullPhiSphere && rho2)  // Check phi intersection
+        {
+          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
+          if (cosPsi >= cosHDPhiIT)
+          {
+            return 0 ;
+          }
+        }
+        else
+        {
+          return 0 ;
+        }
+      }
+
+      // Not entering immediately/travelling through
+
+      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
+      if (t1)
+      { 
         b  = t2/t1 ;
         c  = dist2STheta/t1 ;
@@ -1578 +1529 @@
         {
           d = std::sqrt(d2) ;
-          s = -b + d ;    // Second root
-
-          if ( (s >= 0) && (s < snxt) )
-          {
+          s = -b + d ;
+          if ( (s >= halfCarTolerance) && (s < snxt) && (fSTheta < halfpi) )
+          {  // ^^^^^^^^^^^^^^^^^^^^^  shouldn't it be >=0 instead ?
             xi    = p.x() + s*v.x() ;
             yi    = p.y() + s*v.y() ;
@@ -1592 +1542 @@
               && (zi*(fSTheta - halfpi) <= 0) )
             {
-              if (segPhi && rhoi2)   // Check phi intersection
+              if ( !fFullPhiSphere && rhoi2 )    // Check phi intersection
+              {
+                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+                if ( cosPsi >= cosHDPhiOT )
+                {
+                  snxt = s ;
+                }
+              }
+              else
+              {
+                snxt = s ;
+              }
+            }
+          }
+        }
+      }
+    }   
+    else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
+    {
+
+      // In tolerance of etheta
+      // If entering through solid [r,phi] => 0 to in
+      // else try 2nd root
+
+      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
+
+      if (   ((t2<0) && (eTheta < halfpi)
+          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
+        ||   ((t2>=0) && (eTheta > halfpi)
+          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
+        ||   ((v.z()>0) && (eTheta == halfpi)
+          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))  )
+      {
+        if (!fFullPhiSphere && rho2)   // Check phi intersection
+        {
+          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
+          if (cosPsi >= cosHDPhiIT)
+          {
+            return 0 ;
+          }
+        }
+        else
+        {
+          return 0 ;
+        }
+      }
+
+      // Not entering immediately/travelling through
+
+      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
+      if (t1)
+      { 
+        b  = t2/t1 ;
+        c  = dist2ETheta/t1 ;
+        d2 = b*b - c ;
+
+        if (d2 >= 0)
+        {
+          d = std::sqrt(d2) ;
+          s = -b + d ;
+        
+          if ( (s >= halfCarTolerance)
+            && (s < snxt) && (eTheta > halfpi) )
+          {
+            xi    = p.x() + s*v.x() ;
+            yi    = p.y() + s*v.y() ;
+            zi    = p.z() + s*v.z() ;
+            rhoi2 = xi*xi + yi*yi   ;
+            radi2 = rhoi2 + zi*zi   ;
+
+            if ( (radi2 <= tolORMax2)
+              && (radi2 >= tolORMin2)
+              && (zi*(eTheta - halfpi) <= 0) )
+            {
+              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
               {
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
@@ -1606 +1630 @@
             }
           }
-        }
-      }  
-    }     
-    else if ( (pTheta <tolSTheta + kAngTolerance)
-           && (fSTheta > kAngTolerance) )
-    {
-      // In tolerance of stheta
-      // If entering through solid [r,phi] => 0 to in
-      // else try 2nd root
-
-      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
-      if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<pi*.5)
-        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>pi*.5)
-        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==pi*.5) )
-      {
-        if (segPhi && rho2)  // Check phi intersection
-        {
-          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
-          if (cosPsi >= cosHDPhiIT)
-          {
-            return 0 ;
-          }
-        }
-        else
-        {
-          return 0 ;
-        }
-      }
-
-      // Not entering immediately/travelling through
-
-      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
-      b  = t2/t1 ;
-      c  = dist2STheta/t1 ;
-      d2 = b*b - c ;
-
-      if (d2 >= 0)
-      {
-        d = std::sqrt(d2) ;
-        s = -b + d ;
-        if ( (s >= kCarTolerance*0.5) && (s < snxt) && (fSTheta < pi*0.5) )
-        {
-          xi    = p.x() + s*v.x() ;
-          yi    = p.y() + s*v.y() ;
-          zi    = p.z() + s*v.z() ;
-          rhoi2 = xi*xi + yi*yi   ;
-          radi2 = rhoi2 + zi*zi   ;
-
-          if ( (radi2 <= tolORMax2)
-            && (radi2 >= tolORMin2)
-            && (zi*(fSTheta - halfpi) <= 0) )
-          {
-            if ( segPhi && rhoi2 )    // Check phi intersection
-            {
-              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-              if ( cosPsi >= cosHDPhiOT )
-              {
-                snxt = s ;
-              }
-            }
-            else
-            {
-              snxt = s ;
-            }
-          }
-        }
-      }
-    }   
-    else if ( (pTheta > tolETheta - kAngTolerance)
-           && ((fSTheta + fDTheta) < pi-kAngTolerance) )   
-    {
-
-      // In tolerance of etheta
-      // If entering through solid [r,phi] => 0 to in
-      // else try 2nd root
-
-      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
-
-      if (
-    (t2<0    && (fSTheta+fDTheta) <pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
- || (t2>=0   && (fSTheta+fDTheta) >pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
- || (v.z()>0 && (fSTheta+fDTheta)==pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
-         )
-      {
-        if (segPhi && rho2)   // Check phi intersection
-        {
-          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
-          if (cosPsi >= cosHDPhiIT)
-          {
-            return 0 ;
-          }
-        }
-        else
-        {
-          return 0 ;
-        }
-      }
-
-      // Not entering immediately/travelling through
-
-      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
-      b  = t2/t1 ;
-      c  = dist2ETheta/t1 ;
-      d2 = b*b - c ;
-
-      if (d2 >= 0)
-      {
-        d = std::sqrt(d2) ;
-        s = -b + d ;
-        
-        if ( (s >= kCarTolerance*0.5)
-          && (s < snxt) && ((fSTheta + fDTheta) > pi*0.5) )
-        {
-          xi    = p.x() + s*v.x() ;
-          yi    = p.y() + s*v.y() ;
-          zi    = p.z() + s*v.z() ;
-          rhoi2 = xi*xi + yi*yi   ;
-          radi2 = rhoi2 + zi*zi   ;
-
-          if ( (radi2 <= tolORMax2)
-            && (radi2 >= tolORMin2)
-            && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
-          {
-            if (segPhi && rhoi2)   // Check phi intersection
-            {
-              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-              if (cosPsi>=cosHDPhiOT)
-              {
-                snxt = s ;
-              }
-            }
-            else
-            {
-              snxt = s ;
-            }
-          }
-        }
-      }    
+        } 
+      }   
     }  
     else
@@ -1752 +1640 @@
       t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
       t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
-
-      b  = t2/t1;
-      c  = dist2STheta/t1 ;
-      d2 = b*b - c ;
-
-      if (d2 >= 0)
-      {
-        d = std::sqrt(d2) ;
-        s = -b + d ;    // second root
-
-        if (s >= 0 && s < snxt)
-        {
-          xi    = p.x() + s*v.x() ;
-          yi    = p.y() + s*v.y() ;
-          zi    = p.z() + s*v.z() ;
-          rhoi2 = xi*xi + yi*yi   ;
-          radi2 = rhoi2 + zi*zi   ;
-
-          if ( (radi2 <= tolORMax2)
-            && (radi2 >= tolORMin2)
-            && (zi*(fSTheta - halfpi) <= 0) )
-          {
-            if (segPhi && rhoi2)   // Check phi intersection
-            {
-              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-              if (cosPsi >= cosHDPhiOT)
+      if (t1)
+      { 
+        b  = t2/t1;
+        c  = dist2STheta/t1 ;
+        d2 = b*b - c ;
+
+        if (d2 >= 0)
+        {
+          d = std::sqrt(d2) ;
+          s = -b + d ;    // second root
+
+          if ((s >= 0) && (s < snxt))
+          {
+            xi    = p.x() + s*v.x() ;
+            yi    = p.y() + s*v.y() ;
+            zi    = p.z() + s*v.z() ;
+            rhoi2 = xi*xi + yi*yi   ;
+            radi2 = rhoi2 + zi*zi   ;
+
+            if ( (radi2 <= tolORMax2)
+              && (radi2 >= tolORMin2)
+              && (zi*(fSTheta - halfpi) <= 0) )
+            {
+              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
+              {
+                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+                if (cosPsi >= cosHDPhiOT)
+                {
+                  snxt = s ;
+                }
+              }
+              else
               {
                 snxt = s ;
               }
-            }
-            else
-            {
-              snxt = s ;
             }
           }
@@ -1791 +1681 @@
       t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
       t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
-        
-      b  = t2/t1 ;
-      c  = dist2ETheta/t1 ;
-      d2 = b*b - c ;
-
-      if (d2 >= 0)
-      {
-        d = std::sqrt(d2) ;
-        s = -b + d;    // second root
-
-        if (s >= 0 && s < snxt)
-        {
-          xi    = p.x() + s*v.x() ;
-          yi    = p.y() + s*v.y() ;
-          zi    = p.z() + s*v.z() ;
-          rhoi2 = xi*xi + yi*yi   ;
-          radi2 = rhoi2 + zi*zi   ;
-
-          if ( (radi2 <= tolORMax2)
-            && (radi2 >= tolORMin2)
-            && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
-          {
-            if (segPhi && rhoi2)   // Check phi intersection
-            {
-              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-              if ( cosPsi >= cosHDPhiOT )
-              {
-                snxt=s;
-              }
-            }
-            else
-            {
-              snxt = s ;
+      if (t1)
+      {   
+        b  = t2/t1 ;
+        c  = dist2ETheta/t1 ;
+        d2 = b*b - c ;
+
+        if (d2 >= 0)
+        {
+          d = std::sqrt(d2) ;
+          s = -b + d;    // second root
+
+          if ((s >= 0) && (s < snxt))
+          {
+            xi    = p.x() + s*v.x() ;
+            yi    = p.y() + s*v.y() ;
+            zi    = p.z() + s*v.z() ;
+            rhoi2 = xi*xi + yi*yi   ;
+            radi2 = rhoi2 + zi*zi   ;
+
+            if ( (radi2 <= tolORMax2)
+              && (radi2 >= tolORMin2)
+              && (zi*(eTheta - halfpi) <= 0) )
+            {
+              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
+              {
+                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+                if ( cosPsi >= cosHDPhiOT )
+                {
+                  snxt=s;
+                }
+              }
+              else
+              {
+                snxt = s ;
+              }
             }
           }
@@ -1844 +1736 @@
 {
   G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
-  G4double rho2,rad,rho;
-  G4double phiC,cosPhiC,sinPhiC,cosPsi,ePhi;
+  G4double rho2,rds,rho;
+  G4double cosPsi;
   G4double pTheta,dTheta1,dTheta2;
   rho2=p.x()*p.x()+p.y()*p.y();
-  rad=std::sqrt(rho2+p.z()*p.z());
+  rds=std::sqrt(rho2+p.z()*p.z());
   rho=std::sqrt(rho2);
 
@@ -1856 +1748 @@
   if (fRmin)
   {
-    safeRMin=fRmin-rad;
-    safeRMax=rad-fRmax;
+    safeRMin=fRmin-rds;
+    safeRMax=rds-fRmax;
     if (safeRMin>safeRMax)
     {
@@ -1869 +1761 @@
   else
   {
-    safe=rad-fRmax;
+    safe=rds-fRmax;
   }
 
@@ -1875 +1767 @@
   // Distance to phi extent
   //
-  if (fDPhi<twopi&&rho)
-  {
-    phiC=fSPhi+fDPhi*0.5;
-    cosPhiC=std::cos(phiC);
-    sinPhiC=std::sin(phiC);
-
+  if (!fFullPhiSphere && rho)
+  {
     // Psi=angle from central phi to point
     //
-    cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho;
-    if (cosPsi<std::cos(fDPhi*0.5))
+    cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
+    if (cosPsi<std::cos(hDPhi))
     {
       // Point lies outside phi range
       //
-      if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
-      {
-        safePhi=std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
+      if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
+      {
+        safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
       }
       else
       {
-        ePhi=fSPhi+fDPhi;
-        safePhi=std::fabs(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
-      }
-      if (safePhi>safe) safe=safePhi;
+        safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
+      }
+      if (safePhi>safe)  { safe=safePhi; }
     }
   }
@@ -1903 +1790 @@
   // Distance to Theta extent
   //    
-  if ((rad!=0.0) && (fDTheta<pi))
-  {
-    pTheta=std::acos(p.z()/rad);
-    if (pTheta<0) pTheta+=pi;
+  if ((rds!=0.0) && (!fFullThetaSphere))
+  {
+    pTheta=std::acos(p.z()/rds);
+    if (pTheta<0)  { pTheta+=pi; }
     dTheta1=fSTheta-pTheta;
-    dTheta2=pTheta-(fSTheta+fDTheta);
+    dTheta2=pTheta-eTheta;
     if (dTheta1>dTheta2)
     {
       if (dTheta1>=0)             // WHY ???????????
       {
-        safeTheta=rad*std::sin(dTheta1);
+        safeTheta=rds*std::sin(dTheta1);
         if (safe<=safeTheta)
         {
@@ -1924 +1811 @@
       if (dTheta2>=0)
       {
-        safeTheta=rad*std::sin(dTheta2);
+        safeTheta=rds*std::sin(dTheta2);
         if (safe<=safeTheta)
         {
@@ -1933 +1820 @@
   }
 
-  if (safe<0) safe=0;
+  if (safe<0)  { safe=0; }
   return safe;
 }
@@ -1939 +1826 @@
 /////////////////////////////////////////////////////////////////////
 //
-// Calculate distance to surface of shape from `inside', allowing for tolerance
+// Calculate distance to surface of shape from 'inside', allowing for tolerance
 // - Only Calc rmax intersection if no valid rmin intersection
 
@@ -1952 +1839 @@
   ESide side=kNull,sidephi=kNull,sidetheta=kNull;  
 
+  static const G4double halfCarTolerance = kCarTolerance*0.5;
+  static const G4double halfAngTolerance = kAngTolerance*0.5;
+  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
+  const G4double halfRminTolerance = fRminTolerance*0.5;
+  const G4double Rmax_plus  = fRmax + halfRmaxTolerance;
+  const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
   G4double t1,t2;
   G4double b,c,d;
@@ -1957 +1850 @@
   // Variables for phi intersection:
 
-  G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi;
-  G4double cPhi,sinCPhi,cosCPhi;
   G4double pDistS,compS,pDistE,compE,sphi2,vphi;
     
@@ -1966 +1857 @@
   G4double xi,yi,zi;      // Intersection point
 
-  // G4double Comp; // Phi intersection
-
-  G4bool segPhi;        // Phi flag and precalcs
-  G4double hDPhi,hDPhiOT,hDPhiIT; 
-  G4double cosHDPhiOT,cosHDPhiIT;
-
-  G4bool segTheta;                             // Theta flag and precals
-  G4double tanSTheta=0.,tanETheta=0., rhoSecTheta;
-  G4double tanSTheta2=0.,tanETheta2=0.;
+  // Theta precals
+  //
+  G4double rhoSecTheta;
   G4double dist2STheta, dist2ETheta, distTheta;
   G4double d2,s;
 
   // General Precalcs
-
+  //
   rho2 = p.x()*p.x()+p.y()*p.y();
   rad2 = rho2+p.z()*p.z();
-  //  G4double rad=std::sqrt(rad2);
 
   pTheta = std::atan2(std::sqrt(rho2),p.z());
@@ -1989 +1873 @@
   pDotV3d = pDotV2d+p.z()*v.z();
 
-  // Set phi divided flag and precalcs
-
-  if( fDPhi < twopi )
-  {
-    segPhi=true;
-    hDPhi=0.5*fDPhi;    // half delta phi
-    cPhi=fSPhi+hDPhi;;
-    hDPhiOT=hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi 
-    hDPhiIT=hDPhi-0.5*kAngTolerance;
-    sinCPhi=std::sin(cPhi);
-    cosCPhi=std::cos(cPhi);
-    cosHDPhiOT=std::cos(hDPhiOT);
-    cosHDPhiIT=std::cos(hDPhiIT);
-  }
-  else
-  {
-    segPhi=false;
-  }
-
   // Theta precalcs
     
-  if ( fDTheta < pi )
-  {
-    segTheta  = true;
-    tolSTheta = fSTheta - kAngTolerance*0.5;
-    tolETheta = fSTheta + fDTheta + kAngTolerance*0.5;
-  }
-  else segTheta = false;
-
+  if ( !fFullThetaSphere )
+  {
+    tolSTheta = fSTheta - halfAngTolerance;
+    tolETheta = eTheta + halfAngTolerance;
+  }
     
   // Radial Intersections from G4Sphere::DistanceToIn
@@ -2034 +1896 @@
   //
   // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
-  //
-  // const G4double  fractionTolerance = 1.0e-12;
-
-  const G4double  flexRadMaxTolerance =            // kRadTolerance;
-    std::max(kRadTolerance, fEpsilon * fRmax);
-
-  const G4double  Rmax_plus = fRmax + flexRadMaxTolerance*0.5;
-
-  const G4double  flexRadMinTolerance = std::max(kRadTolerance, 
-                     fEpsilon * fRmin);
-
-  const G4double  Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ;
-
-  if(rad2 <= Rmax_plus*Rmax_plus && rad2 >= Rmin_minus*Rmin_minus)
-    //  if(rad <= Rmax_plus && rad >= Rmin_minus)
+
+  if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
   {
     c = rad2 - fRmax*fRmax;
 
-    if (c < flexRadMaxTolerance*fRmax) 
+    if (c < fRmaxTolerance*fRmax) 
     {
       // Within tolerant Outer radius 
@@ -2065 +1914 @@
       d2 = pDotV3d*pDotV3d - c;
 
-      if( (c >- flexRadMaxTolerance*fRmax)       // on tolerant surface
-       && ((pDotV3d >=0) || (d2 < 0)) )          // leaving outside from Rmax 
-                                                 // not re-entering
+      if( (c >- fRmaxTolerance*fRmax)       // on tolerant surface
+       && ((pDotV3d >=0) || (d2 < 0)) )     // leaving outside from Rmax 
+                                            // not re-entering
       {
         if(calcNorm)
@@ -2092 +1941 @@
       d2 = pDotV3d*pDotV3d - c;
 
-      if ( c >- flexRadMinTolerance*fRmin ) // 2.0 * (0.5*kRadTolerance) * fRmin
-      {
-        if( c < flexRadMinTolerance*fRmin && 
-            d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin
-        {
-          if(calcNorm)  *validNorm = false ;   // Rmin surface is concave         
-                       return snxt = 0 ;
+      if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
+      {
+        if ( (c < fRminTolerance*fRmin)              // leaving from Rmin
+          && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
+        {
+          if(calcNorm)  { *validNorm = false; }  // Rmin surface is concave         
+          return snxt = 0 ;
         }
         else
@@ -2119 +1968 @@
   // Theta segment intersection
 
-  if (segTheta)
+  if ( !fFullThetaSphere )
   {
     // Intersection with theta surfaces
@@ -2143 +1992 @@
     //
   
-    /* ////////////////////////////////////////////////////////
-
-    tanSTheta=std::tan(fSTheta);
-    tanSTheta2=tanSTheta*tanSTheta;
-    tanETheta=std::tan(fSTheta+fDTheta);
-    tanETheta2=tanETheta*tanETheta;
-      
-    if (fSTheta)
-    {
-      dist2STheta=rho2-p.z()*p.z()*tanSTheta2;
-    }
-    else
-    {
-      dist2STheta = kInfinity;
-    }
-    if (fSTheta + fDTheta < pi)
-    {
-      dist2ETheta = rho2-p.z()*p.z()*tanETheta2;
-    }
-    else
-    {
-      dist2ETheta = kInfinity ;
-    }
-    if (pTheta > tolSTheta && pTheta < tolETheta)   // Inside theta  
-    {
-      // In tolerance of STheta and possible leaving out to small thetas N-
-
-      if(pTheta < tolSTheta + kAngTolerance  && fSTheta > kAngTolerance)  
-      {
-        t2=pDotV2d-p.z()*v.z()*tanSTheta2 ; // =(VdotN+)*rhoSecSTheta
-
-        if( fSTheta < pi*0.5 && t2 < 0)
-        {
-          if(calcNorm) *validNorm = false ;
-          return snxt = 0 ;
-        }
-        else if(fSTheta > pi*0.5 && t2 >= 0)
-        {
-          if(calcNorm)
-          {
-            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)) ;
-            *validNorm = true ;
-            *n = G4ThreeVector(-p.x()/rhoSecTheta,   // N-
-                               -p.y()/rhoSecTheta,
-                               tanSTheta/std::sqrt(1+tanSTheta2) ) ;
-          }
-          return snxt = 0 ;
-        }
-        else if( fSTheta == pi*0.5 && v.z() > 0)
-        {
-          if(calcNorm)
-          {
-            *validNorm = true ;
-            *n = G4ThreeVector(0,0,1) ;
-          }
-          return snxt = 0 ;
-        }
-      }
-
-      // In tolerance of ETheta and possible leaving out to larger thetas N+
-
-      if ( (pTheta  > tolETheta - kAngTolerance)
-        && (( fSTheta + fDTheta) < pi - kAngTolerance) )  
-      {
-        t2=pDotV2d-p.z()*v.z()*tanETheta2 ;
-        if((fSTheta+fDTheta)>pi*0.5 && t2<0)
-        {
-          if(calcNorm) *validNorm = false ;
-          return snxt = 0 ;
-        }
-        else if( (fSTheta+fDTheta) < pi*0.5 && t2 >= 0 )
-        {
-          if(calcNorm)
-          {
-            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)) ;
-            *validNorm = true ;
-            *n = G4ThreeVector( p.x()/rhoSecTheta,  // N+
-                                p.y()/rhoSecTheta,
-                                -tanETheta/std::sqrt(1+tanETheta2)  ) ; 
-          }
-          return snxt = 0 ;
-        }
-        else if( ( fSTheta+fDTheta) == pi*0.5 && v.z() < 0 )
-        {
-          if(calcNorm)
-          {
-            *validNorm = true ;
-            *n = G4ThreeVector(0,0,-1) ;
-          }
-          return snxt = 0 ;
-        }
-      }
-      if( fSTheta > 0 )
-      {       
-        // First root of fSTheta cone, second if first root -ve
-
-        t1 = 1-v.z()*v.z()*(1+tanSTheta2);
-        t2 = pDotV2d-p.z()*v.z()*tanSTheta2;
-        
-        b  = t2/t1;
-        c  = dist2STheta/t1;
-        d2 = b*b - c ;
-
-        if ( d2 >= 0 )
-        {
-          d = std::sqrt(d2) ;
-          s = -b - d ;    // First root
-
-          if ( s < 0 )
-          {
-            s = -b + d ;    // Second root
-          }
-          if (s > flexRadMaxTolerance*0.5 )   // && s<sr)
-          {
-            // check against double cone solution
-            zi=p.z()+s*v.z();
-            if (fSTheta<pi*0.5 && zi<0)
-            {
-              s = kInfinity ;  // wrong cone
-            }
-            if (fSTheta>pi*0.5 && zi>0)
-            {
-              s = kInfinity ;  // wrong cone
-            }
-            stheta = s ;
-            sidetheta = kSTheta ;
-          }
-        }
-      }
-
-      // Possible intersection with ETheta cone  
-      
-      if (fSTheta + fDTheta < pi)
-      {
-        t1 = 1-v.z()*v.z()*(1+tanETheta2);
-        t2 = pDotV2d-p.z()*v.z()*tanETheta2;        
-        b  = t2/t1;
-        c  = dist2ETheta/t1;
-        d2 = b*b-c ;
-
-        if ( d2 >= 0 )
-        {
-          d = std::sqrt(d2);
-          s = -b - d ;          // First root
-
-          if ( s < 0 )
-          {
-            s=-b+d;    // Second root
-          }
-          if (s > flexRadMaxTolerance*0.5 && s < stheta )
-          {
-            // check against double cone solution
-            zi=p.z()+s*v.z();
-            if (fSTheta+fDTheta<pi*0.5 && zi<0)
-            {
-              s = kInfinity ;  // wrong cone
-            }
-            if (fSTheta+fDTheta>pi*0.5 && zi>0)
-            {
-              s = kInfinity ;  // wrong cone
-            }
-          }
-          if (s < stheta)
-          {
-            stheta = s ;
-            sidetheta = kETheta ;
-          }
-        }
-      }
-    }  
-    */  ////////////////////////////////////////////////////////////
-
     if(fSTheta) // intersection with first cons
     {
-
-      tanSTheta = std::tan(fSTheta);
-
       if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
       {
         if( v.z() > 0. ) 
         {
-          if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
+          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
           {
             if(calcNorm)
@@ -2333 +2007 @@
             return snxt = 0 ;
           }  
-          // s = -p.z()/v.z();
           stheta    = -p.z()/v.z();
           sidetheta = kSTheta;
@@ -2340 +2013 @@
       else // kons is not plane 
       {
-        tanSTheta2  = tanSTheta*tanSTheta;
         t1          = 1-v.z()*v.z()*(1+tanSTheta2);
         t2          = pDotV2d-p.z()*v.z()*tanSTheta2;  // ~vDotN if p on cons
-        dist2STheta = rho2-p.z()*p.z()*tanSTheta2;      // t3
-
-        // distTheta = std::sqrt(std::fabs(dist2STheta/(1+tanSTheta2)));
+        dist2STheta = rho2-p.z()*p.z()*tanSTheta2;     // t3
+
         distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
 
-        if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
-        {
+        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
+        {                                      // v parallel to kons
           if( v.z() > 0. )
           {
-            if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface
-            {
-              if( fSTheta < halfpi && p.z() > 0. )
-              {
-                if( calcNorm ) *validNorm = false;
-                              return snxt = 0.;
-              }
-              else if( fSTheta > halfpi && p.z() <= 0)
+            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
+            {
+              if( (fSTheta < halfpi) && (p.z() > 0.) )
+              {
+                if( calcNorm )  { *validNorm = false; }
+                return snxt = 0.;
+              }
+              else if( (fSTheta > halfpi) && (p.z() <= 0) )
               {
                 if( calcNorm ) 
@@ -2377 +2048 @@
               }
             }
-            // s = -0.5*dist2STheta/t2;
-
             stheta    = -0.5*dist2STheta/t2;
             sidetheta = kSTheta;
           }  
-        }
-        else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
-        {
-          if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
-          {
-            if( fSTheta > halfpi && t2 >= 0. ) // leave
+        }      // 2nd order equation, 1st root of fSTheta cone,
+        else   // 2nd if 1st root -ve
+        {
+          if( std::fabs(distTheta) < halfRmaxTolerance )
+          {
+            if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
             {
               if( calcNorm ) 
@@ -2394 +2063 @@
                 if (rho2)
                 {
-                    rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
+                  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
                    
-                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
-                                        p.y()/rhoSecTheta,
-                                        std::sin(fSTheta)  );
+                  *n = G4ThreeVector( p.x()/rhoSecTheta,   
+                                      p.y()/rhoSecTheta,
+                                      std::sin(fSTheta)  );
                 }
-                else *n = G4ThreeVector(0.,0.,1.);
-              }                           
+                else  { *n = G4ThreeVector(0.,0.,1.); }
+              }
               return snxt = 0.;
             }
-            else if( fSTheta < halfpi  && t2 < 0. && p.z() >=0. ) // leave
-            {
-                if( calcNorm )   *validNorm = false;                                                 
-                return snxt = 0.;
+            else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
+            {
+              if( calcNorm )  { *validNorm = false; }
+              return snxt = 0.;
             }                               
           }
@@ -2422 +2091 @@
               s = -b - d;         // First root
 
-              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
-                              s < 0.  ||
-                  ( s > 0. && p.z() + s*v.z() > 0.)                   ) 
+              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
+               ||  (s < 0.)  || ( (s > 0.) && (p.z() + s*v.z() > 0.) )     ) 
               {
                 s = -b + d ; // 2nd root
               }
-              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.)  
+              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() <= 0.) )  
               {
                 stheta    = s;
@@ -2438 +2106 @@
               s = -b - d;         // First root
 
-              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
-                              s < 0.                                   ||
-                  ( s > 0. && p.z() + s*v.z() < 0.)                      )
+              if ( ( (std::fabs(s) < halfRmaxTolerance) && (t2 >= 0.) )
+                || (s < 0.) || ( (s > 0.) && (p.z() + s*v.z() < 0.) )   )
               {
                 s = -b + d ; // 2nd root
               }
-              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() >= 0.)  
+              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() >= 0.) )  
               {
                 stheta    = s;
@@ -2454 +2121 @@
       }
     }
-    if (fSTheta + fDTheta < pi) // intersection with second cons
-    {
-
-      tanETheta = std::tan(fSTheta+fDTheta);
-
+    if (eTheta < pi) // intersection with second cons
+    {
       if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
       {
         if( v.z() < 0. ) 
         {
-          if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
+          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
           {
             if(calcNorm)
@@ -2474 +2138 @@
           s = -p.z()/v.z();
 
-          if( s < stheta)
+          if( s < stheta )
           {
             stheta    = s;
@@ -2483 +2147 @@
       else // kons is not plane 
       {
-        tanETheta2  = tanETheta*tanETheta;
         t1          = 1-v.z()*v.z()*(1+tanETheta2);
         t2          = pDotV2d-p.z()*v.z()*tanETheta2;  // ~vDotN if p on cons
-        dist2ETheta = rho2-p.z()*p.z()*tanETheta2;      // t3
-
-        // distTheta = std::sqrt(std::fabs(dist2ETheta/(1+tanETheta2)));
+        dist2ETheta = rho2-p.z()*p.z()*tanETheta2;     // t3
+
         distTheta = std::sqrt(rho2)-p.z()*tanETheta;
 
-        if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
-        {
+        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
+        {                                      // v parallel to kons
           if( v.z() < 0. )
           {
-            if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface
-            {
-              if( fSTheta+fDTheta > halfpi && p.z() < 0. )
-              {
-                if( calcNorm ) *validNorm = false;
-                              return snxt = 0.;
-              }
-              else if( fSTheta+fDTheta < halfpi && p.z() >= 0)
+            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
+            {
+              if( (eTheta > halfpi) && (p.z() < 0.) )
+              {
+                if( calcNorm )  { *validNorm = false; }
+                return snxt = 0.;
+              }
+              else if ( (eTheta < halfpi) && (p.z() >= 0) )
               {
                 if( calcNorm ) 
@@ -2510 +2172 @@
                   {
                     rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
-                   
                     *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                         p.y()/rhoSecTheta,
-                                        -std::sin(fSTheta+fDTheta)  );
+                                        -sinETheta  );
                   }
-                  else *n = G4ThreeVector(0.,0.,-1.);
+                  else  { *n = G4ThreeVector(0.,0.,-1.); }
                 }
                 return snxt = 0.;               
@@ -2522 +2183 @@
             s = -0.5*dist2ETheta/t2;
 
-            if( s < stheta)
+            if( s < stheta )
             {
               stheta    = s;
@@ -2528 +2189 @@
             }
           }  
-        }
-        else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
-        {
-          if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
-          {
-            if( fSTheta+fDTheta < halfpi && t2 >= 0. ) // leave
+        }      // 2nd order equation, 1st root of fSTheta cone
+        else   // 2nd if 1st root -ve
+        {
+          if ( std::fabs(distTheta) < halfRmaxTolerance )
+          {
+            if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
             {
               if( calcNorm ) 
@@ -2541 +2202 @@
                 {
                     rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
-                   
                     *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                         p.y()/rhoSecTheta,
-                                        -std::sin(fSTheta+fDTheta)  );
+                                        -sinETheta  );
                 }
                 else *n = G4ThreeVector(0.,0.,-1.);
@@ -2550 +2210 @@
               return snxt = 0.;
             }
-            else if( fSTheta+fDTheta > halfpi  && t2 < 0. && p.z() <=0. ) // leave
-            {
-                if( calcNorm )   *validNorm = false;                                                 
-                return snxt = 0.;
+            else if ( (eTheta > halfpi)
+                   && (t2 < 0.) && (p.z() <=0.) ) // leave
+            {
+              if( calcNorm )  { *validNorm = false; }
+              return snxt = 0.;
             }                               
           }
@@ -2564 +2225 @@
             d = std::sqrt(d2);
 
-            if( fSTheta+fDTheta < halfpi )
+            if( eTheta < halfpi )
             {
               s = -b - d;         // First root
 
-              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
-                              s < 0. ) 
+              if( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
+               || (s < 0.) ) 
               {
                 s = -b + d ; // 2nd root
               }
-              if( s >  flexRadMaxTolerance*0.5 )  
+              if( s > halfRmaxTolerance )  
               {
                 if( s < stheta )
@@ -2586 +2247 @@
               s = -b - d;         // First root
 
-              if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
-                              s < 0.                                   ||
-                  ( s > 0. && p.z() + s*v.z() > 0.)                      )
+              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 >= 0.))
+                || (s < 0.) || ( (s > 0.) && (p.z() + s*v.z() > 0.) ) )
               {
                 s = -b + d ; // 2nd root
               }
-              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.)  
+              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() <= 0.) )  
               {
                 if( s < stheta )
@@ -2610 +2270 @@
   // Phi Intersection
     
-  if ( fDPhi < twopi)
-  {
-    sinSPhi=std::sin(fSPhi);
-    cosSPhi=std::cos(fSPhi);
-    ePhi=fSPhi+fDPhi;
-    sinEPhi=std::sin(ePhi);
-    cosEPhi=std::cos(ePhi);
-    cPhi=fSPhi+fDPhi*0.5;
-    sinCPhi=std::sin(cPhi);
-    cosCPhi=std::cos(cPhi);
-
-    if ( p.x()||p.y() ) // Check if on z axis (rho not needed later)
+  if ( !fFullPhiSphere )
+  {
+    if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
     {
       // pDist -ve when inside
@@ -2634 +2285 @@
       sidephi = kNull ;
 
-      if ( pDistS <= 0 && pDistE <= 0 )
+      if ( (pDistS <= 0) && (pDistE <= 0) )
       {
         // Inside both phi *full* planes
@@ -2644 +2295 @@
           yi   = p.y()+sphi*v.y() ;
 
-          // Check intersecting with correct half-plane
-          // (if not -> no intersect)
-
-          if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
+          // Check intersection with correct half-plane (if not -> no intersect)
+          //
+          if( (std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance) )
+          {
+            vphi = std::atan2(v.y(),v.x());
+            sidephi = kSPhi;
+            if ( ( (fSPhi-halfAngTolerance) <= vphi)
+              && ( (ePhi+halfAngTolerance)  >= vphi) )
+            {
+              sphi = kInfinity;
+            }
+          }
+          else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
           {
             sphi=kInfinity;
@@ -2654 +2314 @@
           {
             sidephi = kSPhi ;
-            if ( pDistS > -0.5*kCarTolerance) sphi =0 ; // Leave by sphi 
-          }
-        }
-        else sphi = kInfinity ;
+            if ( pDistS > -halfCarTolerance)  { sphi = 0; } // Leave by sphi 
+          }
+        }
+        else  { sphi = kInfinity; }
 
         if ( compE < 0 )
@@ -2667 +2327 @@
             yi = p.y()+sphi2*v.y() ;
 
-            // Check intersecting with correct half-plane
- 
-            if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
+            // Check intersection with correct half-plane
+            //
+            if ((std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance))
+            {
+              // Leaving via ending phi
+              //
+              vphi = std::atan2(v.y(),v.x()) ;
+               
+              if( !((fSPhi-halfAngTolerance <= vphi)
+                  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
+              { 
+                sidephi = kEPhi;
+                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
+                else                                { sphi = 0.0;   }
+              }
+            }
+            else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
             {
               sidephi = kEPhi ;
-              if ( pDistE <= -0.5*kCarTolerance )
+              if ( pDistE <= -halfCarTolerance )
               {
                 sphi=sphi2;
@@ -2684 +2358 @@
         }        
       }
-      else if ( pDistS >= 0 && pDistE >= 0 ) // Outside both *full* phi planes
+      else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
       {
         if ( pDistS <= pDistE )
@@ -2696 +2370 @@
         if ( fDPhi > pi )
         {
-          if ( compS < 0 && compE < 0 ) sphi = 0 ;
-          else                          sphi = kInfinity ;
+          if ( (compS < 0) && (compE < 0) )  { sphi = 0; }
+          else                               { sphi = kInfinity; }
         }
         else
@@ -2704 +2378 @@
           // will remain inside
 
-          if ( compS >= 0 && compE >= 0 )
-          {
-            sphi=kInfinity;
-          }
-          else
-          {
-            sphi=0;
-          }
+          if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
+          else                                { sphi = 0; }
         }    
       }
-      else if ( pDistS > 0 && pDistE < 0 )
+      else if ( (pDistS > 0) && (pDistE < 0) )
       {
         // Outside full starting plane, inside full ending plane
@@ -2729 +2397 @@
             // (if not -> not leaving phi extent)
             //
-            if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
+            if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
+            {
+              vphi = std::atan2(v.y(),v.x());
+              sidephi = kSPhi;
+              if ( ( (fSPhi-halfAngTolerance) <= vphi)
+                && ( (ePhi+halfAngTolerance)  >= vphi) )
+              {
+                sphi = kInfinity;
+              }
+            }
+            else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
             {
               sphi = kInfinity ;
@@ -2736 +2414 @@
             {
               sidephi = kEPhi ;
-              if ( pDistE > -0.5*kCarTolerance ) sphi = 0. ;
+              if ( pDistE > -halfCarTolerance )  { sphi = 0.; }
             }
           }
@@ -2757 +2435 @@
               // (if not -> remain in extent)
               //
-              if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
+              if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
+              {
+                vphi = std::atan2(v.y(),v.x());
+                sidephi = kSPhi;
+                if ( ( (fSPhi-halfAngTolerance) <= vphi)
+                  && ( (ePhi+halfAngTolerance)  >= vphi) )
+                {
+                  sphi = kInfinity;
+                }
+              }
+              else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
               {
                 sphi=kInfinity;
@@ -2787 +2475 @@
             xi=p.x()+sphi*v.x();
             yi=p.y()+sphi*v.y();
-
+            
             // Check intersection in correct half-plane
             // (if not -> not leaving phi extent)
             //
-            if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
+            if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
+            {
+              vphi = std::atan2(v.y(),v.x()) ;
+              sidephi = kSPhi;
+              if ( ( (fSPhi-halfAngTolerance) <= vphi)
+                && ( (ePhi+halfAngTolerance)  >= vphi) )
+              {
+              sphi = kInfinity;
+              }
+            }
+            else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
             {
               sphi = kInfinity ;
             }
-           else  // Leaving via Starting phi
-           {
+            else  // Leaving via Starting phi
+            {
               sidephi = kSPhi ;
-             if ( pDistS > -0.5*kCarTolerance ) sphi = 0 ;
+              if ( pDistS > -halfCarTolerance )  { sphi = 0; }
             }
           }
@@ -2815 +2513 @@
               xi   = p.x()+sphi*v.x() ;
               yi   = p.y()+sphi*v.y() ;
-
+              
               // Check intersection in correct half-plane
               // (if not -> remain in extent)
               //
-              if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
+              if((std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance))
+              {
+                vphi = std::atan2(v.y(),v.x()) ;
+                sidephi = kSPhi;
+                if ( ( (fSPhi-halfAngTolerance) <= vphi)
+                  && ( (ePhi+halfAngTolerance)  >= vphi) )
+                {
+                  sphi = kInfinity;
+                }
+              }
+              else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
               {
                 sphi = kInfinity ;
@@ -2849 +2557 @@
       {
         vphi = std::atan2(v.y(),v.x()) ;
-        if ( fSPhi < vphi && vphi < fSPhi + fDPhi )
-        {
-          sphi=kInfinity;
+        if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
+        {
+          sphi = kInfinity;
         }
         else
@@ -2859 +2567 @@
         }
       }
-      else  // travel along z - no phi intersaction
+      else  // travel along z - no phi intersection
       {
         sphi = kInfinity ;
@@ -2895 +2603 @@
         if ( fDPhi <= pi )     // Normal to Phi-
         {
-          *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
+          *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
           *validNorm=true;
         }
-        else *validNorm=false;
+        else  { *validNorm=false; }
         break ;
 
@@ -2904 +2612 @@
         if ( fDPhi <= pi )      // Normal to Phi+
         {
-          *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
+          *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
           *validNorm=true;
         }
-        else *validNorm=false;
+        else  { *validNorm=false; }
         break;
 
@@ -2920 +2628 @@
           xi = p.x() + snxt*v.x();
           yi = p.y() + snxt*v.y();
-          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2));
-          *n = G4ThreeVector( xi/rhoSecTheta,   // N-
-                              yi/rhoSecTheta,
-                              -tanSTheta/std::sqrt(1+tanSTheta2));
+          rho2=xi*xi+yi*yi;
+          if (rho2)
+          { 
+            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
+            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
+                               -tanSTheta/std::sqrt(1+tanSTheta2));
+          }
+          else
+          {
+            *n = G4ThreeVector(0.,0.,1.);
+          }
           *validNorm=true;
         }
-        else *validNorm=false;  // Concave STheta cone
+        else  { *validNorm=false; }  // Concave STheta cone
         break;
 
       case kETheta:
-        if( ( fSTheta + fDTheta ) == halfpi )
+        if( eTheta == halfpi )
         {
           *n         = G4ThreeVector(0.,0.,-1.);
           *validNorm = true;
         }
-        else if ( ( fSTheta + fDTheta ) < halfpi)
+        else if ( eTheta < halfpi )
         {
           xi=p.x()+snxt*v.x();
           yi=p.y()+snxt*v.y();
-          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2));
-          *n = G4ThreeVector( xi/rhoSecTheta,   // N+
-                              yi/rhoSecTheta,
-                              -tanETheta/std::sqrt(1+tanETheta2) );
+          rho2=xi*xi+yi*yi;
+          if (rho2)
+          { 
+            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
+            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
+                               -tanETheta/std::sqrt(1+tanETheta2) );
+          }
+          else
+          {
+            *n = G4ThreeVector(0.,0.,-1.);
+          }
           *validNorm=true;
         }
-        else *validNorm=false;   // Concave ETheta cone
+        else  { *validNorm=false; }   // Concave ETheta cone
         break;
 
@@ -2995 +2717 @@
 /////////////////////////////////////////////////////////////////////////
 //
-// Calcluate distance (<=actual) to closest surface of shape from inside
+// Calculate distance (<=actual) to closest surface of shape from inside
 
 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
 {
   G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
-  G4double rho2,rad,rho;
-  G4double phiC,cosPhiC,sinPhiC,ePhi;
+  G4double rho2,rds,rho;
   G4double pTheta,dTheta1,dTheta2;
   rho2=p.x()*p.x()+p.y()*p.y();
-  rad=std::sqrt(rho2+p.z()*p.z());
+  rds=std::sqrt(rho2+p.z()*p.z());
   rho=std::sqrt(rho2);
 
@@ -3027 +2748 @@
   if (fRmin)
   {
-    safeRMin=rad-fRmin;
-    safeRMax=fRmax-rad;
+    safeRMin=rds-fRmin;
+    safeRMax=fRmax-rds;
     if (safeRMin<safeRMax)
     {
@@ -3040 +2761 @@
   else
   {
-    safe=fRmax-rad;
+    safe=fRmax-rds;
   }
 
@@ -3046 +2767 @@
   // Distance to phi extent
   //
-  if (fDPhi<twopi && rho)
-  {
-    phiC=fSPhi+fDPhi*0.5;
-    cosPhiC=std::cos(phiC);
-    sinPhiC=std::sin(phiC);
-    if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
-    {
-      safePhi=-(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
+  if (!fFullPhiSphere && rho)
+  {
+    if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
+    {
+      safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
     }
     else
     {
-      ePhi=fSPhi+fDPhi;
-      safePhi=(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
-    }
-    if (safePhi<safe) safe=safePhi;
+      safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
+    }
+    if (safePhi<safe)  { safe=safePhi; }
   }
 
@@ -3066 +2783 @@
   // Distance to Theta extent
   //    
-  if (rad)
-  {
-    pTheta=std::acos(p.z()/rad);
-    if (pTheta<0) pTheta+=pi;
+  if (rds)
+  {
+    pTheta=std::acos(p.z()/rds);
+    if (pTheta<0)  { pTheta+=pi; }
     dTheta1=pTheta-fSTheta;
-    dTheta2=(fSTheta+fDTheta)-pTheta;
-    if (dTheta1<dTheta2)
-    {
-      safeTheta=rad*std::sin(dTheta1);
-      if (safe>safeTheta)
-      {
-        safe=safeTheta;
-      }
-    }
-    else
-    {
-      safeTheta=rad*std::sin(dTheta2);
-      if (safe>safeTheta)
-      {
-        safe=safeTheta;
-      }
-    }
-  }
-
-  if (safe<0) safe=0;
-    return safe;
+    dTheta2=eTheta-pTheta;
+    if (dTheta1<dTheta2)  { safeTheta=rds*std::sin(dTheta1); }
+    else                  { safeTheta=rds*std::sin(dTheta2); }
+    if (safe>safeTheta)   { safe=safeTheta; }
+  }
+
+  if (safe<0)  { safe=0; }
+  return safe;
 }
 
@@ -3119 +2823 @@
   // Phi cross sections
     
-  noPhiCrossSections=G4int (fDPhi/kMeshAngleDefault)+1;
+  noPhiCrossSections = G4int(fDPhi/kMeshAngleDefault)+1;
     
   if (noPhiCrossSections<kMinMeshSections)
@@ -3134 +2838 @@
   // on the x axis. Will give better extent calculations when not rotated.
     
-  if (fDPhi==pi*2.0 && fSPhi==0)
+  if (fFullPhiSphere)
   {
     sAnglePhi = -meshAnglePhi*0.5;
@@ -3160 +2864 @@
   // on the z axis. Will give better extent calculations when not rotated.
     
-  if (fDTheta==pi && fSTheta==0)
+  if (fFullThetaSphere)
   {
     startTheta = -meshTheta*0.5;
@@ -3223 +2927 @@
   }
 
-  delete[] cosCrossTheta;
-  delete[] sinCrossTheta;
+  delete [] cosCrossTheta;
+  delete [] sinCrossTheta;
 
   return vertices;
@@ -3269 +2973 @@
   G4double height1, height2, slant1, slant2, costheta, sintheta,theta,rRand;
 
-  height1 = (fRmax-fRmin)*std::cos(fSTheta);
-  height2 = (fRmax-fRmin)*std::cos(fSTheta+fDTheta);
-  slant1  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta))
-                      + height1*height1);
-  slant2  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta+fDTheta))
-                      + height2*height2);
+  height1 = (fRmax-fRmin)*cosSTheta;
+  height2 = (fRmax-fRmin)*cosETheta;
+  slant1  = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
+  slant2  = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
   rRand   = RandFlat::shoot(fRmin,fRmax);
   
-  aOne = fRmax*fRmax*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
-  aTwo = fRmin*fRmin*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
-  aThr = fDPhi*((fRmax + fRmin)*std::sin(fSTheta))*slant1;
-  aFou = fDPhi*((fRmax + fRmin)*std::sin(fSTheta+fDTheta))*slant2;
+  aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
+  aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
+  aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
+  aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
   aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
   
-  phi = RandFlat::shoot(fSPhi, fSPhi + fDPhi); 
+  phi = RandFlat::shoot(fSPhi, ePhi); 
   cosphi = std::cos(phi); 
   sinphi = std::sin(phi);
-  theta = RandFlat::shoot(fSTheta,fSTheta+fDTheta);
+  theta = RandFlat::shoot(fSTheta,eTheta);
   costheta = std::cos(theta);
   sintheta = std::sqrt(1.-sqr(costheta));
 
-  if( ((fSPhi==0) && (fDPhi==2.*pi)) || (fDPhi==2.*pi) ) {aFiv = 0;}
-  if(fSTheta == 0)  {aThr=0;}
-  if(fDTheta + fSTheta == pi) {aFou = 0;}
-  if(fSTheta == 0.5*pi) {aThr = pi*(fRmax*fRmax-fRmin*fRmin);}
-  if(fSTheta + fDTheta == 0.5*pi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin);}
+  if(fFullPhiSphere) { aFiv = 0; }
+  if(fSTheta == 0)   { aThr=0; }
+  if(eTheta == pi) { aFou = 0; }
+  if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
+  if(eTheta == halfpi)  { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }
 
   chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
@@ -3309 +3011 @@
   else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
   {
-    if (fSTheta != 0.5*pi)
-    {
-      zRand = RandFlat::shoot(fRmin*std::cos(fSTheta),fRmax*std::cos(fSTheta));
-      return G4ThreeVector(std::tan(fSTheta)*zRand*cosphi,
-                           std::tan(fSTheta)*zRand*sinphi,zRand);
+    if (fSTheta != halfpi)
+    {
+      zRand = RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
+      return G4ThreeVector(tanSTheta*zRand*cosphi,
+                           tanSTheta*zRand*sinphi,zRand);
     }
     else
@@ -3322 +3024 @@
   else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
   {
-    if(fSTheta + fDTheta != 0.5*pi)
-    {
-      zRand = RandFlat::shoot(fRmin*std::cos(fSTheta+fDTheta),
-                              fRmax*std::cos(fSTheta+fDTheta));
-      return G4ThreeVector  (std::tan(fSTheta+fDTheta)*zRand*cosphi,
-                             std::tan(fSTheta+fDTheta)*zRand*sinphi,zRand);
+    if(eTheta != halfpi)
+    {
+      zRand = RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
+      return G4ThreeVector  (tanETheta*zRand*cosphi,
+                             tanETheta*zRand*sinphi,zRand);
     }
     else
@@ -3336 +3037 @@
   else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
   {
-    return G4ThreeVector(rRand*sintheta*std::cos(fSPhi),
-                         rRand*sintheta*std::sin(fSPhi),rRand*costheta);
+    return G4ThreeVector(rRand*sintheta*cosSPhi,
+                         rRand*sintheta*sinSPhi,rRand*costheta);
   }
   else
   {
-    return G4ThreeVector(rRand*sintheta*std::cos(fSPhi+fDPhi),
-                         rRand*sintheta*std::sin(fSPhi+fDPhi),rRand*costheta);
-  }
+    return G4ThreeVector(rRand*sintheta*cosEPhi,
+                         rRand*sintheta*sinEPhi,rRand*costheta);
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//
+// GetSurfaceArea
+
+G4double G4Sphere::GetSurfaceArea()
+{
+  if(fSurfaceArea != 0.) {;}
+  else
+  {   
+    G4double Rsq=fRmax*fRmax;
+    G4double rsq=fRmin*fRmin;
+         
+    fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
+    if(!fFullPhiSphere)
+    {
+      fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
+    }
+    if(fSTheta >0)
+    {
+      G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
+                              + std::pow(cosSTheta,2));
+      if(fDPhi>pi)
+      { 
+        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
+      }
+      else
+      {
+        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
+      }
+    }
+    if(eTheta < pi)
+    {
+      G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
+                              + std::pow(cosETheta,2));
+      if(fDPhi>pi)
+      { 
+        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
+      }
+      else
+      {
+        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
+      }
+    }
+  }
+  return fSurfaceArea;
 }
 

trunk/source/geometry/solids/CSG/src/G4Torus.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4Torus.cc,v 1.63 2007/10/02 09:34:17 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// $Id: G4Torus.cc,v 1.65 2009/11/26 10:31:06 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // 
@@ -284 +284 @@
 
     G4double theta = std::atan2(ptmp.y(),ptmp.x());
-
-    if (theta < 0)  { theta += twopi; }
     
+    if ( fSPhi >= 0 )
+    {
+      if ( theta < - kAngTolerance*0.5 )  { theta += twopi; }
+      if ( (std::abs(theta) < kAngTolerance*0.5)
+        && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
+      { 
+        theta += twopi ; // 0 <= theta < 2pi
+      }
+    }
+    if ((fSPhi <= -pi )&&(theta>kAngTolerance*0.5)) { theta = theta-twopi; }
+       
     // We have to verify if this root is inside the region between
     // fSPhi and fSPhi + fDPhi

trunk/source/geometry/solids/CSG/src/G4Trap.cc

@@ -26 +26 @@
 //
 // $Id: G4Trap.cc,v 1.45 2008/04/23 09:49:57 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // class G4Trap

trunk/source/geometry/solids/CSG/src/G4Trd.cc

@@ -26 +26 @@
 //
 // $Id: G4Trd.cc,v 1.34 2006/10/19 15:33:38 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 //

trunk/source/geometry/solids/CSG/src/G4Tubs.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4Tubs.cc,v 1.74 2008/11/06 15:26:53 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-02-ref-02 $
+// $Id: G4Tubs.cc,v 1.79 2009/06/30 10:10:11 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // 
@@ -90 +90 @@
                       G4double pDz,
                       G4double pSPhi, G4double pDPhi )
-  : G4CSGSolid(pName)
+  : G4CSGSolid(pName), fSPhi(0), fDPhi(0)
 {
 
@@ -102 +102 @@
   else
   {
-    G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
-           << "        Negative Z half-length ! - "
-           << pDz << G4endl;
+    G4cerr << "ERROR - G4Tubs()::G4Tubs()" << G4endl
+           << "        Negative Z half-length (" << pDz << ") in solid: "
+           << GetName() << G4endl;
     G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
                 "Invalid Z half-length");
@@ -115 +115 @@
   else
   {
-    G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
-           << "        Invalid values for radii !" << G4endl
+    G4cerr << "ERROR - G4Tubs()::G4Tubs()" << G4endl
+           << "        Invalid values for radii in solid " << GetName()
+           << G4endl
            << "        pRMin = " << pRMin << ", pRMax = " << pRMax << G4endl;
     G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
@@ -122 +123 @@
   }
 
-  fPhiFullTube = true;
-  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
-  {
-    fDPhi=twopi;
-    fSPhi=0;
-  }
-  else
-  {
-    fPhiFullTube = false;
-    if ( pDPhi > 0 )
-    {
-      fDPhi = pDPhi;
-    }
-    else
-    {
-      G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
-             << "        Negative delta-Phi ! - "
-             << pDPhi << G4endl;
-      G4Exception("G4Tubs::G4Tubs()", "InvalidSetup",
-                  FatalException, "Invalid dphi.");
-    }
-
-    // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
-
-    if ( pSPhi < 0 )
-    {
-      fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi);
-    }
-    else
-    {
-      fSPhi = std::fmod(pSPhi,twopi) ;
-    }
-    if ( fSPhi+fDPhi > twopi )
-    {
-      fSPhi -= twopi ;
-    }
-  }
-  InitializeTrigonometry();
+  // Check angles
+
+  CheckPhiAngles(pSPhi, pDPhi);
 }
 
@@ -203 +169 @@
 {
 
-  if ( !pTransform.IsRotated() && fDPhi == twopi && fRMin == 0 )
+  if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
   {
     // Special case handling for unrotated solid tubes
@@ -210 +176 @@
     // and compute exact x and y limit for x/y case
       
-    G4double xoffset, xMin, xMax ;
-    G4double yoffset, yMin, yMax ;
-    G4double zoffset, zMin, zMax ;
-
-    G4double diff1, diff2, maxDiff, newMin, newMax ;
-    G4double xoff1, xoff2, yoff1, yoff2, delta ;
-
-    xoffset = pTransform.NetTranslation().x() ;
-    xMin = xoffset - fRMax ;
-    xMax = xoffset + fRMax ;
+    G4double xoffset, xMin, xMax;
+    G4double yoffset, yMin, yMax;
+    G4double zoffset, zMin, zMax;
+
+    G4double diff1, diff2, maxDiff, newMin, newMax;
+    G4double xoff1, xoff2, yoff1, yoff2, delta;
+
+    xoffset = pTransform.NetTranslation().x();
+    xMin = xoffset - fRMax;
+    xMax = xoffset + fRMax;
 
     if (pVoxelLimit.IsXLimited())
@@ -230 +196 @@
       else
       {
-        if ( xMin < pVoxelLimit.GetMinXExtent() )
-        {
-          xMin = pVoxelLimit.GetMinXExtent() ;
-        }
-        if (xMax > pVoxelLimit.GetMaxXExtent() )
-        {
-          xMax = pVoxelLimit.GetMaxXExtent() ;
-        }
-      }
-    }
-    yoffset = pTransform.NetTranslation().y() ;
-    yMin    = yoffset - fRMax ;
-    yMax    = yoffset + fRMax ;
+        if (xMin < pVoxelLimit.GetMinXExtent())
+        {
+          xMin = pVoxelLimit.GetMinXExtent();
+        }
+        if (xMax > pVoxelLimit.GetMaxXExtent())
+        {
+          xMax = pVoxelLimit.GetMaxXExtent();
+        }
+      }
+    }
+    yoffset = pTransform.NetTranslation().y();
+    yMin    = yoffset - fRMax;
+    yMax    = yoffset + fRMax;
 
     if ( pVoxelLimit.IsYLimited() )
@@ -249 +215 @@
         || (yMax < pVoxelLimit.GetMinYExtent()) )
       {
-        return false ;
+        return false;
       }
       else
       {
-        if ( yMin < pVoxelLimit.GetMinYExtent() )
-        {
-          yMin = pVoxelLimit.GetMinYExtent() ;
-        }
-        if ( yMax > pVoxelLimit.GetMaxYExtent() )
+        if (yMin < pVoxelLimit.GetMinYExtent())
+        {
+          yMin = pVoxelLimit.GetMinYExtent();
+        }
+        if (yMax > pVoxelLimit.GetMaxYExtent())
         {
           yMax=pVoxelLimit.GetMaxYExtent();
@@ -263 +229 @@
       }
     }
-    zoffset = pTransform.NetTranslation().z() ;
-    zMin    = zoffset - fDz ;
-    zMax    = zoffset + fDz ;
+    zoffset = pTransform.NetTranslation().z();
+    zMin    = zoffset - fDz;
+    zMax    = zoffset + fDz;
 
     if ( pVoxelLimit.IsZLimited() )
@@ -272 +238 @@
         || (zMax < pVoxelLimit.GetMinZExtent()) )
       {
-        return false ;
+        return false;
       }
       else
       {
-        if ( zMin < pVoxelLimit.GetMinZExtent() )
-        {
-          zMin = pVoxelLimit.GetMinZExtent() ;
-        }
-        if ( zMax > pVoxelLimit.GetMaxZExtent() )
+        if (zMin < pVoxelLimit.GetMinZExtent())
+        {
+          zMin = pVoxelLimit.GetMinZExtent();
+        }
+        if (zMax > pVoxelLimit.GetMaxZExtent())
         {
           zMax = pVoxelLimit.GetMaxZExtent();
@@ -290 +256 @@
       case kXAxis :
       {
-        yoff1 = yoffset - yMin ;
-        yoff2 = yMax    - yoffset ;
+        yoff1 = yoffset - yMin;
+        yoff2 = yMax    - yoffset;
 
         if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
         {                                   // => no change
-          pMin = xMin ;
-          pMax = xMax ;
+          pMin = xMin;
+          pMax = xMax;
         }
         else
@@ -317 +283 @@
       case kYAxis :
       {
-        xoff1 = xoffset - xMin ;
-        xoff2 = xMax - xoffset ;
+        xoff1 = xoffset - xMin;
+        xoff2 = xMax - xoffset;
 
         if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
         {                                   // => no change
-          pMin = yMin ;
-          pMax = yMax ;
+          pMin = yMin;
+          pMax = yMax;
         }
         else
@@ -334 +300 @@
           delta   = fRMax*fRMax - xoff2*xoff2;
           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
-          maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
-          newMin  = yoffset - maxDiff ;
-          newMax  = yoffset + maxDiff ;
-          pMin    = (newMin < yMin) ? yMin : newMin ;
-          pMax     =(newMax > yMax) ? yMax : newMax ;
-        }
-        break ;
+          maxDiff = (diff1 > diff2) ? diff1 : diff2;
+          newMin  = yoffset - maxDiff;
+          newMax  = yoffset + maxDiff;
+          pMin    = (newMin < yMin) ? yMin : newMin;
+          pMax    = (newMax > yMax) ? yMax : newMax;
+        }
+        break;
       }
       case kZAxis:
       {
-        pMin = zMin ;
-        pMax = zMax ;
-        break ;
+        pMin = zMin;
+        pMax = zMax;
+        break;
       }
       default:
         break;
     }
-    pMin -= kCarTolerance ;
-    pMax += kCarTolerance ;
-    return true;    
+    pMin -= kCarTolerance;
+    pMax += kCarTolerance;
+    return true;
   }
   else // Calculate rotated vertex coordinates
   {
-    G4int i, noEntries, noBetweenSections4 ;
-    G4bool existsAfterClip = false ;
-    G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
+    G4int i, noEntries, noBetweenSections4;
+    G4bool existsAfterClip = false;
+    G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
+
+    pMin =  kInfinity;
+    pMax = -kInfinity;
+
+    noEntries = vertices->size();
+    noBetweenSections4 = noEntries - 4;
     
-    pMin = +kInfinity ;
-    pMax = -kInfinity ;
-
-    noEntries = vertices->size() ;
-    noBetweenSections4 = noEntries - 4 ;
-    
-    for (i = 0 ; i < noEntries ; i += 4 )
-    {
-      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
-    }
-    for (i = 0 ; i < noBetweenSections4 ; i += 4 )
-    {
-      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
-    }
-    if ((pMin != kInfinity) || (pMax != -kInfinity) )
-    {
-      existsAfterClip = true ;
-      pMin -= kCarTolerance ; // Add 2*tolerance to avoid precision troubles
-      pMax += kCarTolerance ;
+    for ( i = 0 ; i < noEntries ; i += 4 )
+    {
+      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
+    }
+    for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
+    {
+      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
+    }
+    if ( (pMin != kInfinity) || (pMax != -kInfinity) )
+    {
+      existsAfterClip = true;
+      pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
+      pMax += kCarTolerance;
     }
     else
@@ -391 +357 @@
              (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
              (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
-             (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
+             (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
         
       if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
       {
-        existsAfterClip = true ;
-        pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
-        pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
+        existsAfterClip = true;
+        pMin            = pVoxelLimit.GetMinExtent(pAxis);
+        pMax            = pVoxelLimit.GetMaxExtent(pAxis);
       }
     }
@@ -808 +774 @@
   G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
   G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
+  const G4double dRmax = 100.*fRMax;
 
   static const G4double halfCarTolerance = 0.5*kCarTolerance;
@@ -908 +875 @@
         if (s >= 0)  // If 'forwards'
         {
+          if ( s>dRmax ) // Avoid rounding errors due to precision issues on
+          {              // 64 bits systems. Split long distances and recompute
+            G4double fTerm = s-std::fmod(s,dRmax);
+            s = fTerm + DistanceToIn(p+fTerm*v,v);
+          } 
           // Check z intersection
           //
@@ -1022 +994 @@
           //
           if(s < 0.0)  { s = 0.0; }
+          if ( s>dRmax ) // Avoid rounding errors due to precision issues seen
+          {              // 64 bits systems. Split long distances and recompute
+            G4double fTerm = s-std::fmod(s,dRmax);
+            s = fTerm + DistanceToIn(p+fTerm*v,v);
+          } 
           zi = p.z() + s*v.z() ;
           if (std::fabs(zi) <= tolODz)