Changeset 921 for trunk/source/geometry/solids/CSG

Timestamp:

Feb 16, 2009, 10:14:30 AM (15 years ago)

Author:

garnier

Message:

en test de gl2ps. Problemes de libraries

Location:

trunk/source/geometry/solids/CSG

Files:

• History (modified) (2 diffs)
• include/G4Cons.hh (modified) (3 diffs)
• include/G4Cons.icc (modified) (8 diffs)
• include/G4Sphere.hh (modified) (4 diffs)
• include/G4Sphere.icc (modified) (3 diffs)
• include/G4Tubs.hh (modified) (8 diffs)
• include/G4Tubs.icc (modified) (7 diffs)
• src/G4Cons.cc (modified) (107 diffs)
• src/G4Tubs.cc (modified) (78 diffs)

trunk/source/geometry/solids/CSG/History

@@ -1 @@
-$Id: History,v 1.105 2008/07/07 10:43:04 gcosmo Exp $
+$Id: History,v 1.109 2008/11/21 09:50:20 gcosmo Exp $
 -------------------------------------------------------------------
 
@@ -17 +17 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+Nov 21, 2008  G.Cosmo                    geom-csg-V09-01-08
+- G4Sphere: defined Get/SetInnerRadius() accessors to be compliant with
+  other CSG solids and allow consistent treatment in persistency code...
+
+Nov 06, 2008  G.Cosmo                    geom-csg-V09-01-07
+- G4Tubs, G4Cons: implemented caching of trigonometric values, now directly
+  computed inside modifiers for Phi angles and required for parametrised
+  cases. Improvement bringing up to 20% speedup in normal tracking for
+  tube/cone-sections placements.
+
+Nov 05, 2008  G.Cosmo                    geom-csg-V09-01-06
+- G4Cons: implemented first speed improvements and corrections from joint
+  code review of G4Cons class. Cached computation for half-tolerance and
+  use of Boolean flag for identifying if full-cone or section.
+- G4Tubs: more refinements to previous review; corrected implementation
+  of constructor to conform to implementation as in G4Cons. Fixed issue in
+  SetDeltaPhi() method after changes introduced in the previous tag.
+
+Sep 18, 2008  T.Nikitina                 geom-csg-V09-01-05
+- G4Tubs: implemented first speed improvements and corrections from joint
+  code review of G4Tubs class. Cached computation for half-tolerance and
+  use of Boolean flag for identifying if full-tube or section.
 
 Jul 07, 2008  V.Grichine                 geom-csg-V09-01-04

trunk/source/geometry/solids/CSG/include/G4Cons.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4Cons.hh,v 1.18 2007/05/18 07:38:00 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Cons.hh,v 1.21 2008/11/06 11:04:00 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 //
@@ -54 +54 @@
 //  fDPhi  delta angle of the segment in radians
 //
+//  fPhiFullCone   Boolean variable used for indicate the Phi Section
+//
 //   Note:
 //      Internally fSPhi & fDPhi are adjusted so that fDPhi<=2PI,
@@ -73 +75 @@
   public:  // with description
 
-        G4Cons(const G4String& pName,
-                     G4double pRmin1, G4double pRmax1,
-                     G4double pRmin2, G4double pRmax2,
-                     G4double pDz,
-                     G4double pSPhi, G4double pDPhi);
+    G4Cons(const G4String& pName,
+                 G4double pRmin1, G4double pRmax1,
+                 G4double pRmin2, G4double pRmax2,
+                 G4double pDz,
+                 G4double pSPhi, G4double pDPhi);
          
-        virtual ~G4Cons() ;
-
-  // Accessors
-
-        inline G4double    GetInnerRadiusMinusZ() const;
-        inline G4double    GetOuterRadiusMinusZ() const;
-        inline G4double    GetInnerRadiusPlusZ()  const;
-        inline G4double    GetOuterRadiusPlusZ()  const;
-  
-        inline G4double    GetZHalfLength()       const;
-  
-        inline G4double    GetStartPhiAngle () const;
-        inline G4double    GetDeltaPhiAngle () const;
-  
-  // Modifiers
-
-        inline void    SetInnerRadiusMinusZ( G4double Rmin1 );
-        inline void    SetOuterRadiusMinusZ( G4double Rmax1 );
-        inline void    SetInnerRadiusPlusZ ( G4double Rmin2 );
-        inline void    SetOuterRadiusPlusZ ( G4double Rmax2 );
+    virtual ~G4Cons() ;
+
+    // Accessors
+
+    inline G4double    GetInnerRadiusMinusZ() const;
+    inline G4double    GetOuterRadiusMinusZ() const;
+    inline G4double    GetInnerRadiusPlusZ()  const;
+    inline G4double    GetOuterRadiusPlusZ()  const;
+  
+    inline G4double    GetZHalfLength()       const;
+  
+    inline G4double    GetStartPhiAngle () const;
+    inline G4double    GetDeltaPhiAngle () const;
+  
+    // Modifiers
+
+    inline void    SetInnerRadiusMinusZ( G4double Rmin1 );
+    inline void    SetOuterRadiusMinusZ( G4double Rmax1 );
+    inline void    SetInnerRadiusPlusZ ( G4double Rmin2 );
+    inline void    SetOuterRadiusPlusZ ( G4double Rmax2 );
          
