Index: trunk/source/geometry/solids/CSG/src/G4Cons.cc
===================================================================
--- trunk/source/geometry/solids/CSG/src/G4Cons.cc	(revision 850)
+++ trunk/source/geometry/solids/CSG/src/G4Cons.cc	(revision 921)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Cons.cc,v 1.56 2008/02/20 08:56:16 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Cons.cc,v 1.60 2008/11/06 15:26:53 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 //
@@ -87,5 +87,7 @@
 
   if ( pDz > 0 )
-    fDz = pDz ;
+  {
+    fDz = pDz;
+  }
   else
   {
@@ -99,5 +101,5 @@
   // Check radii
 
-  if ( pRmin1 < pRmax1 && pRmin2 < pRmax2 && pRmin1 >= 0 && pRmin2 >= 0 )
+  if ( (pRmin1<pRmax1) && (pRmin2<pRmax2) && (pRmin1>=0) && (pRmin2>=0) )
   {
 
@@ -106,6 +108,6 @@
     fRmin2 = pRmin2 ; 
     fRmax2 = pRmax2 ;
-    if( (pRmin1 == 0.0 && pRmin2 > 0.0) ) fRmin1 = 1e3*kRadTolerance ; 
-    if( (pRmin2 == 0.0 && pRmin1 > 0.0) ) fRmin2 = 1e3*kRadTolerance ; 
+    if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
+    if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
   }
   else
@@ -119,7 +121,6 @@
   }
 
-  // Check angles
-
-  if ( pDPhi >= twopi )
+  fPhiFullCone = true;
+  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
   {
     fDPhi=twopi;
@@ -128,5 +129,9 @@
   else
   {
-    if ( pDPhi > 0 ) fDPhi = pDPhi ;
+    fPhiFullCone = false;
+    if ( pDPhi > 0 )
+    {
+      fDPhi = pDPhi;
+    }
     else
     {
@@ -135,14 +140,23 @@
              << pDPhi << G4endl;
       G4Exception("G4Cons::G4Cons()", "InvalidSetup",
-                  FatalException, "Invalid pDPhi.") ;
-    }
-
-    // Ensure pSPhi 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 ;
-  }
+                  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();
 }
 
@@ -173,8 +187,11 @@
   G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
   EInside in;
-
-  if (std::fabs(p.z()) > fDz + kCarTolerance*0.5 ) return in = kOutside;
-  else if(std::fabs(p.z()) >= fDz - kCarTolerance*0.5 )   in = kSurface;
-  else                                                    in = kInside;
+  static const G4double halfCarTolerance=kCarTolerance*0.5;
+  static const G4double halfRadTolerance=kRadTolerance*0.5;
+  static const G4double halfAngTolerance=kAngTolerance*0.5;
+
+  if (std::fabs(p.z()) > fDz + halfCarTolerance )  { return in = kOutside; }
+  else if(std::fabs(p.z()) >= fDz - halfCarTolerance )    { in = kSurface; }
+  else                                                    { in = kInside;  }
 
   r2 = p.x()*p.x() + p.y()*p.y() ;
@@ -184,37 +201,35 @@
   // rh2 = rh*rh;
 
-  tolRMin = rl - kRadTolerance*0.5 ;
-  if ( tolRMin < 0 ) tolRMin = 0 ;
-  tolRMax = rh + kRadTolerance*0.5 ;
-
-  if ( r2 < tolRMin*tolRMin || r2 > tolRMax*tolRMax ) return in = kOutside;
-
-  if (rl) tolRMin = rl + kRadTolerance*0.5 ;
-  else    tolRMin = 0.0 ;
-  tolRMax         = rh - kRadTolerance*0.5 ;
+  tolRMin = rl - halfRadTolerance;
+  if ( tolRMin < 0 )  { tolRMin = 0; }
+  tolRMax = rh + halfRadTolerance;
+
+  if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
+
+  if (rl) { tolRMin = rl + halfRadTolerance; }
+  else    { tolRMin = 0.0; }
+  tolRMax = rh - halfRadTolerance;
       
   if (in == kInside) // else it's kSurface already
   {
-     if (r2 < tolRMin*tolRMin || r2 >= tolRMax*tolRMax) in = kSurface;
-    //  if (r2 <= tolRMin*tolRMin || r2-rh2 >= -rh*kRadTolerance) in = kSurface;
-  }
-  if ( ( fDPhi < twopi - kAngTolerance ) &&
-       ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
+     if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
+  }
+  if ( !fPhiFullCone && ((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 - halfAngTolerance  )             { pPhi += twopi; }
+    else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
     
-    if ( (pPhi < fSPhi - kAngTolerance*0.5) ||          
-         (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) )   return in = kOutside;
+    if ( (pPhi < fSPhi - halfAngTolerance) ||          
+         (pPhi > fSPhi + fDPhi + 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 ;
-    }
-  }
-  else if( fDPhi < twopi - kAngTolerance )  in = kSurface ;
+       if ( (pPhi < fSPhi + halfAngTolerance) || 
+            (pPhi > fSPhi + fDPhi - halfAngTolerance) )  { in = kSurface; }
+    }
+  }
+  else if ( !fPhiFullCone )  { in = kSurface; }
 
   return in ;
