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