-        inline void    SetZHalfLength      ( G4double newDz );
-        inline void    SetStartPhiAngle    ( G4double newSPhi);
-        inline void    SetDeltaPhiAngle    ( G4double newDPhi);
-
-  // Other methods for solid
-
-        inline G4double    GetCubicVolume();
-        inline G4double    GetSurfaceArea();
-
-        void ComputeDimensions(G4VPVParameterisation* p,
-                               const G4int n,
-                               const G4VPhysicalVolume* pRep);
-
-        G4bool CalculateExtent(const EAxis pAxis,
-                               const G4VoxelLimits& pVoxelLimit,
-                               const G4AffineTransform& pTransform,
-                                     G4double& pmin, G4double& pmax) const;         
-
-        EInside Inside(const G4ThreeVector& p) const;
-
-        G4ThreeVector SurfaceNormal(const G4ThreeVector& p) const;
-
-        G4double DistanceToIn(const G4ThreeVector& p,
-                               const G4ThreeVector& v) const;
-        G4double DistanceToIn(const G4ThreeVector& p) const;
-        G4double DistanceToOut(const G4ThreeVector& p,
-                               const G4ThreeVector& v,
-                               const G4bool calcNorm=G4bool(false),
-                                     G4bool *validNorm=0,
-                                     G4ThreeVector *n=0) const;             
-        G4double DistanceToOut(const G4ThreeVector& p) const;
-
-        G4GeometryType  GetEntityType() const;
+    inline void    SetZHalfLength      ( G4double newDz );
+    inline void    SetStartPhiAngle    ( G4double newSPhi);
+    inline void    SetDeltaPhiAngle    ( G4double newDPhi);
+
+    // Other methods for solid
+
+    inline G4double    GetCubicVolume();
+    inline G4double    GetSurfaceArea();
+
+    void ComputeDimensions(G4VPVParameterisation* p,
+                           const G4int n,
+                           const G4VPhysicalVolume* pRep);
+
+    G4bool CalculateExtent(const EAxis pAxis,
+                           const G4VoxelLimits& pVoxelLimit,
+                           const G4AffineTransform& pTransform,
+                                 G4double& pmin, G4double& pmax) const;         
+
+    EInside Inside(const G4ThreeVector& p) const;
+
+    G4ThreeVector SurfaceNormal(const G4ThreeVector& p) const;
+
+    G4double DistanceToIn (const G4ThreeVector& p,
+                           const G4ThreeVector& v) const;
+    G4double DistanceToIn (const G4ThreeVector& p) const;
+    G4double DistanceToOut(const G4ThreeVector& p,
+                           const G4ThreeVector& v,
+                           const G4bool calcNorm=G4bool(false),
+                                 G4bool *validNorm=0,
+                                 G4ThreeVector *n=0) const;             
+    G4double DistanceToOut(const G4ThreeVector& p) const;
+
+    G4GeometryType  GetEntityType() const;
         
-        G4ThreeVector GetPointOnSurface() const; 
+    G4ThreeVector GetPointOnSurface() const; 
         
-        std::ostream& StreamInfo(std::ostream& os) const;
-
-  // Visualisation functions
-
-        void          DescribeYourselfTo( G4VGraphicsScene& scene ) const;
-        G4Polyhedron* CreatePolyhedron() const;
-        G4NURBS*      CreateNURBS() const;
+    std::ostream& StreamInfo(std::ostream& os) const;
+
+    // Visualisation functions
+
+    void          DescribeYourselfTo( G4VGraphicsScene& scene ) const;
+    G4Polyhedron* CreatePolyhedron() const;
+    G4NURBS*      CreateNURBS() const;
 
   public:  // without description
        
-        G4Cons(__void__&);
-          // Fake default constructor for usage restricted to direct object
-          // persistency for clients requiring preallocation of memory for
-          // persistifiable objects.
-
-        //  Old access functions
-
-        inline G4double    GetRmin1() const;
-        inline G4double    GetRmax1() const;
-        inline G4double    GetRmin2() const;
-        inline G4double    GetRmax2() const;
-  
-        inline G4double    GetDz()    const;
-  
-        inline G4double    GetSPhi() const;
-        inline G4double    GetDPhi() const;
+    G4Cons(__void__&);
+      //
+      // Fake default constructor for usage restricted to direct object
+      // persistency for clients requiring preallocation of memory for
+      // persistifiable objects.
+
+    //  Old access functions
+
+    inline G4double    GetRmin1() const;
+    inline G4double    GetRmax1() const;
+    inline G4double    GetRmin2() const;
+    inline G4double    GetRmax2() const;
+  
+    inline G4double    GetDz()    const;
+  
+    inline G4double    GetSPhi() const;
+    inline G4double    GetDPhi() const;
 
   protected:
  
-        G4ThreeVectorList*
-        CreateRotatedVertices(const G4AffineTransform& pTransform) const;
-  
-        // Used by distanceToOut
-  
-        enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
-  
-        // used by normal
-  
-        enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
+    G4ThreeVectorList*
+    CreateRotatedVertices(const G4AffineTransform& pTransform) const;
+  
+    G4double fRmin1, fRmin2, fRmax1, fRmax2, fDz, fSPhi, fDPhi;
+    G4bool fPhiFullCone;
+
+    // Used by distanceToOut
+  
+    enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
+  
+    // used by normal
+  
+    enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
 
   private:
 
