Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4CutTubs.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/CSG/src/G4CutTubs.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4CutTubs.cc (Version 10.1.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4CutTubs implementation                    << 
 27 //                                                 26 //
                                                   >>  27 // $Id: G4CutTubs.cc 83572 2014-09-01 15:23:27Z gcosmo $
                                                   >>  28 //
                                                   >>  29 // 
                                                   >>  30 // class G4CutTubs
                                                   >>  31 //
                                                   >>  32 // History:
                                                   >>  33 //
                                                   >>  34 // 05.04.12 M.Kelsey   - GetPointOnSurface() throw flat in sqrt(r)
 28 // 01.06.11 T.Nikitina - Derived from G4Tubs       35 // 01.06.11 T.Nikitina - Derived from G4Tubs
 29 // 30.10.16 E.Tcherniaev - reimplemented Calcu <<  36 //
 30 //                         removed CreateRotat <<  37 /////////////////////////////////////////////////////////////////////////
 31 // ------------------------------------------- << 
 32                                                    38 
 33 #include "G4CutTubs.hh"                            39 #include "G4CutTubs.hh"
 34                                                    40 
 35 #if !defined(G4GEOM_USE_UCTUBS)                << 
 36                                                << 
 37 #include "G4GeomTools.hh"                      << 
 38 #include "G4VoxelLimits.hh"                        41 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    42 #include "G4AffineTransform.hh"
 40 #include "G4GeometryTolerance.hh"                  43 #include "G4GeometryTolerance.hh"
 41 #include "G4BoundingEnvelope.hh"               << 
 42                                                    44 
 43 #include "G4VPVParameterisation.hh"                45 #include "G4VPVParameterisation.hh"
 44 #include "G4QuickRand.hh"                      << 
 45                                                    46 
 46 #include "G4VGraphicsScene.hh"                 <<  47 #include "Randomize.hh"
 47 #include "G4Polyhedron.hh"                     << 
 48                                                    48 
 49 #include "G4AutoLock.hh"                       <<  49 #include "meshdefs.hh"
 50                                                    50 
 51 namespace                                      <<  51 #include "G4VGraphicsScene.hh"
 52 {                                              <<  52 #include "G4Polyhedron.hh"
 53   G4Mutex zminmaxMutex = G4MUTEX_INITIALIZER;  << 
 54 }                                              << 
 55                                                    53 
 56 using namespace CLHEP;                             54 using namespace CLHEP;
 57                                                    55 
 58 //////////////////////////////////////////////     56 /////////////////////////////////////////////////////////////////////////
 59 //                                                 57 //
 60 // Constructor - check parameters, convert ang     58 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
 61 //             - note if pdphi>2PI then reset      59 //             - note if pdphi>2PI then reset to 2PI
 62                                                    60 
 63 G4CutTubs::G4CutTubs( const G4String &pName,       61 G4CutTubs::G4CutTubs( const G4String &pName,
 64                       G4double pRMin, G4double     62                       G4double pRMin, G4double pRMax,
 65                       G4double pDz,                63                       G4double pDz,
 66                       G4double pSPhi, G4double     64                       G4double pSPhi, G4double pDPhi,
 67                       G4ThreeVector pLowNorm,G     65                       G4ThreeVector pLowNorm,G4ThreeVector pHighNorm )
 68   : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRM <<  66   : G4OTubs(pName, pRMin, pRMax, pDz, pSPhi, pDPhi),
 69     fSPhi(0.), fDPhi(0.), fZMin(0.), fZMax(0.) <<  67     fPhiFullCutTube(true)
 70 {                                                  68 {
 71   kRadTolerance = G4GeometryTolerance::GetInst     69   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
 72   kAngTolerance = G4GeometryTolerance::GetInst     70   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 73                                                    71 
 74   halfCarTolerance = kCarTolerance*0.5;            72   halfCarTolerance = kCarTolerance*0.5;
 75   halfRadTolerance = kRadTolerance*0.5;            73   halfRadTolerance = kRadTolerance*0.5;
 76   halfAngTolerance = kAngTolerance*0.5;            74   halfAngTolerance = kAngTolerance*0.5;
 77                                                    75 
 78   if (pDz<=0) // Check z-len                   << 
 79   {                                            << 
 80     std::ostringstream message;                << 
 81     message << "Negative Z half-length (" << p << 
 82     G4Exception("G4CutTubs::G4CutTubs()", "Geo << 
 83   }                                            << 
 84   if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch << 
 85   {                                            << 
 86     std::ostringstream message;                << 
 87     message << "Invalid values for radii in so << 
 88             << G4endl                          << 
 89             << "        pRMin = " << pRMin <<  << 
 90     G4Exception("G4CutTubs::G4CutTubs()", "Geo << 
 91   }                                            << 
 92                                                << 
 93   // Check angles                              << 
 94   //                                           << 
 95   CheckPhiAngles(pSPhi, pDPhi);                << 
 96                                                << 
 97   // Check on Cutted Planes Normals                76   // Check on Cutted Planes Normals
 98   // If there is NO CUT, propose to use G4Tubs     77   // If there is NO CUT, propose to use G4Tubs instead
 99   //                                               78   //
100   if ( ( pLowNorm.x() == 0.0) && ( pLowNorm.y( <<  79   if(pDPhi<twopi)  { fPhiFullCutTube=false; }
101     && ( pHighNorm.x() == 0.0) && (pHighNorm.y <<  80  
                                                   >>  81   if ( ( !pLowNorm.x()) && ( !pLowNorm.y())
                                                   >>  82     && ( !pHighNorm.x()) && (!pHighNorm.y()) )
102   {                                                83   {
103     std::ostringstream message;                    84     std::ostringstream message;
104     message << "Inexisting Low/High Normal to      85     message << "Inexisting Low/High Normal to Z plane or Parallel to Z."
105             << G4endl                              86             << G4endl
106             << "Normals to Z plane are " << pL <<  87             << "Normals to Z plane are (" << pLowNorm <<" and "
107             << pHighNorm << " in solid: " << G <<  88             << pHighNorm << ") in solid: " << GetName();
108     G4Exception("G4CutTubs::G4CutTubs()", "Geo     89     G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids1001",
109                 JustWarning, message, "Should      90                 JustWarning, message, "Should use G4Tubs!");
110   }                                                91   }
111                                                    92 
112   // If Normal is (0,0,0),means parallel to R,     93   // If Normal is (0,0,0),means parallel to R, give it value of (0,0,+/-1)
113   //                                           <<  94   // 
114   if (pLowNorm.mag2() == 0.)  { pLowNorm.setZ(     95   if (pLowNorm.mag2() == 0.)  { pLowNorm.setZ(-1.); }
115   if (pHighNorm.mag2()== 0.)  { pHighNorm.setZ     96   if (pHighNorm.mag2()== 0.)  { pHighNorm.setZ(1.); }
116                                                    97 
117   // Given Normals to Cut Planes have to be an <<  98   // Given Normals to Cut Planes have to be an unit vectors. 
118   // Normalize if it is needed.                    99   // Normalize if it is needed.
119   //                                              100   //
120   if (pLowNorm.mag2() != 1.)  { pLowNorm  = pL    101   if (pLowNorm.mag2() != 1.)  { pLowNorm  = pLowNorm.unit();  }
121   if (pHighNorm.mag2()!= 1.)  { pHighNorm = pH    102   if (pHighNorm.mag2()!= 1.)  { pHighNorm = pHighNorm.unit(); }
122                                                   103 
123   // Normals to cutted planes have to point ou    104   // Normals to cutted planes have to point outside Solid
124   //                                              105   //
125   if( (pLowNorm.mag2() != 0.) && (pHighNorm.ma    106   if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) )
126   {                                               107   {
127     if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z    108     if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.))
128     {                                             109     {
129       std::ostringstream message;                 110       std::ostringstream message;
130       message << "Invalid Low or High Normal t    111       message << "Invalid Low or High Normal to Z plane; "
131                  "has to point outside Solid."    112                  "has to point outside Solid." << G4endl
132               << "Invalid Norm to Z plane (" <    113               << "Invalid Norm to Z plane (" << pLowNorm << " or  "
133               << pHighNorm << ") in solid: " <    114               << pHighNorm << ") in solid: " << GetName();
134       G4Exception("G4CutTubs::G4CutTubs()", "G    115       G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
135                   FatalException, message);       116                   FatalException, message);
136     }                                             117     }
137   }                                               118   }
138   fLowNorm  = pLowNorm;                           119   fLowNorm  = pLowNorm;
139   fHighNorm = pHighNorm;                          120   fHighNorm = pHighNorm;
140                                                   121 
141   // Check intersection of cut planes, they MU << 122   // Check Intersection of Cutted planes. They MUST NOT Intersect
142   // each other inside the lateral surface     << 
143   //                                              123   //
144   if(IsCrossingCutPlanes())                       124   if(IsCrossingCutPlanes())
145   {                                               125   {
146     std::ostringstream message;                   126     std::ostringstream message;
147     message << "Invalid normals to Z plane in  << 127     message << "Invalid Low or High Normal to Z plane; "
148             << "Cut planes are crossing inside << 128             << "Crossing Cutted Planes." << G4endl
149             << " Solid type: G4CutTubs\n"      << 129             << "Invalid Norm to Z plane (" << pLowNorm << " and "
150             << " Parameters: \n"               << 130             << pHighNorm << ") in solid: " << GetName();
151             << "    inner radius : " << fRMin/ << 
152             << "    outer radius : " << fRMax/ << 
153             << "    half length Z: " << fDz/mm << 
154             << "    starting phi : " << fSPhi/ << 
155             << "    delta phi    : " << fDPhi/ << 
156             << "    low Norm     : " << fLowNo << 
157             << "    high Norm    : " << fHighN << 
158     G4Exception("G4CutTubs::G4CutTubs()", "Geo    131     G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
159                 FatalException, message);         132                 FatalException, message);
160   }                                               133   }
161 }                                                 134 }
162                                                   135 
163 //////////////////////////////////////////////    136 ///////////////////////////////////////////////////////////////////////
164 //                                                137 //
165 // Fake default constructor - sets only member    138 // Fake default constructor - sets only member data and allocates memory
166 //                            for usage restri    139 //                            for usage restricted to object persistency.
167 //                                                140 //
168 G4CutTubs::G4CutTubs( __void__& a )               141 G4CutTubs::G4CutTubs( __void__& a )
169   : G4CSGSolid(a)                              << 142   : G4OTubs(a), fLowNorm(G4ThreeVector()),
                                                   >> 143     fHighNorm(G4ThreeVector()), fPhiFullCutTube(false),
                                                   >> 144     halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
170 {                                                 145 {
171 }                                                 146 }
172                                                   147 
173 //////////////////////////////////////////////    148 //////////////////////////////////////////////////////////////////////////
174 //                                                149 //
175 // Destructor                                     150 // Destructor
176                                                   151 
177 G4CutTubs::~G4CutTubs() = default;             << 152 G4CutTubs::~G4CutTubs()
                                                   >> 153 {
                                                   >> 154 }
178                                                   155 
179 //////////////////////////////////////////////    156 //////////////////////////////////////////////////////////////////////////
180 //                                                157 //
181 // Copy constructor                               158 // Copy constructor
182                                                   159 
183 G4CutTubs::G4CutTubs(const G4CutTubs&) = defau << 160 G4CutTubs::G4CutTubs(const G4CutTubs& rhs)
                                                   >> 161   : G4OTubs(rhs), fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm),
                                                   >> 162     fPhiFullCutTube(rhs.fPhiFullCutTube),
                                                   >> 163     halfCarTolerance(rhs.halfCarTolerance),
                                                   >> 164     halfRadTolerance(rhs.halfRadTolerance),
                                                   >> 165     halfAngTolerance(rhs.halfAngTolerance)
                                                   >> 166 {
                                                   >> 167 }
184                                                   168 
185 //////////////////////////////////////////////    169 //////////////////////////////////////////////////////////////////////////
186 //                                                170 //
187 // Assignment operator                            171 // Assignment operator
188                                                   172 
189 G4CutTubs& G4CutTubs::operator = (const G4CutT << 173 G4CutTubs& G4CutTubs::operator = (const G4CutTubs& rhs) 
190 {                                                 174 {
191    // Check assignment to self                    175    // Check assignment to self
192    //                                             176    //
193    if (this == &rhs)  { return *this; }           177    if (this == &rhs)  { return *this; }
194                                                   178 
195    // Copy base class data                        179    // Copy base class data
196    //                                             180    //
197    G4CSGSolid::operator=(rhs);                 << 181    G4OTubs::operator=(rhs);
198                                                   182 
199    // Copy data                                   183    // Copy data
200    //                                             184    //
201    kRadTolerance = rhs.kRadTolerance; kAngTole << 185    fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm;
202    fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = << 
203    fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;       << 
204    fZMin = rhs.fZMin; fZMax = rhs.fZMax;       << 
205    sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh << 
206    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r << 
207    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh << 
208    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh << 
209    fPhiFullCutTube = rhs.fPhiFullCutTube;         186    fPhiFullCutTube = rhs.fPhiFullCutTube;
210    halfCarTolerance = rhs.halfCarTolerance;       187    halfCarTolerance = rhs.halfCarTolerance;
211    halfRadTolerance = rhs.halfRadTolerance;       188    halfRadTolerance = rhs.halfRadTolerance;
212    halfAngTolerance = rhs.halfAngTolerance;       189    halfAngTolerance = rhs.halfAngTolerance;
213    fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fH << 
214                                                   190 
215    return *this;                                  191    return *this;
216 }                                                 192 }
217                                                   193 
218 ////////////////////////////////////////////// << 194 ////////////////////////////////////////////////////////////////////////
219 //                                                195 //
220 // Get volume                                  << 196 // Calculate extent under transform and specified limit
221                                                   197 
222 G4double G4CutTubs::GetCubicVolume()           << 198 G4bool G4CutTubs::CalculateExtent( const EAxis              pAxis,
                                                   >> 199                                    const G4VoxelLimits&     pVoxelLimit,
                                                   >> 200                                    const G4AffineTransform& pTransform,
                                                   >> 201                                          G4double&          pMin, 
                                                   >> 202                                          G4double&          pMax    ) const
223 {                                                 203 {
224   constexpr G4int nphi = 200, nrho = 100;      << 204   if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
225                                                << 
226   if (fCubicVolume == 0.)                      << 
227   {                                               205   {
228     // get parameters                          << 206     // Special case handling for unrotated solid tubes
229     G4double rmin = GetInnerRadius();          << 207     // Compute x/y/z mins and maxs fro bounding box respecting limits,
230     G4double rmax = GetOuterRadius();          << 208     // with early returns if outside limits. Then switch() on pAxis,
231     G4double dz   = GetZHalfLength();          << 209     // and compute exact x and y limit for x/y case
232     G4double sphi = GetStartPhiAngle();        << 210       
233     G4double dphi = GetDeltaPhiAngle();        << 211     G4double xoffset, xMin, xMax;
                                                   >> 212     G4double yoffset, yMin, yMax;
                                                   >> 213     G4double zoffset, zMin, zMax;
                                                   >> 214 
                                                   >> 215     G4double diff1, diff2, maxDiff, newMin, newMax;
                                                   >> 216     G4double xoff1, xoff2, yoff1, yoff2, delta;
                                                   >> 217 
                                                   >> 218     xoffset = pTransform.NetTranslation().x();
                                                   >> 219     xMin = xoffset - fRMax;
                                                   >> 220     xMax = xoffset + fRMax;
234                                                   221 
235     // calculate volume                        << 222     if (pVoxelLimit.IsXLimited())
236     G4double volume = dz*dphi*(rmax*rmax - rmi << 
237     if (dphi < twopi) // make recalculation    << 
238     {                                             223     {
239       // set values for calculation of h - dis << 224       if ( (xMin > pVoxelLimit.GetMaxXExtent())
240       // opposite points on bases              << 225         || (xMax < pVoxelLimit.GetMinXExtent()) )
241       G4ThreeVector nbot = GetLowNorm();       << 226       {
242       G4ThreeVector ntop = GetHighNorm();      << 227         return false;
243       G4double nx = nbot.x()/nbot.z() - ntop.x << 228       }
244       G4double ny = nbot.y()/nbot.z() - ntop.y << 229       else
                                                   >> 230       {
                                                   >> 231         if (xMin < pVoxelLimit.GetMinXExtent())
                                                   >> 232         {
                                                   >> 233           xMin = pVoxelLimit.GetMinXExtent();
                                                   >> 234         }
                                                   >> 235         if (xMax > pVoxelLimit.GetMaxXExtent())
                                                   >> 236         {
                                                   >> 237           xMax = pVoxelLimit.GetMaxXExtent();
                                                   >> 238         }
                                                   >> 239       }
                                                   >> 240     }
                                                   >> 241     yoffset = pTransform.NetTranslation().y();
                                                   >> 242     yMin    = yoffset - fRMax;
                                                   >> 243     yMax    = yoffset + fRMax;
245                                                   244 
246       // compute volume by integration         << 245     if ( pVoxelLimit.IsYLimited() )
247       G4double delrho = (rmax - rmin)/nrho;    << 246     {
248       G4double delphi = dphi/nphi;             << 247       if ( (yMin > pVoxelLimit.GetMaxYExtent())
249       volume = 0.;                             << 248         || (yMax < pVoxelLimit.GetMinYExtent()) )
250       for (G4int irho=0; irho<nrho; ++irho)    << 
251       {                                           249       {
252         G4double r1  = rmin + delrho*irho;     << 250         return false;
253         G4double r2  = rmin + delrho*(irho + 1 << 251       }
254         G4double rho = 0.5*(r1 + r2);          << 252       else
255         G4double sector = 0.5*delphi*(r2*r2 -  << 253       {
256         for (G4int iphi=0; iphi<nphi; ++iphi)  << 254         if (yMin < pVoxelLimit.GetMinYExtent())
                                                   >> 255         {
                                                   >> 256           yMin = pVoxelLimit.GetMinYExtent();
                                                   >> 257         }
                                                   >> 258         if (yMax > pVoxelLimit.GetMaxYExtent())
257         {                                         259         {
258           G4double phi = sphi + delphi*(iphi + << 260           yMax=pVoxelLimit.GetMaxYExtent();
259           G4double h = nx*rho*std::cos(phi) +  << 
260           volume += sector*h;                  << 
261         }                                         261         }
262       }                                           262       }
263     }                                             263     }
264     fCubicVolume = volume;                     << 264     zoffset = pTransform.NetTranslation().z();
265   }                                            << 265     GetMaxMinZ(zMin,zMax);
266   return fCubicVolume;                         << 266     zMin    += zoffset;
267 }                                              << 267     zMax    += zoffset;
268                                                   268 
269 ////////////////////////////////////////////// << 269     if ( pVoxelLimit.IsZLimited() )
270 //                                             << 270     {
271 // Get surface area                            << 271       if ( (zMin > pVoxelLimit.GetMaxZExtent())
                                                   >> 272         || (zMax < pVoxelLimit.GetMinZExtent()) )
                                                   >> 273       {
                                                   >> 274         return false;
                                                   >> 275       }
                                                   >> 276       else
                                                   >> 277       {
                                                   >> 278         if (zMin < pVoxelLimit.GetMinZExtent())
                                                   >> 279         {
                                                   >> 280           zMin = pVoxelLimit.GetMinZExtent();
                                                   >> 281         }
                                                   >> 282         if (zMax > pVoxelLimit.GetMaxZExtent())
                                                   >> 283         {
                                                   >> 284           zMax = pVoxelLimit.GetMaxZExtent();
                                                   >> 285         }
                                                   >> 286       }
                                                   >> 287     }
                                                   >> 288     switch ( pAxis )  // Known to cut cylinder
                                                   >> 289     {
                                                   >> 290       case kXAxis :
                                                   >> 291       {
                                                   >> 292         yoff1 = yoffset - yMin;
                                                   >> 293         yoff2 = yMax    - yoffset;
272                                                   294 
273 G4double G4CutTubs::GetSurfaceArea()           << 295         if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
274 {                                              << 296         {                                   // => no change
275   constexpr G4int nphi = 400;                  << 297           pMin = xMin;
                                                   >> 298           pMax = xMax;
                                                   >> 299         }
                                                   >> 300         else
                                                   >> 301         {
                                                   >> 302           // Y limits don't cross max/min x => compute max delta x,
                                                   >> 303           // hence new mins/maxs
276                                                   304 
277   if (fSurfaceArea == 0.)                      << 305           delta   = fRMax*fRMax - yoff1*yoff1;
278   {                                            << 306           diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
279     // get parameters                          << 307           delta   = fRMax*fRMax - yoff2*yoff2;
280     G4double rmin = GetInnerRadius();          << 308           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
281     G4double rmax = GetOuterRadius();          << 309           maxDiff = (diff1 > diff2) ? diff1:diff2;
282     G4double dz   = GetZHalfLength();          << 310           newMin  = xoffset - maxDiff;
283     G4double sphi = GetStartPhiAngle();        << 311           newMax  = xoffset + maxDiff;
284     G4double dphi = GetDeltaPhiAngle();        << 312           pMin    = (newMin < xMin) ? xMin : newMin;
285     G4ThreeVector nbot = GetLowNorm();         << 313           pMax    = (newMax > xMax) ? xMax : newMax;
286     G4ThreeVector ntop = GetHighNorm();        << 314         }    
287                                                << 315         break;
288     // calculate lateral surface area          << 316       }
289     G4double sinner = 2.*dz*dphi*rmin;         << 317       case kYAxis :
290     G4double souter = 2.*dz*dphi*rmax;         << 318       {
291     if (dphi < twopi) // make recalculation    << 319         xoff1 = xoffset - xMin;
292     {                                          << 320         xoff2 = xMax - xoffset;
293       // set values for calculation of h - dis << 
294       // opposite points on bases              << 
295       G4double nx = nbot.x()/nbot.z() - ntop.x << 
296       G4double ny = nbot.y()/nbot.z() - ntop.y << 
297                                                << 
298       // compute lateral surface area by integ << 
299       G4double delphi = dphi/nphi;             << 
300       sinner = 0.;                             << 
301       souter = 0.;                             << 
302       for (G4int iphi=0; iphi<nphi; ++iphi)    << 
303       {                                        << 
304         G4double phi = sphi + delphi*(iphi + 0 << 
305         G4double cosphi = std::cos(phi);       << 
306         G4double sinphi = std::sin(phi);       << 
307         sinner += rmin*(nx*cosphi + ny*sinphi) << 
308         souter += rmax*(nx*cosphi + ny*sinphi) << 
309       }                                        << 
310       sinner *= delphi*rmin;                   << 
311       souter *= delphi*rmax;                   << 
312     }                                          << 
313     // set surface area                        << 
314     G4double scut  = (dphi == twopi) ? 0. : 2. << 
315     G4double szero = 0.5*dphi*(rmax*rmax - rmi << 
316     G4double slow  = szero/std::abs(nbot.z()); << 
317     G4double shigh = szero/std::abs(ntop.z()); << 
318     fSurfaceArea = sinner + souter + 2.*scut + << 
319   }                                            << 
320   return fSurfaceArea;                         << 
321 }                                              << 
322                                                   321 
323 ////////////////////////////////////////////// << 322         if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
324 //                                             << 323         {                                   // => no change
325 // Get bounding box                            << 324           pMin = yMin;
                                                   >> 325           pMax = yMax;
                                                   >> 326         }
                                                   >> 327         else
                                                   >> 328         {
                                                   >> 329           // X limits don't cross max/min y => compute max delta y,
                                                   >> 330           // hence new mins/maxs
326                                                   331 
327 void G4CutTubs::BoundingLimits(G4ThreeVector&  << 332           delta   = fRMax*fRMax - xoff1*xoff1;
328 {                                              << 333           diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
329   G4double rmin = GetInnerRadius();            << 334           delta   = fRMax*fRMax - xoff2*xoff2;
330   G4double rmax = GetOuterRadius();            << 335           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
331   G4double dz   = GetZHalfLength();            << 336           maxDiff = (diff1 > diff2) ? diff1 : diff2;
332   G4double dphi = GetDeltaPhiAngle();          << 337           newMin  = yoffset - maxDiff;
333                                                << 338           newMax  = yoffset + maxDiff;
334   G4double sinSphi = GetSinStartPhi();         << 339           pMin    = (newMin < yMin) ? yMin : newMin;
335   G4double cosSphi = GetCosStartPhi();         << 340           pMax    = (newMax > yMax) ? yMax : newMax;
336   G4double sinEphi = GetSinEndPhi();           << 341         }
337   G4double cosEphi = GetCosEndPhi();           << 342         break;
338                                                << 343       }
339   G4ThreeVector norm;                          << 344       case kZAxis:
340   G4double mag, topx, topy, dists, diste;      << 345       {
341   G4bool iftop;                                << 346         pMin = zMin;
342                                                << 347         pMax = zMax;
343   // Find Zmin                                 << 348         break;
344   //                                           << 349       }
345   G4double zmin;                               << 350       default:
346   norm = GetLowNorm();                         << 351         break;
347   mag  = std::sqrt(norm.x()*norm.x() + norm.y( << 352     }
348   topx = (mag == 0) ? 0 : -rmax*norm.x()/mag;  << 353     pMin -= kCarTolerance;
349   topy = (mag == 0) ? 0 : -rmax*norm.y()/mag;  << 354     pMax += kCarTolerance;
350   dists =  sinSphi*topx - cosSphi*topy;        << 355     return true;
351   diste = -sinEphi*topx + cosEphi*topy;        << 
352   if (dphi > pi)                               << 
353   {                                            << 
354     iftop = true;                              << 
355     if (dists > 0 && diste > 0)iftop = false;  << 
356   }                                            << 
357   else                                         << 
358   {                                            << 
359     iftop = false;                             << 
360     if (dists <= 0 && diste <= 0) iftop = true << 
361   }                                            << 
362   if (iftop)                                   << 
363   {                                            << 
364     zmin = -(norm.x()*topx + norm.y()*topy)/no << 
365   }                                            << 
366   else                                         << 
367   {                                            << 
368     G4double z1 = -rmin*(norm.x()*cosSphi + no << 
369     G4double z2 = -rmin*(norm.x()*cosEphi + no << 
370     G4double z3 = -rmax*(norm.x()*cosSphi + no << 
371     G4double z4 = -rmax*(norm.x()*cosEphi + no << 
372     zmin = std::min(std::min(std::min(z1,z2),z << 
373   }                                            << 
374                                                << 
375   // Find Zmax                                 << 
376   //                                           << 
377   G4double zmax;                               << 
378   norm = GetHighNorm();                        << 
379   mag  = std::sqrt(norm.x()*norm.x() + norm.y( << 
380   topx = (mag == 0) ? 0 : -rmax*norm.x()/mag;  << 
381   topy = (mag == 0) ? 0 : -rmax*norm.y()/mag;  << 
382   dists =  sinSphi*topx - cosSphi*topy;        << 
383   diste = -sinEphi*topx + cosEphi*topy;        << 
384   if (dphi > pi)                               << 
385   {                                            << 
386     iftop = true;                              << 
387     if (dists > 0 && diste > 0) iftop = false; << 
388   }                                            << 
389   else                                         << 
390   {                                            << 
391     iftop = false;                             << 
392     if (dists <= 0 && diste <= 0) iftop = true << 
393   }                                            << 
394   if (iftop)                                   << 
395   {                                            << 
396     zmax = -(norm.x()*topx + norm.y()*topy)/no << 
397   }                                            << 
398   else                                         << 
399   {                                            << 
400     G4double z1 = -rmin*(norm.x()*cosSphi + no << 
401     G4double z2 = -rmin*(norm.x()*cosEphi + no << 
402     G4double z3 = -rmax*(norm.x()*cosSphi + no << 
403     G4double z4 = -rmax*(norm.x()*cosEphi + no << 
404     zmax = std::max(std::max(std::max(z1,z2),z << 
405   }                                            << 
406                                                << 
407   // Find bounding box                         << 
408   //                                           << 
409   if (dphi < twopi)                            << 
410   {                                            << 
411     G4TwoVector vmin,vmax;                     << 
412     G4GeomTools::DiskExtent(rmin,rmax,         << 
413                             GetSinStartPhi(),G << 
414                             GetSinEndPhi(),Get << 
415                             vmin,vmax);        << 
416     pMin.set(vmin.x(),vmin.y(), zmin);         << 
417     pMax.set(vmax.x(),vmax.y(), zmax);         << 
418   }                                            << 
419   else                                         << 
420   {                                            << 
421     pMin.set(-rmax,-rmax, zmin);               << 
422     pMax.set( rmax, rmax, zmax);               << 
423   }                                               356   }
424                                                << 357   else // Calculate rotated vertex coordinates
425   // Check correctness of the bounding box     << 
426   //                                           << 
427   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
428   {                                               358   {
429     std::ostringstream message;                << 359     G4int i, noEntries, noBetweenSections4;
430     message << "Bad bounding box (min >= max)  << 360     G4bool existsAfterClip = false;
431             << GetName() << " !"               << 361     G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
432             << "\npMin = " << pMin             << 
433             << "\npMax = " << pMax;            << 
434     G4Exception("G4CutTubs::BoundingLimits()", << 
435                 JustWarning, message);         << 
436     DumpInfo();                                << 
437   }                                            << 
438 }                                              << 
439                                                << 
440 ////////////////////////////////////////////// << 
441 //                                             << 
442 // Calculate extent under transform and specif << 
443                                                << 
444 G4bool G4CutTubs::CalculateExtent( const EAxis << 
445                                    const G4Vox << 
446                                    const G4Aff << 
447                                          G4dou << 
448                                          G4dou << 
449 {                                              << 
450   G4ThreeVector bmin, bmax;                    << 
451   G4bool exist;                                << 
452                                                   362 
453   // Get bounding box                          << 363     pMin =  kInfinity;
454   BoundingLimits(bmin,bmax);                   << 364     pMax = -kInfinity;
455                                                   365 
456   // Check bounding box                        << 366     noEntries = vertices->size();
457   G4BoundingEnvelope bbox(bmin,bmax);          << 367     noBetweenSections4 = noEntries - 4;
458 #ifdef G4BBOX_EXTENT                           << 368     
459   return bbox.CalculateExtent(pAxis,pVoxelLimi << 369     for ( i = 0 ; i < noEntries ; i += 4 )
460 #endif                                         << 370     {
461   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 371       ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
462   {                                            << 372     }
463     return exist = pMin < pMax;                << 373     for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
464   }                                            << 374     {
465                                                << 375       ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
466   // Get parameters of the solid               << 376     }
467   G4double rmin = GetInnerRadius();            << 377     if ( (pMin != kInfinity) || (pMax != -kInfinity) )
468   G4double rmax = GetOuterRadius();            << 378     {
469   G4double dphi = GetDeltaPhiAngle();          << 379       existsAfterClip = true;
470   G4double zmin = bmin.z();                    << 380       pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
471   G4double zmax = bmax.z();                    << 381       pMax += kCarTolerance;
472                                                << 382     }
473   // Find bounding envelope and calculate exte << 383     else
474   //                                           << 384     {
475   const G4int NSTEPS = 24;            // numbe << 385       // Check for case where completely enveloping clipping volume
476   G4double astep  = twopi/NSTEPS;     // max a << 386       // If point inside then we are confident that the solid completely
477   G4int    ksteps = (dphi <= astep) ? 1 : (G4i << 387       // envelopes the clipping volume. Hence set min/max extents according
478   G4double ang    = dphi/ksteps;               << 388       // to clipping volume extents along the specified axis.
479                                                << 389 
480   G4double sinHalf = std::sin(0.5*ang);        << 390       G4ThreeVector clipCentre(
481   G4double cosHalf = std::cos(0.5*ang);        << 391              (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
482   G4double sinStep = 2.*sinHalf*cosHalf;       << 392              (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
483   G4double cosStep = 1. - 2.*sinHalf*sinHalf;  << 393              (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
484   G4double rext    = rmax/cosHalf;             << 394         
485                                                << 395       if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
486   // bounding envelope for full cylinder consi << 396       {
487   // in other cases it is a sequence of quadri << 397         existsAfterClip = true;
488   if (rmin == 0 && dphi == twopi)              << 398         pMin            = pVoxelLimit.GetMinExtent(pAxis);
489   {                                            << 399         pMax            = pVoxelLimit.GetMaxExtent(pAxis);
490     G4double sinCur = sinHalf;                 << 400       }
491     G4double cosCur = cosHalf;                 << 401     }
492                                                << 402     delete vertices;
493     G4ThreeVectorList baseA(NSTEPS),baseB(NSTE << 403     return existsAfterClip;
494     for (G4int k=0; k<NSTEPS; ++k)             << 
495     {                                          << 
496       baseA[k].set(rext*cosCur,rext*sinCur,zmi << 
497       baseB[k].set(rext*cosCur,rext*sinCur,zma << 
498                                                << 
499       G4double sinTmp = sinCur;                << 
500       sinCur = sinCur*cosStep + cosCur*sinStep << 
501       cosCur = cosCur*cosStep - sinTmp*sinStep << 
502     }                                          << 
503     std::vector<const G4ThreeVectorList *> pol << 
504     polygons[0] = &baseA;                      << 
505     polygons[1] = &baseB;                      << 
506     G4BoundingEnvelope benv(bmin,bmax,polygons << 
507     exist = benv.CalculateExtent(pAxis,pVoxelL << 
508   }                                            << 
509   else                                         << 
510   {                                            << 
511     G4double sinStart = GetSinStartPhi();      << 
512     G4double cosStart = GetCosStartPhi();      << 
513     G4double sinEnd   = GetSinEndPhi();        << 
514     G4double cosEnd   = GetCosEndPhi();        << 
515     G4double sinCur   = sinStart*cosHalf + cos << 
516     G4double cosCur   = cosStart*cosHalf - sin << 
517                                                << 
518     // set quadrilaterals                      << 
519     G4ThreeVectorList pols[NSTEPS+2];          << 
520     for (G4int k=0; k<ksteps+2; ++k) pols[k].r << 
521     pols[0][0].set(rmin*cosStart,rmin*sinStart << 
522     pols[0][1].set(rmin*cosStart,rmin*sinStart << 
523     pols[0][2].set(rmax*cosStart,rmax*sinStart << 
524     pols[0][3].set(rmax*cosStart,rmax*sinStart << 
525     for (G4int k=1; k<ksteps+1; ++k)           << 
526     {                                          << 
527       pols[k][0].set(rmin*cosCur,rmin*sinCur,z << 
528       pols[k][1].set(rmin*cosCur,rmin*sinCur,z << 
529       pols[k][2].set(rext*cosCur,rext*sinCur,z << 
530       pols[k][3].set(rext*cosCur,rext*sinCur,z << 
531                                                << 
532       G4double sinTmp = sinCur;                << 
533       sinCur = sinCur*cosStep + cosCur*sinStep << 
534       cosCur = cosCur*cosStep - sinTmp*sinStep << 
535     }                                          << 
536     pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin << 
537     pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin << 
538     pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin << 
539     pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin << 
540                                                << 
541     // set envelope and calculate extent       << 
542     std::vector<const G4ThreeVectorList *> pol << 
543     polygons.resize(ksteps+2);                 << 
544     for (G4int k=0; k<ksteps+2; ++k) { polygon << 
545     G4BoundingEnvelope benv(bmin,bmax,polygons << 
546     exist = benv.CalculateExtent(pAxis,pVoxelL << 
547   }                                               404   }
548   return exist;                                << 
549 }                                                 405 }
550                                                   406 
551 ////////////////////////////////////////////// << 407 ///////////////////////////////////////////////////////////////////////////
552 //                                                408 //
553 // Return whether point inside/outside/on surf    409 // Return whether point inside/outside/on surface
554                                                   410 
555 EInside G4CutTubs::Inside( const G4ThreeVector    411 EInside G4CutTubs::Inside( const G4ThreeVector& p ) const
556 {                                                 412 {
557   G4ThreeVector vZ = G4ThreeVector(0,0,fDz);   << 413   G4double zinLow,zinHigh,r2,pPhi=0.;
                                                   >> 414   G4double tolRMin,tolRMax;
                                                   >> 415   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
558   EInside in = kInside;                           416   EInside in = kInside;
559                                                   417 
                                                   >> 418   // Check if point is contained in the cut plane in -/+ Z
                                                   >> 419 
560   // Check the lower cut plane                    420   // Check the lower cut plane
561   //                                              421   //
562   G4double zinLow =(p+vZ).dot(fLowNorm);       << 422   zinLow =(p+vZ).dot(fLowNorm);
563   if (zinLow > halfCarTolerance)  { return kOu << 423   if (zinLow>halfCarTolerance)  { return kOutside; }
564                                                   424 
565   // Check the higher cut plane                   425   // Check the higher cut plane
566   //                                              426   //
567   G4double zinHigh = (p-vZ).dot(fHighNorm);    << 427   zinHigh = (p-vZ).dot(fHighNorm);
568   if (zinHigh > halfCarTolerance)  { return kO << 428   if (zinHigh>halfCarTolerance)  { return kOutside; }
569                                                   429 
570   // Check radius                                 430   // Check radius
571   //                                              431   //
572   G4double r2 = p.x()*p.x() + p.y()*p.y() ;    << 432   r2 = p.x()*p.x() + p.y()*p.y() ;
573                                                   433 
574   G4double tolRMin = fRMin - halfRadTolerance; << 434   // First check 'generous' boundaries R+tolerance
575   G4double tolRMax = fRMax + halfRadTolerance; << 435   //
                                                   >> 436   tolRMin = fRMin - halfRadTolerance ;
                                                   >> 437   tolRMax = fRMax + halfRadTolerance ;
576   if ( tolRMin < 0 )  { tolRMin = 0; }            438   if ( tolRMin < 0 )  { tolRMin = 0; }
577                                                << 439     
578   if (r2 < tolRMin*tolRMin || r2 > tolRMax*tol << 440   if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax))
579                                                << 441      && (r2 >=halfRadTolerance*halfRadTolerance) )  { return kOutside; }
580   // Check Phi cut                             << 442    
                                                   >> 443   // Check Phi
581   //                                              444   //
582   if(!fPhiFullCutTube)                            445   if(!fPhiFullCutTube)
583   {                                               446   {
584     if ((tolRMin == 0) && (std::fabs(p.x()) <= << 447     // Try outer tolerant phi boundaries only
585                        && (std::fabs(p.y()) <= << 448 
                                                   >> 449     if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
                                                   >> 450                       && (std::fabs(p.y())<=halfCarTolerance) )
586     {                                             451     {
587       return kSurface;                            452       return kSurface;
588     }                                             453     }
                                                   >> 454  
                                                   >> 455     pPhi = std::atan2(p.y(),p.x()) ;
589                                                   456 
590     G4double phi0 = std::atan2(p.y(),p.x());   << 457     if ( pPhi < -halfAngTolerance)  { pPhi += twopi; } // 0<=pPhi<2pi
591     G4double phi1 = phi0 - twopi;              << 458     if ( fSPhi >= 0 )
592     G4double phi2 = phi0 + twopi;              << 459     {
593                                                << 460       if ( (std::fabs(pPhi) < halfAngTolerance)
594     in = kOutside;                             << 461         && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
595     G4double sphi = fSPhi - halfAngTolerance;  << 462       { 
596     G4double ephi = sphi + fDPhi + kAngToleran << 463         pPhi += twopi ; // 0 <= pPhi < 2pi
597     if ((phi0  >= sphi && phi0  <= ephi) ||    << 464       }
598         (phi1  >= sphi && phi1  <= ephi) ||    << 465       if ( (pPhi <= fSPhi - halfAngTolerance)
599         (phi2  >= sphi && phi2  <= ephi)) in = << 466         || (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
600     if (in == kOutside)  { return kOutside; }  << 467       {
601                                                << 468         in = kOutside ;
602     sphi += kAngTolerance;                     << 469       }
603     ephi -= kAngTolerance;                     << 470       else if ( (pPhi <= fSPhi + halfAngTolerance)
604     if ((phi0  >= sphi && phi0  <= ephi) ||    << 471              || (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
605         (phi1  >= sphi && phi1  <= ephi) ||    << 472       {
606         (phi2  >= sphi && phi2  <= ephi)) in = << 473         in=kSurface;
607     if (in == kSurface)  { return kSurface; }  << 474       }
                                                   >> 475     }
                                                   >> 476     else  // fSPhi < 0
                                                   >> 477     {
                                                   >> 478       if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
                                                   >> 479         && (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
                                                   >> 480       {
                                                   >> 481         in = kOutside;
                                                   >> 482       }
                                                   >> 483       else
                                                   >> 484       {
                                                   >> 485         in = kSurface ;
                                                   >> 486       }
                                                   >> 487     }
608   }                                               488   }
609                                                   489 
610   // Check on the Surface for Z                   490   // Check on the Surface for Z
611   //                                              491   //
612   if ((zinLow >= -halfCarTolerance) || (zinHig << 492   if ((zinLow>=-halfCarTolerance)
                                                   >> 493    || (zinHigh>=-halfCarTolerance))
613   {                                               494   {
614     return kSurface;                           << 495     in=kSurface;
615   }                                               496   }
616                                                   497 
617   // Check on the Surface for R                   498   // Check on the Surface for R
618   //                                              499   //
619   if (fRMin != 0.0) { tolRMin = fRMin + halfRa << 500   if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
620   else       { tolRMin = 0; }                  << 501   else       { tolRMin = 0 ; }
621   tolRMax = fRMax - halfRadTolerance;          << 502   tolRMax = fRMax - halfRadTolerance ;  
622   if (((r2 <= tolRMin*tolRMin) || (r2 >= tolRM << 503   if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&&
623        (r2 >= halfRadTolerance*halfRadToleranc << 504         (r2 >=halfRadTolerance*halfRadTolerance) )
624   {                                               505   {
625     return kSurface;                              506     return kSurface;
626   }                                               507   }
627                                                   508 
628   return in;                                      509   return in;
629 }                                                 510 }
630                                                   511 
631 //////////////////////////////////////////////    512 ///////////////////////////////////////////////////////////////////////////
632 //                                                513 //
633 // Return unit normal of surface closest to p     514 // Return unit normal of surface closest to p
634 // - note if point on z axis, ignore phi divid    515 // - note if point on z axis, ignore phi divided sides
635 // - unsafe if point close to z axis a rmin=0     516 // - unsafe if point close to z axis a rmin=0 - no explicit checks
636                                                   517 
637 G4ThreeVector G4CutTubs::SurfaceNormal( const     518 G4ThreeVector G4CutTubs::SurfaceNormal( const G4ThreeVector& p ) const
638 {                                                 519 {
639   G4int noSurfaces = 0;                           520   G4int noSurfaces = 0;
640   G4double rho, pPhi;                             521   G4double rho, pPhi;
641   G4double distZLow,distZHigh, distRMin, distR    522   G4double distZLow,distZHigh, distRMin, distRMax;
642   G4double distSPhi = kInfinity, distEPhi = kI    523   G4double distSPhi = kInfinity, distEPhi = kInfinity;
643   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);        524   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
644                                                   525 
645   G4ThreeVector norm, sumnorm(0.,0.,0.);          526   G4ThreeVector norm, sumnorm(0.,0.,0.);
646   G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);    527   G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
647   G4ThreeVector nR, nPs, nPe;                     528   G4ThreeVector nR, nPs, nPe;
648                                                   529 
649   rho = std::sqrt(p.x()*p.x() + p.y()*p.y());     530   rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
650                                                   531 
651   distRMin = std::fabs(rho - fRMin);              532   distRMin = std::fabs(rho - fRMin);
652   distRMax = std::fabs(rho - fRMax);              533   distRMax = std::fabs(rho - fRMax);
653                                                   534 
654   // dist to Low Cut                              535   // dist to Low Cut
655   //                                              536   //
656   distZLow =std::fabs((p+vZ).dot(fLowNorm));      537   distZLow =std::fabs((p+vZ).dot(fLowNorm));
657                                                << 538  
658   // dist to High Cut                             539   // dist to High Cut
659   //                                              540   //
660   distZHigh = std::fabs((p-vZ).dot(fHighNorm))    541   distZHigh = std::fabs((p-vZ).dot(fHighNorm));
661                                                   542 
662   if (!fPhiFullCutTube)    // Protected agains << 543   if (!fPhiFullCutTube)    // Protected against (0,0,z) 
663   {                                               544   {
664     if ( rho > halfCarTolerance )                 545     if ( rho > halfCarTolerance )
665     {                                             546     {
666       pPhi = std::atan2(p.y(),p.x());             547       pPhi = std::atan2(p.y(),p.x());
667                                                << 548     
668       if(pPhi  < fSPhi- halfCarTolerance)         549       if(pPhi  < fSPhi- halfCarTolerance)           { pPhi += twopi; }
669       else if(pPhi > fSPhi+fDPhi+ halfCarToler    550       else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
670                                                   551 
671       distSPhi = std::fabs(pPhi - fSPhi);      << 552       distSPhi = std::fabs(pPhi - fSPhi);       
672       distEPhi = std::fabs(pPhi - fSPhi - fDPh << 553       distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 
673     }                                             554     }
674     else if( fRMin == 0.0 )                    << 555     else if( !fRMin )
675     {                                             556     {
676       distSPhi = 0.;                           << 557       distSPhi = 0.; 
677       distEPhi = 0.;                           << 558       distEPhi = 0.; 
678     }                                             559     }
679     nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0  << 560     nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
680     nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0  << 561     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
681   }                                               562   }
682   if ( rho > halfCarTolerance ) { nR = G4Three    563   if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
683                                                   564 
684   if( distRMax <= halfCarTolerance )           << 565   if( distRMax <= halfCarTolerance ) 
685   {                                               566   {
686     ++noSurfaces;                              << 567     noSurfaces ++;
687     sumnorm += nR;                                568     sumnorm += nR;
688   }                                               569   }
689   if( (fRMin != 0.0) && (distRMin <= halfCarTo << 570   if( fRMin && (distRMin <= halfCarTolerance) )
690   {                                               571   {
691     ++noSurfaces;                              << 572     noSurfaces ++;
692     sumnorm -= nR;                                573     sumnorm -= nR;
693   }                                               574   }
694   if( fDPhi < twopi )                          << 575   if( fDPhi < twopi )   
695   {                                               576   {
696     if (distSPhi <= halfAngTolerance)          << 577     if (distSPhi <= halfAngTolerance)  
697     {                                             578     {
698       ++noSurfaces;                            << 579       noSurfaces ++;
699       sumnorm += nPs;                             580       sumnorm += nPs;
700     }                                             581     }
701     if (distEPhi <= halfAngTolerance)          << 582     if (distEPhi <= halfAngTolerance)  
702     {                                             583     {
703       ++noSurfaces;                            << 584       noSurfaces ++;
704       sumnorm += nPe;                             585       sumnorm += nPe;
705     }                                             586     }
706   }                                               587   }
707   if (distZLow <= halfCarTolerance)            << 588   if (distZLow <= halfCarTolerance)  
708   {                                               589   {
709     ++noSurfaces;                              << 590     noSurfaces ++;
710     sumnorm += fLowNorm;                          591     sumnorm += fLowNorm;
711   }                                               592   }
712   if (distZHigh <= halfCarTolerance)           << 593   if (distZHigh <= halfCarTolerance)  
713   {                                               594   {
714     ++noSurfaces;                              << 595     noSurfaces ++;
715     sumnorm += fHighNorm;                         596     sumnorm += fHighNorm;
716   }                                               597   }
717   if ( noSurfaces == 0 )                          598   if ( noSurfaces == 0 )
718   {                                               599   {
719 #ifdef G4CSGDEBUG                                 600 #ifdef G4CSGDEBUG
720     G4Exception("G4CutTubs::SurfaceNormal(p)",    601     G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002",
721                 JustWarning, "Point p is not o    602                 JustWarning, "Point p is not on surface !?" );
722     G4int oldprc = G4cout.precision(20);          603     G4int oldprc = G4cout.precision(20);
723     G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<    604     G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
724           << G4endl << G4endl;                    605           << G4endl << G4endl;
725     G4cout.precision(oldprc) ;                    606     G4cout.precision(oldprc) ;
726 #endif                                         << 607 #endif 
727      norm = ApproxSurfaceNormal(p);               608      norm = ApproxSurfaceNormal(p);
728   }                                               609   }
729   else if ( noSurfaces == 1 )  { norm = sumnor    610   else if ( noSurfaces == 1 )  { norm = sumnorm; }
730   else                         { norm = sumnor    611   else                         { norm = sumnorm.unit(); }
731                                                   612 
732   return norm;                                    613   return norm;
733 }                                                 614 }
734                                                   615 
735 //////////////////////////////////////////////    616 /////////////////////////////////////////////////////////////////////////////
736 //                                                617 //
737 // Algorithm for SurfaceNormal() following the    618 // Algorithm for SurfaceNormal() following the original specification
738 // for points not on the surface                  619 // for points not on the surface
739                                                   620 
740 G4ThreeVector G4CutTubs::ApproxSurfaceNormal(     621 G4ThreeVector G4CutTubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
741 {                                                 622 {
742   enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ} << 
743                                                << 
744   ENorm side ;                                    623   ENorm side ;
745   G4ThreeVector norm ;                            624   G4ThreeVector norm ;
746   G4double rho, phi ;                             625   G4double rho, phi ;
747   G4double distZLow,distZHigh,distZ;              626   G4double distZLow,distZHigh,distZ;
748   G4double distRMin, distRMax, distSPhi, distE    627   G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
749   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);        628   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
750                                                   629 
751   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;    630   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
752                                                   631 
753   distRMin = std::fabs(rho - fRMin) ;             632   distRMin = std::fabs(rho - fRMin) ;
754   distRMax = std::fabs(rho - fRMax) ;             633   distRMax = std::fabs(rho - fRMax) ;
755                                                   634 
756   //dist to Low Cut                               635   //dist to Low Cut
757   //                                              636   //
758   distZLow =std::fabs((p+vZ).dot(fLowNorm));      637   distZLow =std::fabs((p+vZ).dot(fLowNorm));
759                                                   638 
760   //dist to High Cut                              639   //dist to High Cut
761   //                                              640   //
762   distZHigh = std::fabs((p-vZ).dot(fHighNorm))    641   distZHigh = std::fabs((p-vZ).dot(fHighNorm));
763   distZ=std::min(distZLow,distZHigh);             642   distZ=std::min(distZLow,distZHigh);
764                                                   643 
765   if (distRMin < distRMax) // First minimum       644   if (distRMin < distRMax) // First minimum
766   {                                               645   {
767     if ( distZ < distRMin )                       646     if ( distZ < distRMin )
768     {                                             647     {
769        distMin = distZ ;                          648        distMin = distZ ;
770        side    = kNZ ;                            649        side    = kNZ ;
771     }                                             650     }
772     else                                          651     else
773     {                                             652     {
774       distMin = distRMin ;                        653       distMin = distRMin ;
775       side    = kNRMin   ;                        654       side    = kNRMin   ;
776     }                                             655     }
777   }                                               656   }
778   else                                            657   else
779   {                                               658   {
780     if ( distZ < distRMax )                       659     if ( distZ < distRMax )
781     {                                             660     {
782       distMin = distZ ;                           661       distMin = distZ ;
783       side    = kNZ   ;                           662       side    = kNZ   ;
784     }                                             663     }
785     else                                          664     else
786     {                                             665     {
787       distMin = distRMax ;                        666       distMin = distRMax ;
788       side    = kNRMax   ;                        667       side    = kNRMax   ;
789     }                                             668     }
790   }                                            << 669   }   
791   if (!fPhiFullCutTube  &&  (rho != 0.0) ) //  << 670   if (!fPhiFullCutTube  &&  rho ) // Protected against (0,0,z) 
792   {                                               671   {
793     phi = std::atan2(p.y(),p.x()) ;               672     phi = std::atan2(p.y(),p.x()) ;
794                                                   673 
795     if ( phi < 0 )  { phi += twopi; }             674     if ( phi < 0 )  { phi += twopi; }
796                                                   675 
797     if ( fSPhi < 0 )                              676     if ( fSPhi < 0 )
798     {                                             677     {
799       distSPhi = std::fabs(phi - (fSPhi + twop    678       distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
800     }                                             679     }
801     else                                          680     else
802     {                                             681     {
803       distSPhi = std::fabs(phi - fSPhi)*rho ;     682       distSPhi = std::fabs(phi - fSPhi)*rho ;
804     }                                             683     }
805     distEPhi = std::fabs(phi - fSPhi - fDPhi)*    684     distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
806                                                << 685                                       
807     if (distSPhi < distEPhi) // Find new minim    686     if (distSPhi < distEPhi) // Find new minimum
808     {                                             687     {
809       if ( distSPhi < distMin )                   688       if ( distSPhi < distMin )
810       {                                           689       {
811         side = kNSPhi ;                           690         side = kNSPhi ;
812       }                                           691       }
813     }                                             692     }
814     else                                          693     else
815     {                                             694     {
816       if ( distEPhi < distMin )                   695       if ( distEPhi < distMin )
817       {                                           696       {
818         side = kNEPhi ;                           697         side = kNEPhi ;
819       }                                           698       }
820     }                                             699     }
821   }                                            << 700   }    
822   switch ( side )                                 701   switch ( side )
823   {                                               702   {
824     case kNRMin : // Inner radius                 703     case kNRMin : // Inner radius
825     {                                          << 704     {                      
826       norm = G4ThreeVector(-p.x()/rho, -p.y()/    705       norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
827       break ;                                     706       break ;
828     }                                             707     }
829     case kNRMax : // Outer radius                 708     case kNRMax : // Outer radius
830     {                                          << 709     {                  
831       norm = G4ThreeVector(p.x()/rho, p.y()/rh    710       norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
832       break ;                                     711       break ;
833     }                                             712     }
834     case kNZ :    // + or - dz                    713     case kNZ :    // + or - dz
835     {                                          << 714     {                              
836       if ( distZHigh > distZLow )  { norm = fH    715       if ( distZHigh > distZLow )  { norm = fHighNorm ; }
837       else                         { norm = fL    716       else                         { norm = fLowNorm; }
838       break ;                                     717       break ;
839     }                                             718     }
840     case kNSPhi:                                  719     case kNSPhi:
841     {                                             720     {
842       norm = G4ThreeVector(sinSPhi, -cosSPhi,  << 721       norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
843       break ;                                     722       break ;
844     }                                             723     }
845     case kNEPhi:                                  724     case kNEPhi:
846     {                                             725     {
847       norm = G4ThreeVector(-sinEPhi, cosEPhi,  << 726       norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
848       break;                                      727       break;
849     }                                             728     }
850     default:      // Should never reach this c    729     default:      // Should never reach this case ...
851     {                                             730     {
852       DumpInfo();                                 731       DumpInfo();
853       G4Exception("G4CutTubs::ApproxSurfaceNor    732       G4Exception("G4CutTubs::ApproxSurfaceNormal()",
854                   "GeomSolids1002", JustWarnin    733                   "GeomSolids1002", JustWarning,
855                   "Undefined side for valid su    734                   "Undefined side for valid surface normal to solid.");
856       break ;                                     735       break ;
857     }                                          << 736     }    
858   }                                            << 737   }                
859   return norm;                                    738   return norm;
860 }                                                 739 }
861                                                   740 
862 //////////////////////////////////////////////    741 ////////////////////////////////////////////////////////////////////
863 //                                                742 //
864 //                                                743 //
865 // Calculate distance to shape from outside, a    744 // Calculate distance to shape from outside, along normalised vector
866 // - return kInfinity if no intersection, or i    745 // - return kInfinity if no intersection, or intersection distance <= tolerance
867 //                                                746 //
868 // - Compute the intersection with the z plane << 747 // - Compute the intersection with the z planes 
869 //        - if at valid r, phi, return            748 //        - if at valid r, phi, return
870 //                                                749 //
871 // -> If point is outer outer radius, compute     750 // -> If point is outer outer radius, compute intersection with rmax
872 //        - if at valid phi,z return              751 //        - if at valid phi,z return
873 //                                                752 //
874 // -> Compute intersection with inner radius,     753 // -> Compute intersection with inner radius, taking largest +ve root
875 //        - if valid (in z,phi), save intersct    754 //        - if valid (in z,phi), save intersction
876 //                                                755 //
877 //    -> If phi segmented, compute intersectio    756 //    -> If phi segmented, compute intersections with phi half planes
878 //        - return smallest of valid phi inter    757 //        - return smallest of valid phi intersections and
879 //          inner radius intersection             758 //          inner radius intersection
880 //                                                759 //
881 // NOTE:                                          760 // NOTE:
882 // - 'if valid' implies tolerant checking of i    761 // - 'if valid' implies tolerant checking of intersection points
883                                                   762 
884 G4double G4CutTubs::DistanceToIn( const G4Thre    763 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p,
885                                   const G4Thre    764                                   const G4ThreeVector& v  ) const
886 {                                                 765 {
887   G4double snxt = kInfinity ;      // snxt = d    766   G4double snxt = kInfinity ;      // snxt = default return value
888   G4double tolORMin2, tolIRMax2 ;  // 'generou    767   G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
889   G4double tolORMax2, tolIRMin2;                  768   G4double tolORMax2, tolIRMin2;
890   const G4double dRmax = 100.*fRMax;              769   const G4double dRmax = 100.*fRMax;
891   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);        770   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
892                                                << 771   
893   // Intersection point variables                 772   // Intersection point variables
894   //                                              773   //
895   G4double Dist, sd=0, xi, yi, zi, rho2, inum,    774   G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ;
896   G4double t1, t2, t3, b, c, d ;     // Quadra << 775   G4double t1, t2, t3, b, c, d ;     // Quadratic solver variables 
897   G4double distZLow,distZHigh;                    776   G4double distZLow,distZHigh;
898   // Calculate tolerant rmin and rmax             777   // Calculate tolerant rmin and rmax
899                                                   778 
900   if (fRMin > kRadTolerance)                      779   if (fRMin > kRadTolerance)
901   {                                               780   {
902     tolORMin2 = (fRMin - halfRadTolerance)*(fR    781     tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
903     tolIRMin2 = (fRMin + halfRadTolerance)*(fR    782     tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
904   }                                               783   }
905   else                                            784   else
906   {                                               785   {
907     tolORMin2 = 0.0 ;                             786     tolORMin2 = 0.0 ;
908     tolIRMin2 = 0.0 ;                             787     tolIRMin2 = 0.0 ;
909   }                                               788   }
910   tolORMax2 = (fRMax + halfRadTolerance)*(fRMa    789   tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
911   tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa    790   tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
912                                                   791 
913   // Intersection with ZCut surfaces              792   // Intersection with ZCut surfaces
914                                                   793 
915   // dist to Low Cut                              794   // dist to Low Cut
916   //                                              795   //
917   distZLow =(p+vZ).dot(fLowNorm);                 796   distZLow =(p+vZ).dot(fLowNorm);
918                                                   797 
919   // dist to High Cut                             798   // dist to High Cut
920   //                                              799   //
921   distZHigh = (p-vZ).dot(fHighNorm);              800   distZHigh = (p-vZ).dot(fHighNorm);
922                                                   801 
923   if ( distZLow >= -halfCarTolerance )            802   if ( distZLow >= -halfCarTolerance )
924   {                                               803   {
925     calf = v.dot(fLowNorm);                       804     calf = v.dot(fLowNorm);
926     if (calf<0)                                   805     if (calf<0)
927     {                                             806     {
928       sd = -distZLow/calf;                        807       sd = -distZLow/calf;
929       if(sd < 0.0)  { sd = 0.0; }                 808       if(sd < 0.0)  { sd = 0.0; }
930                                                   809 
931       xi   = p.x() + sd*v.x() ;                   810       xi   = p.x() + sd*v.x() ;                // Intersection coords
932       yi   = p.y() + sd*v.y() ;                   811       yi   = p.y() + sd*v.y() ;
933       rho2 = xi*xi + yi*yi ;                      812       rho2 = xi*xi + yi*yi ;
934                                                   813 
935       // Check validity of intersection           814       // Check validity of intersection
936                                                   815 
937       if ((tolIRMin2 <= rho2) && (rho2 <= tolI    816       if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
938       {                                           817       {
939         if (!fPhiFullCutTube && (rho2 != 0.0)) << 818         if (!fPhiFullCutTube && rho2)
940         {                                         819         {
941           // Psi = angle made with central (av    820           // Psi = angle made with central (average) phi of shape
942           //                                      821           //
943           inum   = xi*cosCPhi + yi*sinCPhi ;      822           inum   = xi*cosCPhi + yi*sinCPhi ;
944           iden   = std::sqrt(rho2) ;              823           iden   = std::sqrt(rho2) ;
945           cosPsi = inum/iden ;                    824           cosPsi = inum/iden ;
946           if (cosPsi >= cosHDPhiIT)  { return     825           if (cosPsi >= cosHDPhiIT)  { return sd ; }
947         }                                         826         }
948         else                                      827         else
949         {                                         828         {
950           return sd ;                             829           return sd ;
951         }                                         830         }
952       }                                           831       }
953     }                                             832     }
954     else                                          833     else
955     {                                             834     {
956       if ( sd<halfCarTolerance )                  835       if ( sd<halfCarTolerance )
957       {                                           836       {
958         if(calf>=0) { sd=kInfinity; }             837         if(calf>=0) { sd=kInfinity; }
959         return sd ;  // On/outside extent, and    838         return sd ;  // On/outside extent, and heading away
960       }              // -> cannot intersect       839       }              // -> cannot intersect
961     }                                             840     }
962   }                                               841   }
963                                                   842 
964   if(distZHigh >= -halfCarTolerance )             843   if(distZHigh >= -halfCarTolerance )
965   {                                               844   {
966     calf = v.dot(fHighNorm);                      845     calf = v.dot(fHighNorm);
967     if (calf<0)                                   846     if (calf<0)
968     {                                             847     {
969       sd = -distZHigh/calf;                       848       sd = -distZHigh/calf;
970                                                   849 
971       if(sd < 0.0)  { sd = 0.0; }                 850       if(sd < 0.0)  { sd = 0.0; }
972                                                   851 
973       xi   = p.x() + sd*v.x() ;                   852       xi   = p.x() + sd*v.x() ;                // Intersection coords
974       yi   = p.y() + sd*v.y() ;                   853       yi   = p.y() + sd*v.y() ;
975       rho2 = xi*xi + yi*yi ;                      854       rho2 = xi*xi + yi*yi ;
976                                                   855 
977       // Check validity of intersection           856       // Check validity of intersection
978                                                   857 
979       if ((tolIRMin2 <= rho2) && (rho2 <= tolI    858       if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
980       {                                           859       {
981         if (!fPhiFullCutTube && (rho2 != 0.0)) << 860         if (!fPhiFullCutTube && rho2)
982         {                                         861         {
983           // Psi = angle made with central (av    862           // Psi = angle made with central (average) phi of shape
984           //                                      863           //
985           inum   = xi*cosCPhi + yi*sinCPhi ;      864           inum   = xi*cosCPhi + yi*sinCPhi ;
986           iden   = std::sqrt(rho2) ;              865           iden   = std::sqrt(rho2) ;
987           cosPsi = inum/iden ;                    866           cosPsi = inum/iden ;
988           if (cosPsi >= cosHDPhiIT)  { return     867           if (cosPsi >= cosHDPhiIT)  { return sd ; }
989         }                                         868         }
990         else                                      869         else
991         {                                         870         {
992           return sd ;                             871           return sd ;
993         }                                         872         }
994       }                                           873       }
995     }                                             874     }
996     else                                          875     else
997     {                                             876     {
998       if ( sd<halfCarTolerance )                  877       if ( sd<halfCarTolerance )
999       {                                        << 878       { 
1000         if(calf>=0) { sd=kInfinity; }            879         if(calf>=0) { sd=kInfinity; }
1001         return sd ;  // On/outside extent, an    880         return sd ;  // On/outside extent, and heading away
1002       }              // -> cannot intersect      881       }              // -> cannot intersect
1003     }                                            882     }
1004   }                                              883   }
1005                                                  884 
1006   // -> Can not intersect z surfaces             885   // -> Can not intersect z surfaces
1007   //                                             886   //
1008   // Intersection with rmax (possible return)    887   // Intersection with rmax (possible return) and rmin (must also check phi)
1009   //                                             888   //
1010   // Intersection point (xi,yi,zi) on line x=    889   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1011   //                                             890   //
1012   // Intersects with x^2+y^2=R^2                 891   // Intersects with x^2+y^2=R^2
1013   //                                             892   //
1014   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v    893   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1015   //            t1                t2             894   //            t1                t2                t3
1016                                                  895 
1017   t1 = 1.0 - v.z()*v.z() ;                       896   t1 = 1.0 - v.z()*v.z() ;
1018   t2 = p.x()*v.x() + p.y()*v.y() ;               897   t2 = p.x()*v.x() + p.y()*v.y() ;
1019   t3 = p.x()*p.x() + p.y()*p.y() ;               898   t3 = p.x()*p.x() + p.y()*p.y() ;
1020   if ( t1 > 0 )        // Check not || to z a    899   if ( t1 > 0 )        // Check not || to z axis
1021   {                                              900   {
1022     b = t2/t1 ;                                  901     b = t2/t1 ;
1023     c = t3 - fRMax*fRMax ;                       902     c = t3 - fRMax*fRMax ;
1024                                               << 903     
1025     if ((t3 >= tolORMax2) && (t2<0))   // Thi    904     if ((t3 >= tolORMax2) && (t2<0))   // This also handles the tangent case
1026     {                                            905     {
1027       // Try outer cylinder intersection, c=(    906       // Try outer cylinder intersection, c=(t3-fRMax*fRMax)/t1;
1028                                                  907 
1029       c /= t1 ;                                  908       c /= t1 ;
1030       d = b*b - c ;                              909       d = b*b - c ;
1031                                                  910 
1032       if (d >= 0)  // If real root               911       if (d >= 0)  // If real root
1033       {                                          912       {
1034         sd = c/(-b+std::sqrt(d));                913         sd = c/(-b+std::sqrt(d));
1035         if (sd >= 0)  // If 'forwards'           914         if (sd >= 0)  // If 'forwards'
1036         {                                        915         {
1037           if ( sd>dRmax ) // Avoid rounding e    916           if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1038           {               // 64 bits systems.    917           {               // 64 bits systems. Split long distances and recompute
1039             G4double fTerm = sd-std::fmod(sd,    918             G4double fTerm = sd-std::fmod(sd,dRmax);
1040             sd = fTerm + DistanceToIn(p+fTerm    919             sd = fTerm + DistanceToIn(p+fTerm*v,v);
1041           }                                   << 920           } 
1042           // Check z intersection                921           // Check z intersection
1043           //                                     922           //
1044           zi = p.z() + sd*v.z() ;                923           zi = p.z() + sd*v.z() ;
1045           xi = p.x() + sd*v.x() ;                924           xi = p.x() + sd*v.x() ;
1046           yi = p.y() + sd*v.y() ;                925           yi = p.y() + sd*v.y() ;
1047           if ((-xi*fLowNorm.x()-yi*fLowNorm.y    926           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1048                -(zi+fDz)*fLowNorm.z())>-halfC    927                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1049           {                                      928           {
1050             if ((-xi*fHighNorm.x()-yi*fHighNo    929             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1051                  +(fDz-zi)*fHighNorm.z())>-ha    930                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1052             {                                    931             {
1053               // Z ok. Check phi intersection    932               // Z ok. Check phi intersection if reqd
1054               //                                 933               //
1055               if (fPhiFullCutTube)               934               if (fPhiFullCutTube)
1056               {                                  935               {
1057                 return sd ;                      936                 return sd ;
1058               }                                  937               }
1059               else                               938               else
1060               {                                  939               {
1061                 xi     = p.x() + sd*v.x() ;      940                 xi     = p.x() + sd*v.x() ;
1062                 yi     = p.y() + sd*v.y() ;      941                 yi     = p.y() + sd*v.y() ;
1063                 cosPsi = (xi*cosCPhi + yi*sin    942                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
1064                 if (cosPsi >= cosHDPhiIT)  {     943                 if (cosPsi >= cosHDPhiIT)  { return sd ; }
1065               }                                  944               }
1066             }  //  end if std::fabs(zi)          945             }  //  end if std::fabs(zi)
1067           }                                      946           }
1068         }    //  end if (sd>=0)                  947         }    //  end if (sd>=0)
1069       }      //  end if (d>=0)                   948       }      //  end if (d>=0)
1070     }        //  end if (r>=fRMax)               949     }        //  end if (r>=fRMax)
1071     else                                      << 950     else 
1072     {                                            951     {
1073       // Inside outer radius :                   952       // Inside outer radius :
1074       // check not inside, and heading throug    953       // check not inside, and heading through tubs (-> 0 to in)
1075       if ((t3 > tolIRMin2) && (t2 < 0)           954       if ((t3 > tolIRMin2) && (t2 < 0)
1076        && (std::fabs(p.z()) <= std::fabs(GetC    955        && (std::fabs(p.z()) <= std::fabs(GetCutZ(p))-halfCarTolerance ))
1077       {                                          956       {
1078         // Inside both radii, delta r -ve, in    957         // Inside both radii, delta r -ve, inside z extent
1079                                                  958 
1080         if (!fPhiFullCutTube)                    959         if (!fPhiFullCutTube)
1081         {                                        960         {
1082           inum   = p.x()*cosCPhi + p.y()*sinC    961           inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
1083           iden   = std::sqrt(t3) ;               962           iden   = std::sqrt(t3) ;
1084           cosPsi = inum/iden ;                   963           cosPsi = inum/iden ;
1085           if (cosPsi >= cosHDPhiIT)              964           if (cosPsi >= cosHDPhiIT)
1086           {                                      965           {
1087             // In the old version, the small     966             // In the old version, the small negative tangent for the point
1088             // on surface was not taken in ac    967             // on surface was not taken in account, and returning 0.0 ...
1089             // New version: check the tangent << 968             // New version: check the tangent for the point on surface and 
1090             // if no intersection, return kIn    969             // if no intersection, return kInfinity, if intersection instead
1091             // return sd.                        970             // return sd.
1092             //                                   971             //
1093             c = t3-fRMax*fRMax;               << 972             c = t3-fRMax*fRMax; 
1094             if ( c<=0.0 )                        973             if ( c<=0.0 )
1095             {                                    974             {
1096               return 0.0;                        975               return 0.0;
1097             }                                    976             }
1098             else                                 977             else
1099             {                                    978             {
1100               c = c/t1 ;                         979               c = c/t1 ;
1101               d = b*b-c;                         980               d = b*b-c;
1102               if ( d>=0.0 )                      981               if ( d>=0.0 )
1103               {                                  982               {
1104                 snxt = c/(-b+std::sqrt(d)); /    983                 snxt = c/(-b+std::sqrt(d)); // using safe solution
1105                                             / << 984                                             // for quadratic equation 
1106                 if ( snxt < halfCarTolerance     985                 if ( snxt < halfCarTolerance ) { snxt=0; }
1107                 return snxt ;                    986                 return snxt ;
1108               }                               << 987               }      
1109               else                               988               else
1110               {                                  989               {
1111                 return kInfinity;                990                 return kInfinity;
1112               }                                  991               }
1113             }                                    992             }
1114           }                                   << 993           } 
1115         }                                        994         }
1116         else                                     995         else
1117         {                                     << 996         {   
1118           // In the old version, the small ne    997           // In the old version, the small negative tangent for the point
1119           // on surface was not taken in acco    998           // on surface was not taken in account, and returning 0.0 ...
1120           // New version: check the tangent f << 999           // New version: check the tangent for the point on surface and 
1121           // if no intersection, return kInfi    1000           // if no intersection, return kInfinity, if intersection instead
1122           // return sd.                          1001           // return sd.
1123           //                                     1002           //
1124           c = t3 - fRMax*fRMax;               << 1003           c = t3 - fRMax*fRMax; 
1125           if ( c<=0.0 )                          1004           if ( c<=0.0 )
1126           {                                      1005           {
1127             return 0.0;                          1006             return 0.0;
1128           }                                      1007           }
1129           else                                   1008           else
1130           {                                      1009           {
1131             c = c/t1 ;                           1010             c = c/t1 ;
1132             d = b*b-c;                           1011             d = b*b-c;
1133             if ( d>=0.0 )                        1012             if ( d>=0.0 )
1134             {                                    1013             {
1135               snxt= c/(-b+std::sqrt(d)); // u    1014               snxt= c/(-b+std::sqrt(d)); // using safe solution
1136                                          // f << 1015                                          // for quadratic equation 
1137               if ( snxt < halfCarTolerance )     1016               if ( snxt < halfCarTolerance ) { snxt=0; }
1138               return snxt ;                      1017               return snxt ;
1139             }                                 << 1018             }      
1140             else                                 1019             else
1141             {                                    1020             {
1142               return kInfinity;                  1021               return kInfinity;
1143             }                                    1022             }
1144           }                                      1023           }
1145         } // end if   (!fPhiFullCutTube)         1024         } // end if   (!fPhiFullCutTube)
1146       }   // end if   (t3>tolIRMin2)             1025       }   // end if   (t3>tolIRMin2)
1147     }     // end if   (Inside Outer Radius)   << 1026     }     // end if   (Inside Outer Radius) 
1148                                               << 1027       
1149     if ( fRMin != 0.0 )    // Try inner cylin << 1028     if ( fRMin )    // Try inner cylinder intersection
1150     {                                            1029     {
1151       c = (t3 - fRMin*fRMin)/t1 ;                1030       c = (t3 - fRMin*fRMin)/t1 ;
1152       d = b*b - c ;                              1031       d = b*b - c ;
1153       if ( d >= 0.0 )  // If real root           1032       if ( d >= 0.0 )  // If real root
1154       {                                          1033       {
1155         // Always want 2nd root - we are outs    1034         // Always want 2nd root - we are outside and know rmax Hit was bad
1156         // - If on surface of rmin also need     1035         // - If on surface of rmin also need farthest root
1157                                               << 1036         
1158         sd =( b > 0. )? c/(-b - std::sqrt(d))    1037         sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1159         if (sd >= -10*halfCarTolerance)  // c    1038         if (sd >= -10*halfCarTolerance)  // check forwards
1160         {                                        1039         {
1161           // Check z intersection                1040           // Check z intersection
1162           //                                     1041           //
1163           if (sd < 0.0)  { sd = 0.0; }           1042           if (sd < 0.0)  { sd = 0.0; }
1164           if (sd>dRmax) // Avoid rounding err    1043           if (sd>dRmax) // Avoid rounding errors due to precision issues seen
1165           {             // 64 bits systems. S    1044           {             // 64 bits systems. Split long distances and recompute
1166             G4double fTerm = sd-std::fmod(sd,    1045             G4double fTerm = sd-std::fmod(sd,dRmax);
1167             sd = fTerm + DistanceToIn(p+fTerm    1046             sd = fTerm + DistanceToIn(p+fTerm*v,v);
1168           }                                   << 1047           } 
1169           zi = p.z() + sd*v.z() ;                1048           zi = p.z() + sd*v.z() ;
1170           xi = p.x() + sd*v.x() ;                1049           xi = p.x() + sd*v.x() ;
1171           yi = p.y() + sd*v.y() ;                1050           yi = p.y() + sd*v.y() ;
1172           if ((-xi*fLowNorm.x()-yi*fLowNorm.y    1051           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1173                -(zi+fDz)*fLowNorm.z())>-halfC    1052                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1174           {                                      1053           {
1175             if ((-xi*fHighNorm.x()-yi*fHighNo    1054             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1176                  +(fDz-zi)*fHighNorm.z())>-ha    1055                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1177             {                                    1056             {
1178               // Z ok. Check phi                 1057               // Z ok. Check phi
1179               //                                 1058               //
1180               if ( fPhiFullCutTube )             1059               if ( fPhiFullCutTube )
1181               {                                  1060               {
1182                 return sd ;                   << 1061                 return sd ; 
1183               }                                  1062               }
1184               else                               1063               else
1185               {                                  1064               {
1186                 cosPsi = (xi*cosCPhi + yi*sin    1065                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1187                 if (cosPsi >= cosHDPhiIT)        1066                 if (cosPsi >= cosHDPhiIT)
1188                 {                                1067                 {
1189                   // Good inner radius isect     1068                   // Good inner radius isect
1190                   // - but earlier phi isect     1069                   // - but earlier phi isect still possible
1191                   //                             1070                   //
1192                   snxt = sd ;                    1071                   snxt = sd ;
1193                 }                                1072                 }
1194               }                                  1073               }
1195             }      //    end if std::fabs(zi)    1074             }      //    end if std::fabs(zi)
1196           }                                      1075           }
1197         }          //    end if (sd>=0)          1076         }          //    end if (sd>=0)
1198       }            //    end if (d>=0)           1077       }            //    end if (d>=0)
1199     }              //    end if (fRMin)          1078     }              //    end if (fRMin)
1200   }                                              1079   }
1201                                                  1080 
1202   // Phi segment intersection                    1081   // Phi segment intersection
1203   //                                             1082   //
1204   // o Tolerant of points inside phi planes b    1083   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1205   //                                             1084   //
1206   // o NOTE: Large duplication of code betwee    1085   // o NOTE: Large duplication of code between sphi & ephi checks
1207   //         -> only diffs: sphi -> ephi, Com    1086   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1208   //            intersection check <=0 -> >=0    1087   //            intersection check <=0 -> >=0
1209   //         -> use some form of loop Constru    1088   //         -> use some form of loop Construct ?
1210   //                                             1089   //
1211   if ( !fPhiFullCutTube )                        1090   if ( !fPhiFullCutTube )
1212   {                                              1091   {
1213     // First phi surface (Starting phi)          1092     // First phi surface (Starting phi)
1214     //                                           1093     //
1215     Comp = v.x()*sinSPhi - v.y()*cosSPhi ;       1094     Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1216                                               << 1095                 
1217     if ( Comp < 0 )  // Component in outwards    1096     if ( Comp < 0 )  // Component in outwards normal dirn
1218     {                                            1097     {
1219       Dist = (p.y()*cosSPhi - p.x()*sinSPhi)     1098       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1220                                                  1099 
1221       if ( Dist < halfCarTolerance )             1100       if ( Dist < halfCarTolerance )
1222       {                                          1101       {
1223         sd = Dist/Comp ;                         1102         sd = Dist/Comp ;
1224                                                  1103 
1225         if (sd < snxt)                           1104         if (sd < snxt)
1226         {                                        1105         {
1227           if ( sd < 0 )  { sd = 0.0; }           1106           if ( sd < 0 )  { sd = 0.0; }
1228           zi = p.z() + sd*v.z() ;                1107           zi = p.z() + sd*v.z() ;
1229           xi = p.x() + sd*v.x() ;                1108           xi = p.x() + sd*v.x() ;
1230           yi = p.y() + sd*v.y() ;                1109           yi = p.y() + sd*v.y() ;
1231           if ((-xi*fLowNorm.x()-yi*fLowNorm.y    1110           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1232                -(zi+fDz)*fLowNorm.z())>-halfC    1111                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1233           {                                      1112           {
1234             if ((-xi*fHighNorm.x()-yi*fHighNo    1113             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1235                  +(fDz-zi)*fHighNorm.z())>-ha << 1114                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 
1236             {                                    1115             {
1237               rho2 = xi*xi + yi*yi ;             1116               rho2 = xi*xi + yi*yi ;
1238               if ( ( (rho2 >= tolIRMin2) && (    1117               if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1239                 || ( (rho2 >  tolORMin2) && (    1118                 || ( (rho2 >  tolORMin2) && (rho2 <  tolIRMin2)
1240                   && ( v.y()*cosSPhi - v.x()*    1119                   && ( v.y()*cosSPhi - v.x()*sinSPhi >  0 )
1241                   && ( v.x()*cosSPhi + v.y()*    1120                   && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 )     )
1242                 || ( (rho2 > tolIRMax2) && (r    1121                 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1243                   && (v.y()*cosSPhi - v.x()*s    1122                   && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1244                   && (v.x()*cosSPhi + v.y()*s    1123                   && (v.x()*cosSPhi + v.y()*sinSPhi < 0) )    )
1245               {                                  1124               {
1246                 // z and r intersections good    1125                 // z and r intersections good
1247                 // - check intersecting with     1126                 // - check intersecting with correct half-plane
1248                 //                               1127                 //
1249                 if ((yi*cosCPhi-xi*sinCPhi) <    1128                 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1250               }                                  1129               }
1251             }   //two Z conditions               1130             }   //two Z conditions
1252           }                                      1131           }
1253         }                                        1132         }
1254       }                                       << 1133       }    
1255     }                                            1134     }
1256                                               << 1135       
1257     // Second phi surface (Ending phi)           1136     // Second phi surface (Ending phi)
1258     //                                           1137     //
1259     Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;    1138     Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1260                                               << 1139         
1261     if (Comp < 0 )  // Component in outwards     1140     if (Comp < 0 )  // Component in outwards normal dirn
1262     {                                            1141     {
1263       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi)    1142       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1264                                                  1143 
1265       if ( Dist < halfCarTolerance )             1144       if ( Dist < halfCarTolerance )
1266       {                                          1145       {
1267         sd = Dist/Comp ;                         1146         sd = Dist/Comp ;
1268                                                  1147 
1269         if (sd < snxt)                           1148         if (sd < snxt)
1270         {                                        1149         {
1271           if ( sd < 0 )  { sd = 0; }             1150           if ( sd < 0 )  { sd = 0; }
1272           zi = p.z() + sd*v.z() ;                1151           zi = p.z() + sd*v.z() ;
1273           xi = p.x() + sd*v.x() ;                1152           xi = p.x() + sd*v.x() ;
1274           yi = p.y() + sd*v.y() ;                1153           yi = p.y() + sd*v.y() ;
1275           if ((-xi*fLowNorm.x()-yi*fLowNorm.y    1154           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1276                -(zi+fDz)*fLowNorm.z())>-halfC    1155                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1277           {                                      1156           {
1278             if ((-xi*fHighNorm.x()-yi*fHighNo    1157             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1279                  +(fDz-zi)*fHighNorm.z())>-ha    1158                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1280             {                                    1159             {
1281               xi   = p.x() + sd*v.x() ;          1160               xi   = p.x() + sd*v.x() ;
1282               yi   = p.y() + sd*v.y() ;          1161               yi   = p.y() + sd*v.y() ;
1283               rho2 = xi*xi + yi*yi ;             1162               rho2 = xi*xi + yi*yi ;
1284               if ( ( (rho2 >= tolIRMin2) && (    1163               if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1285                   || ( (rho2 > tolORMin2)  &&    1164                   || ( (rho2 > tolORMin2)  && (rho2 < tolIRMin2)
1286                     && (v.x()*sinEPhi - v.y()    1165                     && (v.x()*sinEPhi - v.y()*cosEPhi >  0)
1287                     && (v.x()*cosEPhi + v.y()    1166                     && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1288                   || ( (rho2 > tolIRMax2) &&     1167                   || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1289                     && (v.x()*sinEPhi - v.y()    1168                     && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1290                     && (v.x()*cosEPhi + v.y()    1169                     && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1291               {                                  1170               {
1292                 // z and r intersections good    1171                 // z and r intersections good
1293                 // - check intersecting with     1172                 // - check intersecting with correct half-plane
1294                 //                               1173                 //
1295                 if ( (yi*cosCPhi-xi*sinCPhi)     1174                 if ( (yi*cosCPhi-xi*sinCPhi) >= -halfCarTolerance )
1296                 {                                1175                 {
1297                   snxt = sd;                     1176                   snxt = sd;
1298                 }                                1177                 }
1299               }    //?? >=-halfCarTolerance      1178               }    //?? >=-halfCarTolerance
1300             }                                    1179             }
1301           }  // two Z conditions                 1180           }  // two Z conditions
1302         }                                        1181         }
1303       }                                          1182       }
1304     }         //  Comp < 0                       1183     }         //  Comp < 0
1305   }           //  !fPhiFullTube               << 1184   }           //  !fPhiFullTube 
1306   if ( snxt<halfCarTolerance )  { snxt=0; }      1185   if ( snxt<halfCarTolerance )  { snxt=0; }
1307                                                  1186 
1308   return snxt ;                                  1187   return snxt ;
1309 }                                                1188 }
1310                                               << 1189  
1311 /////////////////////////////////////////////    1190 //////////////////////////////////////////////////////////////////
1312 //                                               1191 //
1313 // Calculate distance to shape from outside,     1192 // Calculate distance to shape from outside, along normalised vector
1314 // - return kInfinity if no intersection, or     1193 // - return kInfinity if no intersection, or intersection distance <= tolerance
1315 //                                               1194 //
1316 // - Compute the intersection with the z plan << 1195 // - Compute the intersection with the z planes 
1317 //        - if at valid r, phi, return           1196 //        - if at valid r, phi, return
1318 //                                               1197 //
1319 // -> If point is outer outer radius, compute    1198 // -> If point is outer outer radius, compute intersection with rmax
1320 //        - if at valid phi,z return             1199 //        - if at valid phi,z return
1321 //                                               1200 //
1322 // -> Compute intersection with inner radius,    1201 // -> Compute intersection with inner radius, taking largest +ve root
1323 //        - if valid (in z,phi), save intersc    1202 //        - if valid (in z,phi), save intersction
1324 //                                               1203 //
1325 //    -> If phi segmented, compute intersecti    1204 //    -> If phi segmented, compute intersections with phi half planes
1326 //        - return smallest of valid phi inte    1205 //        - return smallest of valid phi intersections and
1327 //          inner radius intersection            1206 //          inner radius intersection
1328 //                                               1207 //
1329 // NOTE:                                         1208 // NOTE:
1330 // - Precalculations for phi trigonometry are    1209 // - Precalculations for phi trigonometry are Done `just in time'
1331 // - `if valid' implies tolerant checking of     1210 // - `if valid' implies tolerant checking of intersection points
1332 //   Calculate distance (<= actual) to closes    1211 //   Calculate distance (<= actual) to closest surface of shape from outside
1333 // - Calculate distance to z, radial planes      1212 // - Calculate distance to z, radial planes
1334 // - Only to phi planes if outside phi extent    1213 // - Only to phi planes if outside phi extent
1335 // - Return 0 if point inside                    1214 // - Return 0 if point inside
1336                                                  1215 
1337 G4double G4CutTubs::DistanceToIn( const G4Thr    1216 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p ) const
1338 {                                                1217 {
1339   G4double safRMin,safRMax,safZLow,safZHigh,s    1218   G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1340   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);       1219   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1341                                                  1220 
1342   // Distance to R                               1221   // Distance to R
1343   //                                             1222   //
1344   rho = std::sqrt(p.x()*p.x() + p.y()*p.y())     1223   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1345                                                  1224 
1346   safRMin = fRMin- rho ;                         1225   safRMin = fRMin- rho ;
1347   safRMax = rho - fRMax ;                        1226   safRMax = rho - fRMax ;
1348                                                  1227 
1349   // Distances to ZCut(Low/High)                 1228   // Distances to ZCut(Low/High)
1350                                                  1229 
1351   // Dist to Low Cut                             1230   // Dist to Low Cut
1352   //                                             1231   //
1353   safZLow = (p+vZ).dot(fLowNorm);                1232   safZLow = (p+vZ).dot(fLowNorm);
1354                                                  1233 
1355   // Dist to High Cut                            1234   // Dist to High Cut
1356   //                                             1235   //
1357   safZHigh = (p-vZ).dot(fHighNorm);              1236   safZHigh = (p-vZ).dot(fHighNorm);
1358                                                  1237 
1359   safe = std::max(safZLow,safZHigh);             1238   safe = std::max(safZLow,safZHigh);
1360                                                  1239 
1361   if ( safRMin > safe ) { safe = safRMin; }      1240   if ( safRMin > safe ) { safe = safRMin; }
1362   if ( safRMax> safe )  { safe = safRMax; }      1241   if ( safRMax> safe )  { safe = safRMax; }
1363                                                  1242 
1364   // Distance to Phi                             1243   // Distance to Phi
1365   //                                             1244   //
1366   if ( (!fPhiFullCutTube) && ((rho) != 0.0) ) << 1245   if ( (!fPhiFullCutTube) && (rho) )
1367    {                                             1246    {
1368      // Psi=angle from central phi to point      1247      // Psi=angle from central phi to point
1369      //                                          1248      //
1370      cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)    1249      cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1371                                               << 1250      
1372      if ( cosPsi < cosHDPhi )                 << 1251      if ( cosPsi < std::cos(fDPhi*0.5) )
1373      {                                           1252      {
1374        // Point lies outside phi range           1253        // Point lies outside phi range
1375                                               << 1254  
1376        if ( (p.y()*cosCPhi - p.x()*sinCPhi) <    1255        if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1377        {                                         1256        {
1378          safePhi = std::fabs(p.x()*sinSPhi -     1257          safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1379        }                                         1258        }
1380        else                                      1259        else
1381        {                                         1260        {
1382          safePhi = std::fabs(p.x()*sinEPhi -     1261          safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1383        }                                         1262        }
1384        if ( safePhi > safe )  { safe = safePh    1263        if ( safePhi > safe )  { safe = safePhi; }
1385      }                                           1264      }
1386    }                                             1265    }
1387    if ( safe < 0 )  { safe = 0; }                1266    if ( safe < 0 )  { safe = 0; }
1388                                                  1267 
1389    return safe ;                                 1268    return safe ;
1390 }                                                1269 }
1391                                                  1270 
1392 /////////////////////////////////////////////    1271 //////////////////////////////////////////////////////////////////////////////
1393 //                                               1272 //
1394 // Calculate distance to surface of shape fro    1273 // Calculate distance to surface of shape from `inside', allowing for tolerance
1395 // - Only Calc rmax intersection if no valid     1274 // - Only Calc rmax intersection if no valid rmin intersection
1396                                                  1275 
1397 G4double G4CutTubs::DistanceToOut( const G4Th    1276 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p,
1398                                    const G4Th    1277                                    const G4ThreeVector& v,
1399                                    const G4bo    1278                                    const G4bool calcNorm,
1400                                          G4bo << 1279                                          G4bool *validNorm,
1401                                          G4Th << 1280                                          G4ThreeVector *n    ) const
1402 {                                             << 1281 {  
1403   enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,k << 
1404                                               << 
1405   ESide side=kNull , sider=kNull, sidephi=kNu    1282   ESide side=kNull , sider=kNull, sidephi=kNull ;
1406   G4double snxt=kInfinity, srd=kInfinity,sz=k    1283   G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ;
1407   G4double deltaR, t1, t2, t3, b, c, d2, roMi    1284   G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1408   G4double distZLow,distZHigh,calfH,calfL;       1285   G4double distZLow,distZHigh,calfH,calfL;
1409   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);       1286   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1410                                               << 1287  
1411   // Vars for phi intersection:                  1288   // Vars for phi intersection:
1412   //                                             1289   //
1413   G4double pDistS, compS, pDistE, compE, sphi    1290   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1414                                               << 1291  
1415   // Z plane intersection                        1292   // Z plane intersection
1416   // Distances to ZCut(Low/High)                 1293   // Distances to ZCut(Low/High)
1417                                                  1294 
1418   // dist to Low Cut                             1295   // dist to Low Cut
1419   //                                             1296   //
1420   distZLow =(p+vZ).dot(fLowNorm);                1297   distZLow =(p+vZ).dot(fLowNorm);
1421                                                  1298 
1422   // dist to High Cut                            1299   // dist to High Cut
1423   //                                             1300   //
1424   distZHigh = (p-vZ).dot(fHighNorm);             1301   distZHigh = (p-vZ).dot(fHighNorm);
1425                                                  1302 
1426   calfH = v.dot(fHighNorm);                      1303   calfH = v.dot(fHighNorm);
1427   calfL = v.dot(fLowNorm);                       1304   calfL = v.dot(fLowNorm);
1428                                                  1305 
1429   if (calfH > 0 )                                1306   if (calfH > 0 )
1430   {                                              1307   {
1431     if ( distZHigh < halfCarTolerance )          1308     if ( distZHigh < halfCarTolerance )
1432     {                                            1309     {
1433       snxt = -distZHigh/calfH ;                  1310       snxt = -distZHigh/calfH ;
1434       side = kPZ ;                               1311       side = kPZ ;
1435     }                                            1312     }
1436     else                                         1313     else
1437     {                                            1314     {
1438       if (calcNorm)                              1315       if (calcNorm)
1439       {                                          1316       {
1440         *n         = G4ThreeVector(0,0,1) ;      1317         *n         = G4ThreeVector(0,0,1) ;
1441         *validNorm = true ;                      1318         *validNorm = true ;
1442       }                                          1319       }
1443       return snxt = 0 ;                          1320       return snxt = 0 ;
1444     }                                            1321     }
1445  }                                               1322  }
1446   if ( calfL>0)                                  1323   if ( calfL>0)
1447   {                                              1324   {
1448                                               << 1325    
1449     if ( distZLow < halfCarTolerance )           1326     if ( distZLow < halfCarTolerance )
1450     {                                            1327     {
1451       sz = -distZLow/calfL ;                     1328       sz = -distZLow/calfL ;
1452       if(sz<snxt){                               1329       if(sz<snxt){
1453       snxt=sz;                                   1330       snxt=sz;
1454       side = kMZ ;                               1331       side = kMZ ;
1455       }                                          1332       }
1456                                               << 1333       
1457     }                                            1334     }
1458     else                                         1335     else
1459     {                                            1336     {
1460       if (calcNorm)                              1337       if (calcNorm)
1461       {                                          1338       {
1462         *n         = G4ThreeVector(0,0,-1) ;     1339         *n         = G4ThreeVector(0,0,-1) ;
1463         *validNorm = true ;                      1340         *validNorm = true ;
1464       }                                          1341       }
1465       return snxt = 0.0 ;                        1342       return snxt = 0.0 ;
1466     }                                            1343     }
1467   }                                              1344   }
1468   if((calfH<=0)&&(calfL<=0))                     1345   if((calfH<=0)&&(calfL<=0))
1469   {                                              1346   {
1470     snxt = kInfinity ;    // Travel perpendic    1347     snxt = kInfinity ;    // Travel perpendicular to z axis
1471     side = kNull;                                1348     side = kNull;
1472   }                                              1349   }
1473   // Radial Intersections                        1350   // Radial Intersections
1474   //                                             1351   //
1475   // Find intersection with cylinders at rmax    1352   // Find intersection with cylinders at rmax/rmin
1476   // Intersection point (xi,yi,zi) on line x=    1353   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1477   //                                             1354   //
1478   // Intersects with x^2+y^2=R^2                 1355   // Intersects with x^2+y^2=R^2
1479   //                                             1356   //
1480   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v    1357   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1481   //                                             1358   //
1482   //            t1                t2             1359   //            t1                t2                    t3
1483                                                  1360 
1484   t1   = 1.0 - v.z()*v.z() ;      // since v     1361   t1   = 1.0 - v.z()*v.z() ;      // since v normalised
1485   t2   = p.x()*v.x() + p.y()*v.y() ;             1362   t2   = p.x()*v.x() + p.y()*v.y() ;
1486   t3   = p.x()*p.x() + p.y()*p.y() ;             1363   t3   = p.x()*p.x() + p.y()*p.y() ;
1487                                                  1364 
1488   if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fR    1365   if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fRMax*fRMax; }
1489   else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t    1366   else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }        // radius^2 on +-fDz
1490                                                  1367 
1491   if ( t1 > 0 ) // Check not parallel            1368   if ( t1 > 0 ) // Check not parallel
1492   {                                              1369   {
1493     // Calculate srd, r exit distance            1370     // Calculate srd, r exit distance
1494                                               << 1371      
1495     if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax     1372     if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1496     {                                            1373     {
1497       // Delta r not negative => leaving via     1374       // Delta r not negative => leaving via rmax
1498                                                  1375 
1499       deltaR = t3 - fRMax*fRMax ;                1376       deltaR = t3 - fRMax*fRMax ;
1500                                                  1377 
1501       // NOTE: Should use rho-fRMax<-kRadTole    1378       // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1502       // - avoid sqrt for efficiency             1379       // - avoid sqrt for efficiency
1503                                                  1380 
1504       if ( deltaR < -kRadTolerance*fRMax )       1381       if ( deltaR < -kRadTolerance*fRMax )
1505       {                                          1382       {
1506         b     = t2/t1 ;                          1383         b     = t2/t1 ;
1507         c     = deltaR/t1 ;                      1384         c     = deltaR/t1 ;
1508         d2    = b*b-c;                           1385         d2    = b*b-c;
1509         if( d2 >= 0 ) { srd = c/( -b - std::s    1386         if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1510         else          { srd = 0.; }              1387         else          { srd = 0.; }
1511         sider = kRMax ;                          1388         sider = kRMax ;
1512       }                                          1389       }
1513       else                                       1390       else
1514       {                                          1391       {
1515         // On tolerant boundary & heading out    1392         // On tolerant boundary & heading outwards (or perpendicular to)
1516         // outer radial surface -> leaving im    1393         // outer radial surface -> leaving immediately
1517                                                  1394 
1518         if ( calcNorm )                       << 1395         if ( calcNorm ) 
1519         {                                        1396         {
1520           *n         = G4ThreeVector(p.x()/fR    1397           *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1521           *validNorm = true ;                    1398           *validNorm = true ;
1522         }                                        1399         }
1523         return snxt = 0 ; // Leaving by rmax     1400         return snxt = 0 ; // Leaving by rmax immediately
1524       }                                          1401       }
1525     }                                         << 1402     }             
1526     else if ( t2 < 0. ) // i.e.  t2 < 0; Poss    1403     else if ( t2 < 0. ) // i.e.  t2 < 0; Possible rmin intersection
1527     {                                            1404     {
1528       roMin2 = t3 - t2*t2/t1 ; // min ro2 of  << 1405       roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement 
1529                                                  1406 
1530       if ( (fRMin != 0.0) && (roMin2 < fRMin* << 1407       if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1531       {                                          1408       {
1532         deltaR = t3 - fRMin*fRMin ;              1409         deltaR = t3 - fRMin*fRMin ;
1533         b      = t2/t1 ;                         1410         b      = t2/t1 ;
1534         c      = deltaR/t1 ;                     1411         c      = deltaR/t1 ;
1535         d2     = b*b - c ;                       1412         d2     = b*b - c ;
1536                                                  1413 
1537         if ( d2 >= 0 )   // Leaving via rmin     1414         if ( d2 >= 0 )   // Leaving via rmin
1538         {                                        1415         {
1539           // NOTE: SHould use rho-rmin>kRadTo    1416           // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1540           // - avoid sqrt for efficiency         1417           // - avoid sqrt for efficiency
1541                                                  1418 
1542           if (deltaR > kRadTolerance*fRMin)      1419           if (deltaR > kRadTolerance*fRMin)
1543           {                                      1420           {
1544             srd = c/(-b+std::sqrt(d2));       << 1421             srd = c/(-b+std::sqrt(d2)); 
1545             sider = kRMin ;                      1422             sider = kRMin ;
1546           }                                      1423           }
1547           else                                   1424           else
1548           {                                      1425           {
1549             if ( calcNorm ) { *validNorm = fa    1426             if ( calcNorm ) { *validNorm = false; }  // Concave side
1550             return snxt = 0.0;                   1427             return snxt = 0.0;
1551           }                                      1428           }
1552         }                                        1429         }
1553         else    // No rmin intersect -> must     1430         else    // No rmin intersect -> must be rmax intersect
1554         {                                        1431         {
1555           deltaR = t3 - fRMax*fRMax ;            1432           deltaR = t3 - fRMax*fRMax ;
1556           c     = deltaR/t1 ;                    1433           c     = deltaR/t1 ;
1557           d2    = b*b-c;                         1434           d2    = b*b-c;
1558           if( d2 >=0. )                          1435           if( d2 >=0. )
1559           {                                      1436           {
1560             srd    = -b + std::sqrt(d2) ;        1437             srd    = -b + std::sqrt(d2) ;
1561             sider  = kRMax ;                     1438             sider  = kRMax ;
1562           }                                      1439           }
1563           else // Case: On the border+t2<kRad    1440           else // Case: On the border+t2<kRadTolerance
1564                //       (v is perpendicular t    1441                //       (v is perpendicular to the surface)
1565           {                                      1442           {
1566             if (calcNorm)                        1443             if (calcNorm)
1567             {                                    1444             {
1568               *n = G4ThreeVector(p.x()/fRMax,    1445               *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1569               *validNorm = true ;                1446               *validNorm = true ;
1570             }                                    1447             }
1571             return snxt = 0.0;                   1448             return snxt = 0.0;
1572           }                                      1449           }
1573         }                                        1450         }
1574       }                                          1451       }
1575       else if ( roi2 > fRMax*(fRMax + kRadTol    1452       else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1576            // No rmin intersect -> must be rm    1453            // No rmin intersect -> must be rmax intersect
1577       {                                          1454       {
1578         deltaR = t3 - fRMax*fRMax ;              1455         deltaR = t3 - fRMax*fRMax ;
1579         b      = t2/t1 ;                         1456         b      = t2/t1 ;
1580         c      = deltaR/t1;                      1457         c      = deltaR/t1;
1581         d2     = b*b-c;                          1458         d2     = b*b-c;
1582         if( d2 >= 0 )                            1459         if( d2 >= 0 )
1583         {                                        1460         {
1584           srd    = -b + std::sqrt(d2) ;          1461           srd    = -b + std::sqrt(d2) ;
1585           sider  = kRMax ;                       1462           sider  = kRMax ;
1586         }                                        1463         }
1587         else // Case: On the border+t2<kRadTo    1464         else // Case: On the border+t2<kRadTolerance
1588              //       (v is perpendicular to     1465              //       (v is perpendicular to the surface)
1589         {                                        1466         {
1590           if (calcNorm)                          1467           if (calcNorm)
1591           {                                      1468           {
1592             *n = G4ThreeVector(p.x()/fRMax,p.    1469             *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1593             *validNorm = true ;                  1470             *validNorm = true ;
1594           }                                      1471           }
1595           return snxt = 0.0;                     1472           return snxt = 0.0;
1596         }                                        1473         }
1597       }                                          1474       }
1598     }                                            1475     }
1599     // Phi Intersection                          1476     // Phi Intersection
1600                                                  1477 
1601     if ( !fPhiFullCutTube )                      1478     if ( !fPhiFullCutTube )
1602     {                                            1479     {
1603       // add angle calculation with correctio << 1480       // add angle calculation with correction 
1604       // of the difference in domain of atan2    1481       // of the difference in domain of atan2 and Sphi
1605       //                                         1482       //
1606       vphi = std::atan2(v.y(),v.x()) ;           1483       vphi = std::atan2(v.y(),v.x()) ;
1607                                               << 1484      
1608       if ( vphi < fSPhi - halfAngTolerance  )    1485       if ( vphi < fSPhi - halfAngTolerance  )             { vphi += twopi; }
1609       else if ( vphi > fSPhi + fDPhi + halfAn    1486       else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1610                                                  1487 
1611                                                  1488 
1612       if ( (p.x() != 0.0) || (p.y() != 0.0) ) << 1489       if ( p.x() || p.y() )  // Check if on z axis (rho not needed later)
1613       {                                          1490       {
1614         // pDist -ve when inside                 1491         // pDist -ve when inside
1615                                                  1492 
1616         pDistS = p.x()*sinSPhi - p.y()*cosSPh    1493         pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1617         pDistE = -p.x()*sinEPhi + p.y()*cosEP    1494         pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1618                                                  1495 
1619         // Comp -ve when in direction of outw    1496         // Comp -ve when in direction of outwards normal
1620                                                  1497 
1621         compS   = -sinSPhi*v.x() + cosSPhi*v.    1498         compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
1622         compE   =  sinEPhi*v.x() - cosEPhi*v.    1499         compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
1623                                               << 1500        
1624         sidephi = kNull;                         1501         sidephi = kNull;
1625                                               << 1502         
1626         if( ( (fDPhi <= pi) && ( (pDistS <= h    1503         if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1627                               && (pDistE <= h    1504                               && (pDistE <= halfCarTolerance) ) )
1628          || ( (fDPhi >  pi) && ((pDistS <=  h << 1505          || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
1629                               || (pDistE <=   << 1506                               && (pDistE >  halfCarTolerance) ) )  )
1630         {                                        1507         {
1631           // Inside both phi *full* planes       1508           // Inside both phi *full* planes
1632                                               << 1509           
1633           if ( compS < 0 )                       1510           if ( compS < 0 )
1634           {                                      1511           {
1635             sphi = pDistS/compS ;                1512             sphi = pDistS/compS ;
1636                                               << 1513             
1637             if (sphi >= -halfCarTolerance)       1514             if (sphi >= -halfCarTolerance)
1638             {                                    1515             {
1639               xi = p.x() + sphi*v.x() ;          1516               xi = p.x() + sphi*v.x() ;
1640               yi = p.y() + sphi*v.y() ;          1517               yi = p.y() + sphi*v.y() ;
1641                                               << 1518               
1642               // Check intersecting with corr    1519               // Check intersecting with correct half-plane
1643               // (if not -> no intersect)        1520               // (if not -> no intersect)
1644               //                                 1521               //
1645               if( (std::fabs(xi)<=kCarToleran    1522               if( (std::fabs(xi)<=kCarTolerance)
1646                && (std::fabs(yi)<=kCarToleran    1523                && (std::fabs(yi)<=kCarTolerance) )
1647               {                                  1524               {
1648                 sidephi = kSPhi;                 1525                 sidephi = kSPhi;
1649                 if (((fSPhi-halfAngTolerance)    1526                 if (((fSPhi-halfAngTolerance)<=vphi)
1650                    &&((fSPhi+fDPhi+halfAngTol    1527                    &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1651                 {                                1528                 {
1652                   sphi = kInfinity;              1529                   sphi = kInfinity;
1653                 }                                1530                 }
1654               }                                  1531               }
1655               else if ( yi*cosCPhi-xi*sinCPhi    1532               else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1656               {                                  1533               {
1657                 sphi = kInfinity ;               1534                 sphi = kInfinity ;
1658               }                                  1535               }
1659               else                               1536               else
1660               {                                  1537               {
1661                 sidephi = kSPhi ;                1538                 sidephi = kSPhi ;
1662                 if ( pDistS > -halfCarToleran    1539                 if ( pDistS > -halfCarTolerance )
1663                 {                                1540                 {
1664                   sphi = 0.0 ; // Leave by sp    1541                   sphi = 0.0 ; // Leave by sphi immediately
1665                 }                             << 1542                 }    
1666               }                               << 1543               }       
1667             }                                    1544             }
1668             else                                 1545             else
1669             {                                    1546             {
1670               sphi = kInfinity ;                 1547               sphi = kInfinity ;
1671             }                                    1548             }
1672           }                                      1549           }
1673           else                                   1550           else
1674           {                                      1551           {
1675             sphi = kInfinity ;                   1552             sphi = kInfinity ;
1676           }                                      1553           }
1677                                                  1554 
1678           if ( compE < 0 )                       1555           if ( compE < 0 )
1679           {                                      1556           {
1680             sphi2 = pDistE/compE ;               1557             sphi2 = pDistE/compE ;
1681                                               << 1558             
1682             // Only check further if < starti    1559             // Only check further if < starting phi intersection
1683             //                                   1560             //
1684             if ( (sphi2 > -halfCarTolerance)     1561             if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1685             {                                    1562             {
1686               xi = p.x() + sphi2*v.x() ;         1563               xi = p.x() + sphi2*v.x() ;
1687               yi = p.y() + sphi2*v.y() ;         1564               yi = p.y() + sphi2*v.y() ;
1688                                               << 1565               
1689               if ((std::fabs(xi)<=kCarToleran    1566               if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1690               {                                  1567               {
1691                 // Leaving via ending phi        1568                 // Leaving via ending phi
1692                 //                               1569                 //
1693                 if( (fSPhi-halfAngTolerance > << 1570                 if( !((fSPhi-halfAngTolerance <= vphi)
1694                      ||(fSPhi+fDPhi+halfAngTo << 1571                      &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1695                 {                                1572                 {
1696                   sidephi = kEPhi ;              1573                   sidephi = kEPhi ;
1697                   if ( pDistE <= -halfCarTole    1574                   if ( pDistE <= -halfCarTolerance )  { sphi = sphi2 ; }
1698                   else                           1575                   else                                { sphi = 0.0 ;   }
1699                 }                                1576                 }
1700               }                               << 1577               } 
1701               else    // Check intersecting w << 1578               else    // Check intersecting with correct half-plane 
1702                                                  1579 
1703               if ( (yi*cosCPhi-xi*sinCPhi) >=    1580               if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1704               {                                  1581               {
1705                 // Leaving via ending phi        1582                 // Leaving via ending phi
1706                 //                               1583                 //
1707                 sidephi = kEPhi ;                1584                 sidephi = kEPhi ;
1708                 if ( pDistE <= -halfCarTolera    1585                 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1709                 else                             1586                 else                               { sphi = 0.0 ;   }
1710               }                                  1587               }
1711             }                                    1588             }
1712           }                                      1589           }
1713         }                                        1590         }
1714         else                                     1591         else
1715         {                                        1592         {
1716           sphi = kInfinity ;                     1593           sphi = kInfinity ;
1717         }                                        1594         }
1718       }                                          1595       }
1719       else                                       1596       else
1720       {                                          1597       {
1721         // On z axis + travel not || to z axi    1598         // On z axis + travel not || to z axis -> if phi of vector direction
1722         // within phi of shape, Step limited     1599         // within phi of shape, Step limited by rmax, else Step =0
1723                                               << 1600                
1724         if ( (fSPhi - halfAngTolerance <= vph    1601         if ( (fSPhi - halfAngTolerance <= vphi)
1725            && (vphi <= fSPhi + fDPhi + halfAn    1602            && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1726         {                                        1603         {
1727           sphi = kInfinity ;                     1604           sphi = kInfinity ;
1728         }                                        1605         }
1729         else                                     1606         else
1730         {                                        1607         {
1731           sidephi = kSPhi ; // arbitrary      << 1608           sidephi = kSPhi ; // arbitrary 
1732           sphi    = 0.0 ;                        1609           sphi    = 0.0 ;
1733         }                                        1610         }
1734       }                                          1611       }
1735       if (sphi < snxt)  // Order intersecttio    1612       if (sphi < snxt)  // Order intersecttions
1736       {                                          1613       {
1737         snxt = sphi ;                            1614         snxt = sphi ;
1738         side = sidephi ;                         1615         side = sidephi ;
1739       }                                          1616       }
1740     }                                            1617     }
1741     if (srd < snxt)  // Order intersections      1618     if (srd < snxt)  // Order intersections
1742     {                                            1619     {
1743       snxt = srd ;                               1620       snxt = srd ;
1744       side = sider ;                             1621       side = sider ;
1745     }                                            1622     }
1746   }                                              1623   }
1747   if (calcNorm)                                  1624   if (calcNorm)
1748   {                                              1625   {
1749     switch(side)                                 1626     switch(side)
1750     {                                            1627     {
1751       case kRMax:                                1628       case kRMax:
1752         // Note: returned vector not normalis    1629         // Note: returned vector not normalised
1753         // (divide by fRMax for unit vector)     1630         // (divide by fRMax for unit vector)
1754         //                                       1631         //
1755         xi = p.x() + snxt*v.x() ;                1632         xi = p.x() + snxt*v.x() ;
1756         yi = p.y() + snxt*v.y() ;                1633         yi = p.y() + snxt*v.y() ;
1757         *n = G4ThreeVector(xi/fRMax,yi/fRMax,    1634         *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1758         *validNorm = true ;                      1635         *validNorm = true ;
1759         break ;                                  1636         break ;
1760                                                  1637 
1761       case kRMin:                                1638       case kRMin:
1762         *validNorm = false ;  // Rmin is inco    1639         *validNorm = false ;  // Rmin is inconvex
1763         break ;                                  1640         break ;
1764                                                  1641 
1765       case kSPhi:                                1642       case kSPhi:
1766         if ( fDPhi <= pi )                       1643         if ( fDPhi <= pi )
1767         {                                        1644         {
1768           *n         = G4ThreeVector(sinSPhi,    1645           *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1769           *validNorm = true ;                    1646           *validNorm = true ;
1770         }                                        1647         }
1771         else                                     1648         else
1772         {                                        1649         {
1773           *validNorm = false ;                   1650           *validNorm = false ;
1774         }                                        1651         }
1775         break ;                                  1652         break ;
1776                                                  1653 
1777       case kEPhi:                                1654       case kEPhi:
1778         if (fDPhi <= pi)                         1655         if (fDPhi <= pi)
1779         {                                        1656         {
1780           *n = G4ThreeVector(-sinEPhi,cosEPhi    1657           *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1781           *validNorm = true ;                    1658           *validNorm = true ;
1782         }                                        1659         }
1783         else                                     1660         else
1784         {                                        1661         {
1785           *validNorm = false ;                   1662           *validNorm = false ;
1786         }                                        1663         }
1787         break ;                                  1664         break ;
1788                                                  1665 
1789       case kPZ:                                  1666       case kPZ:
1790         *n         = fHighNorm ;                 1667         *n         = fHighNorm ;
1791         *validNorm = true ;                      1668         *validNorm = true ;
1792         break ;                                  1669         break ;
1793                                                  1670 
1794       case kMZ:                                  1671       case kMZ:
1795         *n         = fLowNorm ;                  1672         *n         = fLowNorm ;
1796         *validNorm = true ;                      1673         *validNorm = true ;
1797         break ;                                  1674         break ;
1798                                                  1675 
1799       default:                                   1676       default:
1800         G4cout << G4endl ;                       1677         G4cout << G4endl ;
1801         DumpInfo();                              1678         DumpInfo();
1802         std::ostringstream message;              1679         std::ostringstream message;
1803         G4long oldprc = message.precision(16) << 1680         G4int oldprc = message.precision(16);
1804         message << "Undefined side for valid     1681         message << "Undefined side for valid surface normal to solid."
1805                 << G4endl                        1682                 << G4endl
1806                 << "Position:"  << G4endl <<     1683                 << "Position:"  << G4endl << G4endl
1807                 << "p.x() = "   << p.x()/mm <    1684                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
1808                 << "p.y() = "   << p.y()/mm <    1685                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
1809                 << "p.z() = "   << p.z()/mm <    1686                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
1810                 << "Direction:" << G4endl <<     1687                 << "Direction:" << G4endl << G4endl
1811                 << "v.x() = "   << v.x() << G    1688                 << "v.x() = "   << v.x() << G4endl
1812                 << "v.y() = "   << v.y() << G    1689                 << "v.y() = "   << v.y() << G4endl
1813                 << "v.z() = "   << v.z() << G    1690                 << "v.z() = "   << v.z() << G4endl << G4endl
1814                 << "Proposed distance :" << G    1691                 << "Proposed distance :" << G4endl << G4endl
1815                 << "snxt = "    << snxt/mm <<    1692                 << "snxt = "    << snxt/mm << " mm" << G4endl ;
1816         message.precision(oldprc) ;              1693         message.precision(oldprc) ;
1817         G4Exception("G4CutTubs::DistanceToOut    1694         G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1818                     JustWarning, message);       1695                     JustWarning, message);
1819         break ;                                  1696         break ;
1820     }                                            1697     }
1821   }                                              1698   }
1822   if ( snxt<halfCarTolerance )  { snxt=0 ; }     1699   if ( snxt<halfCarTolerance )  { snxt=0 ; }
1823   return snxt ;                                  1700   return snxt ;
1824 }                                                1701 }
1825                                                  1702 
1826 /////////////////////////////////////////////    1703 //////////////////////////////////////////////////////////////////////////
1827 //                                               1704 //
1828 // Calculate distance (<=actual) to closest s    1705 // Calculate distance (<=actual) to closest surface of shape from inside
1829                                                  1706 
1830 G4double G4CutTubs::DistanceToOut( const G4Th    1707 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p ) const
1831 {                                                1708 {
1832   G4double safRMin,safRMax,safZLow,safZHigh,s    1709   G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1833   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);       1710   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1834                                                  1711 
1835   rho = std::sqrt(p.x()*p.x() + p.y()*p.y())     1712   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;  // Distance to R
1836                                                  1713 
1837   safRMin =  rho - fRMin ;                       1714   safRMin =  rho - fRMin ;
1838   safRMax =  fRMax - rho ;                       1715   safRMax =  fRMax - rho ;
1839                                                  1716 
1840   // Distances to ZCut(Low/High)                 1717   // Distances to ZCut(Low/High)
1841                                                  1718 
1842   // Dist to Low Cut                             1719   // Dist to Low Cut
1843   //                                             1720   //
1844   safZLow = std::fabs((p+vZ).dot(fLowNorm));     1721   safZLow = std::fabs((p+vZ).dot(fLowNorm));
1845                                                  1722 
1846   // Dist to High Cut                            1723   // Dist to High Cut
1847   //                                             1724   //
1848   safZHigh = std::fabs((p-vZ).dot(fHighNorm))    1725   safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1849   safe = std::min(safZLow,safZHigh);             1726   safe = std::min(safZLow,safZHigh);
1850                                                  1727 
1851   if ( safRMin < safe ) { safe = safRMin; }      1728   if ( safRMin < safe ) { safe = safRMin; }
1852   if ( safRMax< safe )  { safe = safRMax; }      1729   if ( safRMax< safe )  { safe = safRMax; }
1853                                                  1730 
1854   // Check if phi divided, Calc distances clo    1731   // Check if phi divided, Calc distances closest phi plane
1855   //                                             1732   //
1856   if ( !fPhiFullCutTube )                        1733   if ( !fPhiFullCutTube )
1857   {                                              1734   {
1858     if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )      1735     if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1859     {                                            1736     {
1860       safePhi = -(p.x()*sinSPhi - p.y()*cosSP    1737       safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1861     }                                            1738     }
1862     else                                         1739     else
1863     {                                            1740     {
1864       safePhi = (p.x()*sinEPhi - p.y()*cosEPh    1741       safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1865     }                                            1742     }
1866     if (safePhi < safe)  { safe = safePhi ; }    1743     if (safePhi < safe)  { safe = safePhi ; }
1867   }                                              1744   }
1868   if ( safe < 0 )  { safe = 0; }                 1745   if ( safe < 0 )  { safe = 0; }
1869                                                  1746 
1870   return safe ;                                  1747   return safe ;
1871 }                                                1748 }
1872                                                  1749 
                                                   >> 1750 /////////////////////////////////////////////////////////////////////////
                                                   >> 1751 //
                                                   >> 1752 // Create a List containing the transformed vertices
                                                   >> 1753 // Ordering [0-3] -fDz cross section
                                                   >> 1754 //          [4-7] +fDz cross section such that [0] is below [4],
                                                   >> 1755 //                                             [1] below [5] etc.
                                                   >> 1756 // Note:
                                                   >> 1757 //  Caller has deletion resposibility
                                                   >> 1758 
                                                   >> 1759 G4ThreeVectorList*
                                                   >> 1760 G4CutTubs::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
                                                   >> 1761 {
                                                   >> 1762   G4ThreeVectorList* vertices ;
                                                   >> 1763   G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
                                                   >> 1764   G4double meshAngle, meshRMax, crossAngle,
                                                   >> 1765            cosCrossAngle, sinCrossAngle, sAngle;
                                                   >> 1766   G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
                                                   >> 1767   G4int crossSection, noCrossSections;
                                                   >> 1768 
                                                   >> 1769   // Compute no of cross-sections necessary to mesh tube
                                                   >> 1770   //
                                                   >> 1771   noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
                                                   >> 1772 
                                                   >> 1773   if ( noCrossSections < kMinMeshSections )
                                                   >> 1774   {
                                                   >> 1775     noCrossSections = kMinMeshSections ;
                                                   >> 1776   }
                                                   >> 1777   else if (noCrossSections>kMaxMeshSections)
                                                   >> 1778   {
                                                   >> 1779     noCrossSections = kMaxMeshSections ;
                                                   >> 1780   }
                                                   >> 1781   // noCrossSections = 4 ;
                                                   >> 1782 
                                                   >> 1783   meshAngle = fDPhi/(noCrossSections - 1) ;
                                                   >> 1784   // meshAngle = fDPhi/(noCrossSections) ;
                                                   >> 1785 
                                                   >> 1786   meshRMax  = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
                                                   >> 1787   meshRMin = fRMin - 100*kCarTolerance ; 
                                                   >> 1788  
                                                   >> 1789   // If complete in phi, set start angle such that mesh will be at fRMax
                                                   >> 1790   // on the x axis. Will give better extent calculations when not rotated.
                                                   >> 1791 
                                                   >> 1792   if (fPhiFullCutTube && (fSPhi == 0) )  { sAngle = -meshAngle*0.5 ; }
                                                   >> 1793   else                                { sAngle =  fSPhi ; }
                                                   >> 1794     
                                                   >> 1795   vertices = new G4ThreeVectorList();
                                                   >> 1796     
                                                   >> 1797   if ( vertices )
                                                   >> 1798   {
                                                   >> 1799     vertices->reserve(noCrossSections*4);
                                                   >> 1800     for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
                                                   >> 1801     {
                                                   >> 1802       // Compute coordinates of cross section at section crossSection
                                                   >> 1803 
                                                   >> 1804       crossAngle    = sAngle + crossSection*meshAngle ;
                                                   >> 1805       cosCrossAngle = std::cos(crossAngle) ;
                                                   >> 1806       sinCrossAngle = std::sin(crossAngle) ;
                                                   >> 1807 
                                                   >> 1808       rMaxX = meshRMax*cosCrossAngle ;
                                                   >> 1809       rMaxY = meshRMax*sinCrossAngle ;
                                                   >> 1810 
                                                   >> 1811       if(meshRMin <= 0.0)
                                                   >> 1812       {
                                                   >> 1813         rMinX = 0.0 ;
                                                   >> 1814         rMinY = 0.0 ;
                                                   >> 1815       }
                                                   >> 1816       else
                                                   >> 1817       {
                                                   >> 1818         rMinX = meshRMin*cosCrossAngle ;
                                                   >> 1819         rMinY = meshRMin*sinCrossAngle ;
                                                   >> 1820       }
                                                   >> 1821       vertex0 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,-fDz))) ;
                                                   >> 1822       vertex1 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,-fDz))) ;
                                                   >> 1823       vertex2 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,+fDz))) ;
                                                   >> 1824       vertex3 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,+fDz))) ;
                                                   >> 1825 
                                                   >> 1826       vertices->push_back(pTransform.TransformPoint(vertex0)) ;
                                                   >> 1827       vertices->push_back(pTransform.TransformPoint(vertex1)) ;
                                                   >> 1828       vertices->push_back(pTransform.TransformPoint(vertex2)) ;
                                                   >> 1829       vertices->push_back(pTransform.TransformPoint(vertex3)) ;
                                                   >> 1830     }
                                                   >> 1831   }
                                                   >> 1832   else
                                                   >> 1833   {
                                                   >> 1834     DumpInfo();
                                                   >> 1835     G4Exception("G4CutTubs::CreateRotatedVertices()",
                                                   >> 1836                 "GeomSolids0003", FatalException,
                                                   >> 1837                 "Error in allocation of vertices. Out of memory !");
                                                   >> 1838   }
                                                   >> 1839   return vertices ;
                                                   >> 1840 }
                                                   >> 1841 
1873 /////////////////////////////////////////////    1842 //////////////////////////////////////////////////////////////////////////
1874 //                                               1843 //
1875 // Stream object contents to an output stream    1844 // Stream object contents to an output stream
1876                                                  1845 
1877 G4GeometryType G4CutTubs::GetEntityType() con    1846 G4GeometryType G4CutTubs::GetEntityType() const
1878 {                                                1847 {
1879   return {"G4CutTubs"};                       << 1848   return G4String("G4CutTubs");
1880 }                                                1849 }
1881                                                  1850 
1882 /////////////////////////////////////////////    1851 //////////////////////////////////////////////////////////////////////////
1883 //                                               1852 //
1884 // Make a clone of the object                    1853 // Make a clone of the object
1885 //                                               1854 //
1886 G4VSolid* G4CutTubs::Clone() const               1855 G4VSolid* G4CutTubs::Clone() const
1887 {                                                1856 {
1888   return new G4CutTubs(*this);                   1857   return new G4CutTubs(*this);
1889 }                                                1858 }
1890                                                  1859 
1891 /////////////////////////////////////////////    1860 //////////////////////////////////////////////////////////////////////////
1892 //                                               1861 //
1893 // Stream object contents to an output stream    1862 // Stream object contents to an output stream
1894                                                  1863 
1895 std::ostream& G4CutTubs::StreamInfo( std::ost    1864 std::ostream& G4CutTubs::StreamInfo( std::ostream& os ) const
1896 {                                                1865 {
1897   G4long oldprc = os.precision(16);           << 1866   G4int oldprc = os.precision(16);
1898   os << "------------------------------------    1867   os << "-----------------------------------------------------------\n"
1899      << "    *** Dump for solid - " << GetNam    1868      << "    *** Dump for solid - " << GetName() << " ***\n"
1900      << "    ================================    1869      << "    ===================================================\n"
1901      << " Solid type: G4CutTubs\n"               1870      << " Solid type: G4CutTubs\n"
1902      << " Parameters: \n"                        1871      << " Parameters: \n"
1903      << "    inner radius : " << fRMin/mm <<     1872      << "    inner radius : " << fRMin/mm << " mm \n"
1904      << "    outer radius : " << fRMax/mm <<     1873      << "    outer radius : " << fRMax/mm << " mm \n"
1905      << "    half length Z: " << fDz/mm << "     1874      << "    half length Z: " << fDz/mm << " mm \n"
1906      << "    starting phi : " << fSPhi/degree    1875      << "    starting phi : " << fSPhi/degree << " degrees \n"
1907      << "    delta phi    : " << fDPhi/degree    1876      << "    delta phi    : " << fDPhi/degree << " degrees \n"
1908      << "    low Norm     : " << fLowNorm     << 1877      << "    low Norm     : " << fLowNorm     << "  \n" 
1909      << "    high Norm    : "  <<fHighNorm       1878      << "    high Norm    : "  <<fHighNorm    << "  \n"
1910      << "------------------------------------    1879      << "-----------------------------------------------------------\n";
1911   os.precision(oldprc);                          1880   os.precision(oldprc);
1912                                                  1881 
1913   return os;                                     1882   return os;
1914 }                                                1883 }
1915                                                  1884 
1916 /////////////////////////////////////////////    1885 /////////////////////////////////////////////////////////////////////////
1917 //                                               1886 //
1918 // GetPointOnSurface                             1887 // GetPointOnSurface
1919                                                  1888 
1920 G4ThreeVector G4CutTubs::GetPointOnSurface()     1889 G4ThreeVector G4CutTubs::GetPointOnSurface() const
1921 {                                                1890 {
1922   // Set min and max z                        << 1891   G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1923   if (fZMin == 0. && fZMax == 0.)             << 1892            aOne, aTwo, aThr, aFou;
                                                   >> 1893   G4double rRand;
                                                   >> 1894  
                                                   >> 1895   aOne = 2.*fDz*fDPhi*fRMax;
                                                   >> 1896   aTwo = 2.*fDz*fDPhi*fRMin;
                                                   >> 1897   aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
                                                   >> 1898   aFou = 2.*fDz*(fRMax-fRMin);
                                                   >> 1899 
                                                   >> 1900   phi    = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
                                                   >> 1901   cosphi = std::cos(phi);
                                                   >> 1902   sinphi = std::sin(phi);
                                                   >> 1903 
                                                   >> 1904   rRand  = GetRadiusInRing(fRMin,fRMax);
                                                   >> 1905   
                                                   >> 1906   if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
                                                   >> 1907   
                                                   >> 1908   chose  = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
                                                   >> 1909 
                                                   >> 1910   if( (chose >=0) && (chose < aOne) )
                                                   >> 1911   {
                                                   >> 1912     xRand = fRMax*cosphi;
                                                   >> 1913     yRand = fRMax*sinphi;
                                                   >> 1914     zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
                                                   >> 1915                             GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
                                                   >> 1916     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1917   }
                                                   >> 1918   else if( (chose >= aOne) && (chose < aOne + aTwo) )
                                                   >> 1919   {
                                                   >> 1920     xRand = fRMin*cosphi;
                                                   >> 1921     yRand = fRMin*sinphi;
                                                   >> 1922     zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
                                                   >> 1923                             GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
                                                   >> 1924     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1925   }
                                                   >> 1926   else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
                                                   >> 1927   {
                                                   >> 1928     xRand = rRand*cosphi;
                                                   >> 1929     yRand = rRand*sinphi;
                                                   >> 1930     zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz));
                                                   >> 1931     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1932   }
                                                   >> 1933   else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
                                                   >> 1934   {
                                                   >> 1935     xRand = rRand*cosphi;
                                                   >> 1936     yRand = rRand*sinphi;
                                                   >> 1937     zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz));
                                                   >> 1938     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1939   }
                                                   >> 1940   else if( (chose >= aOne + aTwo + 2.*aThr)
                                                   >> 1941         && (chose < aOne + aTwo + 2.*aThr + aFou) )
                                                   >> 1942   {
                                                   >> 1943     xRand = rRand*std::cos(fSPhi);
                                                   >> 1944     yRand = rRand*std::sin(fSPhi);
                                                   >> 1945     zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
                                                   >> 1946                             GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
                                                   >> 1947     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1948   }
                                                   >> 1949   else
1924   {                                              1950   {
1925     G4AutoLock l(&zminmaxMutex);              << 1951     xRand = rRand*std::cos(fSPhi+fDPhi);
1926     G4ThreeVector bmin, bmax;                 << 1952     yRand = rRand*std::sin(fSPhi+fDPhi);
1927     BoundingLimits(bmin,bmax);                << 1953     zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1928     fZMin = bmin.z();                         << 1954                             GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1929     fZMax = bmax.z();                         << 1955     return G4ThreeVector  (xRand, yRand, zRand);
1930     l.unlock();                               << 1956   }
1931   }                                           << 
1932                                               << 
1933   // Set parameters                           << 
1934   G4double hmax = fZMax - fZMin;              << 
1935   G4double sphi = fSPhi;                      << 
1936   G4double dphi = fDPhi;                      << 
1937   G4double rmin = fRMin;                      << 
1938   G4double rmax = fRMax;                      << 
1939   G4double rrmax = rmax*rmax;                 << 
1940   G4double rrmin = rmin*rmin;                 << 
1941                                               << 
1942   G4ThreeVector nbot = GetLowNorm();          << 
1943   G4ThreeVector ntop = GetHighNorm();         << 
1944                                               << 
1945   // Set array of surface areas               << 
1946   G4double sbase = 0.5*dphi*(rrmax - rrmin);  << 
1947   G4double sbot = sbase/std::abs(nbot.z());   << 
1948   G4double stop = sbase/std::abs(ntop.z());   << 
1949   G4double scut = (dphi == twopi) ? 0. : hmax << 
1950   G4double ssurf[6] = { scut, scut, sbot, sto << 
1951   ssurf[1] += ssurf[0];                       << 
1952   ssurf[2] += ssurf[1];                       << 
1953   ssurf[3] += ssurf[2];                       << 
1954   ssurf[4] += ssurf[3];                       << 
1955   ssurf[5] += ssurf[4];                       << 
1956                                               << 
1957   constexpr G4int ntry = 100000;              << 
1958   for (G4int i=0; i<ntry; ++i)                << 
1959   {                                           << 
1960     // Select surface                         << 
1961     G4double select = ssurf[5]*G4QuickRand(); << 
1962     G4int k = 5;                              << 
1963     k -= (G4int)(select <= ssurf[4]);         << 
1964     k -= (G4int)(select <= ssurf[3]);         << 
1965     k -= (G4int)(select <= ssurf[2]);         << 
1966     k -= (G4int)(select <= ssurf[1]);         << 
1967     k -= (G4int)(select <= ssurf[0]);         << 
1968                                               << 
1969     // Generate point on selected surface (re << 
1970     G4ThreeVector p(0,0,0);                   << 
1971     switch(k)                                 << 
1972     {                                         << 
1973       case 0: // cut at start phi             << 
1974       {                                       << 
1975         G4double r = rmin + (rmax - rmin)*G4Q << 
1976         p.set(r*cosSPhi, r*sinSPhi, fZMin + h << 
1977         break;                                << 
1978       }                                       << 
1979       case 1: // cut at end phi               << 
1980       {                                       << 
1981         G4double r = rmin + (rmax - rmin)*G4Q << 
1982         p.set(r*cosEPhi, r*sinEPhi, fZMin + h << 
1983         break;                                << 
1984       }                                       << 
1985       case 2: // base at low z                << 
1986       {                                       << 
1987         G4double r = std::sqrt(rrmin + (rrmax << 
1988         G4double phi = sphi + dphi*G4QuickRan << 
1989         G4double x = r*std::cos(phi);         << 
1990         G4double y = r*std::sin(phi);         << 
1991         G4double z = -fDz - (x*nbot.x() + y*n << 
1992         return {x, y, z};                     << 
1993       }                                       << 
1994       case 3: // base at high z               << 
1995       {                                       << 
1996         G4double r = std::sqrt(rrmin + (rrmax << 
1997         G4double phi = sphi + dphi*G4QuickRan << 
1998         G4double x = r*std::cos(phi);         << 
1999         G4double y = r*std::sin(phi);         << 
2000         G4double z = fDz - (x*ntop.x() + y*nt << 
2001         return {x, y, z};                     << 
2002       }                                       << 
2003       case 4: // external lateral surface     << 
2004       {                                       << 
2005         G4double phi = sphi + dphi*G4QuickRan << 
2006         G4double z = fZMin + hmax*G4QuickRand << 
2007         G4double x = rmax*std::cos(phi);      << 
2008         G4double y = rmax*std::sin(phi);      << 
2009         p.set(x, y, z);                       << 
2010         break;                                << 
2011       }                                       << 
2012       case 5: // internal lateral surface     << 
2013       {                                       << 
2014         G4double phi = sphi + dphi*G4QuickRan << 
2015         G4double z = fZMin + hmax*G4QuickRand << 
2016         G4double x = rmin*std::cos(phi);      << 
2017         G4double y = rmin*std::sin(phi);      << 
2018         p.set(x, y, z);                       << 
2019         break;                                << 
2020       }                                       << 
2021     }                                         << 
2022     if ((ntop.dot(p) - fDz*ntop.z()) > 0.) co << 
2023     if ((nbot.dot(p) + fDz*nbot.z()) > 0.) co << 
2024     return p;                                 << 
2025   }                                           << 
2026   // Just in case, if all attempts to generat << 
2027   // Normally should never happen             << 
2028   G4double x = rmax*std::cos(sphi + 0.5*dphi) << 
2029   G4double y = rmax*std::sin(sphi + 0.5*dphi) << 
2030   G4double z = fDz - (x*ntop.x() + y*ntop.y() << 
2031   return {x, y, z};                           << 
2032 }                                                1957 }
2033                                                  1958 
2034 /////////////////////////////////////////////    1959 ///////////////////////////////////////////////////////////////////////////
2035 //                                               1960 //
2036 // Methods for visualisation                     1961 // Methods for visualisation
2037                                                  1962 
2038 void G4CutTubs::DescribeYourselfTo ( G4VGraph << 1963 void G4CutTubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
2039 {                                                1964 {
2040   scene.AddSolid (*this) ;                       1965   scene.AddSolid (*this) ;
2041 }                                                1966 }
2042                                                  1967 
2043 G4Polyhedron* G4CutTubs::CreatePolyhedron ()  << 1968 G4Polyhedron* G4CutTubs::CreatePolyhedron () const 
2044 {                                                1969 {
2045   typedef G4double G4double3[3];                 1970   typedef G4double G4double3[3];
2046   typedef G4int G4int4[4];                       1971   typedef G4int G4int4[4];
2047                                                  1972 
2048   auto ph = new G4Polyhedron;                 << 1973   G4Polyhedron *ph  = new G4Polyhedron;
2049   G4Polyhedron *ph1 = new G4PolyhedronTubs (f << 1974   G4Polyhedron *ph1 = G4OTubs::CreatePolyhedron();
2050   G4int nn=ph1->GetNoVertices();                 1975   G4int nn=ph1->GetNoVertices();
2051   G4int nf=ph1->GetNoFacets();                   1976   G4int nf=ph1->GetNoFacets();
2052   auto xyz = new G4double3[nn];  // number of << 1977   G4double3* xyz = new G4double3[nn];  // number of nodes 
2053   auto faces = new G4int4[nf] ;  // number of << 1978   G4int4*  faces = new G4int4[nf] ;    // number of faces
2054                                                  1979 
2055   for(G4int i=0; i<nn; ++i)                   << 1980   for(G4int i=0;i<nn;i++)
2056   {                                              1981   {
2057     xyz[i][0]=ph1->GetVertex(i+1).x();           1982     xyz[i][0]=ph1->GetVertex(i+1).x();
2058     xyz[i][1]=ph1->GetVertex(i+1).y();           1983     xyz[i][1]=ph1->GetVertex(i+1).y();
2059     G4double tmpZ=ph1->GetVertex(i+1).z();       1984     G4double tmpZ=ph1->GetVertex(i+1).z();
2060     if(tmpZ>=fDz-kCarTolerance)                  1985     if(tmpZ>=fDz-kCarTolerance)
2061     {                                            1986     {
2062       xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][    1987       xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
2063     }                                            1988     }
2064     else if(tmpZ<=-fDz+kCarTolerance)            1989     else if(tmpZ<=-fDz+kCarTolerance)
2065     {                                            1990     {
2066       xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][    1991       xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
2067     }                                            1992     }
2068     else                                         1993     else
2069     {                                            1994     {
2070       xyz[i][2]=tmpZ;                            1995       xyz[i][2]=tmpZ;
2071     }                                            1996     }
2072   }                                              1997   }
2073   G4int iNodes[4];                               1998   G4int iNodes[4];
2074   G4int* iEdge = nullptr;                     << 1999   G4int *iEdge=0;
2075   G4int n;                                       2000   G4int n;
2076   for(G4int i=0; i<nf ; ++i)                  << 2001   for(G4int i=0;i<nf;i++)
2077   {                                              2002   {
2078     ph1->GetFacet(i+1,n,iNodes,iEdge);           2003     ph1->GetFacet(i+1,n,iNodes,iEdge);
2079     for(G4int k=0; k<n; ++k)                  << 2004     for(G4int k=0;k<n;k++)
2080     {                                            2005     {
2081       faces[i][k]=iNodes[k];                     2006       faces[i][k]=iNodes[k];
2082     }                                            2007     }
2083     for(G4int k=n; k<4; ++k)                  << 2008     for(G4int k=n;k<4;k++)
2084     {                                            2009     {
2085       faces[i][k]=0;                             2010       faces[i][k]=0;
2086     }                                            2011     }
2087   }                                              2012   }
2088   ph->createPolyhedron(nn,nf,xyz,faces);         2013   ph->createPolyhedron(nn,nf,xyz,faces);
2089                                                  2014 
2090   delete [] xyz;                                 2015   delete [] xyz;
2091   delete [] faces;                               2016   delete [] faces;
2092   delete ph1;                                    2017   delete ph1;
2093                                                  2018 
2094   return ph;                                     2019   return ph;
2095 }                                                2020 }
2096                                                  2021 
2097 // Auxilary Methods for Solid                    2022 // Auxilary Methods for Solid
2098                                               << 2023  
2099 ///////////////////////////////////////////// << 2024 ///////////////////////////////////////////////////////////////////////////
2100 //                                            << 2025 // Return true if Cutted planes are crossing 
2101 // Check set of points on the outer lateral s << 2026 // Check Intersection Points on OX and OY axes
2102 // if the cut planes are crossing inside the  << 
2103 //                                            << 
2104                                                  2027 
2105 G4bool G4CutTubs::IsCrossingCutPlanes() const    2028 G4bool G4CutTubs::IsCrossingCutPlanes() const
2106 {                                                2029 {
2107   constexpr G4int npoints = 30;               << 2030   G4double zXLow1,zXLow2,zYLow1,zYLow2;
                                                   >> 2031   G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
                                                   >> 2032 
                                                   >> 2033   zXLow1  = GetCutZ(G4ThreeVector(-fRMax,     0,-fDz));
                                                   >> 2034   zXLow2  = GetCutZ(G4ThreeVector( fRMax,     0,-fDz));
                                                   >> 2035   zYLow1  = GetCutZ(G4ThreeVector(     0,-fRMax,-fDz));
                                                   >> 2036   zYLow2  = GetCutZ(G4ThreeVector(     0, fRMax,-fDz));
                                                   >> 2037   zXHigh1 = GetCutZ(G4ThreeVector(-fRMax,     0, fDz));
                                                   >> 2038   zXHigh2 = GetCutZ(G4ThreeVector( fRMax,     0, fDz));
                                                   >> 2039   zYHigh1 = GetCutZ(G4ThreeVector(     0,-fRMax, fDz));
                                                   >> 2040   zYHigh2 = GetCutZ(G4ThreeVector(     0, fRMax, fDz));
                                                   >> 2041   if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
                                                   >> 2042     || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2))  { return true; }
2108                                                  2043 
2109   // set values for calculation of h - distan << 
2110   // opposite points on bases                 << 
2111   G4ThreeVector nbot = GetLowNorm();          << 
2112   G4ThreeVector ntop = GetHighNorm();         << 
2113   if (std::abs(nbot.z()) < kCarTolerance) ret << 
2114   if (std::abs(ntop.z()) < kCarTolerance) ret << 
2115   G4double nx = nbot.x()/nbot.z() - ntop.x()/ << 
2116   G4double ny = nbot.y()/nbot.z() - ntop.y()/ << 
2117                                               << 
2118   // check points                             << 
2119   G4double cosphi = GetCosStartPhi();         << 
2120   G4double sinphi = GetSinStartPhi();         << 
2121   G4double delphi = GetDeltaPhiAngle()/npoint << 
2122   G4double cosdel = std::cos(delphi);         << 
2123   G4double sindel = std::sin(delphi);         << 
2124   G4double hzero = 2.*GetZHalfLength()/GetOut << 
2125   for (G4int i=0; i<npoints+1; ++i)           << 
2126   {                                           << 
2127     G4double h = nx*cosphi + ny*sinphi + hzer << 
2128     if (h < 0.) return true;                  << 
2129     G4double sintmp = sinphi;                 << 
2130     sinphi = sintmp*cosdel + cosphi*sindel;   << 
2131     cosphi = cosphi*cosdel - sintmp*sindel;   << 
2132   }                                           << 
2133   return false;                                  2044   return false;
2134 }                                                2045 }
2135                                                  2046 
2136 /////////////////////////////////////////////    2047 ///////////////////////////////////////////////////////////////////////////
2137 //                                               2048 //
2138 // Return real Z coordinate of point on Cutte    2049 // Return real Z coordinate of point on Cutted +/- fDZ plane
2139                                                  2050 
2140 G4double G4CutTubs::GetCutZ(const G4ThreeVect    2051 G4double G4CutTubs::GetCutZ(const G4ThreeVector& p) const
2141 {                                                2052 {
2142   G4double newz = p.z();  // p.z() should be     2053   G4double newz = p.z();  // p.z() should be either +fDz or -fDz
2143   if (p.z()<0)                                   2054   if (p.z()<0)
2144   {                                              2055   {
2145     if(fLowNorm.z()!=0.)                         2056     if(fLowNorm.z()!=0.)
2146     {                                            2057     {
2147        newz = -fDz-(p.x()*fLowNorm.x()+p.y()*    2058        newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
2148     }                                            2059     }
2149   }                                              2060   }
2150   else                                           2061   else
2151   {                                              2062   {
2152     if(fHighNorm.z()!=0.)                        2063     if(fHighNorm.z()!=0.)
2153     {                                            2064     {
2154        newz = fDz-(p.x()*fHighNorm.x()+p.y()*    2065        newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
2155     }                                            2066     }
2156   }                                              2067   }
2157   return newz;                                   2068   return newz;
2158 }                                                2069 }
2159 #endif                                        << 2070 
                                                   >> 2071 ///////////////////////////////////////////////////////////////////////////
                                                   >> 2072 //
                                                   >> 2073 // Calculate Min and Max Z for CutZ
                                                   >> 2074 
                                                   >> 2075 void G4CutTubs::GetMaxMinZ(G4double& zmin,G4double& zmax)const
                                                   >> 2076 
                                                   >> 2077 {
                                                   >> 2078   G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x());
                                                   >> 2079   G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x());
                                                   >> 2080 
                                                   >> 2081   G4double xc=0, yc=0,z1;
                                                   >> 2082   G4double z[8];
                                                   >> 2083   G4bool in_range_low = false;
                                                   >> 2084   G4bool in_range_hi = false;
                                                   >> 2085  
                                                   >> 2086   G4int i;
                                                   >> 2087   for (i=0; i<2; i++)
                                                   >> 2088   {
                                                   >> 2089     if (phiLow<0)  { phiLow+=twopi; }
                                                   >> 2090     G4double ddp = phiLow-fSPhi;
                                                   >> 2091     if (ddp<0)  { ddp += twopi; }
                                                   >> 2092     if (ddp <= fDPhi)
                                                   >> 2093     {
                                                   >> 2094       xc = fRMin*std::cos(phiLow);
                                                   >> 2095       yc = fRMin*std::sin(phiLow);
                                                   >> 2096       z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz));
                                                   >> 2097       xc = fRMax*std::cos(phiLow);
                                                   >> 2098       yc = fRMax*std::sin(phiLow);
                                                   >> 2099       z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz)));
                                                   >> 2100       if (in_range_low)  { zmin = std::min(zmin, z1); }
                                                   >> 2101       else               { zmin = z1; }
                                                   >> 2102       in_range_low = true;
                                                   >> 2103     }
                                                   >> 2104     phiLow += pi;
                                                   >> 2105     if (phiLow>twopi)  { phiLow-=twopi; }
                                                   >> 2106   }
                                                   >> 2107   for (i=0; i<2; i++)
                                                   >> 2108   {
                                                   >> 2109     if (phiHigh<0)  { phiHigh+=twopi; }
                                                   >> 2110     G4double ddp = phiHigh-fSPhi;
                                                   >> 2111     if (ddp<0)  { ddp += twopi; }
                                                   >> 2112     if (ddp <= fDPhi)
                                                   >> 2113     {
                                                   >> 2114       xc = fRMin*std::cos(phiHigh);
                                                   >> 2115       yc = fRMin*std::sin(phiHigh);
                                                   >> 2116       z1 = GetCutZ(G4ThreeVector(xc, yc, fDz));
                                                   >> 2117       xc = fRMax*std::cos(phiHigh);
                                                   >> 2118       yc = fRMax*std::sin(phiHigh);
                                                   >> 2119       z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz)));
                                                   >> 2120       if (in_range_hi)  { zmax = std::min(zmax, z1); }
                                                   >> 2121       else              { zmax = z1; }
                                                   >> 2122       in_range_hi = true;
                                                   >> 2123     }
                                                   >> 2124     phiHigh += pi;
                                                   >> 2125     if (phiLow>twopi)  { phiHigh-=twopi; }
                                                   >> 2126   }
                                                   >> 2127 
                                                   >> 2128   xc = fRMin*std::cos(fSPhi);
                                                   >> 2129   yc = fRMin*std::sin(fSPhi);
                                                   >> 2130   z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
                                                   >> 2131   z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz));
                                                   >> 2132  
                                                   >> 2133   xc = fRMin*std::cos(fSPhi+fDPhi);
                                                   >> 2134   yc = fRMin*std::sin(fSPhi+fDPhi);
                                                   >> 2135   z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
                                                   >> 2136   z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz));
                                                   >> 2137  
                                                   >> 2138   xc = fRMax*std::cos(fSPhi);
                                                   >> 2139   yc = fRMax*std::sin(fSPhi);
                                                   >> 2140   z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
                                                   >> 2141   z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz));
                                                   >> 2142  
                                                   >> 2143   xc = fRMax*std::cos(fSPhi+fDPhi);
                                                   >> 2144   yc = fRMax*std::sin(fSPhi+fDPhi);
                                                   >> 2145   z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
                                                   >> 2146   z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz));
                                                   >> 2147  
                                                   >> 2148   // Find min/max
                                                   >> 2149 
                                                   >> 2150   z1=z[0];
                                                   >> 2151   for (i = 1; i < 4; i++)
                                                   >> 2152   {
                                                   >> 2153     if(z[i] < z[i-1])z1=z[i];
                                                   >> 2154   }
                                                   >> 2155     
                                                   >> 2156   if (in_range_low)
                                                   >> 2157   {
                                                   >> 2158     zmin = std::min(zmin, z1);
                                                   >> 2159   }
                                                   >> 2160   else
                                                   >> 2161   {
                                                   >> 2162     zmin = z1;
                                                   >> 2163   }
                                                   >> 2164   z1=z[4];
                                                   >> 2165   for (i = 1; i < 4; i++)
                                                   >> 2166   {
                                                   >> 2167     if(z[4+i] > z[4+i-1])  { z1=z[4+i]; }
                                                   >> 2168   }
                                                   >> 2169 
                                                   >> 2170   if (in_range_hi)  { zmax = std::max(zmax, z1); }
                                                   >> 2171   else              { zmax = z1; }
                                                   >> 2172 }
2160                                                  2173