Index: trunk/source/geometry/solids/CSG/src/G4Box.cc
===================================================================
--- trunk/source/geometry/solids/CSG/src/G4Box.cc	(revision 1228)
+++ trunk/source/geometry/solids/CSG/src/G4Box.cc	(revision 1315)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4Box.cc,v 1.49 2010/05/25 10:14:41 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // 
@@ -32,5 +32,5 @@
 // Implementation for G4Box class
 //
-//  24.06.98 - V. Grichine: insideEdge in DistanceToIn(p,v)
+//  24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
 //  20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
 //  07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
@@ -67,5 +67,5 @@
   if ( (pX > 2*kCarTolerance)
     && (pY > 2*kCarTolerance)
-    && (pZ > 2*kCarTolerance) )
+    && (pZ > 2*kCarTolerance) )  // limit to thickness of surfaces
   {
     fDx = pX ;
@@ -105,6 +105,8 @@
 void G4Box::SetXHalfLength(G4double dx)
 {
-  if(dx > 2*kCarTolerance)
+  if(dx > 2*kCarTolerance)  // limit to thickness of surfaces
+  {
     fDx = dx;
+  }
   else
   {
@@ -122,6 +124,8 @@
 void G4Box::SetYHalfLength(G4double dy) 
 {
-  if(dy > 2*kCarTolerance)
+  if(dy > 2*kCarTolerance)  // limit to thickness of surfaces
+  {
     fDy = dy;
+  }
   else
   {
@@ -139,6 +143,8 @@
 void G4Box::SetZHalfLength(G4double dz) 
 {
-  if(dz > 2*kCarTolerance)
+  if(dz > 2*kCarTolerance)  // limit to thickness of surfaces
+  {
     fDz = dz;
+  }
   else
   {
@@ -153,6 +159,4 @@
   fpPolyhedron = 0;
 } 
-    
-
 
 ////////////////////////////////////////////////////////////////////////
@@ -193,16 +197,10 @@
     if (pVoxelLimit.IsXLimited())
     {
-      if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance || 
-           xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance    ) return false ;
+      if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) || 
+          (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; }
       else
       {
-        if (xMin < pVoxelLimit.GetMinXExtent())
-        {
-          xMin = pVoxelLimit.GetMinXExtent() ;
-        }
-        if (xMax > pVoxelLimit.GetMaxXExtent())
-        {
-          xMax = pVoxelLimit.GetMaxXExtent() ;
-        }
+        xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
+        xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
       }
     }
@@ -213,16 +211,10 @@
     if (pVoxelLimit.IsYLimited())
     {
-      if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
-           yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance   ) return false ;
+      if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) ||
+          (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; }
       else
       {
-        if (yMin < pVoxelLimit.GetMinYExtent())
-        {
-          yMin = pVoxelLimit.GetMinYExtent() ;
-        }
-        if (yMax > pVoxelLimit.GetMaxYExtent())
-        {
-          yMax = pVoxelLimit.GetMaxYExtent() ;
-        }
+        yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
+        yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
       }
     }
@@ -233,16 +225,10 @@
     if (pVoxelLimit.IsZLimited())
     {
-      if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
-           zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance   ) return false ;
+      if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) ||
+          (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; }
       else
       {
-        if (zMin < pVoxelLimit.GetMinZExtent())
-        {
-          zMin = pVoxelLimit.GetMinZExtent() ;
-        }
-        if (zMax > pVoxelLimit.GetMaxZExtent())
-        {
-          zMax = pVoxelLimit.GetMaxZExtent() ;
-        }
+        zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
+        zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
       }
     }
@@ -286,5 +272,5 @@
     if (pVoxelLimit.IsLimited(pAxis) == false) 
     {  
-      if ( pMin != kInfinity || pMax != -kInfinity ) 
+      if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 
       {
         existsAfterClip = true ;
@@ -292,6 +278,6 @@
         // Add 2*tolerance to avoid precision troubles
 
-        pMin           -= kCarTolerance;
-        pMax           += kCarTolerance;
+        pMin -= kCarTolerance;
+        pMax += kCarTolerance;
       }
     }      
@@ -303,5 +289,5 @@
        ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
 