-        G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector& p) const;
-          // Algorithm for SurfaceNormal() following the original
-          // specification for points not on the surface
+    inline void Initialise();
+      // Reset relevant values to zero
+
+    inline void InitializeTrigonometry();
+      //
+      // Recompute relevant trigonometric values and cache them
+
+    G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector& p) const;
+      //
+      // Algorithm for SurfaceNormal() following the original
+      // specification for points not on the surface
 
   private:
 
-        G4double kRadTolerance, kAngTolerance;
-
-        G4double fRmin1,fRmin2,
-                 fRmax1,fRmax2,
-                 fDz,
-                 fSPhi,fDPhi;
+    G4double kRadTolerance, kAngTolerance;
+      //
+      // Radial and angular tolerances
+
+    G4double sinCPhi, cosCPhi, cosHDPhiOT, cosHDPhiIT,
+             sinSPhi, cosSPhi, sinEPhi, cosEPhi;
+      //
+      // Cached trigonometric values
 };
 

trunk/source/geometry/solids/CSG/include/G4Cons.icc

@@ -25 +25 @@
 //
 //
-// $Id: G4Cons.icc,v 1.6 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Cons.icc,v 1.8 2008/11/06 10:55:40 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 // --------------------------------------------------------------------
@@ -37 +37 @@
 
 inline
-G4double G4Cons::GetInnerRadiusMinusZ() const
-{
-  return fRmin1 ;
-}
-
-inline
-G4double G4Cons::GetOuterRadiusMinusZ() const
-{
-  return fRmax1 ;
-}
-
-inline
-G4double G4Cons::GetInnerRadiusPlusZ() const
-{
-  return fRmin2 ;
-}
-
-inline
-G4double G4Cons::GetOuterRadiusPlusZ() const
-{
-  return fRmax2 ;
-}
-
-inline
-G4double G4Cons::GetZHalfLength() const
-{
-  return fDz ;
-}
-
-inline  
-G4double G4Cons::GetStartPhiAngle() const
-{
-  return fSPhi ;
-}
-
-inline
-G4double G4Cons::GetDeltaPhiAngle() const
-{
-  return fDPhi;
-}
-
-inline
-void G4Cons::SetInnerRadiusMinusZ( G4double Rmin1 )
-{
-  fRmin1= Rmin1 ;
+void G4Cons::Initialise()
+{
   fCubicVolume= 0.;
   fSurfaceArea= 0.;
@@ -87 +44 @@
 }
 
+inline 
+void G4Cons::InitializeTrigonometry()
+{
+  G4double hDPhi = 0.5*fDPhi;                       // half delta phi
+  G4double cPhi  = fSPhi + hDPhi; 
+  G4double ePhi  = fSPhi + fDPhi;
+
+  sinCPhi    = std::sin(cPhi);
+  cosCPhi    = std::cos(cPhi);
+  cosHDPhiIT = std::cos(hDPhi - 0.5*kAngTolerance); // inner/outer tol half dphi
+  cosHDPhiOT = std::cos(hDPhi + 0.5*kAngTolerance);
+  sinSPhi = std::sin(fSPhi);
+  cosSPhi = std::cos(fSPhi);
+  sinEPhi = std::sin(ePhi);
+  cosEPhi = std::cos(ePhi);
+}
+
+inline
+G4double G4Cons::GetInnerRadiusMinusZ() const
+{
+  return fRmin1 ;
+}
+
+inline
+G4double G4Cons::GetOuterRadiusMinusZ() const
+{
+  return fRmax1 ;
+}
+
+inline
+G4double G4Cons::GetInnerRadiusPlusZ() const
+{
+  return fRmin2 ;
+}
+
+inline
+G4double G4Cons::GetOuterRadiusPlusZ() const
+{
+  return fRmax2 ;
+}
+
+inline
+G4double G4Cons::GetZHalfLength() const
+{
+  return fDz ;
+}
+
+inline  
+G4double G4Cons::GetStartPhiAngle() const
+{
+  return fSPhi ;
+}
+
+inline
+G4double G4Cons::GetDeltaPhiAngle() const
+{
+  return fDPhi;
+}
+
+inline
+void G4Cons::SetInnerRadiusMinusZ( G4double Rmin1 )
+{
+  fRmin1= Rmin1 ;
+  Initialise();
+}
+
 inline
 void G4Cons::SetOuterRadiusMinusZ( G4double Rmax1 )
 {
   fRmax1= Rmax1 ;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialise();
 }
 
@@ -100 +121 @@
 {
   fRmin2= Rmin2 ;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialise();
 }
 
@@ -109 +128 @@
 {
   fRmax2= Rmax2 ;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialise();
 }
 
@@ -118 +135 @@
 {
   fDz= newDz ;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialise();
 }
 
@@ -127 +142 @@
 {
   fSPhi= newSPhi;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialise();
+  InitializeTrigonometry();
 }
 
 void G4Cons::SetDeltaPhiAngle ( G4double newDPhi )
 {
+  if ( newDPhi >= twopi-kAngTolerance*0.5 )
+  {
+    fPhiFullCone = true;
+  }
+  else if ( newDPhi > 0 )
+  {
+    fPhiFullCone = false;
+  }
+  else
+  {
+     G4cerr << "ERROR - G4Cons()::SetDeltaPhiAngle() : " << GetName() << G4endl
+            << "        Negative delta-Phi ! - " << newDPhi << G4endl;
+     G4Exception("G4Cons::SetDeltaPhiAngle()", "InvalidSetup",
+                  FatalException, "Invalid dphi.");
+  }
   fDPhi=  newDPhi;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialise();
+  InitializeTrigonometry();
 }
 
@@ -220 +248 @@
                          + 0.5*(fRmax1*fRmax1-fRmin1*fRmin1
                                +fRmax2*fRmax2-fRmin2*fRmin2 ));
