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


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