@@ -244,6 +259,6 @@
                     G4double&          pMax  ) const
 {
-  if ( !pTransform.IsRotated() && 
-        fDPhi == twopi && fRmin1 == 0 && fRmin2 == 0 )
+  if ( !pTransform.IsRotated() && (fDPhi == twopi)
+    && (fRmin1 == 0) && (fRmin2 == 0) )
   {
     // Special case handling for unrotated solid cones
@@ -265,6 +280,6 @@
     if (pVoxelLimit.IsZLimited())
     {
-      if( zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance || 
-          zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance   )
+      if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) || 
+          (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance)  )
       {
         return false ;
@@ -290,6 +305,6 @@
     if (pVoxelLimit.IsXLimited())
     {
-      if ( xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance || 
-           xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance    )
+      if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) || 
+           (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance)  )
       {
         return false ;
@@ -315,6 +330,6 @@
     if (pVoxelLimit.IsYLimited())
     {
-      if ( yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance || 
-           yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance     )
+      if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) || 
+           (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance)  )
       {
         return false ;
@@ -338,6 +353,6 @@
         yoff2 = yMax - yoffset ;
 
-        if (yoff1 >= 0 && yoff2 >= 0) // Y limits cross max/min x => no change
-        {
+        if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
+        {                                 // => no change
           pMin = xMin ;
           pMax = xMax ;
@@ -362,6 +377,6 @@
         xoff2 = xMax - xoffset ;
 
-        if (xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change
-        {
+        if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
+        {                                  // => no change
           pMin = yMin ;
           pMax = yMax ;
@@ -409,11 +424,11 @@
     for ( i = 0 ; i < noEntries ; i += 4 )
     {
-      ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
+      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
     }
     for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
     {
-      ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
+      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
     }    
-    if ( pMin != kInfinity || pMax != -kInfinity )
+    if ( (pMin != kInfinity) || (pMax != -kInfinity) )
     {
       existsAfterClip = true ;
@@ -462,7 +477,9 @@
   G4double tanRMin, secRMin, pRMin, widRMin;
   G4double tanRMax, secRMax, pRMax, widRMax;
-  G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
+
+  static const G4double delta  = 0.5*kCarTolerance;
+  static const G4double dAngle = 0.5*kAngTolerance;
   
-  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.0);
+  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
   G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
 
@@ -482,5 +499,5 @@
   distRMax = std::fabs(pRMax - widRMax)/secRMax;
 
-  if (fDPhi < twopi)   //  &&  rho ) // Protected against (0,0,z) 
+  if (!fPhiFullCone)   // Protected against (0,0,z) 
   {
     if ( rho )
@@ -488,9 +505,9 @@
       pPhi = std::atan2(p.y(),p.x());
 
-      if(pPhi  < fSPhi-delta)           pPhi     += twopi;
-      else if(pPhi > fSPhi+fDPhi+delta) pPhi     -= twopi;
+      if (pPhi  < fSPhi-delta)            { pPhi += twopi; }
+      else if (pPhi > fSPhi+fDPhi+delta)  { pPhi -= twopi; }
 
       distSPhi = std::fabs( pPhi - fSPhi ); 
-      distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 
+      distEPhi = std::fabs( pPhi - fSPhi - fDPhi ); 
     }
     else if( !(fRmin1) || !(fRmin2) )
@@ -499,11 +516,14 @@
       distEPhi = 0.; 
     }
-    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
-    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
+    nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
+    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
   }
   if ( rho > delta )   
   {
-    nR = G4ThreeVector(p.x()/rho/secRMax,p.y()/rho/secRMax,-tanRMax/secRMax);
-    if (fRmin1 || fRmin2) nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
+    nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
+    if (fRmin1 || fRmin2)
+    {
+      nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
+    }
   }
 
@@ -513,10 +533,10 @@
     sumnorm += nR;
   }
-  if( (fRmin1 || fRmin2) && distRMin <= delta )
+  if( (fRmin1 || fRmin2) && (distRMin <= delta) )
   {
     noSurfaces ++;
     sumnorm += nr;
   }
-  if( fDPhi < twopi )   
+  if( !fPhiFullCone )   
   {
     if (distSPhi <= dAngle)
@@ -534,6 +554,6 @@
   {
     noSurfaces ++;
-    if ( p.z() >= 0.)  sumnorm += nZ;
-    else               sumnorm -= nZ; 
+    if ( p.z() >= 0.)  { sumnorm += nZ; }
+    else               { sumnorm -= nZ; }
   }
   if ( noSurfaces == 0 )
@@ -545,10 +565,11 @@
      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 ;
 }
 
-//////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////
 //
 // Algorithm for SurfaceNormal() following the original specification
@@ -605,12 +626,12 @@
     }
   }
-  if ( fDPhi < twopi && rho )  // Protected against (0,0,z) 
+  if ( !fPhiFullCone && rho )  // Protected against (0,0,z) 
   {
     phi = std::atan2(p.y(),p.x()) ;
 
-    if (phi < 0) phi += twopi ;
-
-    if (fSPhi < 0) distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
-    else           distSPhi = std::fabs(phi - fSPhi)*rho ;
+    if (phi < 0)  { phi += twopi; }
+
+    if (fSPhi < 0)  { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
+    else            { distSPhi = std::fabs(phi - fSPhi)*rho; }
 
     distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
@@ -620,9 +641,9 @@
     if (distSPhi < distEPhi)
     {
-      if (distSPhi < distMin) side = kNSPhi ;
+      if (distSPhi < distMin)  { side = kNSPhi; }
     }
     else 
     {
-      if (distEPhi < distMin) side = kNEPhi ;
+      if (distEPhi < distMin)  { side = kNEPhi; }
     }
   }    
@@ -631,19 +652,19 @@
     case kNRMin:      // Inner radius
       rho *= secRMin ;
-      norm = G4ThreeVector(-p.x()/rho,-p.y()/rho,tanRMin/secRMin) ;
+      norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
       break ;
     case kNRMax:      // Outer radius
       rho *= secRMax ;
-      norm = G4ThreeVector(p.x()/rho,p.y()/rho,-tanRMax/secRMax) ;
+      norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
       break ;
     case kNZ:      // +/- dz
-      if (p.z() > 0) norm = G4ThreeVector(0,0,1) ; 
-      else           norm = G4ThreeVector(0,0,-1) ; 
+      if (p.z() > 0)  { norm = G4ThreeVector(0,0,1);  }
+      else            { norm = G4ThreeVector(0,0,-1); }
       break ;
     case kNSPhi:
-      norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
+      norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
       break ;
     case kNEPhi:
-      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
+      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
       break ;
     default:
@@ -677,5 +698,4 @@
 //
 // NOTE:
-// - Precalculations for phi trigonometry are Done `just in time'
 // - `if valid' implies tolerant checking of intersection points
 // - z, phi intersection from Tubs
@@ -686,8 +706,6 @@
   G4double snxt = kInfinity ;      // snxt = default return value
 
-  G4bool seg ;        // true if segmented in phi
-  G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0. ;
-          // half dphi + outer tolerance
-  G4double cPhi,sinCPhi=0.,cosCPhi=0. ;  // central phi
+  static const G4double halfCarTolerance=kCarTolerance*0.5;
+  static const G4double halfRadTolerance=kRadTolerance*0.5;
 
   G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;  // Data for cones
@@ -704,23 +722,4 @@
   G4double nt1,nt2,nt3 ;
   G4double Comp ;
-  G4double cosSPhi,sinSPhi ;    // Trig for phi start intersect
-  G4double ePhi,cosEPhi,sinEPhi ;  // for phi end intersect
-
-  //
-  // Set phi divided flag and precalcs
-  //
-  if (fDPhi < twopi)
-  {
-    seg        = true ;
-    hDPhi      = 0.5*fDPhi ;    // half delta phi
-    cPhi       = fSPhi + hDPhi ; ;
-    hDPhiOT    = hDPhi + 0.5*kAngTolerance ;  // outers tol' 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     seg = false ;
 
   // Cone Precalcs
@@ -730,8 +729,8 @@
   rMinAv  = (fRmin1 + fRmin2)*0.5 ;
 
-  if (rMinAv > kRadTolerance*0.5)
-  {
-    rMinOAv = rMinAv - kRadTolerance*0.5 ;
-    rMinIAv = rMinAv + kRadTolerance*0.5 ;
+  if (rMinAv > halfRadTolerance)
+  {
+    rMinOAv = rMinAv - halfRadTolerance ;
+    rMinIAv = rMinAv + halfRadTolerance ;
   }
   else
@@ -743,10 +742,10 @@
   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
   rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
-  rMaxOAv = rMaxAv + kRadTolerance*0.5 ;
+  rMaxOAv = rMaxAv + halfRadTolerance ;
    
   // Intersection with z-surfaces
 
-  tolIDz = fDz - kCarTolerance*0.5 ;
-  tolODz = fDz + kCarTolerance*0.5 ;
+  tolIDz = fDz - halfCarTolerance ;
+  tolODz = fDz + halfCarTolerance ;
 
   if (std::fabs(p.z()) >= tolIDz)
@@ -756,9 +755,9 @@
       s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;  // Z intersect distance
 
-      if( s < 0.0 ) s = 0.0 ;                  // negative dist -> zero
+      if( s < 0.0 )  { s = 0.0; }                      // negative dist -> zero
 
       xi   = p.x() + s*v.x() ;  // Intersection coords
       yi   = p.y() + s*v.y() ;
-      rhoi2 = xi*xi + yi*yi   ;
+      rhoi2 = xi*xi + yi*yi  ;
 
       // Check validity of intersection
@@ -767,17 +766,17 @@
       if (v.z() > 0)
       {
-        tolORMin  = fRmin1 - 0.5*kRadTolerance*secRMin ;
-        tolIRMin  = fRmin1 + 0.5*kRadTolerance*secRMin ;
-        tolIRMax  = fRmax1 - 0.5*kRadTolerance*secRMin ;
-        tolORMax2 = (fRmax1 + 0.5*kRadTolerance*secRMax)*
-                    (fRmax1 + 0.5*kRadTolerance*secRMax) ;
+        tolORMin  = fRmin1 - halfRadTolerance*secRMin ;
+        tolIRMin  = fRmin1 + halfRadTolerance*secRMin ;
+        tolIRMax  = fRmax1 - halfRadTolerance*secRMin ;
+        tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
+                    (fRmax1 + halfRadTolerance*secRMax) ;
       }
       else
       {
-        tolORMin  = fRmin2 - 0.5*kRadTolerance*secRMin ;
-        tolIRMin  = fRmin2 + 0.5*kRadTolerance*secRMin ;
-        tolIRMax  = fRmax2 - 0.5*kRadTolerance*secRMin ;
-        tolORMax2 = (fRmax2 + 0.5*kRadTolerance*secRMax)*
-                    (fRmax2 + 0.5*kRadTolerance*secRMax) ;
+        tolORMin  = fRmin2 - halfRadTolerance*secRMin ;
+        tolIRMin  = fRmin2 + halfRadTolerance*secRMin ;
+        tolIRMax  = fRmax2 - halfRadTolerance*secRMin ;
+        tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
+                    (fRmax2 + halfRadTolerance*secRMax) ;
       }
       if ( tolORMin > 0 ) 
@@ -791,33 +790,22 @@
         tolIRMin2 = 0.0 ;
       }
-      if ( tolIRMax > 0 )   tolIRMax2 = tolIRMax*tolIRMax ;      
-      else                  tolIRMax2 = 0.0 ;
+      if ( tolIRMax > 0 )  { tolIRMax2 = tolIRMax*tolIRMax; }     
+      else                 { tolIRMax2 = 0.0; }
       
-      if (tolIRMin2 <= rhoi2 && rhoi2 <= tolIRMax2)
-      {
-  if ( seg && rhoi2 )
-  {
-    // Psi = angle made with central (average) phi of shape
-
-    cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-
-    if (cosPsi >= cosHDPhiIT) return s ;
-  }
-  else return s ;
-      }
-      /*
-      else if (tolORMin2 <= rhoi2 && rhoi2 <= tolORMax2)
-      {
-  if ( seg && rhoi2 )
-  {
-    // Psi = angle made with central (average) phi of shape
-
-    cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
-
-    if (cosPsi >= cosHDPhiIT) return s ;
-  }
-  else return s ;
-      }
-      */
+      if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
+      {
+        if ( !fPhiFullCone && rhoi2 )
+        {
+          // Psi = angle made with central (average) phi of shape
+
+          cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
+
+          if (cosPsi >= cosHDPhiIT)  { return s; }
+        }
+        else
+        {
+          return s;
+        }
+      }
     }
     else  // On/outside extent, and heading away  -> cannot intersect
@@ -861,17 +849,16 @@
   if (std::fabs(nt1) > kRadTolerance)  // Equation quadratic => 2 roots
   {
-      b = nt2/nt1 ;
-      c = nt3/nt1 ;
-      d = b*b-c   ;
-    if ( nt3 > rout*kRadTolerance*secRMax || rout < 0 )
+    b = nt2/nt1;
+    c = nt3/nt1;
+    d = b*b-c  ;
+    if ( (nt3 > rout*kRadTolerance*secRMax) || (rout < 0) )
     {
       // If outside real cone (should be rho-rout>kRadTolerance*0.5
       // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
 
-
       if (d >= 0)
       {
           
-        if (rout < 0 && nt3 <= 0 )
+        if ((rout < 0) && (nt3 <= 0))
         {
           // Inside `shadow cone' with -ve radius
@@ -882,5 +869,5 @@
         else
         {
-          if ( b <= 0 && c >= 0 ) // both >=0, try smaller root
+          if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
           {
             s = -b - std::sqrt(d) ;
@@ -906,5 +893,5 @@
             // Z ok. Check phi intersection if reqd
 
-            if ( ! seg )  return s ;
+            if ( fPhiFullCone )  { return s; }
             else
             {
@@ -914,5 +901,5 @@
               cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-              if ( cosPsi >= cosHDPhiIT ) return s ;
+              if ( cosPsi >= cosHDPhiIT )  { return s; }
             }
           }
@@ -925,22 +912,17 @@
       // check not inside, and heading through G4Cons (-> 0 to in)
 
-      if ( t3  > (rin + kRadTolerance*0.5*secRMin)*
-                 (rin + kRadTolerance*0.5*secRMin) && 
-           nt2 < 0                                 && 
-           d >= 0                                  && 
-           //  nt2 < -kCarTolerance*secRMax/2/fDz                  && 
-           // t2 < std::sqrt(t3)*v.z()*tanRMax                          && 
-           // d > kCarTolerance*secRMax*(rout-b*tanRMax*v.z())/nt1 && 
-           std::fabs(p.z()) <= tolIDz )
+      if ( ( t3  > (rin + halfRadTolerance*secRMin)*
+                   (rin + halfRadTolerance*secRMin) )
+        && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
       {
         // Inside cones, delta r -ve, inside z extent
 
-        if (seg)
+        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)  { return 0.0; }
+        }
+        else  { return 0.0; }
       }
     }
@@ -952,14 +934,14 @@
       s = -0.5*nt3/nt2 ;
 
-      if ( s < 0 )   return kInfinity ;   // travel away
+      if ( s < 0 )  { return kInfinity; }   // travel away
       else  // s >= 0,  If 'forwards'. Check z intersection
       {
         zi = p.z() + s*v.z() ;
 
-        if (std::fabs(zi) <= tolODz && nt2 < 0)
+        if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
         {
           // Z ok. Check phi intersection if reqd
 
-          if ( ! seg ) return s ;
+          if ( fPhiFullCone )  { return s; }
           else
           {
@@ -969,5 +951,5 @@
             cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-            if (cosPsi >= cosHDPhiIT) return s ;
+            if (cosPsi >= cosHDPhiIT)  { return s; }
           }
         }
@@ -1015,5 +997,5 @@
             if ( std::fabs(zi) <= tolODz )
             {
-              if ( seg )
+              if ( !fPhiFullCone )
               {
                 xi     = p.x() + s*v.x() ;
@@ -1022,7 +1004,7 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiIT) snxt = s ; 
+                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
               }
-              else return s ;
+              else  { return s; }
             }
           }
@@ -1046,9 +1028,9 @@
           ri = rMinAv + zi*tanRMin ;
 
-          if ( ri >= 0 )
+          if ( ri > 0 )
           {
-            if ( s >= 0 && std::fabs(zi) <= tolODz )  // s > 0
+            if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s > 0
             {
-              if ( seg )
+              if ( !fPhiFullCone )
               {
                 xi     = p.x() + s*v.x() ;
@@ -1056,7 +1038,7 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiOT) snxt = s ; 
+                if (cosPsi >= cosHDPhiOT)  { snxt = s; }
               }
-              else  return s ;
+              else  { return s; }
             }
           }
@@ -1067,7 +1049,7 @@
             ri = rMinAv + zi*tanRMin ;
 
-            if ( s >= 0 && ri >= 0 && std::fabs(zi) <= tolODz ) // s>0
+            if ( (s >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // s>0
             {
-              if ( seg )
+              if ( !fPhiFullCone )
               {
                 xi     = p.x() + s*v.x() ;
@@ -1075,7 +1057,7 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiIT) snxt = s ;
+                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
               }
-              else  return s ;
+              else  { return s; }
             }
           }
@@ -1095,11 +1077,11 @@
             // Inside inner real cone, heading outwards, inside z range
 
-            if ( seg )
+            if ( !fPhiFullCone )
             {
               cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
 
-              if (cosPsi >= cosHDPhiIT)  return 0.0 ;
+              if (cosPsi >= cosHDPhiIT)  { return 0.0; }
             }
-            else  return 0.0 ;
+            else  { return 0.0; }
           }
           else
@@ -1123,7 +1105,7 @@
                 zi = p.z() + s*v.z() ;
 
-                if ( s >= 0 && std::fabs(zi) <= tolODz )  // s>0
+                if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
                 {
-                  if ( seg )
+                  if ( !fPhiFullCone )
                   {
                     xi     = p.x() + s*v.x() ;
@@ -1132,10 +1114,10 @@
                     cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                    if ( cosPsi >= cosHDPhiIT )  snxt = s ; 
+                    if ( cosPsi >= cosHDPhiIT )  { snxt = s; }
                   }
-                  else return s ;
+                  else  { return s; }
                 }
               }
-              else  return kInfinity ;
+              else  { return kInfinity; }
             }
           }