-    if(fDPhi < twopi )
+    if(!fPhiFullCone)
     {
       fSurfaceArea = fSurfaceArea+4*fDz*(mmax-mmin);

trunk/source/geometry/solids/CSG/include/G4Sphere.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4Sphere.hh,v 1.20 2007/05/18 07:38:00 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Sphere.hh,v 1.21 2008/11/21 09:50:05 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 //
@@ -88 +88 @@
     // Accessors
        
-    inline G4double GetInsideRadius   () const;
+    inline G4double GetInnerRadius    () const;
     inline G4double GetOuterRadius    () const;
     inline G4double GetStartPhiAngle  () const;
@@ -97 +97 @@
     // Modifiers
 
-    inline void SetInsideRadius   (G4double newRmin);
+    inline void SetInnerRadius    (G4double newRMin);
     inline void SetOuterRadius    (G4double newRmax);
     inline void SetStartPhiAngle  (G4double newSphi);
@@ -162 +162 @@
     inline G4double  GetSTheta() const;
     inline G4double  GetDTheta() const;
+    inline G4double  GetInsideRadius() const;
+    inline void SetInsideRadius(G4double newRmin);
 
   protected:

trunk/source/geometry/solids/CSG/include/G4Sphere.icc

@@ -25 +25 @@
 //
 //
-// $Id: G4Sphere.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Sphere.icc,v 1.8 2008/11/21 09:50:05 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 // --------------------------------------------------------------------
@@ -43 +43 @@
 
 inline
+G4double G4Sphere::GetInnerRadius() const
+{
+  return fRmin;
+}
+
+inline
 G4double G4Sphere::GetOuterRadius() const
 {
@@ -73 +79 @@
 inline
 void G4Sphere::SetInsideRadius(G4double newRmin)
+{
+  fRmin= newRmin;
+  fCubicVolume= 0.;
+  fSurfaceArea= 0.;
+  fpPolyhedron = 0;
+}
+
+inline
+void G4Sphere::SetInnerRadius(G4double newRmin)
 {
   fRmin= newRmin;

trunk/source/geometry/solids/CSG/include/G4Tubs.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4Tubs.hh,v 1.17 2007/05/18 07:38:00 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Tubs.hh,v 1.21 2008/11/06 10:55:40 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 // 
@@ -39 +39 @@
 //   A tube or tube segment with curved sides parallel to
 //   the z-axis. The tube has a specified half-length along
-//   the z axis, about which it is centred, and a given
+//   the z-axis, about which it is centered, and a given
 //   minimum and maximum radius. A minimum radius of 0
-//   signifies a filled tube /cylinder. The tube segment is
+//   corresponds to filled tube /cylinder. The tube segment is
 //   specified by starting and delta angles for phi, with 0
 //   being the +x axis, PI/2 the +y axis.
@@ -57 +57 @@
 //
 //   fDPhi  Delta angle of the segment.
+//
+//   fPhiFullTube   Boolean variable used for indicate the Phi Section
 
 // History:
@@ -103 +105 @@
     inline void SetStartPhiAngle (G4double newSPhi);
     inline void SetDeltaPhiAngle (G4double newDPhi);
-
+    
     // Methods for solid
 
@@ -144 +146 @@
 
     G4Tubs(__void__&);
+      //
       // Fake default constructor for usage restricted to direct object
       // persistency for clients requiring preallocation of memory for
@@ -164 +167 @@
       // for G4VSolid:: ClipCrossSection and ClipBetweenSections
 
-    G4double fRMin,fRMax,fDz,fSPhi,fDPhi;
-
-    // Used by distanceToOut
+    G4double fRMin, fRMax, fDz, fSPhi, fDPhi;
+    G4bool fPhiFullTube;
+   
+      // Used by distanceToOut
 
     enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
 
-    // used by normal
+      // Used by normal
 
     enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
@@ -176 +180 @@
   private:
 
+    inline void Initialize();
+      //
+      // Reset relevant values to zero
+
+    inline void InitializeTrigonometry();
+      //
+      // Recompute relevant trigonometric values and cache them
+
     G4ThreeVector ApproxSurfaceNormal( const G4ThreeVector& p ) const;
+      //
       // Algorithm for SurfaceNormal() following the original
       // specification for points not on the surface
@@ -183 +196 @@
 
     G4double kRadTolerance, kAngTolerance;
+      //
       // Radial and angular tolerances
+
+    G4double sinCPhi, cosCPhi, cosHDPhiOT, cosHDPhiIT,
+             sinSPhi, cosSPhi, sinEPhi, cosEPhi;
+      //
+      // Cached trigonometric values
 };
 

trunk/source/geometry/solids/CSG/include/G4Tubs.icc

@@ -25 +25 @@
 //
 //
-// $Id: G4Tubs.icc,v 1.7 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: HEAD $
+// $Id: G4Tubs.icc,v 1.11 2008/11/06 10:55:40 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-cand-01 $
 //
 // --------------------------------------------------------------------
@@ -66 +66 @@
 }
 
+inline 
+void G4Tubs::Initialize()
+{
+  fCubicVolume = 0.;
+  fSurfaceArea = 0.;
+  fpPolyhedron = 0;
+}
+
+inline 
+void G4Tubs::InitializeTrigonometry()
+{
+  G4double hDPhi = 0.5*fDPhi;                       // half delta phi
+  G4double cPhi  = fSPhi + hDPhi; 
+  G4double ePhi  = fSPhi + fDPhi;
+
+  sinCPhi    = std::sin(cPhi);
+  cosCPhi    = std::cos(cPhi);
+  cosHDPhiIT = std::cos(hDPhi - 0.5*kAngTolerance); // inner/outer tol half dphi
+  cosHDPhiOT = std::cos(hDPhi + 0.5*kAngTolerance);
+  sinSPhi = std::sin(fSPhi);
+  cosSPhi = std::cos(fSPhi);
+  sinEPhi = std::sin(ePhi);
+  cosEPhi = std::cos(ePhi);
+}
+
 inline
 void G4Tubs::SetInnerRadius (G4double newRMin)
 {
   fRMin= newRMin;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialize();
 }
 
@@ -79 +102 @@
 {
   fRMax= newRMax;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialize();
 }
 
@@ -88 +109 @@
 {
   fDz= newDz;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialize();
 }
 
@@ -97 +116 @@
 {
   fSPhi= newSPhi;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialize();
+  InitializeTrigonometry();
 }
 
@@ -105 +123 @@
 void G4Tubs::SetDeltaPhiAngle (G4double newDPhi)
 {
+  if ( newDPhi >= twopi-kAngTolerance*0.5 )
+  {
+    fPhiFullTube = true;
+  }
+  else if ( newDPhi > 0 )
+  {
+    fPhiFullTube = false;
+  }
+  else
+  {
+     G4cerr << "ERROR - G4Tubs()::SetDeltaPhiAngle() : " << GetName() << G4endl
+            << "        Negative delta-Phi ! - " << newDPhi << G4endl;
+     G4Exception("G4Tubs::SetDeltaPhiAngle()", "InvalidSetup",
+                  FatalException, "Invalid dphi.");
+  }
   fDPhi= newDPhi;
-  fCubicVolume= 0.;
-  fSurfaceArea= 0.;
-  fpPolyhedron = 0;
+  Initialize();
+  InitializeTrigonometry();
 }
 
@@ -158 +190 @@
   {
     fSurfaceArea = fDPhi*(fRMin+fRMax)*(2*fDz+fRMax-fRMin);
-    if (fDPhi < twopi)
+    if (!fPhiFullTube)
     {
       fSurfaceArea = fSurfaceArea + 4*fDz*(fRMax-fRMin);

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

@@ -25 +25 @@
 //
 //
-// $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 +87 @@
 
   if ( pDz > 0 )
-    fDz = pDz ;
+  {
+    fDz = pDz;
+  }
   else
   {
@@ -99 +101 @@
   // Check radii
 
-  if ( pRmin1 < pRmax1 && pRmin2 < pRmax2 && pRmin1 >= 0 && pRmin2 >= 0 )
+  if ( (pRmin1<pRmax1) && (pRmin2<pRmax2) && (pRmin1>=0) && (pRmin2>=0) )
   {
 
@@ -106 +108 @@
     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 +121 @@
   }
 
-  // Check angles
-
-  if ( pDPhi >= twopi )
+  fPhiFullCone = true;
+  if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles
   {
     fDPhi=twopi;
@@ -128 +129 @@
   else
   {
-    if ( pDPhi > 0 ) fDPhi = pDPhi ;
+    fPhiFullCone = false;
+    if ( pDPhi > 0 )
+    {
+      fDPhi = pDPhi;
+    }
     else
     {
@@ -135 +140 @@
              << 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 +187 @@
   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 +201 @@
   // 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 +259 @@
                     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 +280 @@
     if (pVoxelLimit.IsZLimited())
     {
-      if( zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance || 
-          zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance   )
+      if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) || 
+          (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance)  )
       {
         return false ;
@@ -290 +305 @@
     if (pVoxelLimit.IsXLimited())
     {
-      if ( xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance || 
-           xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance    )
+      if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) || 
+           (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance)  )
       {
         return false ;
@@ -315 +330 @@
     if (pVoxelLimit.IsYLimited())
     {
-      if ( yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance || 
-           yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance     )
+      if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) || 
+           (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance)  )
       {
         return false ;
@@ -338 +353 @@
         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 +377 @@
         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 +424 @@
     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 +477 @@
   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 +499 @@
   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 +505 @@
       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 +516 @@
       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 +533 @@
     sumnorm += nR;
   }
-  if( (fRmin1 || fRmin2) && distRMin <= delta )
+  if( (fRmin1 || fRmin2) && (distRMin <= delta) )
   {
     noSurfaces ++;
     sumnorm += nr;
   }
-  if( fDPhi < twopi )   
+  if( !fPhiFullCone )   
   {
     if (distSPhi <= dAngle)
@@ -534 +554 @@
   {
     noSurfaces ++;
-    if ( p.z() >= 0.)  sumnorm += nZ;
-    else               sumnorm -= nZ; 
+    if ( p.z() >= 0.)  { sumnorm += nZ; }
+    else               { sumnorm -= nZ; }
   }
   if ( noSurfaces == 0 )
@@ -545 +565 @@
      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 +626 @@
     }
   }
-  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 +641 @@
     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 +652 @@
     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 +698 @@
 //
 // NOTE:
-// - Precalculations for phi trigonometry are Done `just in time'
 // - `if valid' implies tolerant checking of intersection points
 // - z, phi intersection from Tubs
@@ -686 +706 @@
   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 +722 @@
   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 +729 @@
   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 +742 @@
   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 +755 @@
       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 +766 @@
       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 +790 @@
         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 +849 @@
   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 +869 @@
         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 +893 @@
             // Z ok. Check phi intersection if reqd
 
-            if ( ! seg )  return s ;
+            if ( fPhiFullCone )  { return s; }
             else
             {
@@ -914 +901 @@
               cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-              if ( cosPsi >= cosHDPhiIT ) return s ;
+              if ( cosPsi >= cosHDPhiIT )  { return s; }
             }
           }
@@ -925 +912 @@
       // 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 +934 @@
       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 +951 @@
             cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-            if (cosPsi >= cosHDPhiIT) return s ;
+            if (cosPsi >= cosHDPhiIT)  { return s; }
           }
         }
@@ -1015 +997 @@
             if ( std::fabs(zi) <= tolODz )
             {
-              if ( seg )
+              if ( !fPhiFullCone )
               {
                 xi     = p.x() + s*v.x() ;
@@ -1022 +1004 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiIT) snxt = s ; 
+                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
               }
-              else return s ;
+              else  { return s; }
             }
           }
