Index: trunk/source/geometry/solids/CSG/src/G4Tubs.cc
===================================================================
--- trunk/source/geometry/solids/CSG/src/G4Tubs.cc	(revision 850)
+++ trunk/source/geometry/solids/CSG/src/G4Tubs.cc	(revision 921)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Tubs.cc,v 1.68 2008/06/23 13:37:39 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Tubs.cc,v 1.74 2008/11/06 15:26:53 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 // 
@@ -108,5 +108,5 @@
                 "Invalid Z half-length");
   }
-  if ( pRMin < pRMax && pRMin >= 0 ) // Check radii
+  if ( (pRMin < pRMax) && (pRMin >= 0) ) // Check radii
   {
     fRMin = pRMin ; 
@@ -121,10 +121,14 @@
                 "Invalid radii.");
   }
-  if ( pDPhi >= twopi ) // Check angles
+
+  fPhiFullTube = true;
+  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
   {
     fDPhi=twopi;
+    fSPhi=0;
   }
   else
   {
+    fPhiFullTube = false;
     if ( pDPhi > 0 )
     {
@@ -136,25 +140,24 @@
              << "        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
-
-  fSPhi = pSPhi;
-
-  if ( fSPhi < 0 )
-  {
-    fSPhi = twopi - std::fmod(std::fabs(fSPhi),twopi) ;
-  }
-  else
-  {
-    fSPhi = std::fmod(fSPhi,twopi) ;
-  }
-  if (fSPhi + fDPhi > twopi )
-  {
-    fSPhi -= twopi ;
-  }
+      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();
 }
 
@@ -290,6 +293,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 ;
@@ -317,6 +320,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 ;
@@ -363,29 +366,14 @@
     noEntries = vertices->size() ;
     noBetweenSections4 = noEntries - 4 ;
-    /*
-    G4cout << "vertices = " << noEntries << "\t"
-           << "v-4 = " << noBetweenSections4 << G4endl;
-    G4cout << G4endl;
-    for(i = 0 ; i < noEntries ; i++ )
-    {
-      G4cout << i << "\t" << "v.x = " << ((*vertices)[i]).x() << "\t"
-                          << "v.y = " << ((*vertices)[i]).y() << "\t"
-                          << "v.z = " << ((*vertices)[i]).z() << "\t" << G4endl;
-    }      
-    G4cout << G4endl;
-    G4cout << "ClipCrossSection" << G4endl;
-    */
+    
     for (i = 0 ; i < noEntries ; i += 4 )
     {
-      // G4cout << "section = " << i << G4endl;
-      ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
-    }
-    // G4cout << "ClipBetweenSections" << G4endl;
+      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
+    }
     for (i = 0 ; i < noBetweenSections4 ; i += 4 )
     {
-      // G4cout << "between sections = " << i << G4endl;
-      ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
-    }
-    if (pMin != kInfinity || pMax != -kInfinity )
+      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
+    }
+    if ((pMin != kInfinity) || (pMax != -kInfinity) )
     {
       existsAfterClip = true ;
@@ -426,17 +414,20 @@
   G4double r2,pPhi,tolRMin,tolRMax;
   EInside in = kOutside ;
-  
-  if (std::fabs(p.z()) <= fDz - kCarTolerance*0.5)
+  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)
   {
     r2 = p.x()*p.x() + p.y()*p.y() ;
 
-    if (fRMin) { tolRMin = fRMin + kRadTolerance*0.5 ; }
+    if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
     else       { tolRMin = 0 ; }
 
-    tolRMax = fRMax - kRadTolerance*0.5 ;
+    tolRMax = fRMax - halfRadTolerance ;
       
-    if (r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
-    {
-      if ( fDPhi == twopi )
+    if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
+    {
+      if ( fPhiFullTube )
       {
         in = kInside ;
@@ -447,6 +438,5 @@
         // if not inside, try outer tolerant phi boundaries
 
-        pPhi = std::atan2(p.y(),p.x()) ;
-        if ((tolRMin==0)&&(p.x()==0)&&(p.y()==0))
+        if ((tolRMin==0)&&(p.x()<=halfCarTolerance)&&(p.y()<=halfCarTolerance))
         {
           in=kSurface;
@@ -454,20 +444,21 @@
         else
         {
-          if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; } // 0<=pPhi<2pi
+          pPhi = std::atan2(p.y(),p.x()) ;
+          if ( pPhi < -halfAngTolerance )  { pPhi += twopi; } // 0<=pPhi<2pi
 
           if ( fSPhi >= 0 )
           {
-            if ( (std::abs(pPhi) < kAngTolerance*0.5)
-              && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
+            if ( (std::abs(pPhi) < halfAngTolerance)
+              && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
             { 
               pPhi += twopi ; // 0 <= pPhi < 2pi
             }
-            if ( (pPhi >= fSPhi + kAngTolerance*0.5)
-              && (pPhi <= fSPhi + fDPhi - kAngTolerance*0.5) )
+            if ( (pPhi >= fSPhi + halfAngTolerance)
+              && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
             {
               in = kInside ;
             }
-            else if ( (pPhi >= fSPhi - kAngTolerance*0.5)
-                   && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
+            else if ( (pPhi >= fSPhi - halfAngTolerance)
+                   && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
             {
               in = kSurface ;
@@ -476,8 +467,8 @@
           else  // fSPhi < 0
           {
-            if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
-              && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
-            else if ( (pPhi <= fSPhi + twopi + kAngTolerance*0.5)
-                   && (pPhi >= fSPhi + fDPhi  - kAngTolerance*0.5) )
+            if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
+              && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;} //kOutside
+            else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
+                   && (pPhi >= fSPhi + fDPhi  - halfAngTolerance) )
             {
               in = kSurface ;
@@ -493,6 +484,6 @@
     else  // Try generous boundaries
     {
-      tolRMin = fRMin - kRadTolerance*0.5 ;
-      tolRMax = fRMax + kRadTolerance*0.5 ;
+      tolRMin = fRMin - halfRadTolerance ;
+      tolRMax = fRMax + halfRadTolerance ;
 
       if ( tolRMin < 0 )  { tolRMin = 0; }
@@ -500,6 +491,6 @@
       if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
       {
-        if ( fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis
-        {
+        if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
+        {                        // Continuous in phi or on z-axis
           in = kSurface ;
         }
@@ -508,14 +499,14 @@
           pPhi = std::atan2(p.y(),p.x()) ;
 
-          if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; } // 0<=pPhi<2pi
+          if ( pPhi < -halfAngTolerance)  { pPhi += twopi; } // 0<=pPhi<2pi
           if ( fSPhi >= 0 )
           {
-            if ( (std::abs(pPhi) < kAngTolerance*0.5)
-              && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
+            if ( (std::abs(pPhi) < halfAngTolerance)
+              && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
             { 
               pPhi += twopi ; // 0 <= pPhi < 2pi
             }
-            if ( (pPhi >= fSPhi - kAngTolerance*0.5)
-              && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
+            if ( (pPhi >= fSPhi - halfAngTolerance)
+              && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
             {
               in = kSurface ;
@@ -524,6 +515,6 @@
           else  // fSPhi < 0
           {
-            if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
-              && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
+            if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
+              && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
             else
             {
@@ -535,9 +526,9 @@
     }
   }
-  else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5)
+  else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
   {                                          // Check within tolerant r limits
     r2      = p.x()*p.x() + p.y()*p.y() ;
-    tolRMin = fRMin - kRadTolerance*0.5 ;
-    tolRMax = fRMax + kRadTolerance*0.5 ;
+    tolRMin = fRMin - halfRadTolerance ;
+    tolRMax = fRMax + halfRadTolerance ;
 
     if ( tolRMin < 0 )  { tolRMin = 0; }
@@ -545,6 +536,6 @@
     if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
     {
-      if (fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis
-      {
+      if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
+      {                        // Continuous in phi or on z-axis
         in = kSurface ;
       }
@@ -553,14 +544,14 @@
         pPhi = std::atan2(p.y(),p.x()) ;
 
-        if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; }  // 0<=pPhi<2pi
+        if ( pPhi < -halfAngTolerance )  { pPhi += twopi; }  // 0<=pPhi<2pi
         if ( fSPhi >= 0 )
         {
-          if ( (std::abs(pPhi) < kAngTolerance*0.5)
-            && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
+          if ( (std::abs(pPhi) < halfAngTolerance)
+            && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
           { 
             pPhi += twopi ; // 0 <= pPhi < 2pi
           }
-          if ( (pPhi >= fSPhi - kAngTolerance*0.5)
-            && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
+          if ( (pPhi >= fSPhi - halfAngTolerance)
+            && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
           {
             in = kSurface;
@@ -569,6 +560,6 @@
         else  // fSPhi < 0
         {
-          if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
-            && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
+          if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
+            && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;}
           else
           {
@@ -589,9 +580,13 @@
 
 G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const
-{ G4int noSurfaces = 0;
+{
+  G4int noSurfaces = 0;
   G4double rho, pPhi;
-  G4double delta   = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
   G4double distZ, distRMin, distRMax;
   G4double distSPhi = kInfinity, distEPhi = kInfinity;
+
+  static const G4double halfCarTolerance = 0.5*kCarTolerance;
+  static const G4double halfAngTolerance = 0.5*kAngTolerance;
+
   G4ThreeVector norm, sumnorm(0.,0.,0.);
   G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
@@ -604,14 +599,14 @@
   distZ    = std::fabs(std::fabs(p.z()) - fDz);
 
-  if (fDPhi < twopi)   //  &&  rho ) // Protected against (0,0,z) 
-  {
-    if ( rho )
+  if (!fPhiFullTube)    // Protected against (0,0,z) 
+  {
+    if ( rho > halfCarTolerance )
     {
       pPhi = std::atan2(p.y(),p.x());
     
-      if(pPhi  < fSPhi-delta)           { pPhi += twopi; }
-      else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
-
-      distSPhi = std::fabs( pPhi - fSPhi );       
+      if(pPhi  < fSPhi- halfCarTolerance)           { pPhi += twopi; }
+      else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
+
+      distSPhi = std::fabs(pPhi - fSPhi);       
       distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 
     }
@@ -624,12 +619,12 @@
     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
   }
-  if ( rho > delta )   { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
-
-  if( distRMax <= delta )
+  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
+
+  if( distRMax <= halfCarTolerance )
   {
     noSurfaces ++;
     sumnorm += nR;
   }
-  if( fRMin && distRMin <= delta )
+  if( fRMin && (distRMin <= halfCarTolerance) )
   {
     noSurfaces ++;
@@ -638,10 +633,10 @@
   if( fDPhi < twopi )   
   {
-    if (distSPhi <= dAngle)  // delta)
+    if (distSPhi <= halfAngTolerance)  
     {
       noSurfaces ++;
       sumnorm += nPs;
     }
-    if (distEPhi <= dAngle) // delta) 
+    if (distEPhi <= halfAngTolerance)  
     {
       noSurfaces ++;
@@ -649,5 +644,5 @@
     }
   }
-  if (distZ <= delta)  
+  if (distZ <= halfCarTolerance)  
   {
     noSurfaces ++;
@@ -668,4 +663,5 @@
   else if ( noSurfaces == 1 )  { norm = sumnorm; }
   else                         { norm = sumnorm.unit(); }
+
   return norm;
 }
@@ -714,5 +710,5 @@
     }
   }   
-  if (fDPhi < twopi  &&  rho ) // Protected against (0,0,z) 
+  if (!fPhiFullTube  &&  rho ) // Protected against (0,0,z) 
   {
     phi = std::atan2(p.y(),p.x()) ;
@@ -749,26 +745,26 @@
     case kNRMin : // Inner radius
     {                      
-      norm = G4ThreeVector(-p.x()/rho,-p.y()/rho,0) ;
+      norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
       break ;
     }
     case kNRMax : // Outer radius
     {                  
-      norm = G4ThreeVector(p.x()/rho,p.y()/rho,0) ;
+      norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
       break ;
     }
     case kNZ : //    + or - 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;
     }
@@ -804,62 +800,27 @@
 //
 // NOTE:
-// - Precalculations for phi trigonometry are Done `just in time'
-// - `if valid' implies tolerant checking of intersection points
+// - 'if valid' implies tolerant checking of intersection points
 
 G4double G4Tubs::DistanceToIn( const G4ThreeVector& p,
                                const G4ThreeVector& v  ) const
 {
-  G4double snxt = kInfinity ;  // snxt = default return value
-
-  // Precalculated trig for phi intersections - used by r,z intersections to
-  //                                            check validity
-
-  G4bool seg ;        // true if segmented
-
-  G4double hDPhi, hDPhiOT, hDPhiIT, cosHDPhiOT=0., cosHDPhiIT=0. ;
-          // half dphi + outer tolerance
-
-  G4double cPhi, sinCPhi=0., cosCPhi=0. ;  // central phi
-
-  G4double tolORMin2, tolIRMax2 ;  // `generous' radii squared
-
+  G4double snxt = kInfinity ;      // snxt = default return value
+  G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
   G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
+
+  static const G4double halfCarTolerance = 0.5*kCarTolerance;
+  static const G4double halfRadTolerance = 0.5*kRadTolerance;
 
   // Intersection point variables
   //
-  G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi ; 
-
-  G4double t1, t2, t3, b, c, d ;   // Quadratic solver variables 
-
-  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  ;
-  }
-
+  G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
+  G4double t1, t2, t3, b, c, d ;     // Quadratic solver variables 
+  
   // Calculate tolerant rmin and rmax
 
   if (fRMin > kRadTolerance)
   {
-    tolORMin2 = (fRMin - 0.5*kRadTolerance)*(fRMin - 0.5*kRadTolerance) ;
-    tolIRMin2 = (fRMin + 0.5*kRadTolerance)*(fRMin + 0.5*kRadTolerance) ;
+    tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
+    tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
   }
   else
@@ -868,11 +829,11 @@
     tolIRMin2 = 0.0 ;
   }
-  tolORMax2 = (fRMax + 0.5*kRadTolerance)*(fRMax + 0.5*kRadTolerance) ;
-  tolIRMax2 = (fRMax - 0.5*kRadTolerance)*(fRMax - 0.5*kRadTolerance) ;
+  tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
+  tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - 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)
@@ -880,5 +841,5 @@
     if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
     {
-      s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;     // Z intersect distance
+      s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;   // Z intersect distance
 
       if(s < 0.0)  { s = 0.0; }
@@ -890,7 +851,7 @@
       // Check validity of intersection
 
-      if (tolIRMin2 <= rho2 && rho2 <= tolIRMax2)
-      {
-        if (seg && rho2)
+      if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
+      {
+        if (!fPhiFullTube && rho2)
         {
           // Psi = angle made with central (average) phi of shape
@@ -899,5 +860,5 @@
           iden   = std::sqrt(rho2) ;
           cosPsi = inum/iden ;
-          if (cosPsi >= cosHDPhiIT) return s ;
+          if (cosPsi >= cosHDPhiIT)  { return s ; }
         }
         else
@@ -909,5 +870,5 @@
     else
     {
-      if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
+      if ( snxt<halfCarTolerance )  { snxt=0; }
       return snxt ;  // On/outside extent, and heading away
                      // -> cannot intersect
@@ -934,5 +895,5 @@
     b = t2/t1 ;
     c = t3 - fRMax*fRMax ;
-    if (t3 >= tolORMax2 && t2<0)   // This also handles the tangent case
+    if ((t3 >= tolORMax2) && (t2<0))   // This also handles the tangent case
     {
       // Try outer cylinder intersection
@@ -954,5 +915,5 @@
             // Z ok. Check phi intersection if reqd
             //
-            if (!seg)
+            if (fPhiFullTube)
             {
               return s ;
@@ -963,5 +924,5 @@
               yi     = p.y() + s*v.y() ;
               cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
-              if (cosPsi >= cosHDPhiIT) return s ;
+              if (cosPsi >= cosHDPhiIT)  { return s ; }
             }
           }  //  end if std::fabs(zi)
@@ -974,9 +935,9 @@
       // check not inside, and heading through tubs (-> 0 to in)
 
-      if (t3 > tolIRMin2 && t2 < 0 && std::fabs(p.z()) <= tolIDz)
+      if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
       {
         // Inside both radii, delta r -ve, inside z extent
 
-        if (seg)
+        if (!fPhiFullTube)
         {
           inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
@@ -1004,5 +965,5 @@
                 snxt = c/(-b+std::sqrt(d)); // using safe solution
                                             // for quadratic equation 
-                if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
+                if ( snxt < halfCarTolerance ) { snxt=0; }
                 return snxt ;
               }      
@@ -1035,5 +996,5 @@
               snxt= c/(-b+std::sqrt(d)); // using safe solution
                                          // for quadratic equation 
-              if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
+              if ( snxt < halfCarTolerance ) { snxt=0; }
               return snxt ;
             }      
@@ -1043,5 +1004,5 @@
             }
           }
-        } // end if   (seg)
+        } // end if   (!fPhiFullTube)
       }   // end if   (t3>tolIRMin2)
     }     // end if   (Inside Outer Radius) 
@@ -1056,5 +1017,5 @@
 
         s = -b + std::sqrt(d) ;
-        if (s >= -0.5*kCarTolerance)  // check forwards
+        if (s >= -halfCarTolerance)  // check forwards
         {
           // Check z intersection
@@ -1066,5 +1027,5 @@
             // Z ok. Check phi
             //
-            if ( !seg )
+            if ( fPhiFullTube )
             {
               return s ; 
@@ -1098,10 +1059,8 @@
   //         -> use some form of loop Construct ?
   //
-  if ( seg )
+  if ( !fPhiFullTube )
   {
     // First phi surface (Starting phi)
-
-    sinSPhi = std::sin(fSPhi) ;
-    cosSPhi = std::cos(fSPhi) ;
+    //
     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
                     
@@ -1110,5 +1069,5 @@
       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
 
-      if ( Dist < kCarTolerance*0.5 )
+      if ( Dist < halfCarTolerance )
       {
         s = Dist/Comp ;
@@ -1135,6 +1094,6 @@
               // - check intersecting with correct half-plane
               //
-              if ((yi*cosCPhi-xi*sinCPhi) <= 0) snxt = s ;
-            }    
+              if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = s; }
+            }
           }
         }
@@ -1142,9 +1101,6 @@
     }
       
-    // 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) ;
         
@@ -1153,5 +1109,5 @@
       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
 
-      if ( Dist < kCarTolerance*0.5 )
+      if ( Dist < halfCarTolerance )
       {
         s = Dist/Comp ;
@@ -1177,12 +1133,12 @@
               // - check intersecting with correct half-plane
               //
-              if ( (yi*cosCPhi-xi*sinCPhi) >= 0 )  { snxt = s; }
-            }    
+              if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = s; }
+            }                         //?? >=-halfCarTolerance
           }
         }
       }
     }         //  Comp < 0
-  }           //  seg != 0
-  if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
+  }           //  !fPhiFullTube 
+  if ( snxt<halfCarTolerance )  { snxt=0; }
   return snxt ;
 }
@@ -1217,5 +1173,5 @@
 {
   G4double safe=0.0, rho, safe1, safe2, safe3 ;
-  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
+  G4double safePhi, cosPsi ;
 
   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
@@ -1228,26 +1184,21 @@
   if ( safe3 > safe )  { safe = safe3; }
 
-  if (fDPhi < twopi && rho)
-  {
-    phiC    = fSPhi + fDPhi*0.5 ;
-    cosPhiC = std::cos(phiC) ;
-    sinPhiC = std::sin(phiC) ;
-
+  if ( (!fPhiFullTube) && (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 )
-      {
-        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)) ;
+        safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
       }
       if ( safePhi > safe )  { safe = safePhi; }
@@ -1269,12 +1220,13 @@
                                       G4ThreeVector *n    ) const
 {  
-  ESide side = kNull , sider = kNull, sidephi = kNull ;
-  G4double snxt, sr = kInfinity, sphi = kInfinity, pdist ;
+  ESide side=kNull , sider=kNull, sidephi=kNull ;
+  G4double snxt, sr=kInfinity, sphi=kInfinity, pdist ;
   G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
 
+  static const G4double halfCarTolerance = kCarTolerance*0.5;
+  static const G4double halfAngTolerance = kAngTolerance*0.5;
+ 
   // Vars for phi intersection:
 
-  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
-  G4double cPhi, sinCPhi, cosCPhi ;
   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
  
@@ -1284,5 +1236,5 @@
   {
     pdist = fDz - p.z() ;
-    if ( pdist > kCarTolerance*0.5 )
+    if ( pdist > halfCarTolerance )
     {
       snxt = pdist/v.z() ;
@@ -1303,5 +1255,5 @@
     pdist = fDz + p.z() ;
 
-    if ( pdist > kCarTolerance*0.5 )
+    if ( pdist > halfCarTolerance )
     {
       snxt = -pdist/v.z() ;
@@ -1326,5 +1278,5 @@
   // Radial Intersections
   //
-  // Find intersction with cylinders at rmax/rmin
+  // Find intersection with cylinders at rmax/rmin
   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
   //
@@ -1359,7 +1311,7 @@
         b     = t2/t1 ;
         c     = deltaR/t1 ;
-        d2= b*b-c;
-        if(d2>=0.){sr    = -b + std::sqrt(d2);}
-        else{sr=0.;};
+        d2    = b*b-c;
+        if( d2 >= 0 ) { sr = -b + std::sqrt(d2); }
+        else          { sr = 0.; }
         sider = kRMax ;
       }
@@ -1371,12 +1323,4 @@
         if ( calcNorm ) 
         {
-          // if ( p.x() || p.y() )
-          // {
-          //  *n=G4ThreeVector(p.x(),p.y(),0);
-          // }
-          // else
-          // {
-          //  *n=v;
-          // }
           *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
           *validNorm = true ;
@@ -1417,5 +1361,5 @@
           c     = deltaR/t1 ;
           d2    = b*b-c;
-          if(d2>=0.)
+          if( d2 >=0. )
           {
             sr     = -b + std::sqrt(d2) ;
@@ -1423,5 +1367,5 @@
           }
           else // Case: On the border+t2<kRadTolerance
-               //       (v is perpendiculair to the surface)
+               //       (v is perpendicular to the surface)
           {
             if (calcNorm)
@@ -1441,5 +1385,5 @@
         c      = deltaR/t1;
         d2     = b*b-c;
-        if(d2>=0.)
+        if( d2 >= 0 )
         {
           sr     = -b + std::sqrt(d2) ;
@@ -1447,5 +1391,5 @@
         }
         else // Case: On the border+t2<kRadTolerance
-             //       (v is perpendiculair to the surface)
+             //       (v is perpendicular to the surface)
         {
           if (calcNorm)
@@ -1461,22 +1405,13 @@
     // 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 ( !fPhiFullTube )
+    {
       // 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; }
+     
+      if ( vphi < fSPhi - halfAngTolerance  )             { vphi += twopi; }
+      else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
 
 
@@ -1492,10 +1427,11 @@
         compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
         compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
+       
         sidephi = kNull;
         
-        if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)
-                              && (pDistE <= 0.5*kCarTolerance) ) )
-         || ( (fDPhi >  pi) && !((pDistS >  0.5*kCarTolerance)
-                              && (pDistE >  0.5*kCarTolerance) ) )  )
+        if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
+                              && (pDistE <= halfCarTolerance) ) )
+         || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
+                              && (pDistE >  halfCarTolerance) ) )  )
         {
           // Inside both phi *full* planes
@@ -1505,5 +1441,5 @@
             sphi = pDistS/compS ;
             
-            if (sphi >= -0.5*kCarTolerance)
+            if (sphi >= -halfCarTolerance)
             {
               xi = p.x() + sphi*v.x() ;
@@ -1513,13 +1449,14 @@
               // (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))
+              if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
+              {
+                sidephi = kSPhi;
+                if (((fSPhi-halfAngTolerance)<=vphi)
+                   &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
                 {
                   sphi = kInfinity;
                 }
               }
-              else if ((yi*cosCPhi-xi*sinCPhi)>=0)
+              else if ( yi*cosCPhi-xi*sinCPhi >=0 )
               {
                 sphi = kInfinity ;
@@ -1528,5 +1465,5 @@
               {
                 sidephi = kSPhi ;
-                if ( pDistS > -kCarTolerance*0.5 )
+                if ( pDistS > -halfCarTolerance )
                 {
                   sphi = 0.0 ; // Leave by sphi immediately
@@ -1550,5 +1487,5 @@
             // Only check further if < starting phi intersection
             //
-            if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) )
+            if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
             {
               xi = p.x() + sphi2*v.x() ;
@@ -1559,10 +1496,10 @@
                 // Leaving via ending phi
                 //
-                if(!(((fSPhi-0.5*kAngTolerance)<=vphi)
-                    &&((ePhi+0.5*kAngTolerance)>=vphi)))
+                if( !((fSPhi-halfAngTolerance <= vphi)
+                     &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
                 {
                   sidephi = kEPhi ;
-                  if ( pDistE <= -kCarTolerance*0.5 )  { sphi = sphi2 ; }
-                  else                                 { sphi = 0.0 ; }
+                  if ( pDistE <= -halfCarTolerance )  { sphi = sphi2 ; }
+                  else                                { sphi = 0.0 ;   }
                 }
               } 
@@ -1574,6 +1511,6 @@
                 //
                 sidephi = kEPhi ;
-                if ( pDistE <= -kCarTolerance*0.5 ) { sphi = sphi2 ; }
-                else                                { sphi = 0.0 ; }
+                if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
+                else                               { sphi = 0.0 ;   }
               }
             }
@@ -1589,13 +1526,7 @@
         // On z axis + travel not || to z axis -> if phi of vector direction
         // within phi of shape, Step limited by rmax, else Step =0
-
-        // vphi = std::atan2(v.y(),v.x()) ;//defined previosly
-        // G4cout<<"In axis vphi="<<vphi
-        //       <<" Sphi="<<fSPhi<<" Ephi="<<ePhi<<G4endl;
-        // old  if ( (fSPhi < vphi) && (vphi < fSPhi + fDPhi) )
-        // new : correction for if statement, must be '<='  
                
-        if ( ((fSPhi-0.5*kAngTolerance) <= vphi)
-           && (vphi <=( ePhi+0.5*kAngTolerance) )) 
+        if ( (fSPhi - halfAngTolerance <= vphi)
+           && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
         {
           sphi = kInfinity ;
@@ -1640,5 +1571,5 @@
         if ( fDPhi <= pi )
         {
-          *n         = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
+          *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
           *validNorm = true ;
         }
@@ -1652,5 +1583,5 @@
         if (fDPhi <= pi)
         {
-          *n = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
+          *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
           *validNorm = true ;
         }
@@ -1662,6 +1593,6 @@
 
       case kPZ:
-        *n=G4ThreeVector(0,0,1) ;
-        *validNorm=true ;
+        *n         = G4ThreeVector(0,0,1) ;
+        *validNorm = true ;
         break ;
 
@@ -1690,5 +1621,5 @@
     }
   }
-  if ( snxt<kCarTolerance*0.5 )  { snxt=0 ; }
+  if ( snxt<halfCarTolerance )  { snxt=0 ; }
 
   return snxt ;
@@ -1701,6 +1632,5 @@
 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
 {
-  G4double safe=0.0, rho, safeR1, safeR2, safeZ ;
-  G4double safePhi, phiC, cosPhiC, sinPhiC, ePhi ;
+  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
 
@@ -1738,20 +1668,13 @@
   // Check if phi divided, Calc distances closest phi plane
   //
-  if ( fDPhi < twopi )
-  {
-    // Above/below central phi of Tubs?
-
-    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 ( !fPhiFullTube )
+  {
+    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)) ;
+      safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
     }
     if (safePhi < safe)  { safe = safePhi ; }
@@ -1806,6 +1729,6 @@
   // on the x axis. Will give better extent calculations when not rotated.
 
-  if (fDPhi == pi*2.0 && fSPhi == 0 )  { sAngle = -meshAngle*0.5 ; }
-  else                                 { sAngle =  fSPhi ; }
+  if (fPhiFullTube && (fSPhi == 0) )  { sAngle = -meshAngle*0.5 ; }
+  else                                { sAngle =  fSPhi ; }
     
   vertices = new G4ThreeVectorList();
@@ -1975,5 +1898,5 @@
   if (fRMin != 0) 
   {
-    if (fDPhi >= twopi) 
+    if (fPhiFullTube) 
     {
       pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
@@ -1986,5 +1909,5 @@
   else 
   {
-    if (fDPhi >= twopi) 
+    if (fPhiFullTube) 
     {
       pNURBS = new G4NURBScylinder (fRMax,fDz) ;