@@ -1152,7 +1134,7 @@
             zi = p.z() + s*v.z() ;
 
-            if ( s >= 0 && std::fabs(zi) <= tolODz )  // s>0
+            if ( (s >= 0) && (std::fabs(zi) <= tolODz) )  // s>0
             {
-              if ( seg )
+              if ( !fPhiFullCone )
               {
                 xi     = p.x() + s*v.x();
@@ -1161,7 +1143,7 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
 
-                if (cosPsi >= cosHDPhiIT) snxt = s ; 
+                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
               }
-              else  return s ;
+              else  { return s; }
             }
           }
@@ -1180,10 +1162,8 @@
   //         -> Should use some form of loop Construct
 
-  if ( seg )
-  {
-    // First phi surface (`S'tarting phi)
-
-    sinSPhi = std::sin(fSPhi) ;
-    cosSPhi = std::cos(fSPhi) ;
+  if ( !fPhiFullCone )
+  {
+    // First phi surface (starting phi)
+
     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
                     
@@ -1192,5 +1172,5 @@
       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
 
-      if (Dist < kCarTolerance*0.5)
+      if (Dist < halfCarTolerance)
       {
         s = Dist/Comp ;
@@ -1198,5 +1178,5 @@
         if ( s < snxt )
         {
-          if ( s < 0 ) s = 0.0 ;
+          if ( s < 0 )  { s = 0.0; }
 
           zi = p.z() + s*v.z() ;
@@ -1210,20 +1190,18 @@
             tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
 
-            if ( rhoi2 >= tolORMin2 && rhoi2 <= tolORMax2 )
+            if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
             {
               // z and r intersections good - check intersecting with
               // correct half-plane
 
-              if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) snxt = s ;
-            }    
+              if ((yi*cosCPhi - xi*sinCPhi) <= 0 )  { snxt = s; }
+            }
           }
         }
       }