@@ -1046 +1028 @@
           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 +1038 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiOT) snxt = s ; 
+                if (cosPsi >= cosHDPhiOT)  { snxt = s; }
               }
-              else  return s ;
+              else  { return s; }
             }
           }
@@ -1067 +1049 @@
             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 +1057 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
 
-                if (cosPsi >= cosHDPhiIT) snxt = s ;
+                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
               }
-              else  return s ;
+              else  { return s; }
             }
           }
@@ -1095 +1077 @@
             // 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 +1105 @@
                 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 +1114 @@
                     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 +1134 @@
             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 +1143 @@
                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
 
-                if (cosPsi >= cosHDPhiIT) snxt = s ; 
+                if (cosPsi >= cosHDPhiIT)  { snxt = s; }
               }
-              else  return s ;
+              else  { return s; }
             }
           }
@@ -1180 +1162 @@
   //         -> 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 +1172 @@
       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
 
-      if (Dist < kCarTolerance*0.5)
+      if (Dist < halfCarTolerance)
       {
         s = Dist/Comp ;
@@ -1198 +1178 @@
         if ( s < snxt )
         {
-          if ( s < 0 ) s = 0.0 ;
+          if ( s < 0 )  { s = 0.0; }
 
           zi = p.z() + s*v.z() ;
@@ -1210 +1190 @@
             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 +1209 @@
     {
       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
-      if (Dist < kCarTolerance*0.5)
+      if (Dist < halfCarTolerance)
       {
         s = Dist/Comp ;
@@ -1237 +1215 @@
         if ( s < snxt )
         {
-          if ( s < 0 )  s = 0.0 ;
+          if ( s < 0 )  { s = 0.0; }
 
           zi = p.z() + s*v.z() ;
@@ -1249 +1227 @@
             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 +1239 @@
     }
   }
-  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 +1253 @@
 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 +1272 @@
     safeR2  = (rho - pRMax)/secRMax ;
 
-    if ( safeR1 > safeR2) safe = safeR1 ;
-    else                  safe = safeR2 ;
+    if ( safeR1 > safeR2) { safe = safeR1; }
+    else                  { safe = safeR2; }
   }
   else
@@ -1314 +1282 @@
     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 +1310 @@
 ///////////////////////////////////////////////////////////////
 //
-// 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 +1340 @@
   // 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 +1349 @@
     pdist = fDz - p.z() ;
 
-    if (pdist > kCarTolerance*0.5)
+    if (pdist > halfCarTolerance)
     {
       snxt = pdist/v.z() ;
@@ -1396 +1361 @@
         *validNorm = true ;
       }
-      return snxt = 0.0 ;
+      return  snxt = 0.0;
     }
   }
