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