-    }    
-    // Second phi surface (`E'nding phi)
-
-    ePhi    = fSPhi + fDPhi ;
-    sinEPhi = std::sin(ePhi) ;
-    cosEPhi = std::cos(ePhi) ;
+    }
+
+    // Second phi surface (Ending phi)
+
     Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
         
@@ -1231,5 +1209,5 @@
     {
       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
-      if (Dist < kCarTolerance*0.5)
+      if (Dist < halfCarTolerance)
       {
         s = Dist/Comp ;
@@ -1237,5 +1215,5 @@
         if ( s < snxt )
         {
-          if ( s < 0 )  s = 0.0 ;
+          if ( s < 0 )  { s = 0.0; }
 
           zi = p.z() + s*v.z() ;
@@ -1249,11 +1227,11 @@
             tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
 
-            if ( rhoi2 >= tolORMin2 && rhoi2 <= tolORMax2 )
+            if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
             {
               // z and r intersections good - check intersecting with
               // correct half-plane
 
-              if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  snxt = s ;
-            }    
+              if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  { snxt = s; }
+            }
           }
         }
@@ -1261,13 +1239,5 @@
     }
   }
-  if (snxt < kCarTolerance*0.5) snxt = 0.;
-
-#ifdef consdebug
-  G4cout.precision(24);
-  G4cout<<"G4Cons::DistanceToIn(p,v) "<<G4endl;
-  G4cout<<"position = "<<p<<G4endl;
-  G4cout<<"direction = "<<v<<G4endl;
-  G4cout<<"distance = "<<snxt<<G4endl;
-#endif
+  if (snxt < halfCarTolerance)  { snxt = 0.; }
 
   return snxt ;