@@ -1403 +1368 @@
     pdist = fDz + p.z() ;
 
-    if ( pdist > kCarTolerance*0.5)
+    if ( pdist > halfCarTolerance)
     {
       snxt = -pdist/v.z() ;
@@ -1465 +1430 @@
                 - 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 +1446 @@
     {
       // 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 +1454 @@
           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 +1465 @@
         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 +1473 @@
           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 +1482 @@
           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 +1509 @@
       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 +1522 @@
     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 +1538 @@
   // Check possible intersection within tolerance
 
-  if ( slentol <= kCarTolerance*0.5 )
+  if ( slentol <= halfCarTolerance )
   {
     // An intersection within the tolerance was found.  
@@ -1580 +1549 @@
     // 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 +1564 @@
     }
     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 +1590 @@
       d = b*b - c ;
 
-      if (d >= 0.0 )
+      if ( d >= 0.0 )
       {
         // NOTE: should be rho-rin<kRadTolerance*0.5,
@@ -1630 +1599 @@
           if ( nt2 < 0.0 )
           {
-            if (calcNorm)  *validNorm = false ;
-            return          snxt      = 0.0 ;
+            if (calcNorm)  { *validNorm = false; }
+            return          snxt      = 0.0;
           }
         }
@@ -1640 +1609 @@
           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 +1617 @@
             sidetol = kRMax ;
           }
-          if( (ri<0) || (sr2 < kRadTolerance*0.5) )
+          if( (ri<0) || (sr2 < halfRadTolerance) )
           {
             sr3 = -b + std::sqrt(d) ;
@@ -1656 +1624 @@
             //         distancetoout
 
-            if  ( sr3 > kCarTolerance*0.5 )
+            if  ( sr3 > halfRadTolerance )
             {
               if( sr3 < sr )
@@ -1670 +1638 @@
               } 
             }
-            else if ( sr3 > -kCarTolerance*0.5 )
+            else if ( sr3 > -halfRadTolerance )
             {
               // Intersection in tolerance. Store to check if it's good
@@ -1678 +1646 @@
             }
           }