-      if ( pMin != kInfinity || pMax != -kInfinity )
+      if ( (pMin != kInfinity) || (pMax != -kInfinity) )
       {
         existsAfterClip = true ;
@@ -356,23 +342,25 @@
 EInside G4Box::Inside(const G4ThreeVector& p) const
 {
+  static const G4double delta=0.5*kCarTolerance;
   EInside in = kOutside ;
-
-  if ( std::fabs(p.x()) <= fDx - kCarTolerance*0.5 )
-  {
-    if (std::fabs(p.y()) <= fDy - kCarTolerance*0.5 )
-    {
-      if      (std::fabs(p.z()) <= fDz - kCarTolerance*0.5 ) in = kInside ;
-      else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
-    }
-    else if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
-    {
-      if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
-    }
-  }
-  else if (std::fabs(p.x()) <= fDx + kCarTolerance*0.5 )
-  {
-    if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
-    {
-      if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5) in = kSurface ;
+  G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
+
+  if ( q.x() <= (fDx - delta) )
+  {
+    if (q.y() <= (fDy - delta) )
+    {
+      if      ( q.z() <= (fDz - delta) ) { in = kInside ;  }
+      else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
+    }
+    else if ( q.y() <= (fDy + delta) )
+    {
+      if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
+    }
+  }
+  else if ( q.x() <= (fDx + delta) )
+  {
+    if ( q.y() <= (fDy + delta) )
+    {
+      if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
     }
   }
@@ -389,5 +377,5 @@
 {
   G4double distx, disty, distz ;
-  G4ThreeVector norm ;
+  G4ThreeVector norm(0.,0.,0.);
 
   // Calculate distances as if in 1st octant
@@ -400,5 +388,5 @@
   // normals
 
-  const G4double delta    = 0.5*kCarTolerance;
+  static const G4double delta    = 0.5*kCarTolerance;
   const G4ThreeVector nX  = G4ThreeVector( 1.0, 0,0  );
   const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0  );
@@ -415,9 +403,6 @@
   {
     noSurfaces ++; 
-    if ( p.x() >= 0.){        // on +X surface
-      normX= nX ;    // G4ThreeVector( 1.0, 0., 0. );
-    }else{
-      normX= nmX;    // G4ThreeVector(-1.0, 0., 0. ); 
-    }
+    if ( p.x() >= 0. )  { normX= nX ; }       // on +X surface : (1,0,0)
+    else                { normX= nmX; }       //                 (-1,0,0)
     sumnorm= normX; 
   }
@@ -426,9 +411,6 @@
   {
     noSurfaces ++; 
-    if ( p.y() >= 0.){        // on +Y surface
-      normY= nY;
-    }else{
-      normY = nmY; 
-    }
+    if ( p.y() >= 0. )  { normY= nY;  }       // on +Y surface
+    else                { normY= nmY; }
     sumnorm += normY; 
   }
@@ -437,32 +419,35 @@
   {
     noSurfaces ++; 
-    if ( p.z() >= 0.){        // on +Z surface
-      normZ= nZ;
-    }else{
-      normZ = nmZ; 
-    }
-    sumnorm += normZ; 
-  }
-
-  // sumnorm= normX + normY + normZ; 
-  const G4double invSqrt2 = 1.0 / std::sqrt( 2.0); 
-  const G4double invSqrt3 = 1.0 / std::sqrt( 3.0); 
-
-  norm= G4ThreeVector( 0., 0., 0.); 
+    if ( p.z() >= 0. )  { normZ= nZ;  }       // on +Z surface
+    else                { normZ= nmZ; }
+    sumnorm += normZ;
+  }
+
+  static const G4double invSqrt2 = 1.0 / std::sqrt(2.0); 
+  static const G4double invSqrt3 = 1.0 / std::sqrt(3.0); 
+
   if( noSurfaces > 0 )
   { 
-    if( noSurfaces == 1 ){ 
+    if( noSurfaces == 1 )
+    { 
       norm= sumnorm; 
-    }else{
+    }
+    else
+    {
       // norm = sumnorm . unit(); 
-      if( noSurfaces == 2 ) { 
+      if( noSurfaces == 2 )
+      { 
         // 2 surfaces -> on edge 
         norm = invSqrt2 * sumnorm; 
-      } else { 
+      }
+      else
+      { 
         // 3 surfaces (on corner)
         norm = invSqrt3 * sumnorm; 
       }
     }
-  }else{
+  }
+  else
+  {
 #ifdef G4CSGDEBUG
      G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning, 
@@ -483,5 +468,5 @@
 {
   G4double distx, disty, distz ;
-  G4ThreeVector norm ;
+  G4ThreeVector norm(0.,0.,0.);
 
   // Calculate distances as if in 1st octant
@@ -495,11 +480,11 @@
     if ( distx <= distz )     // Closest to X
     {
-      if ( p.x() < 0 ) norm = G4ThreeVector(-1.0,0,0) ;
-      else             norm = G4ThreeVector( 1.0,0,0) ;
+      if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
+      else             { norm = G4ThreeVector( 1.0,0,0) ; }
     }
     else                      // Closest to Z
     {
-      if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
-      else             norm = G4ThreeVector(0,0, 1.0) ;
+      if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
+      else             { norm = G4ThreeVector(0,0, 1.0) ; }
     }
   }
@@ -508,11 +493,11 @@
     if ( disty <= distz )      // Closest to Y
     {
-      if ( p.y() < 0 ) norm = G4ThreeVector(0,-1.0,0) ;
-      else             norm = G4ThreeVector(0, 1.0,0) ;
+      if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
+      else             { norm = G4ThreeVector(0, 1.0,0) ; }
     }
     else                       // Closest to Z
     {
-      if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
-      else             norm = G4ThreeVector(0,0, 1.0) ;
+      if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
+      else             { norm = G4ThreeVector(0,0, 1.0) ; }
     }
   }
@@ -541,5 +526,6 @@
 // shape.
 
-G4double G4Box::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const
+G4double G4Box::DistanceToIn(const G4ThreeVector& p,
+                             const G4ThreeVector& v) const
 {
   G4double safx, safy, safz ;
@@ -549,4 +535,6 @@
   G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
 
+  static const G4double delta = 0.5*kCarTolerance;
+
   safx = std::fabs(p.x()) - fDx ;     // minimum distance to x surface of shape
   safy = std::fabs(p.y()) - fDy ;
@@ -558,7 +546,7 @@
   // travel is in a direction away from the shape.
 
-  if (    ((p.x()*v.x() >= 0.0) && safx > -kCarTolerance*0.5) 
-       || ((p.y()*v.y() >= 0.0) && safy > -kCarTolerance*0.5)
-       || ((p.z()*v.z() >= 0.0) && safz > -kCarTolerance*0.5)   ) 
+  if (    ((p.x()*v.x() >= 0.0) && (safx > -delta)) 
+       || ((p.y()*v.y() >= 0.0) && (safy > -delta))
+       || ((p.z()*v.z() >= 0.0) && (safz > -delta))   ) 
   {
     return kInfinity ;  // travel away or parallel within tolerance
@@ -568,5 +556,5 @@
   // X Planes
 
-  if ( v.x())
+  if ( v.x() )  // != 0
   {
     stmp = 1.0/std::fabs(v.x()) ;
@@ -579,6 +567,6 @@
     else
     {
-      if (v.x() > 0)  sOut = (fDx - p.x())*stmp ;
-      if (v.x() < 0)  sOut = (fDx + p.x())*stmp ;
+      if (v.x() < 0)  { sOut = (fDx + p.x())*stmp ; }
+      else            { sOut = (fDx - p.x())*stmp ; }
     }
   }
@@ -586,5 +574,5 @@
   // Y Planes
 
-  if ( v.y()) 
+  if ( v.y() )  // != 0
   {
     stmp = 1.0/std::fabs(v.y()) ;
@@ -595,8 +583,8 @@
       smaxy = (fDy+std::fabs(p.y()))*stmp ;
 
-      if (sminy > smin) smin=sminy ;
-      if (smaxy < smax) smax=smaxy ;
-
-      if (smin >= smax-kCarTolerance*0.5)
+      if (sminy > smin) { smin=sminy ; }
+      if (smaxy < smax) { smax=smaxy ; }
+
+      if (smin >= (smax-delta))
       {
         return kInfinity ;  // touch XY corner
@@ -605,7 +593,7 @@
     else
     {
-      if (v.y() > 0)  sOuty = (fDy - p.y())*stmp ;
-      if (v.y() < 0)  sOuty = (fDy + p.y())*stmp ;
-      if( sOuty < sOut ) sOut = sOuty ;
+      if (v.y() < 0)  { sOuty = (fDy + p.y())*stmp ; }
+      else            { sOuty = (fDy - p.y())*stmp ; }
+      if( sOuty < sOut ) { sOut = sOuty ; }
     }     
   }
@@ -613,17 +601,17 @@
   // Z planes
 
-  if ( v.z() )
+  if ( v.z() )  // != 0
   {
     stmp = 1.0/std::fabs(v.z()) ;
 
-    if ( safz >= 0.0)
+    if ( safz >= 0.0 )
     {
       sminz = safz*stmp ;
       smaxz = (fDz+std::fabs(p.z()))*stmp ;
 
-      if (sminz > smin) smin = sminz ;
-      if (smaxz < smax) smax = smaxz ;
-
-      if (smin >= smax-kCarTolerance*0.5)
+      if (sminz > smin) { smin = sminz ; }
+      if (smaxz < smax) { smax = smaxz ; }
+
+      if (smin >= (smax-delta))
       { 
         return kInfinity ;    // touch ZX or ZY corners
@@ -632,15 +620,15 @@
     else
     {
-      if (v.z() > 0)  sOutz = (fDz - p.z())*stmp ;
-      if (v.z() < 0)  sOutz = (fDz + p.z())*stmp ;
-      if( sOutz < sOut ) sOut = sOutz ;
-    }
-  }
-
-  if ( sOut <= smin + 0.5*kCarTolerance) // travel over edge
+      if (v.z() < 0)  { sOutz = (fDz + p.z())*stmp ; }
+      else            { sOutz = (fDz - p.z())*stmp ; }
+      if( sOutz < sOut ) { sOut = sOutz ; }
+    }
+  }
+
+  if (sOut <= (smin + delta)) // travel over edge
   {
     return kInfinity ;
   }
-  if (smin < 0.5*kCarTolerance)  smin = 0.0 ;
+  if (smin < delta)  { smin = 0.0 ; }
 
   return smin ;
@@ -662,7 +650,7 @@
   safez = std::fabs(p.z()) - fDz ;
 
-  if (safex > safe) safe = safex ;
-  if (safey > safe) safe = safey ;
-  if (safez > safe) safe = safez ;
+  if (safex > safe) { safe = safex ; }
+  if (safey > safe) { safe = safey ; }
+  if (safez > safe) { safe = safez ; }
 
   return safe ;
@@ -671,5 +659,5 @@
 /////////////////////////////////////////////////////////////////////////
 //
-// Calcluate distance to surface of box from inside
+// Calculate distance to surface of box from inside
 // by calculating distances to box's x/y/z planes.
 // Smallest distance is exact distance to exiting.
@@ -682,13 +670,15 @@
 {
   ESide side = kUndefined ;
-  G4double pdist,stmp,snxt;
-
-  if (calcNorm) *validNorm = true ; // All normals are valid
-
-  if (v.x() > 0)   // X planes  
+  G4double pdist,stmp,snxt=kInfinity;
+
+  static const G4double delta = 0.5*kCarTolerance;
+
+  if (calcNorm) { *validNorm = true ; }  // All normals are valid
+
+  if (v.x() > 0)   // X planes
   {
     pdist = fDx - p.x() ;
 
-    if (pdist > kCarTolerance*0.5)
+    if (pdist > delta)
     {
       snxt = pdist/v.x() ;
@@ -697,13 +687,13 @@
     else
     {
-      if (calcNorm) *n    = G4ThreeVector(1,0,0) ;
-      return         snxt = 0 ;
-    }
-  }
-  else if (v.x() < 0) 
+      if (calcNorm) { *n   = G4ThreeVector(1,0,0) ; }
+      return        snxt = 0 ;
+    }
+  }
+  else if (v.x() < 0)
   {
     pdist = fDx + p.x() ;
 
-    if (pdist > kCarTolerance*0.5)
+    if (pdist > delta)
     {
       snxt = -pdist/v.x() ;
@@ -712,90 +702,91 @@
     else
     {
-      if (calcNorm) *n   = G4ThreeVector(-1,0,0) ;
+      if (calcNorm) { *n   = G4ThreeVector(-1,0,0) ; }
       return        snxt = 0 ;
     }
   }
-  else snxt = kInfinity ;
-
-  if ( v.y() > 0 )   // Y planes  
-  {
-    pdist=fDy-p.y();
-
-    if (pdist>kCarTolerance*0.5)
-    {
-      stmp=pdist/v.y();
-
-      if (stmp<snxt)
-      {
-        snxt=stmp;
-        side=kPY;
-      }
-    }
-    else
-    {
-      if (calcNorm) *n    = G4ThreeVector(0,1,0) ;
-      return         snxt = 0 ;
-    }
-  }
-  else if ( v.y() < 0 ) 
-  {
-    pdist = fDy + p.y() ;
-
-    if (pdist > kCarTolerance*0.5)
-    {
-      stmp=-pdist/v.y();
-
-      if (stmp<snxt)
-      {
-        snxt=stmp;
-        side=kMY;
-      }
-    }
-    else
-    {
-      if (calcNorm) *n    = G4ThreeVector(0,-1,0) ;
-      return         snxt = 0 ;
-    }
-  }
-  if (v.z()>0)        // Z planes 
-  {
-    pdist=fDz-p.z();
-
-    if (pdist > kCarTolerance*0.5)
-    {
-      stmp=pdist/v.z();
+
+  if (v.y() > 0)   // Y planes
+  {
+    pdist = fDy-p.y();
+
+    if (pdist > delta)
+    {
+      stmp = pdist/v.y();
 
       if (stmp < snxt)
       {
-        snxt=stmp;
-        side=kPZ;
+        snxt = stmp;
+        side = kPY;
       }
     }
     else
     {
-      if (calcNorm) *n    = G4ThreeVector(0,0,1) ;
-      return         snxt = 0 ;
-    }
-  }
-  else if (v.z()<0) 
-  {
-    pdist = fDz + p.z() ;
-
-    if (pdist > kCarTolerance*0.5)
-    {
-      stmp=-pdist/v.z();
-
-      if (stmp < snxt)
+      if (calcNorm) { *n   = G4ThreeVector(0,1,0) ; }
+      return        snxt = 0 ;
+    }
+  }
+  else if (v.y() < 0)
+  {
+    pdist = fDy + p.y() ;
+
+    if (pdist > delta)
+    {
+      stmp = -pdist/v.y();
+
+      if ( stmp < snxt )
       {
-        snxt=stmp;
-        side=kMZ;
+        snxt = stmp;
+        side = kMY;
       }
     }
     else
     {
-      if (calcNorm) *n    = G4ThreeVector(0,0,-1) ;
-      return         snxt = 0 ;
-    }
-  }
+      if (calcNorm) { *n   = G4ThreeVector(0,-1,0) ; }
+      return        snxt = 0 ;
+    }
+  }
+
+  if (v.z() > 0)        // Z planes
+  {
+    pdist = fDz-p.z();
+
+    if ( pdist > delta )
+    {
+      stmp = pdist/v.z();
+
+      if ( stmp < snxt )
+      {
+        snxt = stmp;
+        side = kPZ;
+      }
+    }
+    else
+    {
+      if (calcNorm) { *n   = G4ThreeVector(0,0,1) ; } 
+      return        snxt = 0 ;
+    }
+  }
+  else if (v.z() < 0)
+  {
+    pdist = fDz + p.z();
+
+    if ( pdist > delta )
+    {
+      stmp = -pdist/v.z();
+
+      if ( stmp < snxt )
+      {
+        snxt = stmp;
+        side = kMZ;
+      }
+    }
+    else
+    {
+      if (calcNorm) { *n   = G4ThreeVector(0,0,-1) ; }
+      return        snxt = 0 ;
+    }
+  }
+
   if (calcNorm)
   {      
@@ -875,12 +866,12 @@
   // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
 
-  if (safx2 < safx1) safe = safx2 ;
-  else               safe = safx1 ;
-  if (safy1 < safe)  safe = safy1 ;
-  if (safy2 < safe)  safe = safy2 ;
-  if (safz1 < safe)  safe = safz1 ;
-  if (safz2 < safe)  safe = safz2 ;
-
-  if (safe < 0) safe = 0 ;
+  if (safx2 < safx1) { safe = safx2; }
+  else               { safe = safx1; }
+  if (safy1 < safe)  { safe = safy1; }
+  if (safy2 < safe)  { safe = safy2; }
+  if (safz1 < safe)  { safe = safz1; }
+  if (safz2 < safe)  { safe = safz2; }
+
+  if (safe < 0) { safe = 0 ; }
   return safe ;  
 }
@@ -979,6 +970,6 @@
     py = -fDy +2*fDy*G4UniformRand();
 
-    if(G4UniformRand() > 0.5) pz =  fDz;
-    else                      pz = -fDz;
+    if(G4UniformRand() > 0.5) { pz =  fDz; }
+    else                      { pz = -fDz; }
   }
   else if ( ( select - Sxy ) < Sxz ) 
@@ -987,6 +978,6 @@
     pz = -fDz +2*fDz*G4UniformRand();
 
-    if(G4UniformRand() > 0.5) py =  fDy;
-    else                      py = -fDy;
+    if(G4UniformRand() > 0.5) { py =  fDy; }
+    else                      { py = -fDy; }
   }
   else  
@@ -995,6 +986,6 @@
     pz = -fDz +2*fDz*G4UniformRand();
 
-    if(G4UniformRand() > 0.5) px =  fDx;
-    else                      px = -fDx;
+    if(G4UniformRand() > 0.5) { px =  fDx; }
+    else                      { px = -fDx; }
   } 
   return G4ThreeVector(px,py,pz);