@@ -1283,9 +1253,7 @@
 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
 {
-  G4double safe=0.0, rho, safeR1, safeR2, safeZ ;
+  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
   G4double tanRMin, secRMin, pRMin ;
   G4double tanRMax, secRMax, pRMax ;
-  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi ;
-  G4double cosPsi ;
 
   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
@@ -1304,6 +1272,6 @@
     safeR2  = (rho - pRMax)/secRMax ;
 
-    if ( safeR1 > safeR2) safe = safeR1 ;
-    else                  safe = safeR2 ;
+    if ( safeR1 > safeR2) { safe = safeR1; }
+    else                  { safe = safeR2; }
   }
   else
@@ -1314,31 +1282,26 @@
     safe    = (rho - pRMax)/secRMax ;
   }
-  if ( safeZ > safe ) safe = safeZ ;
-
-  if ( fDPhi < twopi && rho )
-  {
-    phiC    = fSPhi + fDPhi*0.5 ;
-    cosPhiC = std::cos(phiC) ;
-    sinPhiC = std::sin(phiC) ;
-
+  if ( safeZ > safe )  { safe = safeZ; }
+
+  if ( !fPhiFullCone && rho )
+  {
     // Psi=angle from central phi to point
 
-    cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
+    cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
 
     if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
     {
-      if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0.0 )
-      {
-        safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
+      if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
+      {
+        safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
       }
       else
       {
-        ePhi    = fSPhi + fDPhi ;
-        safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
-      }
-      if ( safePhi > safe ) safe = safePhi ;
-    }
-  }
-  if ( safe < 0.0 ) safe = 0.0 ;
+        safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
+      }
+      if ( safePhi > safe )  { safe = safePhi; }
+    }
+  }
+  if ( safe < 0.0 )  { safe = 0.0; }
 
   return safe ;
@@ -1347,15 +1310,19 @@
 ///////////////////////////////////////////////////////////////
 //
-// 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
 
 G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
-              const G4ThreeVector& v,
-              const G4bool calcNorm,
-                    G4bool *validNorm,
-                    G4ThreeVector *n) const
+                                const G4ThreeVector& v,
+                                const G4bool calcNorm,
+                                      G4bool *validNorm,
+                                      G4ThreeVector *n) const
 {
   ESide side = kNull, sider = kNull, sidephi = kNull;
 
+  static const G4double halfCarTolerance=kCarTolerance*0.5;
+  static const G4double halfRadTolerance=kRadTolerance*0.5;
+  static const G4double halfAngTolerance=kAngTolerance*0.5;
+
   G4double snxt,sr,sphi,pdist ;
 
@@ -1373,6 +1340,4 @@
   // Vars for phi intersection:
 
-  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
-  G4double cPhi, sinCPhi, cosCPhi ;
   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
   G4double zi, ri, deltaRoi2 ;
@@ -1384,5 +1349,5 @@
     pdist = fDz - p.z() ;
 
-    if (pdist > kCarTolerance*0.5)
+    if (pdist > halfCarTolerance)
     {
       snxt = pdist/v.z() ;
@@ -1396,5 +1361,5 @@
         *validNorm = true ;
       }
-      return snxt = 0.0 ;
+      return  snxt = 0.0;
     }
   }
@@ -1403,5 +1368,5 @@
     pdist = fDz + p.z() ;
 
-    if ( pdist > kCarTolerance*0.5)
+    if ( pdist > halfCarTolerance)
     {
       snxt = -pdist/v.z() ;
@@ -1465,7 +1430,10 @@
                 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
   }