-          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 +1658 @@
             sidetol = kRMin ;
           }    
-          if( slentol <= kCarTolerance*0.5  )
+          if( slentol <= halfCarTolerance  )
           {
             // An intersection within the tolerance was found. 
@@ -1706 +1674 @@
             if( Normal.dot(v) > 0 )
             {
-              // We will leave the Cone immediatelly
+              // We will leave the cone immediately
+
               if( calcNorm ) 
               {
@@ -1731 +1700 @@
   // 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 +1724 @@
       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 +1769 @@
             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 +1823 @@
       // 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 +1836 @@
     if ( sphi < snxt )  // Order intersecttions
     {
-        snxt=sphi ;
-        side=sidephi ;
+      snxt=sphi ;
+      side=sidephi ;
     }
   }
@@ -1880 +1848 @@
   {
     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 +1857 @@
         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 +1897 @@
         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 +1915 @@
     }
   }
-  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 +1926 @@
 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 +1940 @@
     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 +1962 @@
     safeR1  = (rho - pRMin)/secRMin ;
   }
-  else safeR1 = kInfinity ;
+  else
+  {
+    safeR1 = kInfinity ;
+  }
 
   tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
@@ -2004 +1972 @@
   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 +2032 @@
   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 +2038 @@
   // 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 +2064 @@
       rMaxY2 = meshRMax2*sinCrossAngle ;
         
-      // G4double RMin = (fRmin2 <= fRmin1) ? fRmin2 : fRmin1  ;
-
       rMinX1 = fRmin1*cosCrossAngle ;
       rMinY1 = fRmin1*sinCrossAngle ;
@@ -2127 +2087 @@
                 "Error in allocation of vertices. Out of memory !");
   }
+
   return vertices ;
 }
@@ -2192 +2153 @@
   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 +2186 @@
   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)

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

@@ -25 +25 @@
 //
 //
-// $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 +108 @@
                 "Invalid Z half-length");
   }
-  if ( pRMin < pRMax && pRMin >= 0 ) // Check radii
+  if ( (pRMin < pRMax) && (pRMin >= 0) ) // Check radii
   {
     fRMin = pRMin ; 
@@ -121 +121 @@
                 "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 +140 @@
              << "        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 +293 @@
         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 +320 @@
         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 +366 @@
     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 +414 @@
   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 +438 @@
         // 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 +444 @@
         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 +467 @@
           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 +484 @@
     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 +491 @@
       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 +499 @@
           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 +515 @@
           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 +526 @@
     }
   }
-  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 +536 @@
     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 +544 @@
         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 +560 @@
         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 +580 @@
 
 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 +599 @@
   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 +619 @@
     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 +633 @@
   if( fDPhi < twopi )   
   {
-    if (distSPhi <= dAngle)  // delta)
+    if (distSPhi <= halfAngTolerance)  
     {
       noSurfaces ++;
       sumnorm += nPs;
     }
-    if (distEPhi <= dAngle) // delta) 
+    if (distEPhi <= halfAngTolerance)  
     {
       noSurfaces ++;
@@ -649 +644 @@
     }
   }
-  if (distZ <= delta)  
+  if (distZ <= halfCarTolerance)  
   {
     noSurfaces ++;
@@ -668 +663 @@
   else if ( noSurfaces == 1 )  { norm = sumnorm; }
   else                         { norm = sumnorm.unit(); }
+
   return norm;
 }
@@ -714 +710 @@
     }
   }   
-  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 +745 @@
     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 +800 @@
 //
 // 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 +829 @@
     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 +841 @@
     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 +851 @@
       // 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 +860 @@
           iden   = std::sqrt(rho2) ;
           cosPsi = inum/iden ;
-          if (cosPsi >= cosHDPhiIT) return s ;
+          if (cosPsi >= cosHDPhiIT)  { return s ; }
         }
         else
@@ -909 +870 @@
     else
     {
-      if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
+      if ( snxt<halfCarTolerance )  { snxt=0; }
       return snxt ;  // On/outside extent, and heading away
                      // -> cannot intersect
@@ -934 +895 @@
     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 +915 @@
             // Z ok. Check phi intersection if reqd
             //
-            if (!seg)
+            if (fPhiFullTube)
             {
               return s ;
@@ -963 +924 @@
               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 +935 @@
       // 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 +965 @@
                 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 +996 @@
               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 +1004 @@
             }
           }
