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 9.5)


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