-  else deltaRoi2 = 1.0 ;
-
-  if ( nt1 && deltaRoi2 > 0.0 )  
+  else
+  {
+    deltaRoi2 = 1.0;
+  }
+
+  if ( nt1 && (deltaRoi2 > 0.0) )  
   {
     // Equation quadratic => 2 roots : second root must be leaving
@@ -1478,7 +1446,7 @@
     {
       // Check if on outer cone & heading outwards
-      // NOTE: Should use rho-rout>-kRadtolerance*0.5
+      // NOTE: Should use rho-rout>-kRadTolerance*0.5
         
-      if (nt3 > -kRadTolerance*0.5 && nt2 >= 0 )
+      if (nt3 > -halfRadTolerance && nt2 >= 0 )
       {
         if (calcNorm)
@@ -1486,5 +1454,5 @@
           risec      = std::sqrt(t3)*secRMax ;
           *validNorm = true ;
-          *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
+          *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
         }
         return snxt=0 ;
@@ -1497,6 +1465,5 @@
         ri    = tanRMax*zi + rMaxAv ;
           
-        if ( (ri >= 0) && (-kRadTolerance*0.5 <= sr) && 
-                          ( sr <= kRadTolerance*0.5)     )
+        if ((ri >= 0) && (-halfRadTolerance <= sr) && (sr <= halfRadTolerance))
         {
           // An intersection within the tolerance
@@ -1506,5 +1473,5 @@
           sidetol = kRMax ;
         }            
-        if ( (ri < 0) || (sr < kRadTolerance*0.5) )
+        if ( (ri < 0) || (sr < halfRadTolerance) )
         {
           // Safety: if both roots -ve ensure that sr cannot `win'
@@ -1515,11 +1482,13 @@
           ri  = tanRMax*zi + rMaxAv ;
 
-          if (ri >= 0 && sr2 > kRadTolerance*0.5)  sr = sr2 ;
+          if ((ri >= 0) && (sr2 > halfRadTolerance))
+          {
+            sr = sr2;
+          }
           else
           {
             sr = kInfinity ;
 
-            if(    (-kRadTolerance*0.5 <= sr2)
-                && ( sr2 <= kRadTolerance*0.5)  )
+            if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
             {
               // An intersection within the tolerance.
@@ -1540,12 +1509,12 @@
       if ( calcNorm )
       {
-        risec      = std::sqrt(t3)*secRMax ;
-        *validNorm = true ;
-        *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
+        risec      = std::sqrt(t3)*secRMax;
+        *validNorm = true;
+        *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
       }
       return snxt = 0.0 ;
     }
   }
-  else if ( nt2 && deltaRoi2 > 0.0 )
+  else if ( nt2 && (deltaRoi2 > 0.0) )
   {
     // Linear case (only one intersection) => point outside outer cone
@@ -1553,7 +1522,7 @@
     if ( calcNorm )
     {
-      risec      = std::sqrt(t3)*secRMax ;
-      *validNorm = true ;
-      *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) ;
+      risec      = std::sqrt(t3)*secRMax;
+      *validNorm = true;
+      *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
     }
     return snxt = 0.0 ;
@@ -1569,5 +1538,5 @@
   // Check possible intersection within tolerance
 
-  if ( slentol <= kCarTolerance*0.5 )
+  if ( slentol <= halfCarTolerance )
   {
     // An intersection within the tolerance was found.  
@@ -1580,8 +1549,8 @@
     // Calculate a normal vector,  as below
 
-    xi    = p.x() + slentol*v.x() ;
-    yi    = p.y() + slentol*v.y() ;
-    risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
-    G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
+    xi    = p.x() + slentol*v.x();
+    yi    = p.y() + slentol*v.y();
+    risec = std::sqrt(xi*xi + yi*yi)*secRMax;
+    G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
 
     if ( Normal.dot(v) > 0 )    // We will leave the Cone immediatelly
@@ -1595,5 +1564,5 @@
     }
     else // On the surface, but not heading out so we ignore this intersection
-    {    //                                    (as it is within tolerance).  
+    {    //                                        (as it is within tolerance).
       slentol = kInfinity ;
     }
@@ -1621,5 +1590,5 @@
       d = b*b - c ;
 
-      if (d >= 0.0 )
+      if ( d >= 0.0 )
       {
         // NOTE: should be rho-rin<kRadTolerance*0.5,
@@ -1630,6 +1599,6 @@
           if ( nt2 < 0.0 )
           {
-            if (calcNorm)  *validNorm = false ;
-            return          snxt      = 0.0 ;
+            if (calcNorm)  { *validNorm = false; }
+            return          snxt      = 0.0;
           }
         }
@@ -1640,6 +1609,5 @@
           ri  = tanRMin*zi + rMinAv ;
 
-          if( (ri >= 0.0) && (-kRadTolerance*0.5 <= sr2) && 
-                             ( sr2 <= kRadTolerance*0.5)      )
+          if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
           {
             // An intersection within the tolerance
@@ -1649,5 +1617,5 @@
             sidetol = kRMax ;
           }
-          if( (ri<0) || (sr2 < kRadTolerance*0.5) )
+          if( (ri<0) || (sr2 < halfRadTolerance) )
           {
             sr3 = -b + std::sqrt(d) ;
@@ -1656,5 +1624,5 @@
             //         distancetoout
 
-            if  ( sr3 > kCarTolerance*0.5 )
+            if  ( sr3 > halfRadTolerance )
             {
               if( sr3 < sr )
@@ -1670,5 +1638,5 @@
               } 
             }
-            else if ( sr3 > -kCarTolerance*0.5 )
+            else if ( sr3 > -halfRadTolerance )
             {
               // Intersection in tolerance. Store to check if it's good
@@ -1678,10 +1646,10 @@
             }
           }
-          else if ( sr2 < sr && sr2 > kCarTolerance*0.5 )
+          else if ( (sr2 < sr) && (sr2 > halfCarTolerance) )
           {
             sr    = sr2 ;
             sider = kRMin ;
           }
-          else if (sr2 > -kCarTolerance*0.5)
+          else if (sr2 > -halfCarTolerance)
           {
             // Intersection in tolerance. Store to check if it's good
@@ -1690,5 +1658,5 @@
             sidetol = kRMin ;
           }    
-          if( slentol <= kCarTolerance*0.5  )
+          if( slentol <= halfCarTolerance  )
           {
             // An intersection within the tolerance was found. 
@@ -1706,5 +1674,6 @@
             if( Normal.dot(v) > 0 )
             {
-              // We will leave the Cone immediatelly
+              // We will leave the cone immediately
+
               if( calcNorm ) 
               {
@@ -1731,20 +1700,14 @@
   // 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 ( !fPhiFullCone )
+  {
     // add angle calculation with correction 
-      // of the  difference in domain of atan2 and Sphi
-        vphi = std::atan2(v.y(),v.x()) ;
-
-         if ( vphi < fSPhi - kAngTolerance*0.5  )             vphi += twopi ; 
-         else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5 ) vphi -= twopi; 
+    // of the difference in domain of atan2 and Sphi
+
+    vphi = std::atan2(v.y(),v.x()) ;
+
+    if ( vphi < fSPhi - halfAngTolerance  )              { vphi += twopi; }
+    else if ( vphi > fSPhi + fDPhi + halfAngTolerance )  { vphi -= twopi; }
+
     if ( p.x() || p.y() )   // Check if on z axis (rho not needed later)
     {
@@ -1761,45 +1724,44 @@
       sidephi = kNull ;
 
-      if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)
-                              && (pDistE <= 0.5*kCarTolerance) ) )
-         || ( (fDPhi >  pi) && !((pDistS >  0.5*kCarTolerance)
-                              && (pDistE >  0.5*kCarTolerance) ) )  )
-        {
-          // Inside both phi *full* planes
-          if ( compS < 0 )
+      if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
+                            && (pDistE <= halfCarTolerance) ) )
+         || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
+                              && (pDistE >  halfCarTolerance) ) )  )
+      {
+        // Inside both phi *full* planes
+        if ( compS < 0 )
+        {
+          sphi = pDistS/compS ;
+          if (sphi >= -halfCarTolerance)
           {
-            sphi = pDistS/compS ;
-            if (sphi >= -0.5*kCarTolerance)
+            xi = p.x() + sphi*v.x() ;
+            yi = p.y() + sphi*v.y() ;
+
+            // Check intersecting with correct half-plane
+            // (if not -> no intersect)
+            //
+            if ( (std::abs(xi)<=kCarTolerance)
+              && (std::abs(yi)<=kCarTolerance) )
             {
-              xi = p.x() + sphi*v.x() ;
-              yi = p.y() + sphi*v.y() ;
-
-              // Check intersecting with correct half-plane
-              // (if not -> no intersect)
-              //
-               if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)){
-		 sidephi= kSPhi;
-		if(((fSPhi-0.5*kAngTolerance)<=vphi)&&((ePhi+0.5*kAngTolerance)>=vphi))
-                    { sphi = kInfinity; }
-				       
+              sidephi= kSPhi;
+              if ( ( fSPhi-halfAngTolerance <= vphi )
+                && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
+              {
+                sphi = kInfinity;
               }
-              else
-              if ((yi*cosCPhi-xi*sinCPhi)>=0)
-              {
-                sphi = kInfinity ;
-              }
-              else
-              {
-                sidephi = kSPhi ;
-                if ( pDistS > -kCarTolerance*0.5 )
-                {
-                  sphi = 0.0 ; // Leave by sphi immediately
-                }    
-              }       
+            }
+            else
+            if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
+            {
+              sphi = kInfinity ;
             }
             else
             {
-              sphi = kInfinity ;
-            }
+              sidephi = kSPhi ;
+              if ( pDistS > -halfCarTolerance )
+              {
+                sphi = 0.0 ; // Leave by sphi immediately
+              }    
+            }       
           }
           else
@@ -1807,43 +1769,52 @@
             sphi = kInfinity ;
           }