-        } // end if   (seg)
+        } // end if   (!fPhiFullTube)
       }   // end if   (t3>tolIRMin2)
     }     // end if   (Inside Outer Radius) 
@@ -1056 +1017 @@
 
         s = -b + std::sqrt(d) ;
-        if (s >= -0.5*kCarTolerance)  // check forwards
+        if (s >= -halfCarTolerance)  // check forwards
         {
           // Check z intersection
@@ -1066 +1027 @@
             // Z ok. Check phi
             //
-            if ( !seg )
+            if ( fPhiFullTube )
             {
               return s ; 
@@ -1098 +1059 @@
   //         -> 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 +1069 @@
       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
 
-      if ( Dist < kCarTolerance*0.5 )
+      if ( Dist < halfCarTolerance )
       {
         s = Dist/Comp ;
@@ -1135 +1094 @@
               // - check intersecting with correct half-plane
               //
-              if ((yi*cosCPhi-xi*sinCPhi) <= 0) snxt = s ;
-            }    
+              if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = s; }
+            }
           }
         }
@@ -1142 +1101 @@
     }
       
-    // 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 +1109 @@
       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
 
-      if ( Dist < kCarTolerance*0.5 )
+      if ( Dist < halfCarTolerance )
       {
         s = Dist/Comp ;
@@ -1177 +1133 @@
               // - 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 +1173 @@
 {
   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 +1184 @@
   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 +1220 @@
                                       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 +1236 @@
   {
     pdist = fDz - p.z() ;
-    if ( pdist > kCarTolerance*0.5 )
+    if ( pdist > halfCarTolerance )
     {
       snxt = pdist/v.z() ;
@@ -1303 +1255 @@
     pdist = fDz + p.z() ;
 
-    if ( pdist > kCarTolerance*0.5 )
+    if ( pdist > halfCarTolerance )
     {
       snxt = -pdist/v.z() ;
@@ -1326 +1278 @@
   // 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 +1311 @@
         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 +1323 @@
         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 +1361 @@
           c     = deltaR/t1 ;
           d2    = b*b-c;
-          if(d2>=0.)
+          if( d2 >=0. )
           {
             sr     = -b + std::sqrt(d2) ;
@@ -1423 +1367 @@
           }
           else // Case: On the border+t2<kRadTolerance
-               //       (v is perpendiculair to the surface)
+               //       (v is perpendicular to the surface)
           {
             if (calcNorm)
@@ -1441 +1385 @@
         c      = deltaR/t1;
         d2     = b*b-c;
-        if(d2>=0.)
+        if( d2 >= 0 )
         {
           sr     = -b + std::sqrt(d2) ;
@@ -1447 +1391 @@
         }
         else // Case: On the border+t2<kRadTolerance
-             //       (v is perpendiculair to the surface)
+             //       (v is perpendicular to the surface)
         {
           if (calcNorm)
@@ -1461 +1405 @@
     // 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 +1427 @@
         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 +1441 @@
             sphi = pDistS/compS ;
             
-            if (sphi >= -0.5*kCarTolerance)
+            if (sphi >= -halfCarTolerance)
             {
               xi = p.x() + sphi*v.x() ;
@@ -1513 +1449 @@
               // (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 +1465 @@
               {
                 sidephi = kSPhi ;
-                if ( pDistS > -kCarTolerance*0.5 )
+                if ( pDistS > -halfCarTolerance )
                 {
                   sphi = 0.0 ; // Leave by sphi immediately
@@ -1550 +1487 @@
             // 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 +1496 @@
                 // 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 +1511 @@
                 //
                 sidephi = kEPhi ;
-                if ( pDistE <= -kCarTolerance*0.5 ) { sphi = sphi2 ; }
-                else                                { sphi = 0.0 ; }
+                if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
+                else                               { sphi = 0.0 ;   }
               }
             }
@@ -1589 +1526 @@
         // 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 +1571 @@
         if ( fDPhi <= pi )
         {
-          *n         = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
+          *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
           *validNorm = true ;
         }
@@ -1652 +1583 @@
         if (fDPhi <= pi)
         {
-          *n = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
+          *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
           *validNorm = true ;
         }
@@ -1662 +1593 @@
 
       case kPZ:
-        *n=G4ThreeVector(0,0,1) ;
-        *validNorm=true ;
+        *n         = G4ThreeVector(0,0,1) ;
+        *validNorm = true ;
         break ;
 
@@ -1690 +1621 @@
     }
   }
-  if ( snxt<kCarTolerance*0.5 )  { snxt=0 ; }
+  if ( snxt<halfCarTolerance )  { snxt=0 ; }
 
   return snxt ;
@@ -1701 +1632 @@
 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 +1668 @@
   // 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 +1729 @@
   // 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 +1898 @@
   if (fRMin != 0) 
   {
-    if (fDPhi >= twopi) 
+    if (fPhiFullTube) 
     {
       pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
@@ -1986 +1909 @@
   else 
   {
-    if (fDPhi >= twopi) 
+    if (fPhiFullTube) 
     {
       pNURBS = new G4NURBScylinder (fRMax,fDz) ;