Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Tubs.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/G4Tubs.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Tubs.cc (Version 9.6.p1)


  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 // G4Tubs implementation                       << 
 27 //                                                 26 //
 28 // 1994-95 P.Kent: first implementation        <<  27 // $Id$
 29 // 08.08.00 V.Grichine: more stable roots of 2 <<  28 //
                                                   >>  29 // 
                                                   >>  30 // class G4Tubs
                                                   >>  31 //
                                                   >>  32 // History:
                                                   >>  33 //
                                                   >>  34 // 05.04.12 M.Kelsey:   Use sqrt(r) in GetPointOnSurface() for uniform points
                                                   >>  35 // 02.08.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for negative value under sqrt
                                                   >>  36 //                      for the case: p on the surface and v is tangent to the surface
                                                   >>  37 // 11.05.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for phi < 2pi
                                                   >>  38 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
                                                   >>  39 // 16.03.05 V.Grichine: SurfaceNormal(p) with edges/corners for boolean
                                                   >>  40 // 20.07.01 V.Grichine: bug fixed in Inside(p)
                                                   >>  41 // 20.02.01 V.Grichine: bug fixed in Inside(p) and CalculateExtent was 
                                                   >>  42 //                      simplified base on G4Box::CalculateExtent
 30 // 07.12.00 V.Grichine: phi-section algorithm      43 // 07.12.00 V.Grichine: phi-section algorithm was changed in Inside(p)
 31 // 03.05.05 V.Grichine: SurfaceNormal(p) accor <<  44 // 28.11.00 V.Grichine: bug fixed in Inside(p)
 32 // 24.08.16 E.Tcherniaev: reimplemented Calcul <<  45 // 31.10.00 V.Grichine: assign srd, sphi in Distance ToOut(p,v,...)
 33 // ------------------------------------------- <<  46 // 08.08.00 V.Grichine: more stable roots of 2-equation in DistanceToOut(p,v,..)
                                                   >>  47 // 02.08.00 V.Grichine: point is outside check in Distance ToOut(p)
                                                   >>  48 // 17.05.00 V.Grichine: bugs (#76,#91) fixed in Distance ToOut(p,v,...)
                                                   >>  49 // 31.03.00 V.Grichine: bug fixed in Inside(p)
                                                   >>  50 // 19.11.99 V.Grichine: side = kNull in DistanceToOut(p,v,...)
                                                   >>  51 // 13.10.99 V.Grichine: bugs fixed in DistanceToIn(p,v) 
                                                   >>  52 // 28.05.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
                                                   >>  53 // 25.05.99 V.Grichine: bugs fixed in DistanceToIn(p,v) 
                                                   >>  54 // 23.03.99 V.Grichine: bug fixed in DistanceToIn(p,v) 
                                                   >>  55 // 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
                                                   >>  56 // 18.06.98 V.Grichine: n-normalisation in DistanceToOut(p,v)
                                                   >>  57 // 
                                                   >>  58 // 1994-95  P.Kent:     implementation
                                                   >>  59 //
                                                   >>  60 /////////////////////////////////////////////////////////////////////////
 34                                                    61 
 35 #include "G4Tubs.hh"                               62 #include "G4Tubs.hh"
 36                                                    63 
 37 #if !defined(G4GEOM_USE_UTUBS)                 << 
 38                                                << 
 39 #include "G4GeomTools.hh"                      << 
 40 #include "G4VoxelLimits.hh"                        64 #include "G4VoxelLimits.hh"
 41 #include "G4AffineTransform.hh"                    65 #include "G4AffineTransform.hh"
 42 #include "G4GeometryTolerance.hh"                  66 #include "G4GeometryTolerance.hh"
 43 #include "G4BoundingEnvelope.hh"               << 
 44                                                    67 
 45 #include "G4VPVParameterisation.hh"                68 #include "G4VPVParameterisation.hh"
 46 #include "G4QuickRand.hh"                      <<  69 
                                                   >>  70 #include "Randomize.hh"
                                                   >>  71 
                                                   >>  72 #include "meshdefs.hh"
 47                                                    73 
 48 #include "G4VGraphicsScene.hh"                     74 #include "G4VGraphicsScene.hh"
 49 #include "G4Polyhedron.hh"                         75 #include "G4Polyhedron.hh"
                                                   >>  76 #include "G4NURBS.hh"
                                                   >>  77 #include "G4NURBStube.hh"
                                                   >>  78 #include "G4NURBScylinder.hh"
                                                   >>  79 #include "G4NURBStubesector.hh"
 50                                                    80 
 51 using namespace CLHEP;                             81 using namespace CLHEP;
 52                                                    82 
 53 //////////////////////////////////////////////     83 /////////////////////////////////////////////////////////////////////////
 54 //                                                 84 //
 55 // Constructor - check parameters, convert ang     85 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
 56 //             - note if pdphi>2PI then reset      86 //             - note if pdphi>2PI then reset to 2PI
 57                                                    87 
 58 G4Tubs::G4Tubs( const G4String& pName,         <<  88 G4Tubs::G4Tubs( const G4String &pName,
 59                       G4double pRMin, G4double     89                       G4double pRMin, G4double pRMax,
 60                       G4double pDz,                90                       G4double pDz,
 61                       G4double pSPhi, G4double     91                       G4double pSPhi, G4double pDPhi )
 62    : G4CSGSolid(pName), fRMin(pRMin), fRMax(pR <<  92   : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
 63      fSPhi(0), fDPhi(0),                       << 
 64      fInvRmax( pRMax > 0.0 ? 1.0/pRMax : 0.0 ) << 
 65      fInvRmin( pRMin > 0.0 ? 1.0/pRMin : 0.0 ) << 
 66 {                                                  93 {
                                                   >>  94 
 67   kRadTolerance = G4GeometryTolerance::GetInst     95   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
 68   kAngTolerance = G4GeometryTolerance::GetInst     96   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 69                                                    97 
 70   halfCarTolerance=kCarTolerance*0.5;          << 
 71   halfRadTolerance=kRadTolerance*0.5;          << 
 72   halfAngTolerance=kAngTolerance*0.5;          << 
 73                                                << 
 74   if (pDz<=0) // Check z-len                       98   if (pDz<=0) // Check z-len
 75   {                                                99   {
 76     std::ostringstream message;                   100     std::ostringstream message;
 77     message << "Negative Z half-length (" << p    101     message << "Negative Z half-length (" << pDz << ") in solid: " << GetName();
 78     G4Exception("G4Tubs::G4Tubs()", "GeomSolid    102     G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
 79   }                                               103   }
 80   if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch    104   if ( (pRMin >= pRMax) || (pRMin < 0) ) // Check radii
 81   {                                               105   {
 82     std::ostringstream message;                   106     std::ostringstream message;
 83     message << "Invalid values for radii in so    107     message << "Invalid values for radii in solid: " << GetName()
 84             << G4endl                             108             << G4endl
 85             << "        pRMin = " << pRMin <<     109             << "        pRMin = " << pRMin << ", pRMax = " << pRMax;
 86     G4Exception("G4Tubs::G4Tubs()", "GeomSolid    110     G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
 87   }                                               111   }
 88                                                   112 
 89   // Check angles                                 113   // Check angles
 90   //                                              114   //
 91   CheckPhiAngles(pSPhi, pDPhi);                   115   CheckPhiAngles(pSPhi, pDPhi);
 92 }                                                 116 }
 93                                                   117 
 94 //////////////////////////////////////////////    118 ///////////////////////////////////////////////////////////////////////
 95 //                                                119 //
 96 // Fake default constructor - sets only member    120 // Fake default constructor - sets only member data and allocates memory
 97 //                            for usage restri    121 //                            for usage restricted to object persistency.
 98 //                                                122 //
 99 G4Tubs::G4Tubs( __void__& a )                     123 G4Tubs::G4Tubs( __void__& a )
100   : G4CSGSolid(a)                              << 124   : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
                                                   >> 125     fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
                                                   >> 126     sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
                                                   >> 127     sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
                                                   >> 128     fPhiFullTube(false)
101 {                                                 129 {
102 }                                                 130 }
103                                                   131 
104 //////////////////////////////////////////////    132 //////////////////////////////////////////////////////////////////////////
105 //                                                133 //
106 // Destructor                                     134 // Destructor
107                                                   135 
108 G4Tubs::~G4Tubs() = default;                   << 136 G4Tubs::~G4Tubs()
                                                   >> 137 {
                                                   >> 138 }
109                                                   139 
110 //////////////////////////////////////////////    140 //////////////////////////////////////////////////////////////////////////
111 //                                                141 //
112 // Copy constructor                               142 // Copy constructor
113                                                   143 
114 G4Tubs::G4Tubs(const G4Tubs&) = default;       << 144 G4Tubs::G4Tubs(const G4Tubs& rhs)
                                                   >> 145   : G4CSGSolid(rhs),
                                                   >> 146     kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
                                                   >> 147     fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
                                                   >> 148     fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
                                                   >> 149     sinCPhi(rhs.sinCPhi), cosCPhi(rhs.sinCPhi),
                                                   >> 150     cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiOT),
                                                   >> 151     sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
                                                   >> 152     sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube)
                                                   >> 153 {
                                                   >> 154 }
115                                                   155 
116 //////////////////////////////////////////////    156 //////////////////////////////////////////////////////////////////////////
117 //                                                157 //
118 // Assignment operator                            158 // Assignment operator
119                                                   159 
120 G4Tubs& G4Tubs::operator = (const G4Tubs& rhs) << 160 G4Tubs& G4Tubs::operator = (const G4Tubs& rhs) 
121 {                                                 161 {
122    // Check assignment to self                    162    // Check assignment to self
123    //                                             163    //
124    if (this == &rhs)  { return *this; }           164    if (this == &rhs)  { return *this; }
125                                                   165 
126    // Copy base class data                        166    // Copy base class data
127    //                                             167    //
128    G4CSGSolid::operator=(rhs);                    168    G4CSGSolid::operator=(rhs);
129                                                   169 
130    // Copy data                                   170    // Copy data
131    //                                             171    //
132    kRadTolerance = rhs.kRadTolerance; kAngTole    172    kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance;
133    fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz =    173    fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = rhs.fDz;
134    fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;          174    fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
135    sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh << 175    sinCPhi = rhs.sinCPhi; cosCPhi = rhs.sinCPhi;
136    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r << 176    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiOT;
137    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh    177    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
138    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh    178    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
139    fPhiFullTube = rhs.fPhiFullTube;               179    fPhiFullTube = rhs.fPhiFullTube;
140    fInvRmax = rhs.fInvRmax;                    << 
141    fInvRmin = rhs.fInvRmin;                    << 
142    halfCarTolerance = rhs.halfCarTolerance;    << 
143    halfRadTolerance = rhs.halfRadTolerance;    << 
144    halfAngTolerance = rhs.halfAngTolerance;    << 
145                                                   180 
146    return *this;                                  181    return *this;
147 }                                                 182 }
148                                                   183 
149 //////////////////////////////////////////////    184 /////////////////////////////////////////////////////////////////////////
150 //                                                185 //
151 // Dispatch to parameterisation for replicatio    186 // Dispatch to parameterisation for replication mechanism dimension
152 // computation & modification.                    187 // computation & modification.
153                                                   188 
154 void G4Tubs::ComputeDimensions(       G4VPVPar    189 void G4Tubs::ComputeDimensions(       G4VPVParameterisation* p,
155                                 const G4int n,    190                                 const G4int n,
156                                 const G4VPhysi    191                                 const G4VPhysicalVolume* pRep )
157 {                                                 192 {
158   p->ComputeDimensions(*this,n,pRep) ;            193   p->ComputeDimensions(*this,n,pRep) ;
159 }                                                 194 }
160                                                   195 
161 ////////////////////////////////////////////// << 196 ////////////////////////////////////////////////////////////////////////
162 //                                             << 
163 // Get bounding box                            << 
164                                                << 
165 void G4Tubs::BoundingLimits(G4ThreeVector& pMi << 
166 {                                              << 
167   G4double rmin = GetInnerRadius();            << 
168   G4double rmax = GetOuterRadius();            << 
169   G4double dz   = GetZHalfLength();            << 
170                                                << 
171   // Find bounding box                         << 
172   //                                           << 
173   if (GetDeltaPhiAngle() < twopi)              << 
174   {                                            << 
175     G4TwoVector vmin,vmax;                     << 
176     G4GeomTools::DiskExtent(rmin,rmax,         << 
177                             GetSinStartPhi(),G << 
178                             GetSinEndPhi(),Get << 
179                             vmin,vmax);        << 
180     pMin.set(vmin.x(),vmin.y(),-dz);           << 
181     pMax.set(vmax.x(),vmax.y(), dz);           << 
182   }                                            << 
183   else                                         << 
184   {                                            << 
185     pMin.set(-rmax,-rmax,-dz);                 << 
186     pMax.set( rmax, rmax, dz);                 << 
187   }                                            << 
188                                                << 
189   // Check correctness of the bounding box     << 
190   //                                           << 
191   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
192   {                                            << 
193     std::ostringstream message;                << 
194     message << "Bad bounding box (min >= max)  << 
195             << GetName() << " !"               << 
196             << "\npMin = " << pMin             << 
197             << "\npMax = " << pMax;            << 
198     G4Exception("G4Tubs::BoundingLimits()", "G << 
199                 JustWarning, message);         << 
200     DumpInfo();                                << 
201   }                                            << 
202 }                                              << 
203                                                << 
204 ////////////////////////////////////////////// << 
205 //                                                197 //
206 // Calculate extent under transform and specif    198 // Calculate extent under transform and specified limit
207                                                   199 
208 G4bool G4Tubs::CalculateExtent( const EAxis       200 G4bool G4Tubs::CalculateExtent( const EAxis              pAxis,
209                                 const G4VoxelL    201                                 const G4VoxelLimits&     pVoxelLimit,
210                                 const G4Affine    202                                 const G4AffineTransform& pTransform,
211                                       G4double << 203                                       G4double&          pMin, 
212                                       G4double    204                                       G4double&          pMax    ) const
213 {                                                 205 {
214   G4ThreeVector bmin, bmax;                    << 
215   G4bool exist;                                << 
216                                                << 
217   // Get bounding box                          << 
218   BoundingLimits(bmin,bmax);                   << 
219                                                   206 
220   // Check bounding box                        << 207   if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
221   G4BoundingEnvelope bbox(bmin,bmax);          << 
222 #ifdef G4BBOX_EXTENT                           << 
223   return bbox.CalculateExtent(pAxis,pVoxelLimi << 
224 #endif                                         << 
225   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 
226   {                                               208   {
227     return exist = pMin < pMax;                << 209     // Special case handling for unrotated solid tubes
228   }                                            << 210     // Compute x/y/z mins and maxs fro bounding box respecting limits,
                                                   >> 211     // with early returns if outside limits. Then switch() on pAxis,
                                                   >> 212     // and compute exact x and y limit for x/y case
                                                   >> 213       
                                                   >> 214     G4double xoffset, xMin, xMax;
                                                   >> 215     G4double yoffset, yMin, yMax;
                                                   >> 216     G4double zoffset, zMin, zMax;
                                                   >> 217 
                                                   >> 218     G4double diff1, diff2, maxDiff, newMin, newMax;
                                                   >> 219     G4double xoff1, xoff2, yoff1, yoff2, delta;
                                                   >> 220 
                                                   >> 221     xoffset = pTransform.NetTranslation().x();
                                                   >> 222     xMin = xoffset - fRMax;
                                                   >> 223     xMax = xoffset + fRMax;
229                                                   224 
230   // Get parameters of the solid               << 225     if (pVoxelLimit.IsXLimited())
231   G4double rmin = GetInnerRadius();            << 226     {
232   G4double rmax = GetOuterRadius();            << 227       if ( (xMin > pVoxelLimit.GetMaxXExtent())
233   G4double dz   = GetZHalfLength();            << 228         || (xMax < pVoxelLimit.GetMinXExtent()) )
234   G4double dphi = GetDeltaPhiAngle();          << 229       {
                                                   >> 230         return false;
                                                   >> 231       }
                                                   >> 232       else
                                                   >> 233       {
                                                   >> 234         if (xMin < pVoxelLimit.GetMinXExtent())
                                                   >> 235         {
                                                   >> 236           xMin = pVoxelLimit.GetMinXExtent();
                                                   >> 237         }
                                                   >> 238         if (xMax > pVoxelLimit.GetMaxXExtent())
                                                   >> 239         {
                                                   >> 240           xMax = pVoxelLimit.GetMaxXExtent();
                                                   >> 241         }
                                                   >> 242       }
                                                   >> 243     }
                                                   >> 244     yoffset = pTransform.NetTranslation().y();
                                                   >> 245     yMin    = yoffset - fRMax;
                                                   >> 246     yMax    = yoffset + fRMax;
235                                                   247 
236   // Find bounding envelope and calculate exte << 248     if ( pVoxelLimit.IsYLimited() )
237   //                                           << 249     {
238   const G4int NSTEPS = 24;            // numbe << 250       if ( (yMin > pVoxelLimit.GetMaxYExtent())
239   G4double astep  = twopi/NSTEPS;     // max a << 251         || (yMax < pVoxelLimit.GetMinYExtent()) )
240   G4int    ksteps = (dphi <= astep) ? 1 : (G4i << 252       {
241   G4double ang    = dphi/ksteps;               << 253         return false;
242                                                << 254       }
243   G4double sinHalf = std::sin(0.5*ang);        << 255       else
244   G4double cosHalf = std::cos(0.5*ang);        << 256       {
245   G4double sinStep = 2.*sinHalf*cosHalf;       << 257         if (yMin < pVoxelLimit.GetMinYExtent())
246   G4double cosStep = 1. - 2.*sinHalf*sinHalf;  << 258         {
247   G4double rext    = rmax/cosHalf;             << 259           yMin = pVoxelLimit.GetMinYExtent();
248                                                << 260         }
249   // bounding envelope for full cylinder consi << 261         if (yMax > pVoxelLimit.GetMaxYExtent())
250   // in other cases it is a sequence of quadri << 262         {
251   if (rmin == 0 && dphi == twopi)              << 263           yMax=pVoxelLimit.GetMaxYExtent();
252   {                                            << 264         }
253     G4double sinCur = sinHalf;                 << 265       }
254     G4double cosCur = cosHalf;                 << 266     }
255                                                << 267     zoffset = pTransform.NetTranslation().z();
256     G4ThreeVectorList baseA(NSTEPS),baseB(NSTE << 268     zMin    = zoffset - fDz;
257     for (G4int k=0; k<NSTEPS; ++k)             << 269     zMax    = zoffset + fDz;
258     {                                          << 270 
259       baseA[k].set(rext*cosCur,rext*sinCur,-dz << 271     if ( pVoxelLimit.IsZLimited() )
260       baseB[k].set(rext*cosCur,rext*sinCur, dz << 272     {
261                                                << 273       if ( (zMin > pVoxelLimit.GetMaxZExtent())
262       G4double sinTmp = sinCur;                << 274         || (zMax < pVoxelLimit.GetMinZExtent()) )
263       sinCur = sinCur*cosStep + cosCur*sinStep << 275       {
264       cosCur = cosCur*cosStep - sinTmp*sinStep << 276         return false;
265     }                                          << 277       }
266     std::vector<const G4ThreeVectorList *> pol << 278       else
267     polygons[0] = &baseA;                      << 279       {
268     polygons[1] = &baseB;                      << 280         if (zMin < pVoxelLimit.GetMinZExtent())
269     G4BoundingEnvelope benv(bmin,bmax,polygons << 281         {
270     exist = benv.CalculateExtent(pAxis,pVoxelL << 282           zMin = pVoxelLimit.GetMinZExtent();
                                                   >> 283         }
                                                   >> 284         if (zMax > pVoxelLimit.GetMaxZExtent())
                                                   >> 285         {
                                                   >> 286           zMax = pVoxelLimit.GetMaxZExtent();
                                                   >> 287         }
                                                   >> 288       }
                                                   >> 289     }
                                                   >> 290     switch ( pAxis )  // Known to cut cylinder
                                                   >> 291     {
                                                   >> 292       case kXAxis :
                                                   >> 293       {
                                                   >> 294         yoff1 = yoffset - yMin;
                                                   >> 295         yoff2 = yMax    - yoffset;
                                                   >> 296 
                                                   >> 297         if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
                                                   >> 298         {                                   // => no change
                                                   >> 299           pMin = xMin;
                                                   >> 300           pMax = xMax;
                                                   >> 301         }
                                                   >> 302         else
                                                   >> 303         {
                                                   >> 304           // Y limits don't cross max/min x => compute max delta x,
                                                   >> 305           // hence new mins/maxs
                                                   >> 306 
                                                   >> 307           delta   = fRMax*fRMax - yoff1*yoff1;
                                                   >> 308           diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
                                                   >> 309           delta   = fRMax*fRMax - yoff2*yoff2;
                                                   >> 310           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
                                                   >> 311           maxDiff = (diff1 > diff2) ? diff1:diff2;
                                                   >> 312           newMin  = xoffset - maxDiff;
                                                   >> 313           newMax  = xoffset + maxDiff;
                                                   >> 314           pMin    = (newMin < xMin) ? xMin : newMin;
                                                   >> 315           pMax    = (newMax > xMax) ? xMax : newMax;
                                                   >> 316         }    
                                                   >> 317         break;
                                                   >> 318       }
                                                   >> 319       case kYAxis :
                                                   >> 320       {
                                                   >> 321         xoff1 = xoffset - xMin;
                                                   >> 322         xoff2 = xMax - xoffset;
                                                   >> 323 
                                                   >> 324         if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
                                                   >> 325         {                                   // => no change
                                                   >> 326           pMin = yMin;
                                                   >> 327           pMax = yMax;
                                                   >> 328         }
                                                   >> 329         else
                                                   >> 330         {
                                                   >> 331           // X limits don't cross max/min y => compute max delta y,
                                                   >> 332           // hence new mins/maxs
                                                   >> 333 
                                                   >> 334           delta   = fRMax*fRMax - xoff1*xoff1;
                                                   >> 335           diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
                                                   >> 336           delta   = fRMax*fRMax - xoff2*xoff2;
                                                   >> 337           diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
                                                   >> 338           maxDiff = (diff1 > diff2) ? diff1 : diff2;
                                                   >> 339           newMin  = yoffset - maxDiff;
                                                   >> 340           newMax  = yoffset + maxDiff;
                                                   >> 341           pMin    = (newMin < yMin) ? yMin : newMin;
                                                   >> 342           pMax    = (newMax > yMax) ? yMax : newMax;
                                                   >> 343         }
                                                   >> 344         break;
                                                   >> 345       }
                                                   >> 346       case kZAxis:
                                                   >> 347       {
                                                   >> 348         pMin = zMin;
                                                   >> 349         pMax = zMax;
                                                   >> 350         break;
                                                   >> 351       }
                                                   >> 352       default:
                                                   >> 353         break;
                                                   >> 354     }
                                                   >> 355     pMin -= kCarTolerance;
                                                   >> 356     pMax += kCarTolerance;
                                                   >> 357     return true;
271   }                                               358   }
272   else                                         << 359   else // Calculate rotated vertex coordinates
273   {                                               360   {
274     G4double sinStart = GetSinStartPhi();      << 361     G4int i, noEntries, noBetweenSections4;
275     G4double cosStart = GetCosStartPhi();      << 362     G4bool existsAfterClip = false;
276     G4double sinEnd   = GetSinEndPhi();        << 363     G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
277     G4double cosEnd   = GetCosEndPhi();        << 364 
278     G4double sinCur   = sinStart*cosHalf + cos << 365     pMin =  kInfinity;
279     G4double cosCur   = cosStart*cosHalf - sin << 366     pMax = -kInfinity;
280                                                << 367 
281     // set quadrilaterals                      << 368     noEntries = vertices->size();
282     G4ThreeVectorList pols[NSTEPS+2];          << 369     noBetweenSections4 = noEntries - 4;
283     for (G4int k=0; k<ksteps+2; ++k) pols[k].r << 370     
284     pols[0][0].set(rmin*cosStart,rmin*sinStart << 371     for ( i = 0 ; i < noEntries ; i += 4 )
285     pols[0][1].set(rmin*cosStart,rmin*sinStart << 372     {
286     pols[0][2].set(rmax*cosStart,rmax*sinStart << 373       ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
287     pols[0][3].set(rmax*cosStart,rmax*sinStart << 374     }
288     for (G4int k=1; k<ksteps+1; ++k)           << 375     for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
289     {                                          << 376     {
290       pols[k][0].set(rmin*cosCur,rmin*sinCur,  << 377       ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
291       pols[k][1].set(rmin*cosCur,rmin*sinCur,- << 378     }
292       pols[k][2].set(rext*cosCur,rext*sinCur,- << 379     if ( (pMin != kInfinity) || (pMax != -kInfinity) )
293       pols[k][3].set(rext*cosCur,rext*sinCur,  << 380     {
294                                                << 381       existsAfterClip = true;
295       G4double sinTmp = sinCur;                << 382       pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
296       sinCur = sinCur*cosStep + cosCur*sinStep << 383       pMax += kCarTolerance;
297       cosCur = cosCur*cosStep - sinTmp*sinStep << 384     }
298     }                                          << 385     else
299     pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin << 386     {
300     pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin << 387       // Check for case where completely enveloping clipping volume
301     pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin << 388       // If point inside then we are confident that the solid completely
302     pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin << 389       // envelopes the clipping volume. Hence set min/max extents according
303                                                << 390       // to clipping volume extents along the specified axis.
304     // set envelope and calculate extent       << 391 
305     std::vector<const G4ThreeVectorList *> pol << 392       G4ThreeVector clipCentre(
306     polygons.resize(ksteps+2);                 << 393              (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
307     for (G4int k=0; k<ksteps+2; ++k) polygons[ << 394              (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
308     G4BoundingEnvelope benv(bmin,bmax,polygons << 395              (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
309     exist = benv.CalculateExtent(pAxis,pVoxelL << 396         
                                                   >> 397       if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
                                                   >> 398       {
                                                   >> 399         existsAfterClip = true;
                                                   >> 400         pMin            = pVoxelLimit.GetMinExtent(pAxis);
                                                   >> 401         pMax            = pVoxelLimit.GetMaxExtent(pAxis);
                                                   >> 402       }
                                                   >> 403     }
                                                   >> 404     delete vertices;
                                                   >> 405     return existsAfterClip;
310   }                                               406   }
311   return exist;                                << 
312 }                                                 407 }
313                                                   408 
                                                   >> 409 
314 //////////////////////////////////////////////    410 ///////////////////////////////////////////////////////////////////////////
315 //                                                411 //
316 // Return whether point inside/outside/on surf    412 // Return whether point inside/outside/on surface
317                                                   413 
318 EInside G4Tubs::Inside( const G4ThreeVector& p    414 EInside G4Tubs::Inside( const G4ThreeVector& p ) const
319 {                                                 415 {
320   G4double r2,pPhi,tolRMin,tolRMax;               416   G4double r2,pPhi,tolRMin,tolRMax;
321   EInside in = kOutside ;                         417   EInside in = kOutside ;
                                                   >> 418   static const G4double halfCarTolerance=kCarTolerance*0.5;
                                                   >> 419   static const G4double halfRadTolerance=kRadTolerance*0.5;
                                                   >> 420   static const G4double halfAngTolerance=kAngTolerance*0.5;
322                                                   421 
323   if (std::fabs(p.z()) <= fDz - halfCarToleran    422   if (std::fabs(p.z()) <= fDz - halfCarTolerance)
324   {                                               423   {
325     r2 = p.x()*p.x() + p.y()*p.y() ;              424     r2 = p.x()*p.x() + p.y()*p.y() ;
326                                                   425 
327     if (fRMin != 0.0) { tolRMin = fRMin + half << 426     if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
328     else       { tolRMin = 0 ; }                  427     else       { tolRMin = 0 ; }
329                                                   428 
330     tolRMax = fRMax - halfRadTolerance ;          429     tolRMax = fRMax - halfRadTolerance ;
331                                                << 430       
332     if ((r2 >= tolRMin*tolRMin) && (r2 <= tolR    431     if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
333     {                                             432     {
334       if ( fPhiFullTube )                         433       if ( fPhiFullTube )
335       {                                           434       {
336         in = kInside ;                            435         in = kInside ;
337       }                                           436       }
338       else                                        437       else
339       {                                           438       {
340         // Try inner tolerant phi boundaries (    439         // Try inner tolerant phi boundaries (=>inside)
341         // if not inside, try outer tolerant p    440         // if not inside, try outer tolerant phi boundaries
342                                                   441 
343         if ( (tolRMin==0) && (std::fabs(p.x())    442         if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
344                           && (std::fabs(p.y())    443                           && (std::fabs(p.y())<=halfCarTolerance) )
345         {                                         444         {
346           in=kSurface;                            445           in=kSurface;
347         }                                         446         }
348         else                                      447         else
349         {                                         448         {
350           pPhi = std::atan2(p.y(),p.x()) ;        449           pPhi = std::atan2(p.y(),p.x()) ;
351           if ( pPhi < -halfAngTolerance )  { p    450           if ( pPhi < -halfAngTolerance )  { pPhi += twopi; } // 0<=pPhi<2pi
352                                                   451 
353           if ( fSPhi >= 0 )                       452           if ( fSPhi >= 0 )
354           {                                       453           {
355             if ( (std::fabs(pPhi) < halfAngTol    454             if ( (std::fabs(pPhi) < halfAngTolerance)
356               && (std::fabs(fSPhi + fDPhi - tw    455               && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
357             {                                  << 456             { 
358               pPhi += twopi ; // 0 <= pPhi < 2    457               pPhi += twopi ; // 0 <= pPhi < 2pi
359             }                                     458             }
360             if ( (pPhi >= fSPhi + halfAngToler    459             if ( (pPhi >= fSPhi + halfAngTolerance)
361               && (pPhi <= fSPhi + fDPhi - half    460               && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
362             {                                     461             {
363               in = kInside ;                      462               in = kInside ;
364             }                                     463             }
365             else if ( (pPhi >= fSPhi - halfAng    464             else if ( (pPhi >= fSPhi - halfAngTolerance)
366                    && (pPhi <= fSPhi + fDPhi +    465                    && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
367             {                                     466             {
368               in = kSurface ;                     467               in = kSurface ;
369             }                                     468             }
370           }                                       469           }
371           else  // fSPhi < 0                      470           else  // fSPhi < 0
372           {                                       471           {
373             if ( (pPhi <= fSPhi + twopi - half    472             if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
374               && (pPhi >= fSPhi + fDPhi  + hal    473               && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;} //kOutside
375             else if ( (pPhi <= fSPhi + twopi +    474             else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
376                    && (pPhi >= fSPhi + fDPhi      475                    && (pPhi >= fSPhi + fDPhi  - halfAngTolerance) )
377             {                                     476             {
378               in = kSurface ;                     477               in = kSurface ;
379             }                                     478             }
380             else                                  479             else
381             {                                     480             {
382               in = kInside ;                      481               in = kInside ;
383             }                                     482             }
384           }                                       483           }
385         }                                      << 484         }                    
386       }                                           485       }
387     }                                             486     }
388     else  // Try generous boundaries              487     else  // Try generous boundaries
389     {                                             488     {
390       tolRMin = fRMin - halfRadTolerance ;        489       tolRMin = fRMin - halfRadTolerance ;
391       tolRMax = fRMax + halfRadTolerance ;        490       tolRMax = fRMax + halfRadTolerance ;
392                                                   491 
393       if ( tolRMin < 0 )  { tolRMin = 0; }        492       if ( tolRMin < 0 )  { tolRMin = 0; }
394                                                   493 
395       if ( (r2 >= tolRMin*tolRMin) && (r2 <= t    494       if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
396       {                                           495       {
397         if (fPhiFullTube || (r2 <=halfRadToler    496         if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
398         {                        // Continuous    497         {                        // Continuous in phi or on z-axis
399           in = kSurface ;                         498           in = kSurface ;
400         }                                         499         }
401         else // Try outer tolerant phi boundar    500         else // Try outer tolerant phi boundaries only
402         {                                         501         {
403           pPhi = std::atan2(p.y(),p.x()) ;        502           pPhi = std::atan2(p.y(),p.x()) ;
404                                                   503 
405           if ( pPhi < -halfAngTolerance)  { pP    504           if ( pPhi < -halfAngTolerance)  { pPhi += twopi; } // 0<=pPhi<2pi
406           if ( fSPhi >= 0 )                       505           if ( fSPhi >= 0 )
407           {                                       506           {
408             if ( (std::fabs(pPhi) < halfAngTol    507             if ( (std::fabs(pPhi) < halfAngTolerance)
409               && (std::fabs(fSPhi + fDPhi - tw    508               && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
410             {                                  << 509             { 
411               pPhi += twopi ; // 0 <= pPhi < 2    510               pPhi += twopi ; // 0 <= pPhi < 2pi
412             }                                     511             }
413             if ( (pPhi >= fSPhi - halfAngToler    512             if ( (pPhi >= fSPhi - halfAngTolerance)
414               && (pPhi <= fSPhi + fDPhi + half    513               && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
415             {                                     514             {
416               in = kSurface ;                     515               in = kSurface ;
417             }                                     516             }
418           }                                       517           }
419           else  // fSPhi < 0                      518           else  // fSPhi < 0
420           {                                       519           {
421             if ( (pPhi <= fSPhi + twopi - half    520             if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
422               && (pPhi >= fSPhi + fDPhi + half    521               && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
423             else                                  522             else
424             {                                     523             {
425               in = kSurface ;                     524               in = kSurface ;
426             }                                     525             }
427           }                                       526           }
428         }                                         527         }
429       }                                           528       }
430     }                                             529     }
431   }                                               530   }
432   else if (std::fabs(p.z()) <= fDz + halfCarTo    531   else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
433   {                                          /    532   {                                          // Check within tolerant r limits
434     r2      = p.x()*p.x() + p.y()*p.y() ;         533     r2      = p.x()*p.x() + p.y()*p.y() ;
435     tolRMin = fRMin - halfRadTolerance ;          534     tolRMin = fRMin - halfRadTolerance ;
436     tolRMax = fRMax + halfRadTolerance ;          535     tolRMax = fRMax + halfRadTolerance ;
437                                                   536 
438     if ( tolRMin < 0 )  { tolRMin = 0; }          537     if ( tolRMin < 0 )  { tolRMin = 0; }
439                                                   538 
440     if ( (r2 >= tolRMin*tolRMin) && (r2 <= tol    539     if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
441     {                                             540     {
442       if (fPhiFullTube || (r2 <=halfRadToleran    541       if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
443       {                        // Continuous i    542       {                        // Continuous in phi or on z-axis
444         in = kSurface ;                           543         in = kSurface ;
445       }                                           544       }
446       else // Try outer tolerant phi boundarie    545       else // Try outer tolerant phi boundaries
447       {                                           546       {
448         pPhi = std::atan2(p.y(),p.x()) ;          547         pPhi = std::atan2(p.y(),p.x()) ;
449                                                   548 
450         if ( pPhi < -halfAngTolerance )  { pPh    549         if ( pPhi < -halfAngTolerance )  { pPhi += twopi; }  // 0<=pPhi<2pi
451         if ( fSPhi >= 0 )                         550         if ( fSPhi >= 0 )
452         {                                         551         {
453           if ( (std::fabs(pPhi) < halfAngToler    552           if ( (std::fabs(pPhi) < halfAngTolerance)
454             && (std::fabs(fSPhi + fDPhi - twop    553             && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
455           {                                    << 554           { 
456             pPhi += twopi ; // 0 <= pPhi < 2pi    555             pPhi += twopi ; // 0 <= pPhi < 2pi
457           }                                       556           }
458           if ( (pPhi >= fSPhi - halfAngToleran    557           if ( (pPhi >= fSPhi - halfAngTolerance)
459             && (pPhi <= fSPhi + fDPhi + halfAn    558             && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
460           {                                       559           {
461             in = kSurface;                        560             in = kSurface;
462           }                                       561           }
463         }                                         562         }
464         else  // fSPhi < 0                        563         else  // fSPhi < 0
465         {                                         564         {
466           if ( (pPhi <= fSPhi + twopi - halfAn    565           if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
467             && (pPhi >= fSPhi + fDPhi  + halfA    566             && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) ) {;}
468           else                                    567           else
469           {                                       568           {
470             in = kSurface ;                       569             in = kSurface ;
471           }                                       570           }
472         }                                      << 571         }      
473       }                                           572       }
474     }                                             573     }
475   }                                               574   }
476   return in;                                      575   return in;
477 }                                                 576 }
478                                                   577 
479 //////////////////////////////////////////////    578 ///////////////////////////////////////////////////////////////////////////
480 //                                                579 //
481 // Return unit normal of surface closest to p     580 // Return unit normal of surface closest to p
482 // - note if point on z axis, ignore phi divid    581 // - note if point on z axis, ignore phi divided sides
483 // - unsafe if point close to z axis a rmin=0     582 // - unsafe if point close to z axis a rmin=0 - no explicit checks
484                                                   583 
485 G4ThreeVector G4Tubs::SurfaceNormal( const G4T    584 G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const
486 {                                                 585 {
487   G4int noSurfaces = 0;                           586   G4int noSurfaces = 0;
488   G4double rho, pPhi;                             587   G4double rho, pPhi;
489   G4double distZ, distRMin, distRMax;             588   G4double distZ, distRMin, distRMax;
490   G4double distSPhi = kInfinity, distEPhi = kI    589   G4double distSPhi = kInfinity, distEPhi = kInfinity;
491                                                   590 
                                                   >> 591   static const G4double halfCarTolerance = 0.5*kCarTolerance;
                                                   >> 592   static const G4double halfAngTolerance = 0.5*kAngTolerance;
                                                   >> 593 
492   G4ThreeVector norm, sumnorm(0.,0.,0.);          594   G4ThreeVector norm, sumnorm(0.,0.,0.);
493   G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);    595   G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
494   G4ThreeVector nR, nPs, nPe;                     596   G4ThreeVector nR, nPs, nPe;
495                                                   597 
496   rho = std::sqrt(p.x()*p.x() + p.y()*p.y());     598   rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
497                                                   599 
498   distRMin = std::fabs(rho - fRMin);              600   distRMin = std::fabs(rho - fRMin);
499   distRMax = std::fabs(rho - fRMax);              601   distRMax = std::fabs(rho - fRMax);
500   distZ    = std::fabs(std::fabs(p.z()) - fDz)    602   distZ    = std::fabs(std::fabs(p.z()) - fDz);
501                                                   603 
502   if (!fPhiFullTube)    // Protected against ( << 604   if (!fPhiFullTube)    // Protected against (0,0,z) 
503   {                                               605   {
504     if ( rho > halfCarTolerance )                 606     if ( rho > halfCarTolerance )
505     {                                             607     {
506       pPhi = std::atan2(p.y(),p.x());             608       pPhi = std::atan2(p.y(),p.x());
                                                   >> 609     
                                                   >> 610       if(pPhi  < fSPhi- halfCarTolerance)           { pPhi += twopi; }
                                                   >> 611       else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
507                                                   612 
508       if (pPhi  < fSPhi-halfCarTolerance)      << 613       distSPhi = std::fabs(pPhi - fSPhi);       
509       else if (pPhi > fSPhi+fDPhi+halfCarToler << 614       distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 
510                                                << 
511       distSPhi = std::fabs( pPhi - fSPhi );    << 
512       distEPhi = std::fabs( pPhi - fSPhi - fDP << 
513     }                                             615     }
514     else if ( fRMin == 0.0 )                   << 616     else if( !fRMin )
515     {                                             617     {
516       distSPhi = 0.;                           << 618       distSPhi = 0.; 
517       distEPhi = 0.;                           << 619       distEPhi = 0.; 
518     }                                             620     }
519     nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0  << 621     nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
520     nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0  << 622     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
521   }                                               623   }
522   if ( rho > halfCarTolerance ) { nR = G4Three    624   if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
523                                                   625 
524   if( distRMax <= halfCarTolerance )              626   if( distRMax <= halfCarTolerance )
525   {                                               627   {
526     ++noSurfaces;                              << 628     noSurfaces ++;
527     sumnorm += nR;                                629     sumnorm += nR;
528   }                                               630   }
529   if( (fRMin != 0.0) && (distRMin <= halfCarTo << 631   if( fRMin && (distRMin <= halfCarTolerance) )
530   {                                               632   {
531     ++noSurfaces;                              << 633     noSurfaces ++;
532     sumnorm -= nR;                                634     sumnorm -= nR;
533   }                                               635   }
534   if( fDPhi < twopi )                          << 636   if( fDPhi < twopi )   
535   {                                               637   {
536     if (distSPhi <= halfAngTolerance)          << 638     if (distSPhi <= halfAngTolerance)  
537     {                                             639     {
538       ++noSurfaces;                            << 640       noSurfaces ++;
539       sumnorm += nPs;                             641       sumnorm += nPs;
540     }                                             642     }
541     if (distEPhi <= halfAngTolerance)          << 643     if (distEPhi <= halfAngTolerance)  
542     {                                             644     {
543       ++noSurfaces;                            << 645       noSurfaces ++;
544       sumnorm += nPe;                             646       sumnorm += nPe;
545     }                                             647     }
546   }                                               648   }
547   if (distZ <= halfCarTolerance)               << 649   if (distZ <= halfCarTolerance)  
548   {                                               650   {
549     ++noSurfaces;                              << 651     noSurfaces ++;
550     if ( p.z() >= 0.)  { sumnorm += nZ; }         652     if ( p.z() >= 0.)  { sumnorm += nZ; }
551     else               { sumnorm -= nZ; }         653     else               { sumnorm -= nZ; }
552   }                                               654   }
553   if ( noSurfaces == 0 )                          655   if ( noSurfaces == 0 )
554   {                                               656   {
555 #ifdef G4CSGDEBUG                                 657 #ifdef G4CSGDEBUG
556     G4Exception("G4Tubs::SurfaceNormal(p)", "G    658     G4Exception("G4Tubs::SurfaceNormal(p)", "GeomSolids1002",
557                 JustWarning, "Point p is not o    659                 JustWarning, "Point p is not on surface !?" );
558     G4long oldprc = G4cout.precision(20);      << 660     G4int oldprc = G4cout.precision(20);
559     G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y    661     G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
560           << G4endl << G4endl;                    662           << G4endl << G4endl;
561     G4cout.precision(oldprc) ;                    663     G4cout.precision(oldprc) ;
562 #endif                                         << 664 #endif 
563      norm = ApproxSurfaceNormal(p);               665      norm = ApproxSurfaceNormal(p);
564   }                                               666   }
565   else if ( noSurfaces == 1 )  { norm = sumnor    667   else if ( noSurfaces == 1 )  { norm = sumnorm; }
566   else                         { norm = sumnor    668   else                         { norm = sumnorm.unit(); }
567                                                   669 
568   return norm;                                    670   return norm;
569 }                                                 671 }
570                                                   672 
571 //////////////////////////////////////////////    673 /////////////////////////////////////////////////////////////////////////////
572 //                                                674 //
573 // Algorithm for SurfaceNormal() following the    675 // Algorithm for SurfaceNormal() following the original specification
574 // for points not on the surface                  676 // for points not on the surface
575                                                   677 
576 G4ThreeVector G4Tubs::ApproxSurfaceNormal( con    678 G4ThreeVector G4Tubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
577 {                                                 679 {
578   ENorm side ;                                    680   ENorm side ;
579   G4ThreeVector norm ;                            681   G4ThreeVector norm ;
580   G4double rho, phi ;                             682   G4double rho, phi ;
581   G4double distZ, distRMin, distRMax, distSPhi    683   G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
582                                                   684 
583   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;    685   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
584                                                   686 
585   distRMin = std::fabs(rho - fRMin) ;             687   distRMin = std::fabs(rho - fRMin) ;
586   distRMax = std::fabs(rho - fRMax) ;             688   distRMax = std::fabs(rho - fRMax) ;
587   distZ    = std::fabs(std::fabs(p.z()) - fDz)    689   distZ    = std::fabs(std::fabs(p.z()) - fDz) ;
588                                                   690 
589   if (distRMin < distRMax) // First minimum       691   if (distRMin < distRMax) // First minimum
590   {                                               692   {
591     if ( distZ < distRMin )                       693     if ( distZ < distRMin )
592     {                                             694     {
593        distMin = distZ ;                          695        distMin = distZ ;
594        side    = kNZ ;                            696        side    = kNZ ;
595     }                                             697     }
596     else                                          698     else
597     {                                             699     {
598       distMin = distRMin ;                        700       distMin = distRMin ;
599       side    = kNRMin   ;                        701       side    = kNRMin   ;
600     }                                             702     }
601   }                                               703   }
602   else                                            704   else
603   {                                               705   {
604     if ( distZ < distRMax )                       706     if ( distZ < distRMax )
605     {                                             707     {
606       distMin = distZ ;                           708       distMin = distZ ;
607       side    = kNZ   ;                           709       side    = kNZ   ;
608     }                                             710     }
609     else                                          711     else
610     {                                             712     {
611       distMin = distRMax ;                        713       distMin = distRMax ;
612       side    = kNRMax   ;                        714       side    = kNRMax   ;
613     }                                             715     }
614   }                                            << 716   }   
615   if (!fPhiFullTube  &&  (rho != 0.0) ) // Pro << 717   if (!fPhiFullTube  &&  rho ) // Protected against (0,0,z) 
616   {                                               718   {
617     phi = std::atan2(p.y(),p.x()) ;               719     phi = std::atan2(p.y(),p.x()) ;
618                                                   720 
619     if ( phi < 0 )  { phi += twopi; }             721     if ( phi < 0 )  { phi += twopi; }
620                                                   722 
621     if ( fSPhi < 0 )                              723     if ( fSPhi < 0 )
622     {                                             724     {
623       distSPhi = std::fabs(phi - (fSPhi + twop    725       distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
624     }                                             726     }
625     else                                          727     else
626     {                                             728     {
627       distSPhi = std::fabs(phi - fSPhi)*rho ;     729       distSPhi = std::fabs(phi - fSPhi)*rho ;
628     }                                             730     }
629     distEPhi = std::fabs(phi - fSPhi - fDPhi)*    731     distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
630                                                << 732                                       
631     if (distSPhi < distEPhi) // Find new minim    733     if (distSPhi < distEPhi) // Find new minimum
632     {                                             734     {
633       if ( distSPhi < distMin )                   735       if ( distSPhi < distMin )
634       {                                           736       {
635         side = kNSPhi ;                           737         side = kNSPhi ;
636       }                                           738       }
637     }                                             739     }
638     else                                          740     else
639     {                                             741     {
640       if ( distEPhi < distMin )                   742       if ( distEPhi < distMin )
641       {                                           743       {
642         side = kNEPhi ;                           744         side = kNEPhi ;
643       }                                           745       }
644     }                                             746     }
645   }                                            << 747   }    
646   switch ( side )                                 748   switch ( side )
647   {                                               749   {
648     case kNRMin : // Inner radius                 750     case kNRMin : // Inner radius
649     {                                          << 751     {                      
650       norm = G4ThreeVector(-p.x()/rho, -p.y()/    752       norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
651       break ;                                     753       break ;
652     }                                             754     }
653     case kNRMax : // Outer radius                 755     case kNRMax : // Outer radius
654     {                                          << 756     {                  
655       norm = G4ThreeVector(p.x()/rho, p.y()/rh    757       norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
656       break ;                                     758       break ;
657     }                                             759     }
658     case kNZ :    // + or - dz                    760     case kNZ :    // + or - dz
659     {                                          << 761     {                              
660       if ( p.z() > 0 )  { norm = G4ThreeVector    762       if ( p.z() > 0 )  { norm = G4ThreeVector(0,0,1) ; }
661       else              { norm = G4ThreeVector    763       else              { norm = G4ThreeVector(0,0,-1); }
662       break ;                                     764       break ;
663     }                                             765     }
664     case kNSPhi:                                  766     case kNSPhi:
665     {                                             767     {
666       norm = G4ThreeVector(sinSPhi, -cosSPhi,  << 768       norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
667       break ;                                     769       break ;
668     }                                             770     }
669     case kNEPhi:                                  771     case kNEPhi:
670     {                                             772     {
671       norm = G4ThreeVector(-sinEPhi, cosEPhi,  << 773       norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
672       break;                                      774       break;
673     }                                             775     }
674     default:      // Should never reach this c    776     default:      // Should never reach this case ...
675     {                                             777     {
676       DumpInfo();                                 778       DumpInfo();
677       G4Exception("G4Tubs::ApproxSurfaceNormal    779       G4Exception("G4Tubs::ApproxSurfaceNormal()",
678                   "GeomSolids1002", JustWarnin    780                   "GeomSolids1002", JustWarning,
679                   "Undefined side for valid su    781                   "Undefined side for valid surface normal to solid.");
680       break ;                                     782       break ;
681     }                                          << 783     }    
682   }                                            << 784   }                
683   return norm;                                    785   return norm;
684 }                                                 786 }
685                                                   787 
686 //////////////////////////////////////////////    788 ////////////////////////////////////////////////////////////////////
687 //                                                789 //
688 //                                                790 //
689 // Calculate distance to shape from outside, a    791 // Calculate distance to shape from outside, along normalised vector
690 // - return kInfinity if no intersection, or i    792 // - return kInfinity if no intersection, or intersection distance <= tolerance
691 //                                                793 //
692 // - Compute the intersection with the z plane << 794 // - Compute the intersection with the z planes 
693 //        - if at valid r, phi, return            795 //        - if at valid r, phi, return
694 //                                                796 //
695 // -> If point is outer outer radius, compute     797 // -> If point is outer outer radius, compute intersection with rmax
696 //        - if at valid phi,z return              798 //        - if at valid phi,z return
697 //                                                799 //
698 // -> Compute intersection with inner radius,     800 // -> Compute intersection with inner radius, taking largest +ve root
699 //        - if valid (in z,phi), save intersct    801 //        - if valid (in z,phi), save intersction
700 //                                                802 //
701 //    -> If phi segmented, compute intersectio    803 //    -> If phi segmented, compute intersections with phi half planes
702 //        - return smallest of valid phi inter    804 //        - return smallest of valid phi intersections and
703 //          inner radius intersection             805 //          inner radius intersection
704 //                                                806 //
705 // NOTE:                                          807 // NOTE:
706 // - 'if valid' implies tolerant checking of i    808 // - 'if valid' implies tolerant checking of intersection points
707                                                   809 
708 G4double G4Tubs::DistanceToIn( const G4ThreeVe    810 G4double G4Tubs::DistanceToIn( const G4ThreeVector& p,
709                                const G4ThreeVe    811                                const G4ThreeVector& v  ) const
710 {                                                 812 {
711   G4double snxt = kInfinity ;      // snxt = d    813   G4double snxt = kInfinity ;      // snxt = default return value
712   G4double tolORMin2, tolIRMax2 ;  // 'generou    814   G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
713   G4double tolORMax2, tolIRMin2, tolODz, tolID    815   G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
714   const G4double dRmax = 100.*fRMax;              816   const G4double dRmax = 100.*fRMax;
715                                                   817 
                                                   >> 818   static const G4double halfCarTolerance = 0.5*kCarTolerance;
                                                   >> 819   static const G4double halfRadTolerance = 0.5*kRadTolerance;
                                                   >> 820 
716   // Intersection point variables                 821   // Intersection point variables
717   //                                              822   //
718   G4double Dist, sd, xi, yi, zi, rho2, inum, i    823   G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
719   G4double t1, t2, t3, b, c, d ;     // Quadra << 824   G4double t1, t2, t3, b, c, d ;     // Quadratic solver variables 
720                                                << 825   
721   // Calculate tolerant rmin and rmax             826   // Calculate tolerant rmin and rmax
722                                                   827 
723   if (fRMin > kRadTolerance)                      828   if (fRMin > kRadTolerance)
724   {                                               829   {
725     tolORMin2 = (fRMin - halfRadTolerance)*(fR    830     tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
726     tolIRMin2 = (fRMin + halfRadTolerance)*(fR    831     tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
727   }                                               832   }
728   else                                            833   else
729   {                                               834   {
730     tolORMin2 = 0.0 ;                             835     tolORMin2 = 0.0 ;
731     tolIRMin2 = 0.0 ;                             836     tolIRMin2 = 0.0 ;
732   }                                               837   }
733   tolORMax2 = (fRMax + halfRadTolerance)*(fRMa    838   tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
734   tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa    839   tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
735                                                   840 
736   // Intersection with Z surfaces                 841   // Intersection with Z surfaces
737                                                   842 
738   tolIDz = fDz - halfCarTolerance ;               843   tolIDz = fDz - halfCarTolerance ;
739   tolODz = fDz + halfCarTolerance ;               844   tolODz = fDz + halfCarTolerance ;
740                                                   845 
741   if (std::fabs(p.z()) >= tolIDz)                 846   if (std::fabs(p.z()) >= tolIDz)
742   {                                               847   {
743     if ( p.z()*v.z() < 0 )    // at +Z going i    848     if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
744     {                                             849     {
745       sd = (std::fabs(p.z()) - fDz)/std::fabs(    850       sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;  // Z intersect distance
746                                                   851 
747       if(sd < 0.0)  { sd = 0.0; }                 852       if(sd < 0.0)  { sd = 0.0; }
748                                                   853 
749       xi   = p.x() + sd*v.x() ;                   854       xi   = p.x() + sd*v.x() ;                // Intersection coords
750       yi   = p.y() + sd*v.y() ;                   855       yi   = p.y() + sd*v.y() ;
751       rho2 = xi*xi + yi*yi ;                      856       rho2 = xi*xi + yi*yi ;
752                                                   857 
753       // Check validity of intersection           858       // Check validity of intersection
754                                                   859 
755       if ((tolIRMin2 <= rho2) && (rho2 <= tolI    860       if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
756       {                                           861       {
757         if (!fPhiFullTube && (rho2 != 0.0))    << 862         if (!fPhiFullTube && rho2)
758         {                                         863         {
759           // Psi = angle made with central (av    864           // Psi = angle made with central (average) phi of shape
760           //                                      865           //
761           inum   = xi*cosCPhi + yi*sinCPhi ;      866           inum   = xi*cosCPhi + yi*sinCPhi ;
762           iden   = std::sqrt(rho2) ;              867           iden   = std::sqrt(rho2) ;
763           cosPsi = inum/iden ;                    868           cosPsi = inum/iden ;
764           if (cosPsi >= cosHDPhiIT)  { return     869           if (cosPsi >= cosHDPhiIT)  { return sd ; }
765         }                                         870         }
766         else                                      871         else
767         {                                         872         {
768           return sd ;                             873           return sd ;
769         }                                         874         }
770       }                                           875       }
771     }                                             876     }
772     else                                          877     else
773     {                                             878     {
774       if ( snxt<halfCarTolerance )  { snxt=0;     879       if ( snxt<halfCarTolerance )  { snxt=0; }
775       return snxt ;  // On/outside extent, and    880       return snxt ;  // On/outside extent, and heading away
776                      // -> cannot intersect       881                      // -> cannot intersect
777     }                                             882     }
778   }                                               883   }
779                                                   884 
780   // -> Can not intersect z surfaces              885   // -> Can not intersect z surfaces
781   //                                              886   //
782   // Intersection with rmax (possible return)     887   // Intersection with rmax (possible return) and rmin (must also check phi)
783   //                                              888   //
784   // Intersection point (xi,yi,zi) on line x=p    889   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
785   //                                              890   //
786   // Intersects with x^2+y^2=R^2                  891   // Intersects with x^2+y^2=R^2
787   //                                              892   //
788   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.    893   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
789   //            t1                t2              894   //            t1                t2                t3
790                                                   895 
791   t1 = 1.0 - v.z()*v.z() ;                        896   t1 = 1.0 - v.z()*v.z() ;
792   t2 = p.x()*v.x() + p.y()*v.y() ;                897   t2 = p.x()*v.x() + p.y()*v.y() ;
793   t3 = p.x()*p.x() + p.y()*p.y() ;                898   t3 = p.x()*p.x() + p.y()*p.y() ;
794                                                   899 
795   if ( t1 > 0 )        // Check not || to z ax    900   if ( t1 > 0 )        // Check not || to z axis
796   {                                               901   {
797     b = t2/t1 ;                                   902     b = t2/t1 ;
798     c = t3 - fRMax*fRMax ;                        903     c = t3 - fRMax*fRMax ;
799     if ((t3 >= tolORMax2) && (t2<0))   // This    904     if ((t3 >= tolORMax2) && (t2<0))   // This also handles the tangent case
800     {                                             905     {
801       // Try outer cylinder intersection          906       // Try outer cylinder intersection
802       //          c=(t3-fRMax*fRMax)/t1;          907       //          c=(t3-fRMax*fRMax)/t1;
803                                                   908 
804       c /= t1 ;                                   909       c /= t1 ;
805       d = b*b - c ;                               910       d = b*b - c ;
806                                                   911 
807       if (d >= 0)  // If real root                912       if (d >= 0)  // If real root
808       {                                           913       {
809         sd = c/(-b+std::sqrt(d));                 914         sd = c/(-b+std::sqrt(d));
810         if (sd >= 0)  // If 'forwards'            915         if (sd >= 0)  // If 'forwards'
811         {                                         916         {
812           if ( sd>dRmax ) // Avoid rounding er    917           if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
813           {               // 64 bits systems.     918           {               // 64 bits systems. Split long distances and recompute
814             G4double fTerm = sd-std::fmod(sd,d    919             G4double fTerm = sd-std::fmod(sd,dRmax);
815             sd = fTerm + DistanceToIn(p+fTerm*    920             sd = fTerm + DistanceToIn(p+fTerm*v,v);
816           }                                    << 921           } 
817           // Check z intersection                 922           // Check z intersection
818           //                                      923           //
819           zi = p.z() + sd*v.z() ;                 924           zi = p.z() + sd*v.z() ;
820           if (std::fabs(zi)<=tolODz)              925           if (std::fabs(zi)<=tolODz)
821           {                                       926           {
822             // Z ok. Check phi intersection if    927             // Z ok. Check phi intersection if reqd
823             //                                    928             //
824             if (fPhiFullTube)                     929             if (fPhiFullTube)
825             {                                     930             {
826               return sd ;                         931               return sd ;
827             }                                     932             }
828             else                                  933             else
829             {                                     934             {
830               xi     = p.x() + sd*v.x() ;         935               xi     = p.x() + sd*v.x() ;
831               yi     = p.y() + sd*v.y() ;         936               yi     = p.y() + sd*v.y() ;
832               cosPsi = (xi*cosCPhi + yi*sinCPh    937               cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
833               if (cosPsi >= cosHDPhiIT)  { ret    938               if (cosPsi >= cosHDPhiIT)  { return sd ; }
834             }                                     939             }
835           }  //  end if std::fabs(zi)             940           }  //  end if std::fabs(zi)
836         }    //  end if (sd>=0)                   941         }    //  end if (sd>=0)
837       }      //  end if (d>=0)                    942       }      //  end if (d>=0)
838     }        //  end if (r>=fRMax)                943     }        //  end if (r>=fRMax)
839     else                                       << 944     else 
840     {                                             945     {
841       // Inside outer radius :                    946       // Inside outer radius :
842       // check not inside, and heading through    947       // check not inside, and heading through tubs (-> 0 to in)
843                                                   948 
844       if ((t3 > tolIRMin2) && (t2 < 0) && (std    949       if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
845       {                                           950       {
846         // Inside both radii, delta r -ve, ins    951         // Inside both radii, delta r -ve, inside z extent
847                                                   952 
848         if (!fPhiFullTube)                        953         if (!fPhiFullTube)
849         {                                         954         {
850           inum   = p.x()*cosCPhi + p.y()*sinCP    955           inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
851           iden   = std::sqrt(t3) ;                956           iden   = std::sqrt(t3) ;
852           cosPsi = inum/iden ;                    957           cosPsi = inum/iden ;
853           if (cosPsi >= cosHDPhiIT)               958           if (cosPsi >= cosHDPhiIT)
854           {                                       959           {
855             // In the old version, the small n    960             // In the old version, the small negative tangent for the point
856             // on surface was not taken in acc    961             // on surface was not taken in account, and returning 0.0 ...
857             // New version: check the tangent  << 962             // New version: check the tangent for the point on surface and 
858             // if no intersection, return kInf    963             // if no intersection, return kInfinity, if intersection instead
859             // return sd.                         964             // return sd.
860             //                                    965             //
861             c = t3-fRMax*fRMax;                << 966             c = t3-fRMax*fRMax; 
862             if ( c<=0.0 )                         967             if ( c<=0.0 )
863             {                                     968             {
864               return 0.0;                         969               return 0.0;
865             }                                     970             }
866             else                                  971             else
867             {                                     972             {
868               c = c/t1 ;                          973               c = c/t1 ;
869               d = b*b-c;                          974               d = b*b-c;
870               if ( d>=0.0 )                       975               if ( d>=0.0 )
871               {                                   976               {
872                 snxt = c/(-b+std::sqrt(d)); //    977                 snxt = c/(-b+std::sqrt(d)); // using safe solution
873                                             // << 978                                             // for quadratic equation 
874                 if ( snxt < halfCarTolerance )    979                 if ( snxt < halfCarTolerance ) { snxt=0; }
875                 return snxt ;                     980                 return snxt ;
876               }                                << 981               }      
877               else                                982               else
878               {                                   983               {
879                 return kInfinity;                 984                 return kInfinity;
880               }                                   985               }
881             }                                     986             }
882           }                                    << 987           } 
883         }                                         988         }
884         else                                      989         else
885         {                                      << 990         {   
886           // In the old version, the small neg    991           // In the old version, the small negative tangent for the point
887           // on surface was not taken in accou    992           // on surface was not taken in account, and returning 0.0 ...
888           // New version: check the tangent fo << 993           // New version: check the tangent for the point on surface and 
889           // if no intersection, return kInfin    994           // if no intersection, return kInfinity, if intersection instead
890           // return sd.                           995           // return sd.
891           //                                      996           //
892           c = t3 - fRMax*fRMax;                << 997           c = t3 - fRMax*fRMax; 
893           if ( c<=0.0 )                           998           if ( c<=0.0 )
894           {                                       999           {
895             return 0.0;                           1000             return 0.0;
896           }                                       1001           }
897           else                                    1002           else
898           {                                       1003           {
899             c = c/t1 ;                            1004             c = c/t1 ;
900             d = b*b-c;                            1005             d = b*b-c;
901             if ( d>=0.0 )                         1006             if ( d>=0.0 )
902             {                                     1007             {
903               snxt= c/(-b+std::sqrt(d)); // us    1008               snxt= c/(-b+std::sqrt(d)); // using safe solution
904                                          // fo << 1009                                          // for quadratic equation 
905               if ( snxt < halfCarTolerance ) {    1010               if ( snxt < halfCarTolerance ) { snxt=0; }
906               return snxt ;                       1011               return snxt ;
907             }                                  << 1012             }      
908             else                                  1013             else
909             {                                     1014             {
910               return kInfinity;                   1015               return kInfinity;
911             }                                     1016             }
912           }                                       1017           }
913         } // end if   (!fPhiFullTube)             1018         } // end if   (!fPhiFullTube)
914       }   // end if   (t3>tolIRMin2)              1019       }   // end if   (t3>tolIRMin2)
915     }     // end if   (Inside Outer Radius)    << 1020     }     // end if   (Inside Outer Radius) 
916     if ( fRMin != 0.0 )    // Try inner cylind << 1021     if ( fRMin )    // Try inner cylinder intersection
917     {                                             1022     {
918       c = (t3 - fRMin*fRMin)/t1 ;                 1023       c = (t3 - fRMin*fRMin)/t1 ;
919       d = b*b - c ;                               1024       d = b*b - c ;
920       if ( d >= 0.0 )  // If real root            1025       if ( d >= 0.0 )  // If real root
921       {                                           1026       {
922         // Always want 2nd root - we are outsi    1027         // Always want 2nd root - we are outside and know rmax Hit was bad
923         // - If on surface of rmin also need f    1028         // - If on surface of rmin also need farthest root
924                                                   1029 
925         sd =( b > 0. )? c/(-b - std::sqrt(d))     1030         sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
926         if (sd >= -halfCarTolerance)  // check    1031         if (sd >= -halfCarTolerance)  // check forwards
927         {                                         1032         {
928           // Check z intersection                 1033           // Check z intersection
929           //                                      1034           //
930           if(sd < 0.0)  { sd = 0.0; }             1035           if(sd < 0.0)  { sd = 0.0; }
931           if ( sd>dRmax ) // Avoid rounding er    1036           if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen
932           {               // 64 bits systems.     1037           {               // 64 bits systems. Split long distances and recompute
933             G4double fTerm = sd-std::fmod(sd,d    1038             G4double fTerm = sd-std::fmod(sd,dRmax);
934             sd = fTerm + DistanceToIn(p+fTerm*    1039             sd = fTerm + DistanceToIn(p+fTerm*v,v);
935           }                                    << 1040           } 
936           zi = p.z() + sd*v.z() ;                 1041           zi = p.z() + sd*v.z() ;
937           if (std::fabs(zi) <= tolODz)            1042           if (std::fabs(zi) <= tolODz)
938           {                                       1043           {
939             // Z ok. Check phi                    1044             // Z ok. Check phi
940             //                                    1045             //
941             if ( fPhiFullTube )                   1046             if ( fPhiFullTube )
942             {                                     1047             {
943               return sd ;                      << 1048               return sd ; 
944             }                                     1049             }
945             else                                  1050             else
946             {                                     1051             {
947               xi     = p.x() + sd*v.x() ;         1052               xi     = p.x() + sd*v.x() ;
948               yi     = p.y() + sd*v.y() ;         1053               yi     = p.y() + sd*v.y() ;
949               cosPsi = (xi*cosCPhi + yi*sinCPh << 1054               cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
950               if (cosPsi >= cosHDPhiIT)           1055               if (cosPsi >= cosHDPhiIT)
951               {                                   1056               {
952                 // Good inner radius isect        1057                 // Good inner radius isect
953                 // - but earlier phi isect sti    1058                 // - but earlier phi isect still possible
954                                                   1059 
955                 snxt = sd ;                       1060                 snxt = sd ;
956               }                                   1061               }
957             }                                     1062             }
958           }        //    end if std::fabs(zi)     1063           }        //    end if std::fabs(zi)
959         }          //    end if (sd>=0)           1064         }          //    end if (sd>=0)
960       }            //    end if (d>=0)            1065       }            //    end if (d>=0)
961     }              //    end if (fRMin)           1066     }              //    end if (fRMin)
962   }                                               1067   }
963                                                   1068 
964   // Phi segment intersection                     1069   // Phi segment intersection
965   //                                              1070   //
966   // o Tolerant of points inside phi planes by    1071   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
967   //                                              1072   //
968   // o NOTE: Large duplication of code between    1073   // o NOTE: Large duplication of code between sphi & ephi checks
969   //         -> only diffs: sphi -> ephi, Comp    1074   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
970   //            intersection check <=0 -> >=0     1075   //            intersection check <=0 -> >=0
971   //         -> use some form of loop Construc    1076   //         -> use some form of loop Construct ?
972   //                                              1077   //
973   if ( !fPhiFullTube )                            1078   if ( !fPhiFullTube )
974   {                                               1079   {
975     // First phi surface (Starting phi)           1080     // First phi surface (Starting phi)
976     //                                            1081     //
977     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;     1082     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
978                                                << 1083                     
979     if ( Comp < 0 )  // Component in outwards     1084     if ( Comp < 0 )  // Component in outwards normal dirn
980     {                                             1085     {
981       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;    1086       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
982                                                   1087 
983       if ( Dist < halfCarTolerance )              1088       if ( Dist < halfCarTolerance )
984       {                                           1089       {
985         sd = Dist/Comp ;                          1090         sd = Dist/Comp ;
986                                                   1091 
987         if (sd < snxt)                            1092         if (sd < snxt)
988         {                                         1093         {
989           if ( sd < 0 )  { sd = 0.0; }            1094           if ( sd < 0 )  { sd = 0.0; }
990           zi = p.z() + sd*v.z() ;                 1095           zi = p.z() + sd*v.z() ;
991           if ( std::fabs(zi) <= tolODz )          1096           if ( std::fabs(zi) <= tolODz )
992           {                                       1097           {
993             xi   = p.x() + sd*v.x() ;             1098             xi   = p.x() + sd*v.x() ;
994             yi   = p.y() + sd*v.y() ;             1099             yi   = p.y() + sd*v.y() ;
995             rho2 = xi*xi + yi*yi ;                1100             rho2 = xi*xi + yi*yi ;
996                                                   1101 
997             if ( ( (rho2 >= tolIRMin2) && (rho    1102             if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
998               || ( (rho2 >  tolORMin2) && (rho    1103               || ( (rho2 >  tolORMin2) && (rho2 <  tolIRMin2)
999                 && ( v.y()*cosSPhi - v.x()*sin    1104                 && ( v.y()*cosSPhi - v.x()*sinSPhi >  0 )
1000                 && ( v.x()*cosSPhi + v.y()*si    1105                 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 )     )
1001               || ( (rho2 > tolIRMax2) && (rho    1106               || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1002                 && (v.y()*cosSPhi - v.x()*sin    1107                 && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1003                 && (v.x()*cosSPhi + v.y()*sin    1108                 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) )    )
1004             {                                    1109             {
1005               // z and r intersections good      1110               // z and r intersections good
1006               // - check intersecting with co    1111               // - check intersecting with correct half-plane
1007               //                                 1112               //
1008               if ((yi*cosCPhi-xi*sinCPhi) <=     1113               if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1009             }                                    1114             }
1010           }                                      1115           }
1011         }                                        1116         }
1012       }                                       << 1117       }    
1013     }                                            1118     }
1014                                               << 1119       
1015     // Second phi surface (Ending phi)           1120     // Second phi surface (Ending phi)
1016                                                  1121 
1017     Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi    1122     Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1018                                               << 1123         
1019     if (Comp < 0 )  // Component in outwards     1124     if (Comp < 0 )  // Component in outwards normal dirn
1020     {                                            1125     {
1021       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi)    1126       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1022                                                  1127 
1023       if ( Dist < halfCarTolerance )             1128       if ( Dist < halfCarTolerance )
1024       {                                          1129       {
1025         sd = Dist/Comp ;                         1130         sd = Dist/Comp ;
1026                                                  1131 
1027         if (sd < snxt)                           1132         if (sd < snxt)
1028         {                                        1133         {
1029           if ( sd < 0 )  { sd = 0; }             1134           if ( sd < 0 )  { sd = 0; }
1030           zi = p.z() + sd*v.z() ;                1135           zi = p.z() + sd*v.z() ;
1031           if ( std::fabs(zi) <= tolODz )         1136           if ( std::fabs(zi) <= tolODz )
1032           {                                      1137           {
1033             xi   = p.x() + sd*v.x() ;            1138             xi   = p.x() + sd*v.x() ;
1034             yi   = p.y() + sd*v.y() ;            1139             yi   = p.y() + sd*v.y() ;
1035             rho2 = xi*xi + yi*yi ;               1140             rho2 = xi*xi + yi*yi ;
1036             if ( ( (rho2 >= tolIRMin2) && (rh    1141             if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1037                 || ( (rho2 > tolORMin2)  && (    1142                 || ( (rho2 > tolORMin2)  && (rho2 < tolIRMin2)
1038                   && (v.x()*sinEPhi - v.y()*c    1143                   && (v.x()*sinEPhi - v.y()*cosEPhi >  0)
1039                   && (v.x()*cosEPhi + v.y()*s    1144                   && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1040                 || ( (rho2 > tolIRMax2) && (r    1145                 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1041                   && (v.x()*sinEPhi - v.y()*c    1146                   && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1042                   && (v.x()*cosEPhi + v.y()*s    1147                   && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1043             {                                    1148             {
1044               // z and r intersections good      1149               // z and r intersections good
1045               // - check intersecting with co    1150               // - check intersecting with correct half-plane
1046               //                                 1151               //
1047               if ( (yi*cosCPhi-xi*sinCPhi) >=    1152               if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = sd; }
1048             }                         //?? >=    1153             }                         //?? >=-halfCarTolerance
1049           }                                      1154           }
1050         }                                        1155         }
1051       }                                          1156       }
1052     }         //  Comp < 0                       1157     }         //  Comp < 0
1053   }           //  !fPhiFullTube               << 1158   }           //  !fPhiFullTube 
1054   if ( snxt<halfCarTolerance )  { snxt=0; }      1159   if ( snxt<halfCarTolerance )  { snxt=0; }
1055   return snxt ;                                  1160   return snxt ;
1056 }                                                1161 }
1057                                               << 1162  
1058 /////////////////////////////////////////////    1163 //////////////////////////////////////////////////////////////////
1059 //                                               1164 //
1060 // Calculate distance to shape from outside,     1165 // Calculate distance to shape from outside, along normalised vector
1061 // - return kInfinity if no intersection, or     1166 // - return kInfinity if no intersection, or intersection distance <= tolerance
1062 //                                               1167 //
1063 // - Compute the intersection with the z plan << 1168 // - Compute the intersection with the z planes 
1064 //        - if at valid r, phi, return           1169 //        - if at valid r, phi, return
1065 //                                               1170 //
1066 // -> If point is outer outer radius, compute    1171 // -> If point is outer outer radius, compute intersection with rmax
1067 //        - if at valid phi,z return             1172 //        - if at valid phi,z return
1068 //                                               1173 //
1069 // -> Compute intersection with inner radius,    1174 // -> Compute intersection with inner radius, taking largest +ve root
1070 //        - if valid (in z,phi), save intersc    1175 //        - if valid (in z,phi), save intersction
1071 //                                               1176 //
1072 //    -> If phi segmented, compute intersecti    1177 //    -> If phi segmented, compute intersections with phi half planes
1073 //        - return smallest of valid phi inte    1178 //        - return smallest of valid phi intersections and
1074 //          inner radius intersection            1179 //          inner radius intersection
1075 //                                               1180 //
1076 // NOTE:                                         1181 // NOTE:
1077 // - Precalculations for phi trigonometry are    1182 // - Precalculations for phi trigonometry are Done `just in time'
1078 // - `if valid' implies tolerant checking of     1183 // - `if valid' implies tolerant checking of intersection points
1079 //   Calculate distance (<= actual) to closes    1184 //   Calculate distance (<= actual) to closest surface of shape from outside
1080 // - Calculate distance to z, radial planes      1185 // - Calculate distance to z, radial planes
1081 // - Only to phi planes if outside phi extent    1186 // - Only to phi planes if outside phi extent
1082 // - Return 0 if point inside                    1187 // - Return 0 if point inside
1083                                                  1188 
1084 G4double G4Tubs::DistanceToIn( const G4ThreeV    1189 G4double G4Tubs::DistanceToIn( const G4ThreeVector& p ) const
1085 {                                                1190 {
1086   G4double safe=0.0, rho, safe1, safe2, safe3    1191   G4double safe=0.0, rho, safe1, safe2, safe3 ;
1087   G4double safePhi, cosPsi ;                     1192   G4double safePhi, cosPsi ;
1088                                                  1193 
1089   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()    1194   rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1090   safe1 = fRMin - rho ;                          1195   safe1 = fRMin - rho ;
1091   safe2 = rho - fRMax ;                          1196   safe2 = rho - fRMax ;
1092   safe3 = std::fabs(p.z()) - fDz ;               1197   safe3 = std::fabs(p.z()) - fDz ;
1093                                                  1198 
1094   if ( safe1 > safe2 ) { safe = safe1; }         1199   if ( safe1 > safe2 ) { safe = safe1; }
1095   else                 { safe = safe2; }         1200   else                 { safe = safe2; }
1096   if ( safe3 > safe )  { safe = safe3; }         1201   if ( safe3 > safe )  { safe = safe3; }
1097                                                  1202 
1098   if ( (!fPhiFullTube) && ((rho) != 0.0) )    << 1203   if ( (!fPhiFullTube) && (rho) )
1099   {                                              1204   {
1100     // Psi=angle from central phi to point       1205     // Psi=angle from central phi to point
1101     //                                           1206     //
1102     cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/    1207     cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1103                                               << 1208     
1104     if ( cosPsi < cosHDPhi )                  << 1209     if ( cosPsi < std::cos(fDPhi*0.5) )
1105     {                                            1210     {
1106       // Point lies outside phi range            1211       // Point lies outside phi range
1107                                                  1212 
1108       if ( (p.y()*cosCPhi - p.x()*sinCPhi) <=    1213       if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1109       {                                          1214       {
1110         safePhi = std::fabs(p.x()*sinSPhi - p    1215         safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1111       }                                          1216       }
1112       else                                       1217       else
1113       {                                          1218       {
1114         safePhi = std::fabs(p.x()*sinEPhi - p    1219         safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1115       }                                          1220       }
1116       if ( safePhi > safe )  { safe = safePhi    1221       if ( safePhi > safe )  { safe = safePhi; }
1117     }                                            1222     }
1118   }                                              1223   }
1119   if ( safe < 0 )  { safe = 0; }                 1224   if ( safe < 0 )  { safe = 0; }
1120   return safe ;                                  1225   return safe ;
1121 }                                                1226 }
1122                                                  1227 
1123 /////////////////////////////////////////////    1228 //////////////////////////////////////////////////////////////////////////////
1124 //                                               1229 //
1125 // Calculate distance to surface of shape fro    1230 // Calculate distance to surface of shape from `inside', allowing for tolerance
1126 // - Only Calc rmax intersection if no valid     1231 // - Only Calc rmax intersection if no valid rmin intersection
1127                                                  1232 
1128 G4double G4Tubs::DistanceToOut( const G4Three    1233 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p,
1129                                 const G4Three    1234                                 const G4ThreeVector& v,
1130                                 const G4bool     1235                                 const G4bool calcNorm,
1131                                       G4bool* << 1236                                       G4bool *validNorm,
1132                                       G4Three << 1237                                       G4ThreeVector *n    ) const
1133 {                                             << 1238 {  
1134   ESide side=kNull , sider=kNull, sidephi=kNu    1239   ESide side=kNull , sider=kNull, sidephi=kNull ;
1135   G4double snxt, srd=kInfinity, sphi=kInfinit    1240   G4double snxt, srd=kInfinity, sphi=kInfinity, pdist ;
1136   G4double deltaR, t1, t2, t3, b, c, d2, roMi    1241   G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1137                                                  1242 
                                                   >> 1243   static const G4double halfCarTolerance = kCarTolerance*0.5;
                                                   >> 1244   static const G4double halfAngTolerance = kAngTolerance*0.5;
                                                   >> 1245  
1138   // Vars for phi intersection:                  1246   // Vars for phi intersection:
1139                                                  1247 
1140   G4double pDistS, compS, pDistE, compE, sphi    1248   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1141                                               << 1249  
1142   // Z plane intersection                        1250   // Z plane intersection
1143                                                  1251 
1144   if (v.z() > 0 )                                1252   if (v.z() > 0 )
1145   {                                              1253   {
1146     pdist = fDz - p.z() ;                        1254     pdist = fDz - p.z() ;
1147     if ( pdist > halfCarTolerance )              1255     if ( pdist > halfCarTolerance )
1148     {                                            1256     {
1149       snxt = pdist/v.z() ;                       1257       snxt = pdist/v.z() ;
1150       side = kPZ ;                               1258       side = kPZ ;
1151     }                                            1259     }
1152     else                                         1260     else
1153     {                                            1261     {
1154       if (calcNorm)                              1262       if (calcNorm)
1155       {                                          1263       {
1156         *n         = G4ThreeVector(0,0,1) ;      1264         *n         = G4ThreeVector(0,0,1) ;
1157         *validNorm = true ;                      1265         *validNorm = true ;
1158       }                                          1266       }
1159       return snxt = 0 ;                          1267       return snxt = 0 ;
1160     }                                            1268     }
1161   }                                              1269   }
1162   else if ( v.z() < 0 )                          1270   else if ( v.z() < 0 )
1163   {                                              1271   {
1164     pdist = fDz + p.z() ;                        1272     pdist = fDz + p.z() ;
1165                                                  1273 
1166     if ( pdist > halfCarTolerance )              1274     if ( pdist > halfCarTolerance )
1167     {                                            1275     {
1168       snxt = -pdist/v.z() ;                      1276       snxt = -pdist/v.z() ;
1169       side = kMZ ;                               1277       side = kMZ ;
1170     }                                            1278     }
1171     else                                         1279     else
1172     {                                            1280     {
1173       if (calcNorm)                              1281       if (calcNorm)
1174       {                                          1282       {
1175         *n         = G4ThreeVector(0,0,-1) ;     1283         *n         = G4ThreeVector(0,0,-1) ;
1176         *validNorm = true ;                      1284         *validNorm = true ;
1177       }                                          1285       }
1178       return snxt = 0.0 ;                        1286       return snxt = 0.0 ;
1179     }                                            1287     }
1180   }                                              1288   }
1181   else                                           1289   else
1182   {                                              1290   {
1183     snxt = kInfinity ;    // Travel perpendic    1291     snxt = kInfinity ;    // Travel perpendicular to z axis
1184     side = kNull;                                1292     side = kNull;
1185   }                                              1293   }
1186                                                  1294 
1187   // Radial Intersections                        1295   // Radial Intersections
1188   //                                             1296   //
1189   // Find intersection with cylinders at rmax    1297   // Find intersection with cylinders at rmax/rmin
1190   // Intersection point (xi,yi,zi) on line x=    1298   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1191   //                                             1299   //
1192   // Intersects with x^2+y^2=R^2                 1300   // Intersects with x^2+y^2=R^2
1193   //                                             1301   //
1194   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v    1302   // 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
1195   //                                             1303   //
1196   //            t1                t2             1304   //            t1                t2                    t3
1197                                                  1305 
1198   t1   = 1.0 - v.z()*v.z() ;      // since v     1306   t1   = 1.0 - v.z()*v.z() ;      // since v normalised
1199   t2   = p.x()*v.x() + p.y()*v.y() ;             1307   t2   = p.x()*v.x() + p.y()*v.y() ;
1200   t3   = p.x()*p.x() + p.y()*p.y() ;             1308   t3   = p.x()*p.x() + p.y()*p.y() ;
1201                                                  1309 
1202   if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fR    1310   if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fRMax*fRMax; }
1203   else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t    1311   else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }        // radius^2 on +-fDz
1204                                                  1312 
1205   if ( t1 > 0 ) // Check not parallel            1313   if ( t1 > 0 ) // Check not parallel
1206   {                                              1314   {
1207     // Calculate srd, r exit distance            1315     // Calculate srd, r exit distance
1208                                               << 1316      
1209     if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax     1317     if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1210     {                                            1318     {
1211       // Delta r not negative => leaving via     1319       // Delta r not negative => leaving via rmax
1212                                                  1320 
1213       deltaR = t3 - fRMax*fRMax ;                1321       deltaR = t3 - fRMax*fRMax ;
1214                                                  1322 
1215       // NOTE: Should use rho-fRMax<-kRadTole    1323       // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1216       // - avoid sqrt for efficiency             1324       // - avoid sqrt for efficiency
1217                                                  1325 
1218       if ( deltaR < -kRadTolerance*fRMax )       1326       if ( deltaR < -kRadTolerance*fRMax )
1219       {                                          1327       {
1220         b     = t2/t1 ;                          1328         b     = t2/t1 ;
1221         c     = deltaR/t1 ;                      1329         c     = deltaR/t1 ;
1222         d2    = b*b-c;                           1330         d2    = b*b-c;
1223         if( d2 >= 0 ) { srd = c/( -b - std::s    1331         if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1224         else          { srd = 0.; }              1332         else          { srd = 0.; }
1225         sider = kRMax ;                          1333         sider = kRMax ;
1226       }                                          1334       }
1227       else                                       1335       else
1228       {                                          1336       {
1229         // On tolerant boundary & heading out    1337         // On tolerant boundary & heading outwards (or perpendicular to)
1230         // outer radial surface -> leaving im    1338         // outer radial surface -> leaving immediately
1231                                                  1339 
1232         if ( calcNorm )                       << 1340         if ( calcNorm ) 
1233         {                                        1341         {
1234           G4double invRho = FastInverseRxy( p << 1342           *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1235           *n         = G4ThreeVector(p.x()*in << 
1236           *validNorm = true ;                    1343           *validNorm = true ;
1237         }                                        1344         }
1238         return snxt = 0 ; // Leaving by rmax     1345         return snxt = 0 ; // Leaving by rmax immediately
1239       }                                          1346       }
1240     }                                         << 1347     }             
1241     else if ( t2 < 0. ) // i.e.  t2 < 0; Poss    1348     else if ( t2 < 0. ) // i.e.  t2 < 0; Possible rmin intersection
1242     {                                            1349     {
1243       roMin2 = t3 - t2*t2/t1 ; // min ro2 of  << 1350       roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement 
1244                                                  1351 
1245       if ( (fRMin != 0.0) && (roMin2 < fRMin* << 1352       if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1246       {                                          1353       {
1247         deltaR = t3 - fRMin*fRMin ;              1354         deltaR = t3 - fRMin*fRMin ;
1248         b      = t2/t1 ;                         1355         b      = t2/t1 ;
1249         c      = deltaR/t1 ;                     1356         c      = deltaR/t1 ;
1250         d2     = b*b - c ;                       1357         d2     = b*b - c ;
1251                                                  1358 
1252         if ( d2 >= 0 )   // Leaving via rmin     1359         if ( d2 >= 0 )   // Leaving via rmin
1253         {                                        1360         {
1254           // NOTE: SHould use rho-rmin>kRadTo    1361           // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1255           // - avoid sqrt for efficiency         1362           // - avoid sqrt for efficiency
1256                                                  1363 
1257           if (deltaR > kRadTolerance*fRMin)      1364           if (deltaR > kRadTolerance*fRMin)
1258           {                                      1365           {
1259             srd = c/(-b+std::sqrt(d2));       << 1366             srd = c/(-b+std::sqrt(d2)); 
1260             sider = kRMin ;                      1367             sider = kRMin ;
1261           }                                      1368           }
1262           else                                   1369           else
1263           {                                      1370           {
1264             if ( calcNorm ) {                 << 1371             if ( calcNorm ) { *validNorm = false; }  // Concave side
1265                *validNorm = false;            << 
1266             }  // Concave side                << 
1267             return snxt = 0.0;                   1372             return snxt = 0.0;
1268           }                                      1373           }
1269         }                                        1374         }
1270         else    // No rmin intersect -> must     1375         else    // No rmin intersect -> must be rmax intersect
1271         {                                        1376         {
1272           deltaR = t3 - fRMax*fRMax ;            1377           deltaR = t3 - fRMax*fRMax ;
1273           c     = deltaR/t1 ;                    1378           c     = deltaR/t1 ;
1274           d2    = b*b-c;                         1379           d2    = b*b-c;
1275           if( d2 >=0. )                          1380           if( d2 >=0. )
1276           {                                      1381           {
1277             srd     = -b + std::sqrt(d2) ;       1382             srd     = -b + std::sqrt(d2) ;
1278             sider  = kRMax ;                     1383             sider  = kRMax ;
1279           }                                      1384           }
1280           else // Case: On the border+t2<kRad    1385           else // Case: On the border+t2<kRadTolerance
1281                //       (v is perpendicular t    1386                //       (v is perpendicular to the surface)
1282           {                                      1387           {
1283             if (calcNorm)                        1388             if (calcNorm)
1284             {                                    1389             {
1285               G4double invRho = FastInverseRx << 1390               *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1286               *n = G4ThreeVector(p.x()*invRho << 
1287               *validNorm = true ;                1391               *validNorm = true ;
1288             }                                    1392             }
1289             return snxt = 0.0;                   1393             return snxt = 0.0;
1290           }                                      1394           }
1291         }                                        1395         }
1292       }                                          1396       }
1293       else if ( roi2 > fRMax*(fRMax + kRadTol    1397       else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1294            // No rmin intersect -> must be rm    1398            // No rmin intersect -> must be rmax intersect
1295       {                                          1399       {
1296         deltaR = t3 - fRMax*fRMax ;              1400         deltaR = t3 - fRMax*fRMax ;
1297         b      = t2/t1 ;                         1401         b      = t2/t1 ;
1298         c      = deltaR/t1;                      1402         c      = deltaR/t1;
1299         d2     = b*b-c;                          1403         d2     = b*b-c;
1300         if( d2 >= 0 )                            1404         if( d2 >= 0 )
1301         {                                        1405         {
1302           srd     = -b + std::sqrt(d2) ;         1406           srd     = -b + std::sqrt(d2) ;
1303           sider  = kRMax ;                       1407           sider  = kRMax ;
1304         }                                        1408         }
1305         else // Case: On the border+t2<kRadTo    1409         else // Case: On the border+t2<kRadTolerance
1306              //       (v is perpendicular to     1410              //       (v is perpendicular to the surface)
1307         {                                        1411         {
1308           if (calcNorm)                          1412           if (calcNorm)
1309           {                                      1413           {
1310             G4double invRho = FastInverseRxy( << 1414             *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1311             *n = G4ThreeVector(p.x()*invRho,p << 
1312             *validNorm = true ;                  1415             *validNorm = true ;
1313           }                                      1416           }
1314           return snxt = 0.0;                     1417           return snxt = 0.0;
1315         }                                        1418         }
1316       }                                          1419       }
1317     }                                            1420     }
1318                                               << 1421     
1319     // Phi Intersection                          1422     // Phi Intersection
1320                                                  1423 
1321     if ( !fPhiFullTube )                         1424     if ( !fPhiFullTube )
1322     {                                            1425     {
1323       // add angle calculation with correctio << 1426       // add angle calculation with correction 
1324       // of the difference in domain of atan2    1427       // of the difference in domain of atan2 and Sphi
1325       //                                         1428       //
1326       vphi = std::atan2(v.y(),v.x()) ;           1429       vphi = std::atan2(v.y(),v.x()) ;
1327                                               << 1430      
1328       if ( vphi < fSPhi - halfAngTolerance  )    1431       if ( vphi < fSPhi - halfAngTolerance  )             { vphi += twopi; }
1329       else if ( vphi > fSPhi + fDPhi + halfAn    1432       else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1330                                                  1433 
1331                                                  1434 
1332       if ( (p.x() != 0.0) || (p.y() != 0.0) ) << 1435       if ( p.x() || p.y() )  // Check if on z axis (rho not needed later)
1333       {                                          1436       {
1334         // pDist -ve when inside                 1437         // pDist -ve when inside
1335                                                  1438 
1336         pDistS = p.x()*sinSPhi - p.y()*cosSPh    1439         pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1337         pDistE = -p.x()*sinEPhi + p.y()*cosEP    1440         pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1338                                                  1441 
1339         // Comp -ve when in direction of outw    1442         // Comp -ve when in direction of outwards normal
1340                                                  1443 
1341         compS = -sinSPhi*v.x() + cosSPhi*v.y( << 1444         compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
1342         compE =  sinEPhi*v.x() - cosEPhi*v.y( << 1445         compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
1343                                               << 1446        
1344         sidephi = kNull;                         1447         sidephi = kNull;
1345                                               << 1448         
1346         if( ( (fDPhi <= pi) && ( (pDistS <= h    1449         if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1347                               && (pDistE <= h    1450                               && (pDistE <= halfCarTolerance) ) )
1348          || ( (fDPhi >  pi) && ((pDistS <=  h << 1451          || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
1349                               || (pDistE <=   << 1452                               && (pDistE >  halfCarTolerance) ) )  )
1350         {                                        1453         {
1351           // Inside both phi *full* planes       1454           // Inside both phi *full* planes
1352                                               << 1455           
1353           if ( compS < 0 )                       1456           if ( compS < 0 )
1354           {                                      1457           {
1355             sphi = pDistS/compS ;                1458             sphi = pDistS/compS ;
1356                                               << 1459             
1357             if (sphi >= -halfCarTolerance)       1460             if (sphi >= -halfCarTolerance)
1358             {                                    1461             {
1359               xi = p.x() + sphi*v.x() ;          1462               xi = p.x() + sphi*v.x() ;
1360               yi = p.y() + sphi*v.y() ;          1463               yi = p.y() + sphi*v.y() ;
1361                                               << 1464               
1362               // Check intersecting with corr    1465               // Check intersecting with correct half-plane
1363               // (if not -> no intersect)        1466               // (if not -> no intersect)
1364               //                                 1467               //
1365               if((std::fabs(xi)<=kCarToleranc << 1468               if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
1366               {                                  1469               {
1367                 sidephi = kSPhi;                 1470                 sidephi = kSPhi;
1368                 if (((fSPhi-halfAngTolerance)    1471                 if (((fSPhi-halfAngTolerance)<=vphi)
1369                    &&((fSPhi+fDPhi+halfAngTol    1472                    &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1370                 {                                1473                 {
1371                   sphi = kInfinity;              1474                   sphi = kInfinity;
1372                 }                                1475                 }
1373               }                                  1476               }
1374               else if ( yi*cosCPhi-xi*sinCPhi    1477               else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1375               {                                  1478               {
1376                 sphi = kInfinity ;               1479                 sphi = kInfinity ;
1377               }                                  1480               }
1378               else                               1481               else
1379               {                                  1482               {
1380                 sidephi = kSPhi ;                1483                 sidephi = kSPhi ;
1381                 if ( pDistS > -halfCarToleran    1484                 if ( pDistS > -halfCarTolerance )
1382                 {                                1485                 {
1383                   sphi = 0.0 ; // Leave by sp    1486                   sphi = 0.0 ; // Leave by sphi immediately
1384                 }                             << 1487                 }    
1385               }                               << 1488               }       
1386             }                                    1489             }
1387             else                                 1490             else
1388             {                                    1491             {
1389               sphi = kInfinity ;                 1492               sphi = kInfinity ;
1390             }                                    1493             }
1391           }                                      1494           }
1392           else                                   1495           else
1393           {                                      1496           {
1394             sphi = kInfinity ;                   1497             sphi = kInfinity ;
1395           }                                      1498           }
1396                                                  1499 
1397           if ( compE < 0 )                       1500           if ( compE < 0 )
1398           {                                      1501           {
1399             sphi2 = pDistE/compE ;               1502             sphi2 = pDistE/compE ;
1400                                               << 1503             
1401             // Only check further if < starti    1504             // Only check further if < starting phi intersection
1402             //                                   1505             //
1403             if ( (sphi2 > -halfCarTolerance)     1506             if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1404             {                                    1507             {
1405               xi = p.x() + sphi2*v.x() ;         1508               xi = p.x() + sphi2*v.x() ;
1406               yi = p.y() + sphi2*v.y() ;         1509               yi = p.y() + sphi2*v.y() ;
1407                                               << 1510               
1408               if((std::fabs(xi)<=kCarToleranc << 1511               if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1409               {                                  1512               {
1410                 // Leaving via ending phi        1513                 // Leaving via ending phi
1411                 //                               1514                 //
1412                 if( (fSPhi-halfAngTolerance > << 1515                 if( !((fSPhi-halfAngTolerance <= vphi)
1413                      ||(fSPhi+fDPhi+halfAngTo << 1516                      &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1414                 {                                1517                 {
1415                   sidephi = kEPhi ;              1518                   sidephi = kEPhi ;
1416                   if ( pDistE <= -halfCarTole    1519                   if ( pDistE <= -halfCarTolerance )  { sphi = sphi2 ; }
1417                   else                           1520                   else                                { sphi = 0.0 ;   }
1418                 }                                1521                 }
1419               }                               << 1522               } 
1420               else    // Check intersecting w << 1523               else    // Check intersecting with correct half-plane 
1421                                                  1524 
1422               if ( (yi*cosCPhi-xi*sinCPhi) >=    1525               if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1423               {                                  1526               {
1424                 // Leaving via ending phi        1527                 // Leaving via ending phi
1425                 //                               1528                 //
1426                 sidephi = kEPhi ;                1529                 sidephi = kEPhi ;
1427                 if ( pDistE <= -halfCarTolera    1530                 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1428                 else                             1531                 else                               { sphi = 0.0 ;   }
1429               }                                  1532               }
1430             }                                    1533             }
1431           }                                      1534           }
1432         }                                        1535         }
1433         else                                     1536         else
1434         {                                        1537         {
1435           sphi = kInfinity ;                     1538           sphi = kInfinity ;
1436         }                                        1539         }
1437       }                                          1540       }
1438       else                                       1541       else
1439       {                                          1542       {
1440         // On z axis + travel not || to z axi    1543         // On z axis + travel not || to z axis -> if phi of vector direction
1441         // within phi of shape, Step limited     1544         // within phi of shape, Step limited by rmax, else Step =0
1442                                               << 1545                
1443         if ( (fSPhi - halfAngTolerance <= vph    1546         if ( (fSPhi - halfAngTolerance <= vphi)
1444            && (vphi <= fSPhi + fDPhi + halfAn    1547            && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1445         {                                        1548         {
1446           sphi = kInfinity ;                     1549           sphi = kInfinity ;
1447         }                                        1550         }
1448         else                                     1551         else
1449         {                                        1552         {
1450           sidephi = kSPhi ; // arbitrary      << 1553           sidephi = kSPhi ; // arbitrary 
1451           sphi    = 0.0 ;                        1554           sphi    = 0.0 ;
1452         }                                        1555         }
1453       }                                          1556       }
1454       if (sphi < snxt)  // Order intersecttio    1557       if (sphi < snxt)  // Order intersecttions
1455       {                                          1558       {
1456         snxt = sphi ;                            1559         snxt = sphi ;
1457         side = sidephi ;                         1560         side = sidephi ;
1458       }                                          1561       }
1459     }                                            1562     }
1460     if (srd < snxt)  // Order intersections      1563     if (srd < snxt)  // Order intersections
1461     {                                            1564     {
1462       snxt = srd ;                               1565       snxt = srd ;
1463       side = sider ;                             1566       side = sider ;
1464     }                                            1567     }
1465   }                                              1568   }
1466   if (calcNorm)                                  1569   if (calcNorm)
1467   {                                              1570   {
1468     switch(side)                                 1571     switch(side)
1469     {                                            1572     {
1470       case kRMax:                                1573       case kRMax:
1471         // Note: returned vector not normalis    1574         // Note: returned vector not normalised
1472         // (divide by fRMax for unit vector)     1575         // (divide by fRMax for unit vector)
1473         //                                       1576         //
1474         xi = p.x() + snxt*v.x() ;                1577         xi = p.x() + snxt*v.x() ;
1475         yi = p.y() + snxt*v.y() ;                1578         yi = p.y() + snxt*v.y() ;
1476         *n = G4ThreeVector(xi/fRMax,yi/fRMax,    1579         *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1477         *validNorm = true ;                      1580         *validNorm = true ;
1478         break ;                                  1581         break ;
1479                                                  1582 
1480       case kRMin:                                1583       case kRMin:
1481         *validNorm = false ;  // Rmin is inco    1584         *validNorm = false ;  // Rmin is inconvex
1482         break ;                                  1585         break ;
1483                                                  1586 
1484       case kSPhi:                                1587       case kSPhi:
1485         if ( fDPhi <= pi )                       1588         if ( fDPhi <= pi )
1486         {                                        1589         {
1487           *n         = G4ThreeVector(sinSPhi,    1590           *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1488           *validNorm = true ;                    1591           *validNorm = true ;
1489         }                                        1592         }
1490         else                                     1593         else
1491         {                                        1594         {
1492           *validNorm = false ;                   1595           *validNorm = false ;
1493         }                                        1596         }
1494         break ;                                  1597         break ;
1495                                                  1598 
1496       case kEPhi:                                1599       case kEPhi:
1497         if (fDPhi <= pi)                         1600         if (fDPhi <= pi)
1498         {                                        1601         {
1499           *n = G4ThreeVector(-sinEPhi,cosEPhi    1602           *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1500           *validNorm = true ;                    1603           *validNorm = true ;
1501         }                                        1604         }
1502         else                                     1605         else
1503         {                                        1606         {
1504           *validNorm = false ;                   1607           *validNorm = false ;
1505         }                                        1608         }
1506         break ;                                  1609         break ;
1507                                                  1610 
1508       case kPZ:                                  1611       case kPZ:
1509         *n         = G4ThreeVector(0,0,1) ;      1612         *n         = G4ThreeVector(0,0,1) ;
1510         *validNorm = true ;                      1613         *validNorm = true ;
1511         break ;                                  1614         break ;
1512                                                  1615 
1513       case kMZ:                                  1616       case kMZ:
1514         *n         = G4ThreeVector(0,0,-1) ;     1617         *n         = G4ThreeVector(0,0,-1) ;
1515         *validNorm = true ;                      1618         *validNorm = true ;
1516         break ;                                  1619         break ;
1517                                                  1620 
1518       default:                                   1621       default:
1519         G4cout << G4endl ;                       1622         G4cout << G4endl ;
1520         DumpInfo();                              1623         DumpInfo();
1521         std::ostringstream message;              1624         std::ostringstream message;
1522         G4long oldprc = message.precision(16) << 1625         G4int oldprc = message.precision(16);
1523         message << "Undefined side for valid     1626         message << "Undefined side for valid surface normal to solid."
1524                 << G4endl                        1627                 << G4endl
1525                 << "Position:"  << G4endl <<     1628                 << "Position:"  << G4endl << G4endl
1526                 << "p.x() = "   << p.x()/mm <    1629                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
1527                 << "p.y() = "   << p.y()/mm <    1630                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
1528                 << "p.z() = "   << p.z()/mm <    1631                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
1529                 << "Direction:" << G4endl <<     1632                 << "Direction:" << G4endl << G4endl
1530                 << "v.x() = "   << v.x() << G    1633                 << "v.x() = "   << v.x() << G4endl
1531                 << "v.y() = "   << v.y() << G    1634                 << "v.y() = "   << v.y() << G4endl
1532                 << "v.z() = "   << v.z() << G    1635                 << "v.z() = "   << v.z() << G4endl << G4endl
1533                 << "Proposed distance :" << G    1636                 << "Proposed distance :" << G4endl << G4endl
1534                 << "snxt = "    << snxt/mm <<    1637                 << "snxt = "    << snxt/mm << " mm" << G4endl ;
1535         message.precision(oldprc) ;              1638         message.precision(oldprc) ;
1536         G4Exception("G4Tubs::DistanceToOut(p,    1639         G4Exception("G4Tubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1537                     JustWarning, message);       1640                     JustWarning, message);
1538         break ;                                  1641         break ;
1539     }                                            1642     }
1540   }                                              1643   }
1541   if ( snxt<halfCarTolerance )  { snxt=0 ; }     1644   if ( snxt<halfCarTolerance )  { snxt=0 ; }
1542                                                  1645 
1543   return snxt ;                                  1646   return snxt ;
1544 }                                                1647 }
1545                                                  1648 
1546 /////////////////////////////////////////////    1649 //////////////////////////////////////////////////////////////////////////
1547 //                                               1650 //
1548 // Calculate distance (<=actual) to closest s    1651 // Calculate distance (<=actual) to closest surface of shape from inside
1549                                                  1652 
1550 G4double G4Tubs::DistanceToOut( const G4Three    1653 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
1551 {                                                1654 {
1552   G4double safe=0.0, rho, safeR1, safeR2, saf    1655   G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1553   rho = std::sqrt(p.x()*p.x() + p.y()*p.y())     1656   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1554                                                  1657 
1555 #ifdef G4CSGDEBUG                                1658 #ifdef G4CSGDEBUG
1556   if( Inside(p) == kOutside )                    1659   if( Inside(p) == kOutside )
1557   {                                              1660   {
1558     G4long oldprc = G4cout.precision(16) ;    << 1661     G4int oldprc = G4cout.precision(16) ;
1559     G4cout << G4endl ;                           1662     G4cout << G4endl ;
1560     DumpInfo();                                  1663     DumpInfo();
1561     G4cout << "Position:"  << G4endl << G4end    1664     G4cout << "Position:"  << G4endl << G4endl ;
1562     G4cout << "p.x() = "   << p.x()/mm << " m    1665     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1563     G4cout << "p.y() = "   << p.y()/mm << " m    1666     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1564     G4cout << "p.z() = "   << p.z()/mm << " m    1667     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1565     G4cout.precision(oldprc) ;                   1668     G4cout.precision(oldprc) ;
1566     G4Exception("G4Tubs::DistanceToOut(p)", "    1669     G4Exception("G4Tubs::DistanceToOut(p)", "GeomSolids1002",
1567                 JustWarning, "Point p is outs    1670                 JustWarning, "Point p is outside !?");
1568   }                                              1671   }
1569 #endif                                           1672 #endif
1570                                                  1673 
1571   if ( fRMin != 0.0 )                         << 1674   if ( fRMin )
1572   {                                              1675   {
1573     safeR1 = rho   - fRMin ;                     1676     safeR1 = rho   - fRMin ;
1574     safeR2 = fRMax - rho ;                       1677     safeR2 = fRMax - rho ;
1575                                               << 1678  
1576     if ( safeR1 < safeR2 ) { safe = safeR1 ;     1679     if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1577     else                   { safe = safeR2 ;     1680     else                   { safe = safeR2 ; }
1578   }                                              1681   }
1579   else                                           1682   else
1580   {                                              1683   {
1581     safe = fRMax - rho ;                         1684     safe = fRMax - rho ;
1582   }                                              1685   }
1583   safeZ = fDz - std::fabs(p.z()) ;               1686   safeZ = fDz - std::fabs(p.z()) ;
1584                                                  1687 
1585   if ( safeZ < safe )  { safe = safeZ ; }        1688   if ( safeZ < safe )  { safe = safeZ ; }
1586                                                  1689 
1587   // Check if phi divided, Calc distances clo    1690   // Check if phi divided, Calc distances closest phi plane
1588   //                                             1691   //
1589   if ( !fPhiFullTube )                           1692   if ( !fPhiFullTube )
1590   {                                              1693   {
1591     if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )      1694     if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1592     {                                            1695     {
1593       safePhi = -(p.x()*sinSPhi - p.y()*cosSP    1696       safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1594     }                                            1697     }
1595     else                                         1698     else
1596     {                                            1699     {
1597       safePhi = (p.x()*sinEPhi - p.y()*cosEPh    1700       safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1598     }                                            1701     }
1599     if (safePhi < safe)  { safe = safePhi ; }    1702     if (safePhi < safe)  { safe = safePhi ; }
1600   }                                              1703   }
1601   if ( safe < 0 )  { safe = 0 ; }                1704   if ( safe < 0 )  { safe = 0 ; }
1602                                                  1705 
1603   return safe ;                               << 1706   return safe ;  
                                                   >> 1707 }
                                                   >> 1708 
                                                   >> 1709 /////////////////////////////////////////////////////////////////////////
                                                   >> 1710 //
                                                   >> 1711 // Create a List containing the transformed vertices
                                                   >> 1712 // Ordering [0-3] -fDz cross section
                                                   >> 1713 //          [4-7] +fDz cross section such that [0] is below [4],
                                                   >> 1714 //                                             [1] below [5] etc.
                                                   >> 1715 // Note:
                                                   >> 1716 //  Caller has deletion resposibility
                                                   >> 1717 //  Potential improvement: For last slice, use actual ending angle
                                                   >> 1718 //                         to avoid rounding error problems.
                                                   >> 1719 
                                                   >> 1720 G4ThreeVectorList*
                                                   >> 1721 G4Tubs::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
                                                   >> 1722 {
                                                   >> 1723   G4ThreeVectorList* vertices ;
                                                   >> 1724   G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
                                                   >> 1725   G4double meshAngle, meshRMax, crossAngle,
                                                   >> 1726            cosCrossAngle, sinCrossAngle, sAngle;
                                                   >> 1727   G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
                                                   >> 1728   G4int crossSection, noCrossSections;
                                                   >> 1729 
                                                   >> 1730   // Compute no of cross-sections necessary to mesh tube
                                                   >> 1731   //
                                                   >> 1732   noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
                                                   >> 1733 
                                                   >> 1734   if ( noCrossSections < kMinMeshSections )
                                                   >> 1735   {
                                                   >> 1736     noCrossSections = kMinMeshSections ;
                                                   >> 1737   }
                                                   >> 1738   else if (noCrossSections>kMaxMeshSections)
                                                   >> 1739   {
                                                   >> 1740     noCrossSections = kMaxMeshSections ;
                                                   >> 1741   }
                                                   >> 1742   // noCrossSections = 4 ;
                                                   >> 1743 
                                                   >> 1744   meshAngle = fDPhi/(noCrossSections - 1) ;
                                                   >> 1745   // meshAngle = fDPhi/(noCrossSections) ;
                                                   >> 1746 
                                                   >> 1747   meshRMax  = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
                                                   >> 1748   meshRMin = fRMin - 100*kCarTolerance ; 
                                                   >> 1749  
                                                   >> 1750   // If complete in phi, set start angle such that mesh will be at fRMax
                                                   >> 1751   // on the x axis. Will give better extent calculations when not rotated.
                                                   >> 1752 
                                                   >> 1753   if (fPhiFullTube && (fSPhi == 0) )  { sAngle = -meshAngle*0.5 ; }
                                                   >> 1754   else                                { sAngle =  fSPhi ; }
                                                   >> 1755     
                                                   >> 1756   vertices = new G4ThreeVectorList();
                                                   >> 1757     
                                                   >> 1758   if ( vertices )
                                                   >> 1759   {
                                                   >> 1760     vertices->reserve(noCrossSections*4);
                                                   >> 1761     for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
                                                   >> 1762     {
                                                   >> 1763       // Compute coordinates of cross section at section crossSection
                                                   >> 1764 
                                                   >> 1765       crossAngle    = sAngle + crossSection*meshAngle ;
                                                   >> 1766       cosCrossAngle = std::cos(crossAngle) ;
                                                   >> 1767       sinCrossAngle = std::sin(crossAngle) ;
                                                   >> 1768 
                                                   >> 1769       rMaxX = meshRMax*cosCrossAngle ;
                                                   >> 1770       rMaxY = meshRMax*sinCrossAngle ;
                                                   >> 1771 
                                                   >> 1772       if(meshRMin <= 0.0)
                                                   >> 1773       {
                                                   >> 1774         rMinX = 0.0 ;
                                                   >> 1775         rMinY = 0.0 ;
                                                   >> 1776       }
                                                   >> 1777       else
                                                   >> 1778       {
                                                   >> 1779         rMinX = meshRMin*cosCrossAngle ;
                                                   >> 1780         rMinY = meshRMin*sinCrossAngle ;
                                                   >> 1781       }
                                                   >> 1782       vertex0 = G4ThreeVector(rMinX,rMinY,-fDz) ;
                                                   >> 1783       vertex1 = G4ThreeVector(rMaxX,rMaxY,-fDz) ;
                                                   >> 1784       vertex2 = G4ThreeVector(rMaxX,rMaxY,+fDz) ;
                                                   >> 1785       vertex3 = G4ThreeVector(rMinX,rMinY,+fDz) ;
                                                   >> 1786 
                                                   >> 1787       vertices->push_back(pTransform.TransformPoint(vertex0)) ;
                                                   >> 1788       vertices->push_back(pTransform.TransformPoint(vertex1)) ;
                                                   >> 1789       vertices->push_back(pTransform.TransformPoint(vertex2)) ;
                                                   >> 1790       vertices->push_back(pTransform.TransformPoint(vertex3)) ;
                                                   >> 1791     }
                                                   >> 1792   }
                                                   >> 1793   else
                                                   >> 1794   {
                                                   >> 1795     DumpInfo();
                                                   >> 1796     G4Exception("G4Tubs::CreateRotatedVertices()",
                                                   >> 1797                 "GeomSolids0003", FatalException,
                                                   >> 1798                 "Error in allocation of vertices. Out of memory !");
                                                   >> 1799   }
                                                   >> 1800   return vertices ;
1604 }                                                1801 }
1605                                                  1802 
1606 /////////////////////////////////////////////    1803 //////////////////////////////////////////////////////////////////////////
1607 //                                               1804 //
1608 // Stream object contents to an output stream    1805 // Stream object contents to an output stream
1609                                                  1806 
1610 G4GeometryType G4Tubs::GetEntityType() const     1807 G4GeometryType G4Tubs::GetEntityType() const
1611 {                                                1808 {
1612   return {"G4Tubs"};                          << 1809   return G4String("G4Tubs");
1613 }                                                1810 }
1614                                                  1811 
1615 /////////////////////////////////////////////    1812 //////////////////////////////////////////////////////////////////////////
1616 //                                               1813 //
1617 // Make a clone of the object                    1814 // Make a clone of the object
1618 //                                               1815 //
1619 G4VSolid* G4Tubs::Clone() const                  1816 G4VSolid* G4Tubs::Clone() const
1620 {                                                1817 {
1621   return new G4Tubs(*this);                      1818   return new G4Tubs(*this);
1622 }                                                1819 }
1623                                                  1820 
1624 /////////////////////////////////////////////    1821 //////////////////////////////////////////////////////////////////////////
1625 //                                               1822 //
1626 // Stream object contents to an output stream    1823 // Stream object contents to an output stream
1627                                                  1824 
1628 std::ostream& G4Tubs::StreamInfo( std::ostrea    1825 std::ostream& G4Tubs::StreamInfo( std::ostream& os ) const
1629 {                                                1826 {
1630   G4long oldprc = os.precision(16);           << 1827   G4int oldprc = os.precision(16);
1631   os << "------------------------------------    1828   os << "-----------------------------------------------------------\n"
1632      << "    *** Dump for solid - " << GetNam    1829      << "    *** Dump for solid - " << GetName() << " ***\n"
1633      << "    ================================    1830      << "    ===================================================\n"
1634      << " Solid type: G4Tubs\n"                  1831      << " Solid type: G4Tubs\n"
1635      << " Parameters: \n"                        1832      << " Parameters: \n"
1636      << "    inner radius : " << fRMin/mm <<     1833      << "    inner radius : " << fRMin/mm << " mm \n"
1637      << "    outer radius : " << fRMax/mm <<     1834      << "    outer radius : " << fRMax/mm << " mm \n"
1638      << "    half length Z: " << fDz/mm << "     1835      << "    half length Z: " << fDz/mm << " mm \n"
1639      << "    starting phi : " << fSPhi/degree    1836      << "    starting phi : " << fSPhi/degree << " degrees \n"
1640      << "    delta phi    : " << fDPhi/degree    1837      << "    delta phi    : " << fDPhi/degree << " degrees \n"
1641      << "------------------------------------    1838      << "-----------------------------------------------------------\n";
1642   os.precision(oldprc);                          1839   os.precision(oldprc);
1643                                                  1840 
1644   return os;                                     1841   return os;
1645 }                                                1842 }
1646                                                  1843 
1647 /////////////////////////////////////////////    1844 /////////////////////////////////////////////////////////////////////////
1648 //                                               1845 //
1649 // GetPointOnSurface                             1846 // GetPointOnSurface
1650                                                  1847 
1651 G4ThreeVector G4Tubs::GetPointOnSurface() con    1848 G4ThreeVector G4Tubs::GetPointOnSurface() const
1652 {                                                1849 {
1653   G4double Rmax = fRMax;                      << 1850   G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1654   G4double Rmin = fRMin;                      << 1851            aOne, aTwo, aThr, aFou;
1655   G4double hz = 2.*fDz;       // height       << 1852   G4double rRand;
1656   G4double lext = fDPhi*Rmax; // length of ex << 1853 
1657   G4double lint = fDPhi*Rmin; // length of in << 1854   aOne = 2.*fDz*fDPhi*fRMax;
1658                                               << 1855   aTwo = 2.*fDz*fDPhi*fRMin;
1659   // Set array of surface areas               << 1856   aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1660   //                                          << 1857   aFou = 2.*fDz*(fRMax-fRMin);
1661   G4double RRmax = Rmax * Rmax;               << 1858 
1662   G4double RRmin = Rmin * Rmin;               << 1859   phi    = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1663   G4double sbase = 0.5*fDPhi*(RRmax - RRmin); << 1860   cosphi = std::cos(phi);
1664   G4double scut = (fDPhi == twopi) ? 0. : hz* << 1861   sinphi = std::sin(phi);
1665   G4double ssurf[6] = { scut, scut, sbase, sb << 1862 
1666   ssurf[1] += ssurf[0];                       << 1863   rRand  = GetRadiusInRing(fRMin,fRMax);
1667   ssurf[2] += ssurf[1];                       << 1864   
1668   ssurf[3] += ssurf[2];                       << 1865   if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1669   ssurf[4] += ssurf[3];                       << 1866   
1670   ssurf[5] += ssurf[4];                       << 1867   chose  = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1671                                               << 1868 
1672   // Select surface                           << 1869   if( (chose >=0) && (chose < aOne) )
1673   //                                          << 1870   {
1674   G4double select = ssurf[5]*G4QuickRand();   << 1871     xRand = fRMax*cosphi;
1675   G4int k = 5;                                << 1872     yRand = fRMax*sinphi;
1676   k -= (G4int)(select <= ssurf[4]);           << 1873     zRand = RandFlat::shoot(-1.*fDz,fDz);
1677   k -= (G4int)(select <= ssurf[3]);           << 1874     return G4ThreeVector  (xRand, yRand, zRand);
1678   k -= (G4int)(select <= ssurf[2]);           << 1875   }
1679   k -= (G4int)(select <= ssurf[1]);           << 1876   else if( (chose >= aOne) && (chose < aOne + aTwo) )
1680   k -= (G4int)(select <= ssurf[0]);           << 1877   {
1681                                               << 1878     xRand = fRMin*cosphi;
1682   // Generate point on selected surface       << 1879     yRand = fRMin*sinphi;
1683   //                                          << 1880     zRand = RandFlat::shoot(-1.*fDz,fDz);
1684   switch(k)                                   << 1881     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1882   }
                                                   >> 1883   else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
                                                   >> 1884   {
                                                   >> 1885     xRand = rRand*cosphi;
                                                   >> 1886     yRand = rRand*sinphi;
                                                   >> 1887     zRand = fDz;
                                                   >> 1888     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1889   }
                                                   >> 1890   else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
                                                   >> 1891   {
                                                   >> 1892     xRand = rRand*cosphi;
                                                   >> 1893     yRand = rRand*sinphi;
                                                   >> 1894     zRand = -1.*fDz;
                                                   >> 1895     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1896   }
                                                   >> 1897   else if( (chose >= aOne + aTwo + 2.*aThr)
                                                   >> 1898         && (chose < aOne + aTwo + 2.*aThr + aFou) )
                                                   >> 1899   {
                                                   >> 1900     xRand = rRand*std::cos(fSPhi);
                                                   >> 1901     yRand = rRand*std::sin(fSPhi);
                                                   >> 1902     zRand = RandFlat::shoot(-1.*fDz,fDz);
                                                   >> 1903     return G4ThreeVector  (xRand, yRand, zRand);
                                                   >> 1904   }
                                                   >> 1905   else
1685   {                                              1906   {
1686     case 0: // start phi cut                  << 1907     xRand = rRand*std::cos(fSPhi+fDPhi);
1687     {                                         << 1908     yRand = rRand*std::sin(fSPhi+fDPhi);
1688       G4double r = Rmin + (Rmax - Rmin)*G4Qui << 1909     zRand = RandFlat::shoot(-1.*fDz,fDz);
1689       return { r*cosSPhi, r*sinSPhi, hz*G4Qui << 1910     return G4ThreeVector  (xRand, yRand, zRand);
1690     }                                         << 
1691     case 1: // end phi cut                    << 
1692     {                                         << 
1693       G4double r = Rmin + (Rmax - Rmin)*G4Qui << 
1694       return { r*cosEPhi, r*sinEPhi, hz*G4Qui << 
1695     }                                         << 
1696     case 2: // base at -dz                    << 
1697     {                                         << 
1698       G4double r = std::sqrt(RRmin + (RRmax - << 
1699       G4double phi = fSPhi + fDPhi*G4QuickRan << 
1700       return { r*std::cos(phi), r*std::sin(ph << 
1701     }                                         << 
1702     case 3: // base at +dz                    << 
1703     {                                         << 
1704       G4double r = std::sqrt(RRmin + (RRmax - << 
1705       G4double phi = fSPhi + fDPhi*G4QuickRan << 
1706       return { r*std::cos(phi), r*std::sin(ph << 
1707     }                                         << 
1708     case 4: // external lateral surface       << 
1709     {                                         << 
1710       G4double phi = fSPhi + fDPhi*G4QuickRan << 
1711       G4double z = hz*G4QuickRand() - fDz;    << 
1712       G4double x = Rmax*std::cos(phi);        << 
1713       G4double y = Rmax*std::sin(phi);        << 
1714       return { x,y,z };                       << 
1715     }                                         << 
1716     case 5: // internal lateral surface       << 
1717     {                                         << 
1718       G4double phi = fSPhi + fDPhi*G4QuickRan << 
1719       G4double z = hz*G4QuickRand() - fDz;    << 
1720       G4double x = Rmin*std::cos(phi);        << 
1721       G4double y = Rmin*std::sin(phi);        << 
1722       return { x,y,z };                       << 
1723     }                                         << 
1724   }                                              1911   }
1725   return {0., 0., 0.};                        << 
1726 }                                                1912 }
1727                                                  1913 
1728 /////////////////////////////////////////////    1914 ///////////////////////////////////////////////////////////////////////////
1729 //                                               1915 //
1730 // Methods for visualisation                     1916 // Methods for visualisation
1731                                                  1917 
1732 void G4Tubs::DescribeYourselfTo ( G4VGraphics << 1918 void G4Tubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
1733 {                                                1919 {
1734   scene.AddSolid (*this) ;                       1920   scene.AddSolid (*this) ;
1735 }                                                1921 }
1736                                                  1922 
1737 G4Polyhedron* G4Tubs::CreatePolyhedron () con << 1923 G4Polyhedron* G4Tubs::CreatePolyhedron () const 
1738 {                                                1924 {
1739   return new G4PolyhedronTubs (fRMin, fRMax,     1925   return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1740 }                                                1926 }
1741                                                  1927 
1742 #endif                                        << 1928 G4NURBS* G4Tubs::CreateNURBS () const 
                                                   >> 1929 {
                                                   >> 1930   G4NURBS* pNURBS ;
                                                   >> 1931   if (fRMin != 0) 
                                                   >> 1932   {
                                                   >> 1933     if (fPhiFullTube) 
                                                   >> 1934     {
                                                   >> 1935       pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
                                                   >> 1936     }
                                                   >> 1937     else 
                                                   >> 1938     {
                                                   >> 1939       pNURBS = new G4NURBStubesector (fRMin,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
                                                   >> 1940     }
                                                   >> 1941   }
                                                   >> 1942   else 
                                                   >> 1943   {
                                                   >> 1944     if (fPhiFullTube) 
                                                   >> 1945     {
                                                   >> 1946       pNURBS = new G4NURBScylinder (fRMax,fDz) ;
                                                   >> 1947     }
                                                   >> 1948     else 
                                                   >> 1949     {
                                                   >> 1950       const G4double epsilon = 1.e-4 ; // Cylinder sector not yet available!
                                                   >> 1951       pNURBS = new G4NURBStubesector (epsilon,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
                                                   >> 1952     }
                                                   >> 1953   }
                                                   >> 1954   return pNURBS ;
                                                   >> 1955 }
1743                                                  1956