-
-          if ( compE < 0 )
+        }
+        else
+        {
+          sphi = kInfinity ;
+        }
+
+        if ( compE < 0 )
+        {
+          sphi2 = pDistE/compE ;
+
+          // Only check further if < starting phi intersection
+          //
+          if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
           {
-            sphi2 = pDistE/compE ;
-
-            // Only check further if < starting phi intersection
-            //
-            if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) )
+            xi = p.x() + sphi2*v.x() ;
+            yi = p.y() + sphi2*v.y() ;
+
+            // Check intersecting with correct half-plane
+
+            if ( (std::abs(xi)<=kCarTolerance)
+              && (std::abs(yi)<=kCarTolerance) )
             {
-              xi = p.x() + sphi2*v.x() ;
-              yi = p.y() + sphi2*v.y() ;
-
-              // Check intersecting with correct half-plane 
-               if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)){
-               // Leaving via ending phi
-                if(!(((fSPhi-0.5*kAngTolerance)<=vphi)&&((ePhi+0.5*kAngTolerance)>=vphi))){
-                  sidephi = kEPhi ;
-                  if ( pDistE <= -kCarTolerance*0.5 ) sphi = sphi2 ;
-                  else                                sphi = 0.0 ;
-                  }
-                }
-	     else // Check intersecting with correct half-plane
-              if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
+              // Leaving via ending phi
+
+              if(!( (fSPhi-halfAngTolerance <= vphi)
+                 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
               {
-                // Leaving via ending phi
-
                 sidephi = kEPhi ;
-                if ( pDistE <= -kCarTolerance*0.5 ) sphi = sphi2 ;
-                else                                sphi = 0.0 ;
+                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
+                else                                { sphi = 0.0; }
               }
             }
+            else // Check intersecting with correct half-plane
+            if ( yi*cosCPhi-xi*sinCPhi >= 0 )
+            {
+              // Leaving via ending phi
+
+              sidephi = kEPhi ;
+              if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
+              else                                { sphi = 0.0; }
+            }
           }
         }
-        else
-        {
-          sphi = kInfinity ;
-        }
-
-
+      }
+      else
+      {
+        sphi = kInfinity ;
+      }
     }
     else
@@ -1852,12 +1823,9 @@
       // within phi of shape, Step limited by rmax, else Step =0
 
