Geant4 Cross Reference

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


  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 // G4Torus implementation                      << 
 27 //                                                 26 //
 28 // 30.10.96 V.Grichine: first implementation w <<  27 // $Id: G4Torus.cc 104316 2017-05-24 13:04:23Z gcosmo $
 29 // 26.05.00 V.Grichine: added new fuctions dev <<  28 //
 30 // 31.08.00 E.Medernach: numerical computation <<  29 // 
 31 // 11.01.01 E.Medernach: Use G4PolynomialSolve <<  30 // class G4Torus
 32 // 03.05.05 V.Grichine: SurfaceNormal(p) accor <<  31 //
 33 // 25.08.05 O.Link: new methods for DistanceTo <<  32 // Implementation
 34 // 28.10.16 E.Tcherniaev: new CalculateExtent( <<  33 //
 35 // 16.12.16 H.Burkhardt: use radius difference     34 // 16.12.16 H.Burkhardt: use radius differences and hypot to improve precision
 36 // ------------------------------------------- <<  35 // 28.10.16 E.Tcherniaev: reimplemented CalculateExtent(),
                                                   >>  36 //                      removed CreateRotatedVertices()
                                                   >>  37 // 05.04.12 M.Kelsey:   Use sqrt(r) in GetPointOnSurface() for uniform points
                                                   >>  38 // 02.10.07 T.Nikitina: Bug fixed in SolveNumericJT(), b.969:segmentation fault.
                                                   >>  39 //                      rootsrefined is used only if the number of refined roots
                                                   >>  40 //                      is the same as for primary roots. 
                                                   >>  41 // 02.10.07 T.Nikitina: Bug fixed in CalculateExtent() for case of non-rotated
                                                   >>  42 //                      full-phi torus:protect against negative value for sqrt,
                                                   >>  43 //                      correct  formula for delta.  
                                                   >>  44 // 20.11.05 V.Grichine: Bug fixed in Inside(p) for phi sections, b.810 
                                                   >>  45 // 25.08.05 O.Link: new methods for DistanceToIn/Out using JTPolynomialSolver
                                                   >>  46 // 07.06.05 V.Grichine: SurfaceNormal(p) for rho=0, Constructor as G4Cons 
                                                   >>  47 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
                                                   >>  48 // 18.03.04 V.Grichine: bug fixed in DistanceToIn(p)
                                                   >>  49 // 11.01.01 E.Medernach: Use G4PolynomialSolver to find roots
                                                   >>  50 // 03.10.00 E.Medernach: SafeNewton added
                                                   >>  51 // 31.08.00 E.Medernach: numerical computation of roots wuth bounding
                                                   >>  52 //                       volume technique
                                                   >>  53 // 26.05.00 V.Grichine: new fuctions developed by O.Cremonesi were added
                                                   >>  54 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
                                                   >>  55 // 19.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
                                                   >>  56 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
                                                   >>  57 // 30.10.96 V.Grichine: first implementation with G4Tubs elements in Fs
                                                   >>  58 //
 37                                                    59 
 38 #include "G4Torus.hh"                              60 #include "G4Torus.hh"
 39                                                    61 
 40 #if !(defined(G4GEOM_USE_UTORUS) && defined(G4     62 #if !(defined(G4GEOM_USE_UTORUS) && defined(G4GEOM_USE_SYS_USOLIDS))
 41                                                    63 
 42 #include "G4GeomTools.hh"                          64 #include "G4GeomTools.hh"
 43 #include "G4VoxelLimits.hh"                        65 #include "G4VoxelLimits.hh"
 44 #include "G4AffineTransform.hh"                    66 #include "G4AffineTransform.hh"
 45 #include "G4BoundingEnvelope.hh"                   67 #include "G4BoundingEnvelope.hh"
 46 #include "G4GeometryTolerance.hh"                  68 #include "G4GeometryTolerance.hh"
 47 #include "G4JTPolynomialSolver.hh"                 69 #include "G4JTPolynomialSolver.hh"
 48                                                    70 
 49 #include "G4VPVParameterisation.hh"                71 #include "G4VPVParameterisation.hh"
 50                                                    72 
 51 #include "meshdefs.hh"                             73 #include "meshdefs.hh"
 52                                                    74 
 53 #include "Randomize.hh"                            75 #include "Randomize.hh"
 54                                                    76 
 55 #include "G4VGraphicsScene.hh"                     77 #include "G4VGraphicsScene.hh"
 56 #include "G4Polyhedron.hh"                         78 #include "G4Polyhedron.hh"
 57                                                    79 
 58 using namespace CLHEP;                             80 using namespace CLHEP;
 59                                                    81 
 60 //////////////////////////////////////////////     82 ///////////////////////////////////////////////////////////////
 61 //                                                 83 //
 62 // Constructor - check parameters, convert ang     84 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
 63 //             - note if pdphi>2PI then reset      85 //             - note if pdphi>2PI then reset to 2PI
 64                                                    86 
 65 G4Torus::G4Torus( const G4String& pName,       <<  87 G4Torus::G4Torus( const G4String &pName,
 66                         G4double pRmin,            88                         G4double pRmin,
 67                         G4double pRmax,            89                         G4double pRmax,
 68                         G4double pRtor,            90                         G4double pRtor,
 69                         G4double pSPhi,            91                         G4double pSPhi,
 70                         G4double pDPhi )       <<  92                         G4double pDPhi)
 71   : G4CSGSolid(pName)                              93   : G4CSGSolid(pName)
 72 {                                                  94 {
 73   SetAllParameters(pRmin, pRmax, pRtor, pSPhi,     95   SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
 74 }                                                  96 }
 75                                                    97 
 76 //////////////////////////////////////////////     98 ////////////////////////////////////////////////////////////////////////////
 77 //                                                 99 //
 78 //                                                100 //
 79                                                   101 
 80 void                                              102 void
 81 G4Torus::SetAllParameters( G4double pRmin,        103 G4Torus::SetAllParameters( G4double pRmin,
 82                            G4double pRmax,        104                            G4double pRmax,
 83                            G4double pRtor,        105                            G4double pRtor,
 84                            G4double pSPhi,        106                            G4double pSPhi,
 85                            G4double pDPhi )       107                            G4double pDPhi )
 86 {                                                 108 {
 87   const G4double fEpsilon = 4.e-11;  // relati    109   const G4double fEpsilon = 4.e-11;  // relative tolerance of radii
 88                                                   110 
 89   fCubicVolume = 0.;                              111   fCubicVolume = 0.;
 90   fSurfaceArea = 0.;                              112   fSurfaceArea = 0.;
 91   fRebuildPolyhedron = true;                      113   fRebuildPolyhedron = true;
 92                                                   114 
 93   kRadTolerance = G4GeometryTolerance::GetInst    115   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
 94   kAngTolerance = G4GeometryTolerance::GetInst    116   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 95                                                   117 
 96   halfCarTolerance = 0.5*kCarTolerance;           118   halfCarTolerance = 0.5*kCarTolerance;
 97   halfAngTolerance = 0.5*kAngTolerance;           119   halfAngTolerance = 0.5*kAngTolerance;
 98                                                   120 
 99   if ( pRtor >= pRmax+1.e3*kCarTolerance )  //    121   if ( pRtor >= pRmax+1.e3*kCarTolerance )  // Check swept radius, as in G4Cons
100   {                                               122   {
101     fRtor = pRtor ;                               123     fRtor = pRtor ;
102   }                                               124   }
103   else                                            125   else
104   {                                               126   {
105     std::ostringstream message;                   127     std::ostringstream message;
106     message << "Invalid swept radius for Solid    128     message << "Invalid swept radius for Solid: " << GetName() << G4endl
107             << "        pRtor = " << pRtor <<     129             << "        pRtor = " << pRtor << ", pRmax = " << pRmax;
108     G4Exception("G4Torus::SetAllParameters()",    130     G4Exception("G4Torus::SetAllParameters()",
109                 "GeomSolids0002", FatalExcepti    131                 "GeomSolids0002", FatalException, message);
110   }                                               132   }
111                                                   133 
112   // Check radii, as in G4Cons                    134   // Check radii, as in G4Cons
113   //                                              135   //
114   if ( pRmin < pRmax - 1.e2*kCarTolerance && p    136   if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
115   {                                               137   {
116     if (pRmin >= 1.e2*kCarTolerance) { fRmin =    138     if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; }
117     else                             { fRmin =    139     else                             { fRmin = 0.0   ; }
118     fRmax = pRmax ;                               140     fRmax = pRmax ;
119   }                                               141   }
120   else                                            142   else
121   {                                               143   {
122     std::ostringstream message;                   144     std::ostringstream message;
123     message << "Invalid values of radii for So    145     message << "Invalid values of radii for Solid: " << GetName() << G4endl
124             << "        pRmin = " << pRmin <<     146             << "        pRmin = " << pRmin << ", pRmax = " << pRmax;
125     G4Exception("G4Torus::SetAllParameters()",    147     G4Exception("G4Torus::SetAllParameters()",
126                 "GeomSolids0002", FatalExcepti    148                 "GeomSolids0002", FatalException, message);
127   }                                               149   }
128                                                   150 
129   // Relative tolerances                          151   // Relative tolerances
130   //                                              152   //
131   fRminTolerance = (fRmin) != 0.0              << 153   fRminTolerance = (fRmin)
132                  ? 0.5*std::max( kRadTolerance    154                  ? 0.5*std::max( kRadTolerance, fEpsilon*(fRtor-fRmin )) : 0;
133   fRmaxTolerance = 0.5*std::max( kRadTolerance    155   fRmaxTolerance = 0.5*std::max( kRadTolerance, fEpsilon*(fRtor+fRmax) );
134                                                   156 
135   // Check angles                                 157   // Check angles
136   //                                              158   //
137   if ( pDPhi >= twopi )  { fDPhi = twopi ; }      159   if ( pDPhi >= twopi )  { fDPhi = twopi ; }
138   else                                            160   else
139   {                                               161   {
140     if (pDPhi > 0)       { fDPhi = pDPhi ; }      162     if (pDPhi > 0)       { fDPhi = pDPhi ; }
141     else                                          163     else
142     {                                             164     {
143       std::ostringstream message;                 165       std::ostringstream message;
144       message << "Invalid Z delta-Phi for Soli    166       message << "Invalid Z delta-Phi for Solid: " << GetName() << G4endl
145               << "        pDPhi = " << pDPhi;     167               << "        pDPhi = " << pDPhi;
146       G4Exception("G4Torus::SetAllParameters()    168       G4Exception("G4Torus::SetAllParameters()",
147                   "GeomSolids0002", FatalExcep    169                   "GeomSolids0002", FatalException, message);
148     }                                             170     }
149   }                                               171   }
150                                                   172   
151   // Ensure psphi in 0-2PI or -2PI-0 range if     173   // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
152   //                                              174   //
153   fSPhi = pSPhi;                                  175   fSPhi = pSPhi;
154                                                   176 
155   if (fSPhi < 0)  { fSPhi = twopi-std::fmod(st    177   if (fSPhi < 0)  { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; }
156   else            { fSPhi = std::fmod(fSPhi,tw    178   else            { fSPhi = std::fmod(fSPhi,twopi) ; }
157                                                   179 
158   if (fSPhi+fDPhi > twopi)  { fSPhi-=twopi ; }    180   if (fSPhi+fDPhi > twopi)  { fSPhi-=twopi ; }
159 }                                                 181 }
160                                                   182 
161 //////////////////////////////////////////////    183 ///////////////////////////////////////////////////////////////////////
162 //                                                184 //
163 // Fake default constructor - sets only member    185 // Fake default constructor - sets only member data and allocates memory
164 //                            for usage restri    186 //                            for usage restricted to object persistency.
165 //                                                187 //
166 G4Torus::G4Torus( __void__& a )                   188 G4Torus::G4Torus( __void__& a )
167   : G4CSGSolid(a)                              << 189   : G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
                                                   >> 190     fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
                                                   >> 191     kRadTolerance(0.), kAngTolerance(0.),
                                                   >> 192     halfCarTolerance(0.), halfAngTolerance(0.)
168 {                                                 193 {
169 }                                                 194 }
170                                                   195 
171 //////////////////////////////////////////////    196 //////////////////////////////////////////////////////////////////////
172 //                                                197 //
173 // Destructor                                     198 // Destructor
174                                                   199 
175 G4Torus::~G4Torus() = default;                 << 200 G4Torus::~G4Torus()
                                                   >> 201 {}
176                                                   202 
177 //////////////////////////////////////////////    203 //////////////////////////////////////////////////////////////////////////
178 //                                                204 //
179 // Copy constructor                               205 // Copy constructor
180                                                   206 
181 G4Torus::G4Torus(const G4Torus&) = default;    << 207 G4Torus::G4Torus(const G4Torus& rhs)
                                                   >> 208   : G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
                                                   >> 209     fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
                                                   >> 210     fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
                                                   >> 211     kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
                                                   >> 212     halfCarTolerance(rhs.halfCarTolerance),
                                                   >> 213     halfAngTolerance(rhs.halfAngTolerance)
                                                   >> 214 {
                                                   >> 215 }
182                                                   216 
183 //////////////////////////////////////////////    217 //////////////////////////////////////////////////////////////////////////
184 //                                                218 //
185 // Assignment operator                            219 // Assignment operator
186                                                   220 
187 G4Torus& G4Torus::operator = (const G4Torus& r    221 G4Torus& G4Torus::operator = (const G4Torus& rhs) 
188 {                                                 222 {
189    // Check assignment to self                    223    // Check assignment to self
190    //                                             224    //
191    if (this == &rhs)  { return *this; }           225    if (this == &rhs)  { return *this; }
192                                                   226 
193    // Copy base class data                        227    // Copy base class data
194    //                                             228    //
195    G4CSGSolid::operator=(rhs);                    229    G4CSGSolid::operator=(rhs);
196                                                   230 
197    // Copy data                                   231    // Copy data
198    //                                             232    //
199    fRmin = rhs.fRmin; fRmax = rhs.fRmax;          233    fRmin = rhs.fRmin; fRmax = rhs.fRmax;
200    fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi    234    fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
201    fRminTolerance = rhs.fRminTolerance; fRmaxT    235    fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
202    kRadTolerance = rhs.kRadTolerance; kAngTole    236    kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance;
203    halfCarTolerance = rhs.halfCarTolerance;       237    halfCarTolerance = rhs.halfCarTolerance;
204    halfAngTolerance = rhs.halfAngTolerance;       238    halfAngTolerance = rhs.halfAngTolerance;
205                                                   239 
206    return *this;                                  240    return *this;
207 }                                                 241 }
208                                                   242 
209 //////////////////////////////////////////////    243 //////////////////////////////////////////////////////////////////////
210 //                                                244 //
211 // Dispatch to parameterisation for replicatio    245 // Dispatch to parameterisation for replication mechanism dimension
212 // computation & modification.                    246 // computation & modification.
213                                                   247 
214 void G4Torus::ComputeDimensions(       G4VPVPa    248 void G4Torus::ComputeDimensions(       G4VPVParameterisation* p,
215                                  const G4int n    249                                  const G4int n,
216                                  const G4VPhys    250                                  const G4VPhysicalVolume* pRep )
217 {                                                 251 {
218   p->ComputeDimensions(*this,n,pRep);             252   p->ComputeDimensions(*this,n,pRep);
219 }                                                 253 }
220                                                   254 
221                                                   255 
222                                                   256 
223 //////////////////////////////////////////////    257 ////////////////////////////////////////////////////////////////////////////////
224 //                                                258 //
225 // Calculate the real roots to torus surface.     259 // Calculate the real roots to torus surface. 
226 // Returns negative solutions as well.            260 // Returns negative solutions as well.
227                                                   261 
228 void G4Torus::TorusRootsJT( const G4ThreeVecto    262 void G4Torus::TorusRootsJT( const G4ThreeVector& p,
229                             const G4ThreeVecto    263                             const G4ThreeVector& v,
230                                   G4double r,     264                                   G4double r,
231                                   std::vector<    265                                   std::vector<G4double>& roots ) const
232 {                                                 266 {
233                                                   267 
234   G4int i, num ;                                  268   G4int i, num ;
235   G4double c[5], srd[4], si[4] ;                  269   G4double c[5], srd[4], si[4] ;
236                                                   270 
237   G4double Rtor2 = fRtor*fRtor, r2 = r*r  ;       271   G4double Rtor2 = fRtor*fRtor, r2 = r*r  ;
238                                                   272 
239   G4double pDotV = p.x()*v.x() + p.y()*v.y() +    273   G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
240   G4double pRad2 = p.x()*p.x() + p.y()*p.y() +    274   G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
241                                                   275 
242   G4double d=pRad2 - Rtor2;                       276   G4double d=pRad2 - Rtor2;
243   c[0] = 1.0 ;                                    277   c[0] = 1.0 ;
244   c[1] = 4*pDotV ;                                278   c[1] = 4*pDotV ;
245   c[2] = 2*( (d + 2*pDotV*pDotV  - r2) + 2*Rto    279   c[2] = 2*( (d + 2*pDotV*pDotV  - r2) + 2*Rtor2*v.z()*v.z());
246   c[3] = 4*(pDotV*(d - r2) + 2*Rtor2*p.z()*v.z    280   c[3] = 4*(pDotV*(d - r2) + 2*Rtor2*p.z()*v.z()) ;
247   c[4] = (d-r2)*(d-r2) +4*Rtor2*(p.z()*p.z()-r    281   c[4] = (d-r2)*(d-r2) +4*Rtor2*(p.z()*p.z()-r2);
248                                                   282 
249   G4JTPolynomialSolver  torusEq;                  283   G4JTPolynomialSolver  torusEq;
250                                                   284 
251   num = torusEq.FindRoots( c, 4, srd, si );       285   num = torusEq.FindRoots( c, 4, srd, si );
252                                                   286   
253   for ( i = 0; i < num; ++i )                  << 287   for ( i = 0; i < num; i++ ) 
254   {                                               288   {
255     if( si[i] == 0. )  { roots.push_back(srd[i    289     if( si[i] == 0. )  { roots.push_back(srd[i]) ; }  // store real roots
256   }                                               290   }  
257                                                   291 
258   std::sort(roots.begin() , roots.end() ) ;  /    292   std::sort(roots.begin() , roots.end() ) ;  // sorting  with <
259 }                                                 293 }
260                                                   294 
261 //////////////////////////////////////////////    295 //////////////////////////////////////////////////////////////////////////////
262 //                                                296 //
263 // Interface for DistanceToIn and DistanceToOu    297 // Interface for DistanceToIn and DistanceToOut.
264 // Calls TorusRootsJT and returns the smalles     298 // Calls TorusRootsJT and returns the smalles possible distance to 
265 // the surface.                                   299 // the surface.
266 // Attention: Difference in DistanceToIn/Out f    300 // Attention: Difference in DistanceToIn/Out for points p on the surface.
267                                                   301 
268 G4double G4Torus::SolveNumericJT( const G4Thre    302 G4double G4Torus::SolveNumericJT( const G4ThreeVector& p,
269                                   const G4Thre    303                                   const G4ThreeVector& v,
270                                         G4doub    304                                         G4double r,
271                                         G4bool    305                                         G4bool IsDistanceToIn ) const
272 {                                                 306 {
273   G4double bigdist = 10*mm ;                      307   G4double bigdist = 10*mm ;
274   G4double tmin = kInfinity ;                     308   G4double tmin = kInfinity ;
275   G4double t, scal ;                              309   G4double t, scal ;
276                                                   310 
277   // calculate the distances to the intersecti    311   // calculate the distances to the intersections with the Torus
278   // from a given point p and direction v.        312   // from a given point p and direction v.
279   //                                              313   //
280   std::vector<G4double> roots ;                   314   std::vector<G4double> roots ;
281   std::vector<G4double> rootsrefined ;            315   std::vector<G4double> rootsrefined ;
282   TorusRootsJT(p,v,r,roots) ;                     316   TorusRootsJT(p,v,r,roots) ;
283                                                   317 
284   G4ThreeVector ptmp ;                            318   G4ThreeVector ptmp ;
285                                                   319 
286   // determine the smallest non-negative solut    320   // determine the smallest non-negative solution
287   //                                              321   //
288   for ( std::size_t k = 0 ; k<roots.size() ; + << 322   for ( size_t k = 0 ; k<roots.size() ; k++ )
289   {                                               323   {
290     t = roots[k] ;                                324     t = roots[k] ;
291                                                   325 
292     if ( t < -halfCarTolerance )  { continue ;    326     if ( t < -halfCarTolerance )  { continue ; }  // skip negative roots
293                                                   327 
294     if ( t > bigdist && t<kInfinity )    // pr    328     if ( t > bigdist && t<kInfinity )    // problem with big distances
295     {                                             329     {
296       ptmp = p + t*v ;                            330       ptmp = p + t*v ;
297       TorusRootsJT(ptmp,v,r,rootsrefined) ;       331       TorusRootsJT(ptmp,v,r,rootsrefined) ;
298       if ( rootsrefined.size()==roots.size() )    332       if ( rootsrefined.size()==roots.size() )
299       {                                           333       {
300         t = t + rootsrefined[k] ;                 334         t = t + rootsrefined[k] ;
301       }                                           335       }
302     }                                             336     }
303                                                   337 
304     ptmp = p + t*v ;   // calculate the positi    338     ptmp = p + t*v ;   // calculate the position of the proposed intersection
305                                                   339 
306     G4double theta = std::atan2(ptmp.y(),ptmp.    340     G4double theta = std::atan2(ptmp.y(),ptmp.x());
307                                                   341     
308     if ( fSPhi >= 0 )                             342     if ( fSPhi >= 0 )
309     {                                             343     {
310       if ( theta < - halfAngTolerance )  { the    344       if ( theta < - halfAngTolerance )  { theta += twopi; }
311       if ( (std::fabs(theta) < halfAngToleranc    345       if ( (std::fabs(theta) < halfAngTolerance)
312         && (std::fabs(fSPhi + fDPhi - twopi) <    346         && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
313       {                                           347       { 
314         theta += twopi ; // 0 <= theta < 2pi      348         theta += twopi ; // 0 <= theta < 2pi
315       }                                           349       }
316     }                                             350     }
317     if ((fSPhi <= -pi )&&(theta>halfAngToleran    351     if ((fSPhi <= -pi )&&(theta>halfAngTolerance)) { theta = theta-twopi; }
318                                                   352        
319     // We have to verify if this root is insid    353     // We have to verify if this root is inside the region between
320     // fSPhi and fSPhi + fDPhi                    354     // fSPhi and fSPhi + fDPhi
321     //                                            355     //
322     if ( (theta - fSPhi >= - halfAngTolerance)    356     if ( (theta - fSPhi >= - halfAngTolerance)
323       && (theta - (fSPhi + fDPhi) <=  halfAngT    357       && (theta - (fSPhi + fDPhi) <=  halfAngTolerance) )
324     {                                             358     {
325       // check if P is on the surface, and cal    359       // check if P is on the surface, and called from DistanceToIn
326       // DistanceToIn has to return 0.0 if par    360       // DistanceToIn has to return 0.0 if particle is going inside the solid
327                                                   361 
328       if ( IsDistanceToIn )                    << 362       if ( IsDistanceToIn == true )
329       {                                           363       {
330         if (std::fabs(t) < halfCarTolerance )     364         if (std::fabs(t) < halfCarTolerance )
331         {                                         365         {
332           // compute scalar product at positio    366           // compute scalar product at position p : v.n
333           // ( n taken from SurfaceNormal, not    367           // ( n taken from SurfaceNormal, not normalized )
334                                                   368 
335           scal = v* G4ThreeVector( p.x()*(1-fR    369           scal = v* G4ThreeVector( p.x()*(1-fRtor/std::hypot(p.x(),p.y())),
336                                    p.y()*(1-fR    370                                    p.y()*(1-fRtor/std::hypot(p.x(),p.y())),
337                                    p.z() );       371                                    p.z() );
338                                                   372 
339           // change sign in case of inner radi    373           // change sign in case of inner radius
340           //                                      374           //
341           if ( r == GetRmin() )  { scal = -sca    375           if ( r == GetRmin() )  { scal = -scal ; }
342           if ( scal < 0 )  { return 0.0  ; }      376           if ( scal < 0 )  { return 0.0  ; }
343         }                                         377         }
344       }                                           378       }
345                                                   379 
346       // check if P is on the surface, and cal    380       // check if P is on the surface, and called from DistanceToOut
347       // DistanceToIn has to return 0.0 if par    381       // DistanceToIn has to return 0.0 if particle is leaving the solid
348                                                   382 
349       if ( !IsDistanceToIn )                   << 383       if ( IsDistanceToIn == false )
350       {                                           384       {
351         if (std::fabs(t) < halfCarTolerance )     385         if (std::fabs(t) < halfCarTolerance )
352         {                                         386         {
353           // compute scalar product at positio    387           // compute scalar product at position p : v.n   
354           //                                      388           //
355           scal = v* G4ThreeVector( p.x()*(1-fR    389           scal = v* G4ThreeVector( p.x()*(1-fRtor/std::hypot(p.x(),p.y())),
356                                    p.y()*(1-fR    390                                    p.y()*(1-fRtor/std::hypot(p.x(),p.y())),
357                                    p.z() );       391                                    p.z() );
358                                                   392 
359           // change sign in case of inner radi    393           // change sign in case of inner radius
360           //                                      394           //
361           if ( r == GetRmin() )  { scal = -sca    395           if ( r == GetRmin() )  { scal = -scal ; }
362           if ( scal > 0 )  { return 0.0  ; }      396           if ( scal > 0 )  { return 0.0  ; }
363         }                                         397         }
364       }                                           398       }
365                                                   399 
366       // check if distance is larger than 1/2     400       // check if distance is larger than 1/2 kCarTolerance
367       //                                          401       //
368       if(  t > halfCarTolerance  )                402       if(  t > halfCarTolerance  )
369       {                                           403       {
370         tmin = t  ;                               404         tmin = t  ;
371         return tmin  ;                            405         return tmin  ;
372       }                                           406       }
373     }                                             407     }
374   }                                               408   }
375                                                   409 
376   return tmin;                                    410   return tmin;
377 }                                                 411 }
378                                                   412 
379 //////////////////////////////////////////////    413 /////////////////////////////////////////////////////////////////////////////
380 //                                                414 //
381 // Get bounding box                               415 // Get bounding box
382                                                   416 
383 void G4Torus::BoundingLimits(G4ThreeVector& pM    417 void G4Torus::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
384 {                                                 418 {
385   G4double rmax = GetRmax();                      419   G4double rmax = GetRmax();
386   G4double rtor = GetRtor();                      420   G4double rtor = GetRtor();
387   G4double rint = rtor - rmax;                    421   G4double rint = rtor - rmax;
388   G4double rext = rtor + rmax;                    422   G4double rext = rtor + rmax;
389   G4double dz   = rmax;                           423   G4double dz   = rmax;
390                                                   424 
391   // Find bounding box                            425   // Find bounding box
392   //                                              426   //
393   if (GetDPhi() >= twopi)                         427   if (GetDPhi() >= twopi)
394   {                                               428   {
395     pMin.set(-rext,-rext,-dz);                    429     pMin.set(-rext,-rext,-dz);
396     pMax.set( rext, rext, dz);                    430     pMax.set( rext, rext, dz);
397   }                                               431   }
398   else                                            432   else
399   {                                               433   {
400     G4TwoVector vmin,vmax;                        434     G4TwoVector vmin,vmax;
401     G4GeomTools::DiskExtent(rint,rext,            435     G4GeomTools::DiskExtent(rint,rext,
402                             GetSinStartPhi(),G    436                             GetSinStartPhi(),GetCosStartPhi(),
403                             GetSinEndPhi(),Get    437                             GetSinEndPhi(),GetCosEndPhi(),
404                             vmin,vmax);           438                             vmin,vmax);
405     pMin.set(vmin.x(),vmin.y(),-dz);              439     pMin.set(vmin.x(),vmin.y(),-dz);
406     pMax.set(vmax.x(),vmax.y(), dz);              440     pMax.set(vmax.x(),vmax.y(), dz);
407   }                                               441   }
408                                                   442 
409   // Check correctness of the bounding box        443   // Check correctness of the bounding box
410   //                                              444   //
411   if (pMin.x() >= pMax.x() || pMin.y() >= pMax    445   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
412   {                                               446   {
413     std::ostringstream message;                   447     std::ostringstream message;
414     message << "Bad bounding box (min >= max)     448     message << "Bad bounding box (min >= max) for solid: "
415             << GetName() << " !"                  449             << GetName() << " !"
416             << "\npMin = " << pMin                450             << "\npMin = " << pMin
417             << "\npMax = " << pMax;               451             << "\npMax = " << pMax;
418     G4Exception("G4Torus::BoundingLimits()", "    452     G4Exception("G4Torus::BoundingLimits()", "GeomMgt0001",
419                 JustWarning, message);            453                 JustWarning, message);
420     DumpInfo();                                   454     DumpInfo();
421   }                                               455   }
422 }                                                 456 }
423                                                   457 
424 //////////////////////////////////////////////    458 /////////////////////////////////////////////////////////////////////////////
425 //                                                459 //
426 // Calculate extent under transform and specif    460 // Calculate extent under transform and specified limit
427                                                   461 
428 G4bool G4Torus::CalculateExtent( const EAxis p    462 G4bool G4Torus::CalculateExtent( const EAxis pAxis,
429                                  const G4Voxel    463                                  const G4VoxelLimits& pVoxelLimit,
430                                  const G4Affin    464                                  const G4AffineTransform& pTransform,
431                                        G4doubl    465                                        G4double& pMin, G4double& pMax) const
432 {                                                 466 {
433   G4ThreeVector bmin, bmax;                       467   G4ThreeVector bmin, bmax;
434   G4bool exist;                                   468   G4bool exist;
435                                                   469 
436   // Get bounding box                             470   // Get bounding box
437   BoundingLimits(bmin,bmax);                      471   BoundingLimits(bmin,bmax);
438                                                   472 
439   // Check bounding box                           473   // Check bounding box
440   G4BoundingEnvelope bbox(bmin,bmax);             474   G4BoundingEnvelope bbox(bmin,bmax);
441 #ifdef G4BBOX_EXTENT                              475 #ifdef G4BBOX_EXTENT
442   return bbox.CalculateExtent(pAxis,pVoxelLimi << 476   if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
443 #endif                                            477 #endif
444   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox    478   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
445   {                                               479   {
446     return exist = pMin < pMax;                << 480     return exist = (pMin < pMax) ? true : false;
447   }                                               481   }
448                                                   482 
449   // Get parameters of the solid                  483   // Get parameters of the solid
450   G4double rmin = GetRmin();                      484   G4double rmin = GetRmin();
451   G4double rmax = GetRmax();                      485   G4double rmax = GetRmax();
452   G4double rtor = GetRtor();                      486   G4double rtor = GetRtor();
453   G4double dphi = GetDPhi();                      487   G4double dphi = GetDPhi();
454   G4double sinStart = GetSinStartPhi();           488   G4double sinStart = GetSinStartPhi();
455   G4double cosStart = GetCosStartPhi();           489   G4double cosStart = GetCosStartPhi();
456   G4double sinEnd   = GetSinEndPhi();             490   G4double sinEnd   = GetSinEndPhi();
457   G4double cosEnd   = GetCosEndPhi();             491   G4double cosEnd   = GetCosEndPhi();
458   G4double rint = rtor - rmax;                    492   G4double rint = rtor - rmax;
459   G4double rext = rtor + rmax;                    493   G4double rext = rtor + rmax;
460                                                   494 
461   // Find bounding envelope and calculate exte    495   // Find bounding envelope and calculate extent
462   //                                              496   //
463   static const G4int NPHI  = 24; // number of     497   static const G4int NPHI  = 24; // number of steps for whole torus
464   static const G4int NDISK = 16; // number of     498   static const G4int NDISK = 16; // number of steps for disk
465   static const G4double sinHalfDisk = std::sin    499   static const G4double sinHalfDisk = std::sin(pi/NDISK);
466   static const G4double cosHalfDisk = std::cos    500   static const G4double cosHalfDisk = std::cos(pi/NDISK);
467   static const G4double sinStepDisk = 2.*sinHa    501   static const G4double sinStepDisk = 2.*sinHalfDisk*cosHalfDisk;
468   static const G4double cosStepDisk = 1. - 2.*    502   static const G4double cosStepDisk = 1. - 2.*sinHalfDisk*sinHalfDisk;
469                                                   503 
470   G4double astep = (360/NPHI)*deg; // max angl    504   G4double astep = (360/NPHI)*deg; // max angle for one slice in phi
471   G4int    kphi  = (dphi <= astep) ? 1 : (G4in    505   G4int    kphi  = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
472   G4double ang   = dphi/kphi;                     506   G4double ang   = dphi/kphi;
473                                                   507 
474   G4double sinHalf = std::sin(0.5*ang);           508   G4double sinHalf = std::sin(0.5*ang);
475   G4double cosHalf = std::cos(0.5*ang);           509   G4double cosHalf = std::cos(0.5*ang);
476   G4double sinStep = 2.*sinHalf*cosHalf;          510   G4double sinStep = 2.*sinHalf*cosHalf;
477   G4double cosStep = 1. - 2.*sinHalf*sinHalf;     511   G4double cosStep = 1. - 2.*sinHalf*sinHalf;
478                                                   512 
479   // define vectors for bounding envelope         513   // define vectors for bounding envelope
480   G4ThreeVectorList pols[NDISK+1];                514   G4ThreeVectorList pols[NDISK+1];
481   for (auto & pol : pols) pol.resize(4);       << 515   for (G4int k=0; k<NDISK+1; ++k) pols[k].resize(4);
482                                                   516 
483   std::vector<const G4ThreeVectorList *> polyg    517   std::vector<const G4ThreeVectorList *> polygons;
484   polygons.resize(NDISK+1);                       518   polygons.resize(NDISK+1);
485   for (G4int k=0; k<NDISK+1; ++k) polygons[k]     519   for (G4int k=0; k<NDISK+1; ++k) polygons[k] = &pols[k];
486                                                   520 
487   // set internal and external reference circl    521   // set internal and external reference circles
488   G4TwoVector rzmin[NDISK];                       522   G4TwoVector rzmin[NDISK];
489   G4TwoVector rzmax[NDISK];                       523   G4TwoVector rzmax[NDISK];
490                                                   524 
491   if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+    525   if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+rmin*sinHalfDisk)) rmin = 0;
492   rmax /= cosHalfDisk;                            526   rmax /= cosHalfDisk;
493   G4double sinCurDisk = sinHalfDisk;              527   G4double sinCurDisk = sinHalfDisk;
494   G4double cosCurDisk = cosHalfDisk;              528   G4double cosCurDisk = cosHalfDisk;
495   for (G4int k=0; k<NDISK; ++k)                   529   for (G4int k=0; k<NDISK; ++k)
496   {                                               530   {
497     G4double rmincur = rtor + rmin*cosCurDisk;    531     G4double rmincur = rtor + rmin*cosCurDisk;
498     if (cosCurDisk < 0 && rmin > 0) rmincur /=    532     if (cosCurDisk < 0 && rmin > 0) rmincur /= cosHalf;
499     rzmin[k].set(rmincur,rmin*sinCurDisk);        533     rzmin[k].set(rmincur,rmin*sinCurDisk);
500                                                   534 
501     G4double rmaxcur = rtor + rmax*cosCurDisk;    535     G4double rmaxcur = rtor + rmax*cosCurDisk;
502     if (cosCurDisk > 0) rmaxcur /= cosHalf;       536     if (cosCurDisk > 0) rmaxcur /= cosHalf;
503     rzmax[k].set(rmaxcur,rmax*sinCurDisk);        537     rzmax[k].set(rmaxcur,rmax*sinCurDisk);
504                                                   538 
505     G4double sinTmpDisk = sinCurDisk;             539     G4double sinTmpDisk = sinCurDisk;
506     sinCurDisk = sinCurDisk*cosStepDisk + cosC    540     sinCurDisk = sinCurDisk*cosStepDisk + cosCurDisk*sinStepDisk;
507     cosCurDisk = cosCurDisk*cosStepDisk - sinT    541     cosCurDisk = cosCurDisk*cosStepDisk - sinTmpDisk*sinStepDisk;
508   }                                               542   }
509                                                   543 
510   // Loop along slices in Phi. The extent is c    544   // Loop along slices in Phi. The extent is calculated as cumulative
511   // extent of the slices                         545   // extent of the slices
512   pMin =  kInfinity;                              546   pMin =  kInfinity;
513   pMax = -kInfinity;                              547   pMax = -kInfinity;
514   G4double eminlim = pVoxelLimit.GetMinExtent(    548   G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
515   G4double emaxlim = pVoxelLimit.GetMaxExtent(    549   G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
516   G4double sinCur1 = 0, cosCur1 = 0, sinCur2 =    550   G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = 0, cosCur2 = 0;
517   for (G4int i=0; i<kphi+1; ++i)                  551   for (G4int i=0; i<kphi+1; ++i)
518   {                                               552   {
519     if (i == 0)                                   553     if (i == 0)
520     {                                             554     {
521       sinCur1 = sinStart;                         555       sinCur1 = sinStart;
522       cosCur1 = cosStart;                         556       cosCur1 = cosStart;
523       sinCur2 = sinCur1*cosHalf + cosCur1*sinH    557       sinCur2 = sinCur1*cosHalf + cosCur1*sinHalf;
524       cosCur2 = cosCur1*cosHalf - sinCur1*sinH    558       cosCur2 = cosCur1*cosHalf - sinCur1*sinHalf;
525     }                                             559     }
526     else                                          560     else
527     {                                             561     {
528       sinCur1 = sinCur2;                          562       sinCur1 = sinCur2;
529       cosCur1 = cosCur2;                          563       cosCur1 = cosCur2;
530       sinCur2 = (i == kphi) ? sinEnd : sinCur1    564       sinCur2 = (i == kphi) ? sinEnd : sinCur1*cosStep + cosCur1*sinStep;
531       cosCur2 = (i == kphi) ? cosEnd : cosCur1    565       cosCur2 = (i == kphi) ? cosEnd : cosCur1*cosStep - sinCur1*sinStep;
532     }                                             566     }
533     for (G4int k=0; k<NDISK; ++k)                 567     for (G4int k=0; k<NDISK; ++k)
534     {                                             568     {
535       G4double r1 = rzmin[k].x(), r2 = rzmax[k    569       G4double r1 = rzmin[k].x(), r2 = rzmax[k].x();
536       G4double z1 = rzmin[k].y(), z2 = rzmax[k    570       G4double z1 = rzmin[k].y(), z2 = rzmax[k].y();
537       pols[k][0].set(r1*cosCur1,r1*sinCur1,z1)    571       pols[k][0].set(r1*cosCur1,r1*sinCur1,z1);
538       pols[k][1].set(r2*cosCur1,r2*sinCur1,z2)    572       pols[k][1].set(r2*cosCur1,r2*sinCur1,z2);
539       pols[k][2].set(r2*cosCur2,r2*sinCur2,z2)    573       pols[k][2].set(r2*cosCur2,r2*sinCur2,z2);
540       pols[k][3].set(r1*cosCur2,r1*sinCur2,z1)    574       pols[k][3].set(r1*cosCur2,r1*sinCur2,z1);
541     }                                             575     }
542     pols[NDISK] = pols[0];                        576     pols[NDISK] = pols[0];
543                                                   577 
544     // get bounding box of current slice          578     // get bounding box of current slice
545     G4TwoVector vmin,vmax;                        579     G4TwoVector vmin,vmax;
546     G4GeomTools::                                 580     G4GeomTools::
547       DiskExtent(rint,rext,sinCur1,cosCur1,sin    581       DiskExtent(rint,rext,sinCur1,cosCur1,sinCur2,cosCur2,vmin,vmax);
548     bmin.setX(vmin.x()); bmin.setY(vmin.y());     582     bmin.setX(vmin.x()); bmin.setY(vmin.y());
549     bmax.setX(vmax.x()); bmax.setY(vmax.y());     583     bmax.setX(vmax.x()); bmax.setY(vmax.y());
550                                                   584 
551     // set bounding envelope for current slice    585     // set bounding envelope for current slice and adjust extent
552     G4double emin,emax;                           586     G4double emin,emax;
553     G4BoundingEnvelope benv(bmin,bmax,polygons    587     G4BoundingEnvelope benv(bmin,bmax,polygons);
554     if (!benv.CalculateExtent(pAxis,pVoxelLimi    588     if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
555     if (emin < pMin) pMin = emin;                 589     if (emin < pMin) pMin = emin;
556     if (emax > pMax) pMax = emax;                 590     if (emax > pMax) pMax = emax;
557     if (eminlim > pMin && emaxlim < pMax) brea    591     if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
558   }                                               592   }
559   return (pMin < pMax);                           593   return (pMin < pMax);
560 }                                                 594 }
561                                                   595 
562 //////////////////////////////////////////////    596 //////////////////////////////////////////////////////////////////////////////
563 //                                                597 //
564 // Return whether point inside/outside/on surf    598 // Return whether point inside/outside/on surface
565                                                   599 
566 EInside G4Torus::Inside( const G4ThreeVector&     600 EInside G4Torus::Inside( const G4ThreeVector& p ) const
567 {                                                 601 {
568   G4double r, pt2, pPhi, tolRMin, tolRMax ;       602   G4double r, pt2, pPhi, tolRMin, tolRMax ;
569                                                   603 
570   EInside in = kOutside ;                         604   EInside in = kOutside ;
571                                                   605 
572   // General precals                              606   // General precals
573   //                                              607   //
574   r   = std::hypot(p.x(),p.y());                  608   r   = std::hypot(p.x(),p.y());
575   pt2 = p.z()*p.z() + (r-fRtor)*(r-fRtor);        609   pt2 = p.z()*p.z() + (r-fRtor)*(r-fRtor);
576                                                   610 
577   if (fRmin != 0.0) tolRMin = fRmin + fRminTol << 611   if (fRmin) tolRMin = fRmin + fRminTolerance ;
578   else       tolRMin = 0 ;                        612   else       tolRMin = 0 ;
579                                                   613 
580   tolRMax = fRmax - fRmaxTolerance;               614   tolRMax = fRmax - fRmaxTolerance;
581                                                   615       
582   if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax    616   if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
583   {                                               617   {
584     if ( fDPhi == twopi || pt2 == 0 )  // on t    618     if ( fDPhi == twopi || pt2 == 0 )  // on torus swept axis
585     {                                             619     {
586       in = kInside ;                              620       in = kInside ;
587     }                                             621     }
588     else                                          622     else
589     {                                             623     {
590       // Try inner tolerant phi boundaries (=>    624       // Try inner tolerant phi boundaries (=>inside)
591       // if not inside, try outer tolerant phi    625       // if not inside, try outer tolerant phi boundaries
592                                                   626 
593       pPhi = std::atan2(p.y(),p.x()) ;            627       pPhi = std::atan2(p.y(),p.x()) ;
594                                                   628 
595       if ( pPhi < -halfAngTolerance )  { pPhi     629       if ( pPhi < -halfAngTolerance )  { pPhi += twopi ; }  // 0<=pPhi<2pi
596       if ( fSPhi >= 0 )                           630       if ( fSPhi >= 0 )
597       {                                           631       {
598         if ( (std::fabs(pPhi) < halfAngToleran    632         if ( (std::fabs(pPhi) < halfAngTolerance)
599             && (std::fabs(fSPhi + fDPhi - twop    633             && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
600         {                                         634         { 
601             pPhi += twopi ; // 0 <= pPhi < 2pi    635             pPhi += twopi ; // 0 <= pPhi < 2pi
602         }                                         636         }
603         if ( (pPhi >= fSPhi + halfAngTolerance    637         if ( (pPhi >= fSPhi + halfAngTolerance)
604             && (pPhi <= fSPhi + fDPhi - halfAn    638             && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
605         {                                         639         {
606           in = kInside ;                          640           in = kInside ;
607         }                                         641         }
608           else if ( (pPhi >= fSPhi - halfAngTo    642           else if ( (pPhi >= fSPhi - halfAngTolerance)
609                  && (pPhi <= fSPhi + fDPhi + h    643                  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
610         {                                         644         {
611           in = kSurface ;                         645           in = kSurface ;
612         }                                         646         }
613       }                                           647       }
614       else  // fSPhi < 0                          648       else  // fSPhi < 0
615       {                                           649       {
616           if ( (pPhi <= fSPhi + twopi - halfAn    650           if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
617             && (pPhi >= fSPhi + fDPhi  + halfA    651             && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) )  {;}
618           else                                    652           else
619           {                                       653           {
620             in = kSurface ;                       654             in = kSurface ;
621           }                                       655           }
622       }                                           656       }
623     }                                             657     }
624   }                                               658   }
625   else   // Try generous boundaries               659   else   // Try generous boundaries
626   {                                               660   {
627     tolRMin = fRmin - fRminTolerance ;            661     tolRMin = fRmin - fRminTolerance ;
628     tolRMax = fRmax + fRmaxTolerance ;            662     tolRMax = fRmax + fRmaxTolerance ;
629                                                   663 
630     if (tolRMin < 0 )  { tolRMin = 0 ; }          664     if (tolRMin < 0 )  { tolRMin = 0 ; }
631                                                   665 
632     if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= t    666     if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
633     {                                             667     {
634       if ( (fDPhi == twopi) || (pt2 == 0) ) //    668       if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis
635       {                                           669       {
636         in = kSurface ;                           670         in = kSurface ;
637       }                                           671       }
638       else // Try outer tolerant phi boundarie    672       else // Try outer tolerant phi boundaries only
639       {                                           673       {
640         pPhi = std::atan2(p.y(),p.x()) ;          674         pPhi = std::atan2(p.y(),p.x()) ;
641                                                   675 
642         if ( pPhi < -halfAngTolerance )  { pPh    676         if ( pPhi < -halfAngTolerance )  { pPhi += twopi ; }  // 0<=pPhi<2pi
643         if ( fSPhi >= 0 )                         677         if ( fSPhi >= 0 )
644         {                                         678         {
645           if ( (std::fabs(pPhi) < halfAngToler    679           if ( (std::fabs(pPhi) < halfAngTolerance)
646             && (std::fabs(fSPhi + fDPhi - twop    680             && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
647           {                                       681           { 
648             pPhi += twopi ; // 0 <= pPhi < 2pi    682             pPhi += twopi ; // 0 <= pPhi < 2pi
649           }                                       683           }
650           if ( (pPhi >= fSPhi - halfAngToleran    684           if ( (pPhi >= fSPhi - halfAngTolerance)
651             && (pPhi <= fSPhi + fDPhi + halfAn    685             && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
652           {                                       686           {
653             in = kSurface;                        687             in = kSurface;
654           }                                       688           }
655         }                                         689         }
656         else  // fSPhi < 0                        690         else  // fSPhi < 0
657         {                                         691         {
658           if ( (pPhi <= fSPhi + twopi - halfAn    692           if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
659             && (pPhi >= fSPhi + fDPhi  + halfA    693             && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) )  {;}
660           else                                    694           else
661           {                                       695           {
662             in = kSurface ;                       696             in = kSurface ;
663           }                                       697           }
664         }                                         698         }
665       }                                           699       }
666     }                                             700     }
667   }                                               701   }
668   return in ;                                     702   return in ;
669 }                                                 703 }
670                                                   704 
671 //////////////////////////////////////////////    705 /////////////////////////////////////////////////////////////////////////////
672 //                                                706 //
673 // Return unit normal of surface closest to p     707 // Return unit normal of surface closest to p
674 // - note if point on z axis, ignore phi divid    708 // - note if point on z axis, ignore phi divided sides
675 // - unsafe if point close to z axis a rmin=0     709 // - unsafe if point close to z axis a rmin=0 - no explicit checks
676                                                   710 
677 G4ThreeVector G4Torus::SurfaceNormal( const G4    711 G4ThreeVector G4Torus::SurfaceNormal( const G4ThreeVector& p ) const
678 {                                                 712 {
679   G4int noSurfaces = 0;                           713   G4int noSurfaces = 0;  
680   G4double rho, pt, pPhi;                         714   G4double rho, pt, pPhi;
681   G4double distRMin = kInfinity;                  715   G4double distRMin = kInfinity;
682   G4double distSPhi = kInfinity, distEPhi = kI    716   G4double distSPhi = kInfinity, distEPhi = kInfinity;
683                                                   717 
684   // To cope with precision loss                  718   // To cope with precision loss
685   //                                              719   //
686   const G4double delta = std::max(10.0*kCarTol    720   const G4double delta = std::max(10.0*kCarTolerance,
687                                   1.0e-8*(fRto    721                                   1.0e-8*(fRtor+fRmax));
688   const G4double dAngle = 10.0*kAngTolerance;     722   const G4double dAngle = 10.0*kAngTolerance;
689                                                   723 
690   G4ThreeVector nR, nPs, nPe;                     724   G4ThreeVector nR, nPs, nPe;
691   G4ThreeVector norm, sumnorm(0.,0.,0.);          725   G4ThreeVector norm, sumnorm(0.,0.,0.);
692                                                   726 
693   rho = std::hypot(p.x(),p.y());                  727   rho = std::hypot(p.x(),p.y());
694   pt  = std::hypot(p.z(),rho-fRtor);              728   pt  = std::hypot(p.z(),rho-fRtor);
695                                                   729 
696   G4double  distRMax = std::fabs(pt - fRmax);     730   G4double  distRMax = std::fabs(pt - fRmax);
697   if(fRmin != 0.0) distRMin = std::fabs(pt - f << 731   if(fRmin) distRMin = std::fabs(pt - fRmin);
698                                                   732 
699   if( rho > delta && pt != 0.0 )                  733   if( rho > delta && pt != 0.0 )
700   {                                               734   {
701     G4double redFactor= (rho-fRtor)/rho;          735     G4double redFactor= (rho-fRtor)/rho;
702     nR = G4ThreeVector( p.x()*redFactor,  // p    736     nR = G4ThreeVector( p.x()*redFactor,  // p.x()*(1.-fRtor/rho),
703                         p.y()*redFactor,  // p    737                         p.y()*redFactor,  // p.y()*(1.-fRtor/rho),
704                         p.z()          );         738                         p.z()          );
705     nR *= 1.0/pt;                                 739     nR *= 1.0/pt;
706   }                                               740   }
707                                                   741 
708   if ( fDPhi < twopi ) // && rho ) // old limi    742   if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
709   {                                               743   {
710     if ( rho != 0.0 )                          << 744     if ( rho )
711     {                                             745     {
712       pPhi = std::atan2(p.y(),p.x());             746       pPhi = std::atan2(p.y(),p.x());
713                                                   747 
714       if(pPhi < fSPhi-delta)            { pPhi    748       if(pPhi < fSPhi-delta)            { pPhi += twopi; }
715       else if(pPhi > fSPhi+fDPhi+delta) { pPhi    749       else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
716                                                   750 
717       distSPhi = std::fabs( pPhi - fSPhi );       751       distSPhi = std::fabs( pPhi - fSPhi );
718       distEPhi = std::fabs(pPhi-fSPhi-fDPhi);     752       distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
719     }                                             753     }
720     nPs = G4ThreeVector(std::sin(fSPhi),-std::    754     nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
721     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi)    755     nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
722   }                                               756   } 
723   if( distRMax <= delta )                         757   if( distRMax <= delta )
724   {                                               758   {
725     ++noSurfaces;                              << 759     noSurfaces ++;
726     sumnorm += nR;                                760     sumnorm += nR;
727   }                                               761   }
728   else if( (fRmin != 0.0) && (distRMin <= delt << 762   else if( fRmin && (distRMin <= delta) ) // Must not be on both Outer and Inner
729   {                                               763   {
730     ++noSurfaces;                              << 764     noSurfaces ++;
731     sumnorm -= nR;                                765     sumnorm -= nR;
732   }                                               766   }
733                                                   767 
734   //  To be on one of the 'phi' surfaces,         768   //  To be on one of the 'phi' surfaces,
735   //  it must be within the 'tube' - with tole    769   //  it must be within the 'tube' - with tolerance
736                                                   770 
737   if( (fDPhi < twopi) && (fRmin-delta <= pt) &    771   if( (fDPhi < twopi) && (fRmin-delta <= pt) && (pt <= (fRmax+delta)) )
738   {                                               772   {
739     if (distSPhi <= dAngle)                       773     if (distSPhi <= dAngle)
740     {                                             774     {
741       ++noSurfaces;                            << 775       noSurfaces ++;
742       sumnorm += nPs;                             776       sumnorm += nPs;
743     }                                             777     }
744     if (distEPhi <= dAngle)                       778     if (distEPhi <= dAngle) 
745     {                                             779     {
746       ++noSurfaces;                            << 780       noSurfaces ++;
747       sumnorm += nPe;                             781       sumnorm += nPe;
748     }                                             782     }
749   }                                               783   }
750   if ( noSurfaces == 0 )                          784   if ( noSurfaces == 0 )
751   {                                               785   {
752 #ifdef G4CSGDEBUG                                 786 #ifdef G4CSGDEBUG
753      G4ExceptionDescription ed;                   787      G4ExceptionDescription ed;
754      ed.precision(16);                            788      ed.precision(16);
755                                                   789 
756      EInside  inIt= Inside( p );                  790      EInside  inIt= Inside( p );
757                                                   791      
758      if( inIt != kSurface )                       792      if( inIt != kSurface )
759      {                                            793      {
760         ed << " ERROR>  Surface Normal was cal    794         ed << " ERROR>  Surface Normal was called for Torus,"
761            << " with point not on surface." <<    795            << " with point not on surface." << G4endl;
762      }                                            796      }
763      else                                         797      else
764      {                                            798      {
765         ed << " ERROR>  Surface Normal has not    799         ed << " ERROR>  Surface Normal has not found a surface, "
766            << " despite the point being on the    800            << " despite the point being on the surface. " <<G4endl;
767      }                                            801      }
768                                                   802 
769      if( inIt != kInside)                         803      if( inIt != kInside)
770      {                                            804      {
771          ed << " Safety (Dist To In)  = " << D    805          ed << " Safety (Dist To In)  = " << DistanceToIn(p) << G4endl;
772      }                                            806      }
773      if( inIt != kOutside)                        807      if( inIt != kOutside)
774      {                                            808      {
775          ed << " Safety (Dist to Out) = " << D    809          ed << " Safety (Dist to Out) = " << DistanceToOut(p) << G4endl;
776      }                                            810      }
777      ed << " Coordinates of point : " << p <<     811      ed << " Coordinates of point : " << p << G4endl;
778      ed << " Parameters  of solid : " << G4end    812      ed << " Parameters  of solid : " << G4endl << *this << G4endl;
779                                                   813 
780      if( inIt == kSurface )                       814      if( inIt == kSurface )
781      {                                            815      {
782         G4Exception("G4Torus::SurfaceNormal(p)    816         G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
783                     JustWarning, ed,              817                     JustWarning, ed,
784                     "Failing to find normal, e    818                     "Failing to find normal, even though point is on surface!");
785      }                                            819      }
786      else                                         820      else
787      {                                            821      {
788         static const char* NameInside[3]= { "I    822         static const char* NameInside[3]= { "Inside", "Surface", "Outside" };
789         ed << "  The point is " << NameInside[    823         ed << "  The point is " << NameInside[inIt] << " the solid. "<< G4endl;
790         G4Exception("G4Torus::SurfaceNormal(p)    824         G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
791                     JustWarning, ed, "Point p     825                     JustWarning, ed, "Point p is not on surface !?" );
792      }                                            826      }
793 #endif                                            827 #endif
794      norm = ApproxSurfaceNormal(p);               828      norm = ApproxSurfaceNormal(p);
795   }                                               829   }
796   else if ( noSurfaces == 1 )  { norm = sumnor    830   else if ( noSurfaces == 1 )  { norm = sumnorm; }
797   else                         { norm = sumnor    831   else                         { norm = sumnorm.unit(); }
798                                                   832 
799   return norm ;                                   833   return norm ;
800 }                                                 834 }
801                                                   835 
802 //////////////////////////////////////////////    836 //////////////////////////////////////////////////////////////////////////////
803 //                                                837 //
804 // Algorithm for SurfaceNormal() following the    838 // Algorithm for SurfaceNormal() following the original specification
805 // for points not on the surface                  839 // for points not on the surface
806                                                   840 
807 G4ThreeVector G4Torus::ApproxSurfaceNormal( co    841 G4ThreeVector G4Torus::ApproxSurfaceNormal( const G4ThreeVector& p ) const
808 {                                                 842 {
809   ENorm side ;                                    843   ENorm side ;
810   G4ThreeVector norm;                             844   G4ThreeVector norm;
811   G4double rho,pt,phi;                            845   G4double rho,pt,phi;
812   G4double distRMin,distRMax,distSPhi,distEPhi    846   G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
813                                                   847 
814   rho = std::hypot(p.x(),p.y());                  848   rho = std::hypot(p.x(),p.y());
815   pt  = std::hypot(p.z(),rho-fRtor);              849   pt  = std::hypot(p.z(),rho-fRtor);
816                                                   850 
817 #ifdef G4CSGDEBUG                                 851 #ifdef G4CSGDEBUG
818   G4cout << " G4Torus::ApproximateSurfaceNorma    852   G4cout << " G4Torus::ApproximateSurfaceNormal called for point " << p
819          << G4endl;                               853          << G4endl;
820 #endif                                            854 #endif
821                                                   855    
822   distRMax = std::fabs(pt - fRmax) ;              856   distRMax = std::fabs(pt - fRmax) ;
823                                                   857 
824   if(fRmin != 0.0)  // First minimum radius    << 858   if(fRmin)  // First minimum radius
825   {                                               859   {
826     distRMin = std::fabs(pt - fRmin) ;            860     distRMin = std::fabs(pt - fRmin) ;
827                                                   861 
828     if (distRMin < distRMax)                      862     if (distRMin < distRMax)
829     {                                             863     {
830       distMin = distRMin ;                        864       distMin = distRMin ;
831       side    = kNRMin ;                          865       side    = kNRMin ;
832     }                                             866     }
833     else                                          867     else
834     {                                             868     {
835       distMin = distRMax ;                        869       distMin = distRMax ;
836       side    = kNRMax ;                          870       side    = kNRMax ;
837     }                                             871     }
838   }                                               872   }
839   else                                            873   else
840   {                                               874   {
841     distMin = distRMax ;                          875     distMin = distRMax ;
842     side    = kNRMax ;                            876     side    = kNRMax ;
843   }                                               877   }    
844   if ( (fDPhi < twopi) && (rho != 0.0) )       << 878   if ( (fDPhi < twopi) && rho )
845   {                                               879   {
846     phi = std::atan2(p.y(),p.x()) ; // Protect    880     phi = std::atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho!=0)
847                                                   881 
848     if (phi < 0)  { phi += twopi ; }              882     if (phi < 0)  { phi += twopi ; }
849                                                   883 
850     if (fSPhi < 0 )  { distSPhi = std::fabs(ph    884     if (fSPhi < 0 )  { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; }
851     else             { distSPhi = std::fabs(ph    885     else             { distSPhi = std::fabs(phi-fSPhi)*rho ; }
852                                                   886 
853     distEPhi = std::fabs(phi - fSPhi - fDPhi)*    887     distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
854                                                   888 
855     if (distSPhi < distEPhi) // Find new minim    889     if (distSPhi < distEPhi) // Find new minimum
856     {                                             890     {
857       if (distSPhi<distMin) side = kNSPhi ;       891       if (distSPhi<distMin) side = kNSPhi ;
858     }                                             892     }
859     else                                          893     else
860     {                                             894     {
861       if (distEPhi < distMin)  { side = kNEPhi    895       if (distEPhi < distMin)  { side = kNEPhi ; }
862     }                                             896     }
863   }                                               897   }  
864   switch (side)                                   898   switch (side)
865   {                                               899   {
866     case kNRMin:      // Inner radius             900     case kNRMin:      // Inner radius
867       norm = G4ThreeVector( -p.x()*(1-fRtor/rh    901       norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
868                             -p.y()*(1-fRtor/rh    902                             -p.y()*(1-fRtor/rho)/pt,
869                             -p.z()/pt             903                             -p.z()/pt                 ) ;
870       break ;                                     904       break ;
871     case kNRMax:      // Outer radius             905     case kNRMax:      // Outer radius
872       norm = G4ThreeVector( p.x()*(1-fRtor/rho    906       norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
873                             p.y()*(1-fRtor/rho    907                             p.y()*(1-fRtor/rho)/pt,
874                             p.z()/pt              908                             p.z()/pt                  ) ;
875       break;                                      909       break;
876     case kNSPhi:                                  910     case kNSPhi:
877       norm = G4ThreeVector(std::sin(fSPhi),-st    911       norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
878       break;                                      912       break;
879     case kNEPhi:                                  913     case kNEPhi:
880       norm = G4ThreeVector(-std::sin(fSPhi+fDP    914       norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
881       break;                                      915       break;
882     default:          // Should never reach th    916     default:          // Should never reach this case ...
883       DumpInfo();                                 917       DumpInfo();
884       G4Exception("G4Torus::ApproxSurfaceNorma    918       G4Exception("G4Torus::ApproxSurfaceNormal()",
885                   "GeomSolids1002", JustWarnin    919                   "GeomSolids1002", JustWarning,
886                   "Undefined side for valid su    920                   "Undefined side for valid surface normal to solid.");
887       break ;                                     921       break ;
888   }                                               922   } 
889   return norm ;                                   923   return norm ;
890 }                                                 924 }
891                                                   925 
892 //////////////////////////////////////////////    926 ///////////////////////////////////////////////////////////////////////
893 //                                                927 //
894 // Calculate distance to shape from outside, a    928 // Calculate distance to shape from outside, along normalised vector
895 // - return kInfinity if no intersection, or i    929 // - return kInfinity if no intersection, or intersection distance <= tolerance
896 //                                                930 //
897 // - Compute the intersection with the z plane    931 // - Compute the intersection with the z planes 
898 //        - if at valid r, phi, return            932 //        - if at valid r, phi, return
899 //                                                933 //
900 // -> If point is outer outer radius, compute     934 // -> If point is outer outer radius, compute intersection with rmax
901 //        - if at valid phi,z return              935 //        - if at valid phi,z return
902 //                                                936 //
903 // -> Compute intersection with inner radius,     937 // -> Compute intersection with inner radius, taking largest +ve root
904 //        - if valid (phi), save intersction      938 //        - if valid (phi), save intersction
905 //                                                939 //
906 //    -> If phi segmented, compute intersectio    940 //    -> If phi segmented, compute intersections with phi half planes
907 //        - return smallest of valid phi inter    941 //        - return smallest of valid phi intersections and
908 //          inner radius intersection             942 //          inner radius intersection
909 //                                                943 //
910 // NOTE:                                          944 // NOTE:
911 // - Precalculations for phi trigonometry are     945 // - Precalculations for phi trigonometry are Done `just in time'
912 // - `if valid' implies tolerant checking of i    946 // - `if valid' implies tolerant checking of intersection points
913                                                   947 
914 G4double G4Torus::DistanceToIn( const G4ThreeV    948 G4double G4Torus::DistanceToIn( const G4ThreeVector& p,
915                                 const G4ThreeV    949                                 const G4ThreeVector& v ) const
916 {                                                 950 {
917   // Get bounding box of full torus            << 
918   //                                           << 
919   G4double boxDx  = fRtor + fRmax;             << 
920   G4double boxDy  = boxDx;                     << 
921   G4double boxDz  = fRmax;                     << 
922   G4double boxMax = boxDx;                     << 
923   G4double boxMin = boxDz;                     << 
924                                                   951 
925   // Check if point is traveling away          << 
926   //                                           << 
927   G4double distX = std::abs(p.x()) - boxDx;    << 
928   G4double distY = std::abs(p.y()) - boxDy;    << 
929   G4double distZ = std::abs(p.z()) - boxDz;    << 
930   if (distX >= -halfCarTolerance && p.x()*v.x( << 
931   if (distY >= -halfCarTolerance && p.y()*v.y( << 
932   if (distZ >= -halfCarTolerance && p.z()*v.z( << 
933                                                << 
934   // Calculate safety distance to bounding box << 
935   // If point is too far, move it closer and c << 
936   //                                           << 
937   G4double Dmax = 32*boxMax;                   << 
938   G4double safe = std::max(std::max(distX,dist << 
939   if (safe > Dmax)                             << 
940   {                                            << 
941     G4double dist = safe - 1.e-8*safe - boxMin << 
942     dist += DistanceToIn(p + dist*v, v);       << 
943     return (dist >= kInfinity) ? kInfinity : d << 
944   }                                            << 
945                                                << 
946   // Find intersection with torus              << 
947   //                                           << 
948   G4double snxt=kInfinity, sphi=kInfinity; //     952   G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value
949                                                   953 
950   G4double  sd[4] ;                               954   G4double  sd[4] ;
951                                                   955 
952   // Precalculated trig for phi intersections     956   // Precalculated trig for phi intersections - used by r,z intersections to
953   //                                              957   //                                            check validity
954                                                   958 
955   G4bool seg;        // true if segmented         959   G4bool seg;        // true if segmented
956   G4double hDPhi;    // half dphi                 960   G4double hDPhi;    // half dphi
957   G4double cPhi,sinCPhi=0.,cosCPhi=0.;  // cen    961   G4double cPhi,sinCPhi=0.,cosCPhi=0.;  // central phi
958                                                   962 
959   G4double tolORMin2;  // `generous' radii squ    963   G4double tolORMin2;  // `generous' radii squared
960   G4double tolORMax2;                             964   G4double tolORMax2;
961                                                   965 
962   G4double Dist,xi,yi,zi,rhoi,it2; // Intersec    966   G4double Dist,xi,yi,zi,rhoi,it2; // Intersection point variables
963                                                   967 
964   G4double Comp;                                  968   G4double Comp;
965   G4double cosSPhi,sinSPhi;       // Trig for     969   G4double cosSPhi,sinSPhi;       // Trig for phi start intersect
966   G4double ePhi,cosEPhi,sinEPhi;  // for phi e    970   G4double ePhi,cosEPhi,sinEPhi;  // for phi end intersect
967                                                   971 
968   // Set phi divided flag and precalcs            972   // Set phi divided flag and precalcs
969   //                                              973   //
970   if ( fDPhi < twopi )                            974   if ( fDPhi < twopi )
971   {                                               975   {
972     seg        = true ;                           976     seg        = true ;
973     hDPhi      = 0.5*fDPhi ;    // half delta     977     hDPhi      = 0.5*fDPhi ;    // half delta phi
974     cPhi       = fSPhi + hDPhi ;                  978     cPhi       = fSPhi + hDPhi ;
975     sinCPhi    = std::sin(cPhi) ;                 979     sinCPhi    = std::sin(cPhi) ;
976     cosCPhi    = std::cos(cPhi) ;                 980     cosCPhi    = std::cos(cPhi) ;
977   }                                               981   }
978   else                                            982   else
979   {                                               983   {
980     seg = false ;                                 984     seg = false ;
981   }                                               985   }
982                                                   986 
983   if (fRmin > fRminTolerance) // Calculate tol    987   if (fRmin > fRminTolerance) // Calculate tolerant rmin and rmax
984   {                                               988   {
985     tolORMin2 = (fRmin - fRminTolerance)*(fRmi    989     tolORMin2 = (fRmin - fRminTolerance)*(fRmin - fRminTolerance) ;
986   }                                               990   }
987   else                                            991   else
988   {                                               992   {
989     tolORMin2 = 0 ;                               993     tolORMin2 = 0 ;
990   }                                               994   }
991   tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax     995   tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax + fRmaxTolerance) ;
992                                                   996 
993   // Intersection with Rmax (possible return)     997   // Intersection with Rmax (possible return) and Rmin (must also check phi)
994                                                   998 
995   snxt = SolveNumericJT(p,v,fRmax,true);          999   snxt = SolveNumericJT(p,v,fRmax,true);
996                                                   1000 
997   if (fRmin != 0.0)  // Possible Rmin intersec << 1001   if (fRmin)  // Possible Rmin intersection
998   {                                               1002   {
999     sd[0] = SolveNumericJT(p,v,fRmin,true);       1003     sd[0] = SolveNumericJT(p,v,fRmin,true);
1000     if ( sd[0] < snxt )  { snxt = sd[0] ; }      1004     if ( sd[0] < snxt )  { snxt = sd[0] ; }
1001   }                                              1005   }
1002                                                  1006 
1003   //                                             1007   //
1004   // Phi segment intersection                    1008   // Phi segment intersection
1005   //                                             1009   //
1006   // o Tolerant of points inside phi planes b    1010   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1007   //                                             1011   //
1008   // o NOTE: Large duplication of code betwee    1012   // o NOTE: Large duplication of code between sphi & ephi checks
1009   //         -> only diffs: sphi -> ephi, Com    1013   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1010   //            intersection check <=0 -> >=0    1014   //            intersection check <=0 -> >=0
1011   //         -> use some form of loop Constru    1015   //         -> use some form of loop Construct ?
1012                                                  1016 
1013   if (seg)                                       1017   if (seg)
1014   {                                              1018   {
1015     sinSPhi = std::sin(fSPhi) ; // First phi     1019     sinSPhi = std::sin(fSPhi) ; // First phi surface ('S'tarting phi)
1016     cosSPhi = std::cos(fSPhi) ;                  1020     cosSPhi = std::cos(fSPhi) ;
1017     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;    1021     Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;  // Component in outwards
1018                                                  1022                                                // normal direction
1019     if (Comp < 0 )                               1023     if (Comp < 0 )
1020     {                                            1024     {
1021       Dist = (p.y()*cosSPhi - p.x()*sinSPhi)     1025       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1022                                                  1026 
1023       if (Dist < halfCarTolerance)               1027       if (Dist < halfCarTolerance)
1024       {                                          1028       {
1025         sphi = Dist/Comp ;                       1029         sphi = Dist/Comp ;
1026         if (sphi < snxt)                         1030         if (sphi < snxt)
1027         {                                        1031         {
1028           if ( sphi < 0 )  { sphi = 0 ; }        1032           if ( sphi < 0 )  { sphi = 0 ; }
1029                                                  1033 
1030           xi    = p.x() + sphi*v.x() ;           1034           xi    = p.x() + sphi*v.x() ;
1031           yi    = p.y() + sphi*v.y() ;           1035           yi    = p.y() + sphi*v.y() ;
1032           zi    = p.z() + sphi*v.z() ;           1036           zi    = p.z() + sphi*v.z() ;
1033           rhoi = std::hypot(xi,yi);              1037           rhoi = std::hypot(xi,yi);
1034           it2 = zi*zi + (rhoi-fRtor)*(rhoi-fR    1038           it2 = zi*zi + (rhoi-fRtor)*(rhoi-fRtor);
1035                                                  1039 
1036           if ( it2 >= tolORMin2 && it2 <= tol    1040           if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
1037           {                                      1041           {
1038             // r intersection is good - check    1042             // r intersection is good - check intersecting
1039             // with correct half-plane           1043             // with correct half-plane
1040             //                                   1044             //
1041             if ((yi*cosCPhi-xi*sinCPhi)<=0)      1045             if ((yi*cosCPhi-xi*sinCPhi)<=0)  { snxt=sphi; }
1042           }                                      1046           }
1043         }                                        1047         }
1044       }                                          1048       }
1045     }                                            1049     }
1046     ePhi=fSPhi+fDPhi;    // Second phi surfac    1050     ePhi=fSPhi+fDPhi;    // Second phi surface ('E'nding phi)
1047     sinEPhi=std::sin(ePhi);                      1051     sinEPhi=std::sin(ePhi);
1048     cosEPhi=std::cos(ePhi);                      1052     cosEPhi=std::cos(ePhi);
1049     Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);         1053     Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
1050                                                  1054         
1051     if ( Comp < 0 )   // Component in outward    1055     if ( Comp < 0 )   // Component in outwards normal dirn
1052     {                                            1056     {
1053       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi)    1057       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1054                                                  1058 
1055       if (Dist < halfCarTolerance )              1059       if (Dist < halfCarTolerance )
1056       {                                          1060       {
1057         sphi = Dist/Comp ;                       1061         sphi = Dist/Comp ;
1058                                                  1062 
1059         if (sphi < snxt )                        1063         if (sphi < snxt )
1060         {                                        1064         {
1061           if (sphi < 0 )  { sphi = 0 ; }         1065           if (sphi < 0 )  { sphi = 0 ; }
1062                                                  1066        
1063           xi    = p.x() + sphi*v.x() ;           1067           xi    = p.x() + sphi*v.x() ;
1064           yi    = p.y() + sphi*v.y() ;           1068           yi    = p.y() + sphi*v.y() ;
1065           zi    = p.z() + sphi*v.z() ;           1069           zi    = p.z() + sphi*v.z() ;
1066           rhoi = std::hypot(xi,yi);              1070           rhoi = std::hypot(xi,yi);
1067           it2 = zi*zi + (rhoi-fRtor)*(rhoi-fR    1071           it2 = zi*zi + (rhoi-fRtor)*(rhoi-fRtor);
1068                                                  1072 
1069           if (it2 >= tolORMin2 && it2 <= tolO    1073           if (it2 >= tolORMin2 && it2 <= tolORMax2)
1070           {                                      1074           {
1071             // z and r intersections good - c    1075             // z and r intersections good - check intersecting
1072             // with correct half-plane           1076             // with correct half-plane
1073             //                                   1077             //
1074             if ((yi*cosCPhi-xi*sinCPhi)>=0)      1078             if ((yi*cosCPhi-xi*sinCPhi)>=0)  { snxt=sphi; }
1075           }                                      1079           }    
1076         }                                        1080         }
1077       }                                          1081       }
1078     }                                            1082     }
1079   }                                              1083   }
1080   if(snxt < halfCarTolerance)  { snxt = 0.0 ;    1084   if(snxt < halfCarTolerance)  { snxt = 0.0 ; }
1081                                                  1085 
1082   return snxt ;                                  1086   return snxt ;
1083 }                                                1087 }
1084                                                  1088 
1085 /////////////////////////////////////////////    1089 /////////////////////////////////////////////////////////////////////////////
1086 //                                               1090 //
1087 // Calculate distance (<= actual) to closest     1091 // Calculate distance (<= actual) to closest surface of shape from outside
1088 // - Calculate distance to z, radial planes      1092 // - Calculate distance to z, radial planes
1089 // - Only to phi planes if outside phi extent    1093 // - Only to phi planes if outside phi extent
1090 // - Return 0 if point inside                    1094 // - Return 0 if point inside
1091                                                  1095 
1092 G4double G4Torus::DistanceToIn( const G4Three    1096 G4double G4Torus::DistanceToIn( const G4ThreeVector& p ) const
1093 {                                                1097 {
1094   G4double safe=0.0, safe1, safe2 ;              1098   G4double safe=0.0, safe1, safe2 ;
1095   G4double phiC, cosPhiC, sinPhiC, safePhi, e    1099   G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1096   G4double rho, pt ;                             1100   G4double rho, pt ;
1097                                                  1101   
1098   rho = std::hypot(p.x(),p.y());                 1102   rho = std::hypot(p.x(),p.y());
1099   pt  = std::hypot(p.z(),rho-fRtor);             1103   pt  = std::hypot(p.z(),rho-fRtor);
1100   safe1 = fRmin - pt ;                           1104   safe1 = fRmin - pt ;
1101   safe2 = pt - fRmax ;                           1105   safe2 = pt - fRmax ;
1102                                                  1106 
1103   if (safe1 > safe2)  { safe = safe1; }          1107   if (safe1 > safe2)  { safe = safe1; }
1104   else                { safe = safe2; }          1108   else                { safe = safe2; }
1105                                                  1109 
1106   if ( fDPhi < twopi && (rho != 0.0) )        << 1110   if ( fDPhi < twopi && rho )
1107   {                                              1111   {
1108     phiC    = fSPhi + fDPhi*0.5 ;                1112     phiC    = fSPhi + fDPhi*0.5 ;
1109     cosPhiC = std::cos(phiC) ;                   1113     cosPhiC = std::cos(phiC) ;
1110     sinPhiC = std::sin(phiC) ;                   1114     sinPhiC = std::sin(phiC) ;
1111     cosPsi  = (p.x()*cosPhiC + p.y()*sinPhiC)    1115     cosPsi  = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1112                                                  1116 
1113     if (cosPsi < std::cos(fDPhi*0.5) ) // Psi    1117     if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point
1114     {                                  // Poi    1118     {                                  // Point lies outside phi range
1115       if ((p.y()*cosPhiC - p.x()*sinPhiC) <=     1119       if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1116       {                                          1120       {
1117         safePhi = std::fabs(p.x()*std::sin(fS    1121         safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1118       }                                          1122       }
1119       else                                       1123       else
1120       {                                          1124       {
1121         ePhi    = fSPhi + fDPhi ;                1125         ePhi    = fSPhi + fDPhi ;
1122         safePhi = std::fabs(p.x()*std::sin(eP    1126         safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1123       }                                          1127       }
1124       if (safePhi > safe)  { safe = safePhi ;    1128       if (safePhi > safe)  { safe = safePhi ; }
1125     }                                            1129     }
1126   }                                              1130   }
1127   if (safe < 0 )  { safe = 0 ; }                 1131   if (safe < 0 )  { safe = 0 ; }
1128   return safe;                                   1132   return safe;
1129 }                                                1133 }
1130                                                  1134 
1131 /////////////////////////////////////////////    1135 ///////////////////////////////////////////////////////////////////////////
1132 //                                               1136 //
1133 // Calculate distance to surface of shape fro    1137 // Calculate distance to surface of shape from `inside', allowing for tolerance
1134 // - Only Calc rmax intersection if no valid     1138 // - Only Calc rmax intersection if no valid rmin intersection
1135 //                                               1139 //
1136                                                  1140 
1137 G4double G4Torus::DistanceToOut( const G4Thre    1141 G4double G4Torus::DistanceToOut( const G4ThreeVector& p,
1138                                  const G4Thre    1142                                  const G4ThreeVector& v,
1139                                  const G4bool    1143                                  const G4bool calcNorm,
1140                                        G4bool << 1144                                        G4bool *validNorm,
1141                                        G4Thre << 1145                                        G4ThreeVector  *n  ) const
1142 {                                                1146 {
1143   ESide    side = kNull, sidephi = kNull ;       1147   ESide    side = kNull, sidephi = kNull ;
1144   G4double snxt = kInfinity, sphi, sd[4] ;       1148   G4double snxt = kInfinity, sphi, sd[4] ;
1145                                                  1149 
1146   // Vars for phi intersection                   1150   // Vars for phi intersection
1147   //                                             1151   //
1148   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, c    1152   G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1149   G4double cPhi, sinCPhi, cosCPhi ;              1153   G4double cPhi, sinCPhi, cosCPhi ;
1150   G4double pDistS, compS, pDistE, compE, sphi    1154   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1151                                                  1155 
1152   // Radial Intersections Defenitions & Gener    1156   // Radial Intersections Defenitions & General Precals
1153                                                  1157 
1154   //////////////////////// new calculation //    1158   //////////////////////// new calculation //////////////////////
1155                                                  1159 
1156 #if 1                                            1160 #if 1
1157                                                  1161 
1158   // This is the version with the calculation    1162   // This is the version with the calculation of CalcNorm = true 
1159   // To be done: Check the precision of this     1163   // To be done: Check the precision of this calculation.
1160   // If you want return always validNorm = fa    1164   // If you want return always validNorm = false, then take the version below
1161                                                  1165   
1162                                                  1166   
1163   G4double rho = std::hypot(p.x(),p.y());        1167   G4double rho = std::hypot(p.x(),p.y());
1164   G4double pt = hypot(p.z(),rho-fRtor);          1168   G4double pt = hypot(p.z(),rho-fRtor);
1165                                                  1169 
1166   G4double pDotV = p.x()*v.x() + p.y()*v.y()     1170   G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
1167                                                  1171 
1168   G4double tolRMax = fRmax - fRmaxTolerance ;    1172   G4double tolRMax = fRmax - fRmaxTolerance ;
1169                                                  1173    
1170   G4double vDotNmax   = pDotV - fRtor*(v.x()*    1174   G4double vDotNmax   = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
1171   G4double pDotxyNmax = (1 - fRtor/rho) ;        1175   G4double pDotxyNmax = (1 - fRtor/rho) ;
1172                                                  1176 
1173   if( (pt*pt > tolRMax*tolRMax) && (vDotNmax     1177   if( (pt*pt > tolRMax*tolRMax) && (vDotNmax >= 0) )
1174   {                                              1178   {
1175     // On tolerant boundary & heading outward    1179     // On tolerant boundary & heading outwards (or perpendicular to) outer
1176     // radial surface -> leaving immediately     1180     // radial surface -> leaving immediately with *n for really convex part
1177     // only                                      1181     // only
1178                                                  1182       
1179     if ( calcNorm && (pDotxyNmax >= -2.*fRmax    1183     if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) ) 
1180     {                                            1184     {
1181       *n = G4ThreeVector( p.x()*(1 - fRtor/rh    1185       *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
1182                           p.y()*(1 - fRtor/rh    1186                           p.y()*(1 - fRtor/rho)/pt,
1183                           p.z()/pt               1187                           p.z()/pt                  ) ;
1184       *validNorm = true ;                        1188       *validNorm = true ;
1185     }                                            1189     }
1186                                                  1190      
1187     return snxt = 0 ; // Leaving by Rmax imme    1191     return snxt = 0 ; // Leaving by Rmax immediately
1188   }                                              1192   }
1189                                                  1193   
1190   snxt = SolveNumericJT(p,v,fRmax,false);        1194   snxt = SolveNumericJT(p,v,fRmax,false);  
1191   side = kRMax ;                                 1195   side = kRMax ;
1192                                                  1196 
1193   // rmin                                        1197   // rmin
1194                                                  1198 
1195   if ( fRmin != 0.0 )                         << 1199   if ( fRmin )
1196   {                                              1200   {
1197     G4double tolRMin = fRmin + fRminTolerance    1201     G4double tolRMin = fRmin + fRminTolerance ;
1198                                                  1202 
1199     if ( (pt*pt < tolRMin*tolRMin) && (vDotNm    1203     if ( (pt*pt < tolRMin*tolRMin) && (vDotNmax < 0) )
1200     {                                            1204     {
1201       if (calcNorm)  { *validNorm = false ; }    1205       if (calcNorm)  { *validNorm = false ; } // Concave surface of the torus
1202       return  snxt = 0 ;                         1206       return  snxt = 0 ;                      // Leaving by Rmin immediately
1203     }                                            1207     }
1204                                                  1208     
1205     sd[0] = SolveNumericJT(p,v,fRmin,false);     1209     sd[0] = SolveNumericJT(p,v,fRmin,false);
1206     if ( sd[0] < snxt )                          1210     if ( sd[0] < snxt )
1207     {                                            1211     {
1208       snxt = sd[0] ;                             1212       snxt = sd[0] ;
1209       side = kRMin ;                             1213       side = kRMin ;
1210     }                                            1214     }
1211   }                                              1215   }
1212                                                  1216 
1213 #else                                            1217 #else
1214                                                  1218 
1215   // this is the "conservative" version which    1219   // this is the "conservative" version which return always validnorm = false
1216   // NOTE: using this version the unit test t    1220   // NOTE: using this version the unit test testG4Torus will break
1217                                                  1221 
1218   snxt = SolveNumericJT(p,v,fRmax,false);        1222   snxt = SolveNumericJT(p,v,fRmax,false);  
1219   side = kRMax ;                                 1223   side = kRMax ;
1220                                                  1224 
1221   if ( fRmin )                                   1225   if ( fRmin )
1222   {                                              1226   {
1223     sd[0] = SolveNumericJT(p,v,fRmin,false);     1227     sd[0] = SolveNumericJT(p,v,fRmin,false);
1224     if ( sd[0] < snxt )                          1228     if ( sd[0] < snxt )
1225     {                                            1229     {
1226       snxt = sd[0] ;                             1230       snxt = sd[0] ;
1227       side = kRMin ;                             1231       side = kRMin ;
1228     }                                            1232     }
1229   }                                              1233   }
1230                                                  1234 
1231   if ( calcNorm && (snxt == 0.0) )               1235   if ( calcNorm && (snxt == 0.0) )
1232   {                                              1236   {
1233     *validNorm = false ;    // Leaving solid,    1237     *validNorm = false ;    // Leaving solid, but possible re-intersection
1234     return snxt  ;                               1238     return snxt  ;
1235   }                                              1239   }
1236                                                  1240 
1237 #endif                                           1241 #endif
1238                                                  1242   
1239   if (fDPhi < twopi)  // Phi Intersections       1243   if (fDPhi < twopi)  // Phi Intersections
1240   {                                              1244   {
1241     sinSPhi = std::sin(fSPhi) ;                  1245     sinSPhi = std::sin(fSPhi) ;
1242     cosSPhi = std::cos(fSPhi) ;                  1246     cosSPhi = std::cos(fSPhi) ;
1243     ePhi    = fSPhi + fDPhi ;                    1247     ePhi    = fSPhi + fDPhi ;
1244     sinEPhi = std::sin(ePhi) ;                   1248     sinEPhi = std::sin(ePhi) ;
1245     cosEPhi = std::cos(ePhi) ;                   1249     cosEPhi = std::cos(ePhi) ;
1246     cPhi    = fSPhi + fDPhi*0.5 ;                1250     cPhi    = fSPhi + fDPhi*0.5 ;
1247     sinCPhi = std::sin(cPhi) ;                   1251     sinCPhi = std::sin(cPhi) ;
1248     cosCPhi = std::cos(cPhi) ;                   1252     cosCPhi = std::cos(cPhi) ;
1249                                                  1253     
1250     // angle calculation with correction         1254     // angle calculation with correction 
1251     // of difference in domain of atan2 and S    1255     // of difference in domain of atan2 and Sphi
1252     //                                           1256     //
1253     vphi = std::atan2(v.y(),v.x()) ;             1257     vphi = std::atan2(v.y(),v.x()) ;
1254                                                  1258      
1255     if ( vphi < fSPhi - halfAngTolerance  )      1259     if ( vphi < fSPhi - halfAngTolerance  )    { vphi += twopi; }
1256     else if ( vphi > ePhi + halfAngTolerance     1260     else if ( vphi > ePhi + halfAngTolerance ) { vphi -= twopi; }
1257                                                  1261 
1258     if ( (p.x() != 0.0) || (p.y() != 0.0) ) / << 1262     if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1259     {                                            1263     {
1260       pDistS = p.x()*sinSPhi - p.y()*cosSPhi     1264       pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
1261       pDistE = -p.x()*sinEPhi + p.y()*cosEPhi    1265       pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1262                                                  1266 
1263       // Comp -ve when in direction of outwar    1267       // Comp -ve when in direction of outwards normal
1264       //                                         1268       //
1265       compS   = -sinSPhi*v.x() + cosSPhi*v.y(    1269       compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
1266       compE   = sinEPhi*v.x() - cosEPhi*v.y()    1270       compE   = sinEPhi*v.x() - cosEPhi*v.y() ;
1267       sidephi = kNull ;                          1271       sidephi = kNull ;
1268                                                  1272      
1269       if( ( (fDPhi <= pi) && ( (pDistS <= hal    1273       if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1270                             && (pDistE <= hal    1274                             && (pDistE <= halfCarTolerance) ) )
1271        || ( (fDPhi >  pi) && ((pDistS <=  hal << 1275        || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
1272                             || (pDistE <=  ha << 1276                             && (pDistE >  halfCarTolerance) ) )  )
1273       {                                          1277       {
1274         // Inside both phi *full* planes         1278         // Inside both phi *full* planes
1275                                                  1279 
1276         if ( compS < 0 )                         1280         if ( compS < 0 )
1277         {                                        1281         {
1278           sphi = pDistS/compS ;                  1282           sphi = pDistS/compS ;
1279                                                  1283             
1280           if (sphi >= -halfCarTolerance)         1284           if (sphi >= -halfCarTolerance)
1281           {                                      1285           {
1282             xi = p.x() + sphi*v.x() ;            1286             xi = p.x() + sphi*v.x() ;
1283             yi = p.y() + sphi*v.y() ;            1287             yi = p.y() + sphi*v.y() ;
1284                                                  1288               
1285             // Check intersecting with correc    1289             // Check intersecting with correct half-plane
1286             // (if not -> no intersect)          1290             // (if not -> no intersect)
1287             //                                   1291             //
1288             if ( (std::fabs(xi)<=kCarToleranc    1292             if ( (std::fabs(xi)<=kCarTolerance)
1289               && (std::fabs(yi)<=kCarToleranc    1293               && (std::fabs(yi)<=kCarTolerance) )
1290             {                                    1294             {
1291               sidephi = kSPhi;                   1295               sidephi = kSPhi;
1292               if ( ((fSPhi-halfAngTolerance)<    1296               if ( ((fSPhi-halfAngTolerance)<=vphi)
1293                 && ((ePhi+halfAngTolerance)>=    1297                 && ((ePhi+halfAngTolerance)>=vphi) )
1294               {                                  1298               {
1295                 sphi = kInfinity;                1299                 sphi = kInfinity;
1296               }                                  1300               }
1297             }                                    1301             }
1298             else if ( yi*cosCPhi-xi*sinCPhi >    1302             else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1299             {                                    1303             {
1300               sphi = kInfinity ;                 1304               sphi = kInfinity ;
1301             }                                    1305             }
1302             else                                 1306             else
1303             {                                    1307             {
1304               sidephi = kSPhi ;                  1308               sidephi = kSPhi ;
1305             }                                    1309             }       
1306           }                                      1310           }
1307           else                                   1311           else
1308           {                                      1312           {
1309             sphi = kInfinity ;                   1313             sphi = kInfinity ;
1310           }                                      1314           }
1311         }                                        1315         }
1312         else                                     1316         else
1313         {                                        1317         {
1314           sphi = kInfinity ;                     1318           sphi = kInfinity ;
1315         }                                        1319         }
1316                                                  1320 
1317         if ( compE < 0 )                         1321         if ( compE < 0 )
1318         {                                        1322         {
1319           sphi2 = pDistE/compE ;                 1323           sphi2 = pDistE/compE ;
1320                                                  1324             
1321           // Only check further if < starting    1325           // Only check further if < starting phi intersection
1322           //                                     1326           //
1323           if ( (sphi2 > -kCarTolerance) && (s    1327           if ( (sphi2 > -kCarTolerance) && (sphi2 < sphi) )
1324           {                                      1328           {
1325             xi = p.x() + sphi2*v.x() ;           1329             xi = p.x() + sphi2*v.x() ;
1326             yi = p.y() + sphi2*v.y() ;           1330             yi = p.y() + sphi2*v.y() ;
1327                                                  1331               
1328             if ( (std::fabs(xi)<=kCarToleranc    1332             if ( (std::fabs(xi)<=kCarTolerance)
1329               && (std::fabs(yi)<=kCarToleranc    1333               && (std::fabs(yi)<=kCarTolerance) )
1330             {                                    1334             {
1331               // Leaving via ending phi          1335               // Leaving via ending phi
1332               //                                 1336               //
1333               if( (fSPhi-halfAngTolerance > v << 1337               if( !( (fSPhi-halfAngTolerance <= vphi)
1334                   || (ePhi+halfAngTolerance < << 1338                   && (ePhi+halfAngTolerance >= vphi) ) )
1335               {                                  1339               {
1336                 sidephi = kEPhi ;                1340                 sidephi = kEPhi ;
1337                 sphi = sphi2;                    1341                 sphi = sphi2;
1338               }                                  1342               }
1339             }                                    1343             } 
1340             else    // Check intersecting wit    1344             else    // Check intersecting with correct half-plane 
1341             {                                    1345             {
1342               if ( (yi*cosCPhi-xi*sinCPhi) >=    1346               if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1343               {                                  1347               {
1344                 // Leaving via ending phi        1348                 // Leaving via ending phi
1345                 //                               1349                 //
1346                 sidephi = kEPhi ;                1350                 sidephi = kEPhi ;
1347                 sphi = sphi2;                    1351                 sphi = sphi2;
1348                                                  1352                
1349               }                                  1353               }
1350             }                                    1354             }
1351           }                                      1355           }
1352         }                                        1356         }
1353       }                                          1357       }
1354       else                                       1358       else
1355       {                                          1359       {
1356         sphi = kInfinity ;                       1360         sphi = kInfinity ;
1357       }                                          1361       }
1358     }                                            1362     } 
1359     else                                         1363     else
1360     {                                            1364     {
1361       // On z axis + travel not || to z axis     1365       // On z axis + travel not || to z axis -> if phi of vector direction
1362       // within phi of shape, Step limited by    1366       // within phi of shape, Step limited by rmax, else Step =0
1363                                                  1367 
1364       vphi = std::atan2(v.y(),v.x());            1368       vphi = std::atan2(v.y(),v.x());
1365                                                  1369  
1366       if ( ( fSPhi-halfAngTolerance <= vphi )    1370       if ( ( fSPhi-halfAngTolerance <= vphi ) && 
1367            ( vphi <= ( ePhi+halfAngTolerance     1371            ( vphi <= ( ePhi+halfAngTolerance ) ) )
1368       {                                          1372       {
1369         sphi = kInfinity;                        1373         sphi = kInfinity;
1370       }                                          1374       }
1371       else                                       1375       else
1372       {                                          1376       {
1373         sidephi = kSPhi ; // arbitrary           1377         sidephi = kSPhi ; // arbitrary 
1374         sphi=0;                                  1378         sphi=0;
1375       }                                          1379       }
1376     }                                            1380     }
1377                                                  1381 
1378     // Order intersections                       1382     // Order intersections
1379                                                  1383 
1380     if (sphi<snxt)                               1384     if (sphi<snxt)
1381     {                                            1385     {
1382       snxt=sphi;                                 1386       snxt=sphi;
1383       side=sidephi;                              1387       side=sidephi;
1384     }                                            1388     }     
1385   }                                              1389   }
1386                                                  1390 
1387   G4double rhoi,it,iDotxyNmax ;                  1391   G4double rhoi,it,iDotxyNmax ;
1388   // Note: by numerical computation we know w    1392   // Note: by numerical computation we know where the ray hits the torus
1389   // So I propose to return the side where th    1393   // So I propose to return the side where the ray hits
1390                                                  1394 
1391   if (calcNorm)                                  1395   if (calcNorm)
1392   {                                              1396   {
1393     switch(side)                                 1397     switch(side)
1394     {                                            1398     {
1395       case kRMax:                     // n is    1399       case kRMax:                     // n is unit vector 
1396         xi    = p.x() + snxt*v.x() ;             1400         xi    = p.x() + snxt*v.x() ;
1397         yi    = p.y() + snxt*v.y() ;             1401         yi    = p.y() + snxt*v.y() ;
1398         zi    = p.z() + snxt*v.z() ;             1402         zi    = p.z() + snxt*v.z() ;
1399         rhoi = std::hypot(xi,yi);                1403         rhoi = std::hypot(xi,yi);
1400         it = hypot(zi,rhoi-fRtor);               1404         it = hypot(zi,rhoi-fRtor);
1401                                                  1405 
1402         iDotxyNmax = (1-fRtor/rhoi) ;            1406         iDotxyNmax = (1-fRtor/rhoi) ;
1403         if(iDotxyNmax >= -2.*fRmaxTolerance)     1407         if(iDotxyNmax >= -2.*fRmaxTolerance) // really convex part of Rmax
1404         {                                        1408         {                       
1405           *n = G4ThreeVector( xi*(1-fRtor/rho    1409           *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
1406                               yi*(1-fRtor/rho    1410                               yi*(1-fRtor/rhoi)/it,
1407                               zi/it              1411                               zi/it                 ) ;
1408           *validNorm = true ;                    1412           *validNorm = true ;
1409         }                                        1413         }
1410         else                                     1414         else
1411         {                                        1415         {
1412           *validNorm = false ; // concave-con    1416           *validNorm = false ; // concave-convex part of Rmax
1413         }                                        1417         }
1414         break ;                                  1418         break ;
1415                                                  1419 
1416       case kRMin:                                1420       case kRMin:
1417         *validNorm = false ;  // Rmin is conc    1421         *validNorm = false ;  // Rmin is concave or concave-convex
1418         break;                                   1422         break;
1419                                                  1423 
1420       case kSPhi:                                1424       case kSPhi:
1421         if (fDPhi <= pi )                        1425         if (fDPhi <= pi )
1422         {                                        1426         {
1423           *n=G4ThreeVector(std::sin(fSPhi),-s    1427           *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
1424           *validNorm=true;                       1428           *validNorm=true;
1425         }                                        1429         }
1426         else                                     1430         else
1427         {                                        1431         {
1428           *validNorm = false ;                   1432           *validNorm = false ;
1429         }                                        1433         }
1430         break ;                                  1434         break ;
1431                                                  1435 
1432       case kEPhi:                                1436       case kEPhi:
1433         if (fDPhi <= pi)                         1437         if (fDPhi <= pi)
1434         {                                        1438         {
1435           *n=G4ThreeVector(-std::sin(fSPhi+fD    1439           *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
1436           *validNorm=true;                       1440           *validNorm=true;
1437         }                                        1441         }
1438         else                                     1442         else
1439         {                                        1443         {
1440           *validNorm = false ;                   1444           *validNorm = false ;
1441         }                                        1445         }
1442         break;                                   1446         break;
1443                                                  1447 
1444       default:                                   1448       default:
1445                                                  1449 
1446         // It seems we go here from time to t    1450         // It seems we go here from time to time ...
1447                                                  1451 
1448         G4cout << G4endl;                        1452         G4cout << G4endl;
1449         DumpInfo();                              1453         DumpInfo();
1450         std::ostringstream message;              1454         std::ostringstream message;
1451         G4long oldprc = message.precision(16) << 1455         G4int oldprc = message.precision(16);
1452         message << "Undefined side for valid     1456         message << "Undefined side for valid surface normal to solid."
1453                 << G4endl                        1457                 << G4endl
1454                 << "Position:"  << G4endl <<     1458                 << "Position:"  << G4endl << G4endl
1455                 << "p.x() = "   << p.x()/mm <    1459                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
1456                 << "p.y() = "   << p.y()/mm <    1460                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
1457                 << "p.z() = "   << p.z()/mm <    1461                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
1458                 << "Direction:" << G4endl <<     1462                 << "Direction:" << G4endl << G4endl
1459                 << "v.x() = "   << v.x() << G    1463                 << "v.x() = "   << v.x() << G4endl
1460                 << "v.y() = "   << v.y() << G    1464                 << "v.y() = "   << v.y() << G4endl
1461                 << "v.z() = "   << v.z() << G    1465                 << "v.z() = "   << v.z() << G4endl << G4endl
1462                 << "Proposed distance :" << G    1466                 << "Proposed distance :" << G4endl << G4endl
1463                 << "snxt = " << snxt/mm << "     1467                 << "snxt = " << snxt/mm << " mm" << G4endl;
1464         message.precision(oldprc);               1468         message.precision(oldprc);
1465         G4Exception("G4Torus::DistanceToOut(p    1469         G4Exception("G4Torus::DistanceToOut(p,v,..)",
1466                     "GeomSolids1002",JustWarn    1470                     "GeomSolids1002",JustWarning, message);
1467         break;                                   1471         break;
1468     }                                            1472     }
1469   }                                              1473   }
1470   if ( snxt<halfCarTolerance )  { snxt=0 ; }     1474   if ( snxt<halfCarTolerance )  { snxt=0 ; }
1471                                                  1475 
1472   return snxt;                                   1476   return snxt;
1473 }                                                1477 }
1474                                                  1478 
1475 /////////////////////////////////////////////    1479 /////////////////////////////////////////////////////////////////////////
1476 //                                               1480 //
1477 // Calculate distance (<=actual) to closest s    1481 // Calculate distance (<=actual) to closest surface of shape from inside
1478                                                  1482 
1479 G4double G4Torus::DistanceToOut( const G4Thre    1483 G4double G4Torus::DistanceToOut( const G4ThreeVector& p ) const
1480 {                                                1484 {
1481   G4double safe=0.0,safeR1,safeR2;               1485   G4double safe=0.0,safeR1,safeR2;
1482   G4double rho,pt ;                              1486   G4double rho,pt ;
1483   G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;    1487   G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1484                                                  1488   
1485   rho = std::hypot(p.x(),p.y());                 1489   rho = std::hypot(p.x(),p.y());
1486   pt  = std::hypot(p.z(),rho-fRtor);             1490   pt  = std::hypot(p.z(),rho-fRtor);
1487                                                  1491   
1488 #ifdef G4CSGDEBUG                                1492 #ifdef G4CSGDEBUG
1489   if( Inside(p) == kOutside )                    1493   if( Inside(p) == kOutside )
1490   {                                              1494   {
1491      G4long oldprc = G4cout.precision(16) ;   << 1495      G4int oldprc = G4cout.precision(16) ;
1492      G4cout << G4endl ;                          1496      G4cout << G4endl ;
1493      DumpInfo();                                 1497      DumpInfo();
1494      G4cout << "Position:"  << G4endl << G4en    1498      G4cout << "Position:"  << G4endl << G4endl ;
1495      G4cout << "p.x() = "   << p.x()/mm << "     1499      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1496      G4cout << "p.y() = "   << p.y()/mm << "     1500      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1497      G4cout << "p.z() = "   << p.z()/mm << "     1501      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1498      G4cout.precision(oldprc);                   1502      G4cout.precision(oldprc);
1499      G4Exception("G4Torus::DistanceToOut(p)",    1503      G4Exception("G4Torus::DistanceToOut(p)", "GeomSolids1002",
1500                  JustWarning, "Point p is out    1504                  JustWarning, "Point p is outside !?" );
1501   }                                              1505   }
1502 #endif                                           1506 #endif
1503                                                  1507 
1504   if (fRmin != 0.0)                           << 1508   if (fRmin)
1505   {                                              1509   {
1506     safeR1 = pt - fRmin ;                        1510     safeR1 = pt - fRmin ;
1507     safeR2 = fRmax - pt ;                        1511     safeR2 = fRmax - pt ;
1508                                                  1512 
1509     if (safeR1 < safeR2)  { safe = safeR1 ; }    1513     if (safeR1 < safeR2)  { safe = safeR1 ; }
1510     else                  { safe = safeR2 ; }    1514     else                  { safe = safeR2 ; }
1511   }                                              1515   }
1512   else                                           1516   else
1513   {                                              1517   {
1514     safe = fRmax - pt ;                          1518     safe = fRmax - pt ;
1515   }                                              1519   }  
1516                                                  1520 
1517   // Check if phi divided, Calc distances clo    1521   // Check if phi divided, Calc distances closest phi plane
1518   //                                             1522   //
1519   if (fDPhi < twopi) // Above/below central p << 1523   if (fDPhi<twopi) // Above/below central phi of Torus?
1520   {                                              1524   {
1521     phiC    = fSPhi + fDPhi*0.5 ;                1525     phiC    = fSPhi + fDPhi*0.5 ;
1522     cosPhiC = std::cos(phiC) ;                   1526     cosPhiC = std::cos(phiC) ;
1523     sinPhiC = std::sin(phiC) ;                   1527     sinPhiC = std::sin(phiC) ;
1524                                                  1528 
1525     if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)        1529     if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1526     {                                            1530     {
1527       safePhi = -(p.x()*std::sin(fSPhi) - p.y    1531       safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1528     }                                            1532     }
1529     else                                         1533     else
1530     {                                            1534     {
1531       ePhi    = fSPhi + fDPhi ;                  1535       ePhi    = fSPhi + fDPhi ;
1532       safePhi = (p.x()*std::sin(ePhi) - p.y()    1536       safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1533     }                                            1537     }
1534     if (safePhi < safe)  { safe = safePhi ; }    1538     if (safePhi < safe)  { safe = safePhi ; }
1535   }                                              1539   }
1536   if (safe < 0)  { safe = 0 ; }                  1540   if (safe < 0)  { safe = 0 ; }
1537   return safe ;                                  1541   return safe ;  
1538 }                                                1542 }
1539                                                  1543 
1540 /////////////////////////////////////////////    1544 //////////////////////////////////////////////////////////////////////////
1541 //                                               1545 //
1542 // Stream object contents to an output stream    1546 // Stream object contents to an output stream
1543                                                  1547 
1544 G4GeometryType G4Torus::GetEntityType() const    1548 G4GeometryType G4Torus::GetEntityType() const
1545 {                                                1549 {
1546   return {"G4Torus"};                         << 1550   return G4String("G4Torus");
1547 }                                                1551 }
1548                                                  1552 
1549 /////////////////////////////////////////////    1553 //////////////////////////////////////////////////////////////////////////
1550 //                                               1554 //
1551 // Make a clone of the object                    1555 // Make a clone of the object
1552 //                                               1556 //
1553 G4VSolid* G4Torus::Clone() const                 1557 G4VSolid* G4Torus::Clone() const
1554 {                                                1558 {
1555   return new G4Torus(*this);                     1559   return new G4Torus(*this);
1556 }                                                1560 }
1557                                                  1561 
1558 /////////////////////////////////////////////    1562 //////////////////////////////////////////////////////////////////////////
1559 //                                               1563 //
1560 // Stream object contents to an output stream    1564 // Stream object contents to an output stream
1561                                                  1565 
1562 std::ostream& G4Torus::StreamInfo( std::ostre    1566 std::ostream& G4Torus::StreamInfo( std::ostream& os ) const
1563 {                                                1567 {
1564   G4long oldprc = os.precision(16);           << 1568   G4int oldprc = os.precision(16);
1565   os << "------------------------------------    1569   os << "-----------------------------------------------------------\n"
1566      << "    *** Dump for solid - " << GetNam    1570      << "    *** Dump for solid - " << GetName() << " ***\n"
1567      << "    ================================    1571      << "    ===================================================\n"
1568      << " Solid type: G4Torus\n"                 1572      << " Solid type: G4Torus\n"
1569      << " Parameters: \n"                        1573      << " Parameters: \n"
1570      << "    inner radius: " << fRmin/mm << "    1574      << "    inner radius: " << fRmin/mm << " mm \n"
1571      << "    outer radius: " << fRmax/mm << "    1575      << "    outer radius: " << fRmax/mm << " mm \n"
1572      << "    swept radius: " << fRtor/mm << "    1576      << "    swept radius: " << fRtor/mm << " mm \n"
1573      << "    starting phi: " << fSPhi/degree     1577      << "    starting phi: " << fSPhi/degree << " degrees \n"
1574      << "    delta phi   : " << fDPhi/degree     1578      << "    delta phi   : " << fDPhi/degree << " degrees \n"
1575      << "------------------------------------    1579      << "-----------------------------------------------------------\n";
1576   os.precision(oldprc);                          1580   os.precision(oldprc);
1577                                                  1581 
1578   return os;                                     1582   return os;
1579 }                                                1583 }
1580                                                  1584 
1581 /////////////////////////////////////////////    1585 ////////////////////////////////////////////////////////////////////////////
1582 //                                               1586 //
1583 // GetPointOnSurface                             1587 // GetPointOnSurface
1584                                                  1588 
1585 G4ThreeVector G4Torus::GetPointOnSurface() co    1589 G4ThreeVector G4Torus::GetPointOnSurface() const
1586 {                                                1590 {
1587   G4double cosu, sinu,cosv, sinv, aOut, aIn,     1591   G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
1588                                                  1592    
1589   phi   = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi    1593   phi   = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi);
1590   theta = G4RandFlat::shoot(0.,twopi);           1594   theta = G4RandFlat::shoot(0.,twopi);
1591                                                  1595   
1592   cosu   = std::cos(phi);    sinu = std::sin(    1596   cosu   = std::cos(phi);    sinu = std::sin(phi);
1593   cosv   = std::cos(theta);  sinv = std::sin(    1597   cosv   = std::cos(theta);  sinv = std::sin(theta); 
1594                                                  1598 
1595   // compute the areas                           1599   // compute the areas
1596                                                  1600 
1597   aOut   = (fDPhi)*twopi*fRtor*fRmax;            1601   aOut   = (fDPhi)*twopi*fRtor*fRmax;
1598   aIn    = (fDPhi)*twopi*fRtor*fRmin;            1602   aIn    = (fDPhi)*twopi*fRtor*fRmin;
1599   aSide  = pi*(fRmax*fRmax-fRmin*fRmin);         1603   aSide  = pi*(fRmax*fRmax-fRmin*fRmin);
1600                                                  1604   
1601   if ((fSPhi == 0) && (fDPhi == twopi)){ aSid    1605   if ((fSPhi == 0) && (fDPhi == twopi)){ aSide = 0; }
1602   chose = G4RandFlat::shoot(0.,aOut + aIn + 2    1606   chose = G4RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
1603                                                  1607 
1604   if(chose < aOut)                               1608   if(chose < aOut)
1605   {                                              1609   {
1606     return { (fRtor+fRmax*cosv)*cosu, (fRtor+ << 1610     return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
                                                   >> 1611                           (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
1607   }                                              1612   }
1608   else if( (chose >= aOut) && (chose < aOut +    1613   else if( (chose >= aOut) && (chose < aOut + aIn) )
1609   {                                              1614   {
1610     return { (fRtor+fRmin*cosv)*cosu, (fRtor+ << 1615     return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
                                                   >> 1616                           (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
1611   }                                              1617   }
1612   else if( (chose >= aOut + aIn) && (chose <     1618   else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
1613   {                                              1619   {
1614     rRand = GetRadiusInRing(fRmin,fRmax);        1620     rRand = GetRadiusInRing(fRmin,fRmax);
1615     return { (fRtor+rRand*cosv)*std::cos(fSPh << 1621     return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
1616              (fRtor+rRand*cosv)*std::sin(fSPh << 1622                           (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
1617   }                                              1623   }
1618   else                                           1624   else
1619   {                                              1625   {   
1620     rRand = GetRadiusInRing(fRmin,fRmax);        1626     rRand = GetRadiusInRing(fRmin,fRmax);
1621     return { (fRtor+rRand*cosv)*std::cos(fSPh << 1627     return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
1622              (fRtor+rRand*cosv)*std::sin(fSPh << 1628                           (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi), 
                                                   >> 1629                           rRand*sinv);
1623    }                                             1630    }
1624 }                                                1631 }
1625                                                  1632 
1626 /////////////////////////////////////////////    1633 ///////////////////////////////////////////////////////////////////////
1627 //                                               1634 //
1628 // Visualisation Functions                       1635 // Visualisation Functions
1629                                                  1636 
1630 void G4Torus::DescribeYourselfTo ( G4VGraphic    1637 void G4Torus::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
1631 {                                                1638 {
1632   scene.AddSolid (*this);                        1639   scene.AddSolid (*this);
1633 }                                                1640 }
1634                                                  1641 
1635 G4Polyhedron* G4Torus::CreatePolyhedron () co    1642 G4Polyhedron* G4Torus::CreatePolyhedron () const 
1636 {                                                1643 {
1637   return new G4PolyhedronTorus (fRmin, fRmax,    1644   return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
1638 }                                                1645 }
1639                                                  1646 
1640 #endif // !defined(G4GEOM_USE_TORUS) || !defi    1647 #endif // !defined(G4GEOM_USE_TORUS) || !defined(G4GEOM_USE_SYS_USOLIDS)
1641                                                  1648