-      // vphi = std::atan2(v.y(),v.x()) ;
-
-      // if ( fSPhi < vphi && vphi < fSPhi + fDPhi ) sphi = kInfinity ;
-
-       if ( ((fSPhi-0.5*kAngTolerance) <= vphi) && (vphi <=( fSPhi + fDPhi)+0.5*kAngTolerance) )
-         {
-          sphi = kInfinity ;
-        }
+      if ( (fSPhi-halfAngTolerance <= vphi)
+        && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
+      {
+        sphi = kInfinity ;
+      }
       else
       {
@@ -1868,6 +1836,6 @@
     if ( sphi < snxt )  // Order intersecttions
     {
-        snxt=sphi ;
-        side=sidephi ;
+      snxt=sphi ;
+      side=sidephi ;
     }
   }
@@ -1880,8 +1848,6 @@
   {
     switch(side)
-    {
-      case kRMax:
-          // Note: returned vector not normalised
-          // (divide by frmax for unit vector)
+    {                     // Note: returned vector not normalised
+      case kRMax:         // (divide by frmax for unit vector)
         xi         = p.x() + snxt*v.x() ;
         yi         = p.y() + snxt*v.y() ;
@@ -1891,21 +1857,27 @@
         break ;
       case kRMin:
-        *validNorm=false ;  // Rmin is inconvex
+        *validNorm = false ;  // Rmin is inconvex
         break ;
       case kSPhi:
         if ( fDPhi <= pi )
         {
-          *n         = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
+          *n         = G4ThreeVector(sinSPhi, -cosSPhi, 0);
           *validNorm = true ;
         }
-        else   *validNorm = false ;
+        else
+        {
+          *validNorm = false ;
+        }
         break ;
       case kEPhi:
         if ( fDPhi <= pi )
         {
-          *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 ;
       case kPZ:
@@ -1925,10 +1897,10 @@
         G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
         G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
-        G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm << " mm"
-               << G4endl << G4endl ;
+        G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
+               << " mm" << G4endl << G4endl ;
         if( p.x() != 0. || p.x() != 0.)
         {
-           G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree << " degree" 
-                  << G4endl << G4endl ; 
+           G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree
+                  << " degree" << G4endl << G4endl ; 
         }
         G4cout << "Direction:" << G4endl << G4endl ;
@@ -1943,12 +1915,6 @@
     }
   }
-  if (snxt < kCarTolerance*0.5) snxt = 0.;
-#ifdef consdebug
-  G4cout.precision(24);
-  G4cout<<"G4Cons::DistanceToOut(p,v,...) "<<G4endl;
-  G4cout<<"position = "<<p<<G4endl;
-  G4cout<<"direction = "<<v<<G4endl;
-  G4cout<<"distance = "<<snxt<<G4endl;
-#endif
+  if (snxt < halfCarTolerance)  { snxt = 0.; }
+
   return snxt ;
 }
@@ -1960,8 +1926,7 @@
 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
 {
-  G4double safe=0.0,rho,safeR1,safeR2,safeZ ;
-  G4double tanRMin,secRMin,pRMin ;
-  G4double tanRMax,secRMax,pRMax ;
-  G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi ;
+  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
+  G4double tanRMin, secRMin, pRMin;
+  G4double tanRMax, secRMax, pRMax;
 
 #ifdef G4CSGDEBUG
@@ -1975,13 +1940,13 @@
     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
-    G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm << " mm" 
-           << G4endl << G4endl ;
-    if( p.x() != 0. || p.x() != 0.)
-    {
-      G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree << " degree" 
-             << G4endl << G4endl ; 
-    }
-    G4Exception("G4Cons::DistanceToOut(p)", "Notification", JustWarning, 
-                "Point p is outside !?" );
+    G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
+           << " mm" << G4endl << G4endl ;
+    if( (p.x() != 0.) || (p.x() != 0.) )
+    {
+      G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree
+             << " degree" << G4endl << G4endl ; 
+    }
+    G4Exception("G4Cons::DistanceToOut(p)", "Notification",
+                JustWarning, "Point p is outside !?" );
   }
 #endif
@@ -1997,5 +1962,8 @@
     safeR1  = (rho - pRMin)/secRMin ;
   }
-  else safeR1 = kInfinity ;
+  else
+  {
+    safeR1 = kInfinity ;
+  }
 
   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
@@ -2004,31 +1972,27 @@
   safeR2  = (pRMax - rho)/secRMax ;
 
-  if (safeR1 < safeR2) safe = safeR1 ;
-  else                 safe = safeR2 ;
-  if (safeZ < safe)    safe = safeZ  ;
+  if (safeR1 < safeR2)  { safe = safeR1; }
+  else                  { safe = safeR2; }
+  if (safeZ < safe)     { safe = safeZ ; }
 
   // Check if phi divided, Calc distances closest phi plane
 
-  if (fDPhi < twopi)
+  if (!fPhiFullCone)
   {
     // Above/below central phi of G4Cons?
 
-    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 ( (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 ;
-  }
-  if ( safe < 0 ) safe = 0 ;
-  return safe ;  
+      safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
+    }
+    if (safePhi < safe)  { safe = safePhi; }
+  }
+  if ( safe < 0 )  { safe = 0; }
+
+  return safe ;
 }
 
@@ -2068,6 +2032,4 @@
   meshAngle = fDPhi/(noCrossSections - 1) ;
 
-  // G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1  ;
-
   meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
   meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;
@@ -2076,5 +2038,5 @@
   // on the x axis. Will give better extent calculations when not rotated.
 
-  if (fDPhi == twopi && fSPhi == 0.0 )
+  if ( fPhiFullCone && (fSPhi == 0.0) )
   {
     sAngle = -meshAngle*0.5 ;
@@ -2102,6 +2064,4 @@
       rMaxY2 = meshRMax2*sinCrossAngle ;
         
-      // G4double RMin = (fRmin2 <= fRmin1) ? fRmin2 : fRmin1  ;
-
       rMinX1 = fRmin1*cosCrossAngle ;
       rMinY1 = fRmin1*sinCrossAngle ;
@@ -2127,4 +2087,5 @@
                 "Error in allocation of vertices. Out of memory !");
   }
+
   return vertices ;
 }
@@ -2192,5 +2153,5 @@
   rRand2 = RandFlat::shoot(fRmin2,fRmax2);
   
-  if(fSPhi == 0. && fDPhi == twopi){ Afive = 0.; }
+  if ( (fSPhi == 0.) && fPhiFullCone )  { Afive = 0.; }
   chose  = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
  
@@ -2225,5 +2186,5 @@
   else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
   {
-    return G4ThreeVector (rRand1*cosu,rRand1*sinu,-1*fDz);
+    return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
   }
   else if( (chose >= Aone + Atwo + Athree)
