| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4Torus implementation << 23 //
>> 24 // $Id: G4Torus.cc,v 1.28 2002/10/28 15:18:17 gcosmo Exp $
>> 25 // GEANT4 tag $Name: geant4-05-00 $
>> 26 //
>> 27 //
>> 28 // class G4Torus
>> 29 //
>> 30 // Implementation
27 // 31 //
28 // 30.10.96 V.Grichine: first implementation w 32 // 30.10.96 V.Grichine: first implementation with G4Tubs elements in Fs
29 // 26.05.00 V.Grichine: added new fuctions dev << 33 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
30 // 31.08.00 E.Medernach: numerical computation << 34 // 19.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
>> 35 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
>> 36 // 26.05.00 V.Grichine: new fuctions developed by O.Cremonesi were added
>> 37 // 31.08.00 E.Medernach: numerical computation of roots wuth bounding
>> 38 // volume technique
>> 39 // 03.10.00 E.Medernach: SafeNewton added
31 // 11.01.01 E.Medernach: Use G4PolynomialSolve 40 // 11.01.01 E.Medernach: Use G4PolynomialSolver to find roots
32 // 03.05.05 V.Grichine: SurfaceNormal(p) accor << 41 //
33 // 25.08.05 O.Link: new methods for DistanceTo << 42 // ********************************************************************
34 // 28.10.16 E.Tcherniaev: new CalculateExtent( <<
35 // 16.12.16 H.Burkhardt: use radius difference <<
36 // ------------------------------------------- <<
37 43
38 #include "G4Torus.hh" 44 #include "G4Torus.hh"
39 45
40 #if !(defined(G4GEOM_USE_UTORUS) && defined(G4 <<
41 <<
42 #include "G4GeomTools.hh" <<
43 #include "G4VoxelLimits.hh" 46 #include "G4VoxelLimits.hh"
44 #include "G4AffineTransform.hh" 47 #include "G4AffineTransform.hh"
45 #include "G4BoundingEnvelope.hh" <<
46 #include "G4GeometryTolerance.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" <<
54 <<
55 #include "G4VGraphicsScene.hh" 53 #include "G4VGraphicsScene.hh"
56 #include "G4Polyhedron.hh" 54 #include "G4Polyhedron.hh"
>> 55 #include "G4NURBS.hh"
>> 56 #include "G4NURBStube.hh"
>> 57 #include "G4NURBScylinder.hh"
>> 58 #include "G4NURBStubesector.hh"
>> 59 #include "G4PolynomialSolver.hh"
57 60
58 using namespace CLHEP; << 61 // #define DEBUGTORUS 1
59 62
60 ////////////////////////////////////////////// 63 ///////////////////////////////////////////////////////////////
61 // 64 //
62 // Constructor - check parameters, convert ang 65 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
63 // - note if pdphi>2PI then reset 66 // - note if pdphi>2PI then reset to 2PI
64 67
65 G4Torus::G4Torus( const G4String& pName, << 68 G4Torus::G4Torus( const G4String &pName,
66 G4double pRmin, 69 G4double pRmin,
67 G4double pRmax, 70 G4double pRmax,
68 G4double pRtor, 71 G4double pRtor,
69 G4double pSPhi, 72 G4double pSPhi,
70 G4double pDPhi ) << 73 G4double pDPhi)
71 : G4CSGSolid(pName) 74 : G4CSGSolid(pName)
72 { 75 {
73 SetAllParameters(pRmin, pRmax, pRtor, pSPhi, 76 SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
74 } 77 }
75 78
76 ////////////////////////////////////////////// <<
77 // <<
78 // <<
79 <<
80 void 79 void
81 G4Torus::SetAllParameters( G4double pRmin, 80 G4Torus::SetAllParameters( G4double pRmin,
82 G4double pRmax, 81 G4double pRmax,
83 G4double pRtor, 82 G4double pRtor,
84 G4double pSPhi, 83 G4double pSPhi,
85 G4double pDPhi ) 84 G4double pDPhi )
86 { 85 {
87 const G4double fEpsilon = 4.e-11; // relati << 86 if ( pRtor >= pRmax + kCarTolerance ) // Check swept radius
88 <<
89 fCubicVolume = 0.; <<
90 fSurfaceArea = 0.; <<
91 fRebuildPolyhedron = true; <<
92 <<
93 kRadTolerance = G4GeometryTolerance::GetInst <<
94 kAngTolerance = G4GeometryTolerance::GetInst <<
95 <<
96 halfCarTolerance = 0.5*kCarTolerance; <<
97 halfAngTolerance = 0.5*kAngTolerance; <<
98 <<
99 if ( pRtor >= pRmax+1.e3*kCarTolerance ) // <<
100 { 87 {
101 fRtor = pRtor ; 88 fRtor = pRtor ;
102 } 89 }
103 else 90 else
104 { 91 {
105 std::ostringstream message; << 92 G4cout << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl
106 message << "Invalid swept radius for Solid << 93 << " Invalid swept radius !" << G4endl
107 << " pRtor = " << pRtor << << 94 << "pRtor = " << pRtor << ", pRmax = " << pRmax << G4endl;
108 G4Exception("G4Torus::SetAllParameters()", << 95 G4cerr << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl
109 "GeomSolids0002", FatalExcepti << 96 << " Invalid swept radius !" << G4endl
>> 97 << "pRtor = " << pRtor << ", pRmax = " << pRmax << G4endl;
>> 98 G4Exception("G4Torus::SetAllParameters() - invalid swept radius");
110 } 99 }
111 100
112 // Check radii, as in G4Cons << 101 // Check radii
113 // << 102
114 if ( pRmin < pRmax - 1.e2*kCarTolerance && p << 103 if ( pRmin < pRmax - 2*kCarTolerance && pRmin >= 0 )
115 { 104 {
116 if (pRmin >= 1.e2*kCarTolerance) { fRmin = << 105 if (pRmin >= kCarTolerance) fRmin = pRmin ;
117 else { fRmin = << 106 else fRmin = 0.0 ;
118 fRmax = pRmax ; 107 fRmax = pRmax ;
119 } 108 }
120 else 109 else
121 { 110 {
122 std::ostringstream message; << 111 G4cout << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl
123 message << "Invalid values of radii for So << 112 << " Invalid values for radii !" << G4endl
124 << " pRmin = " << pRmin << << 113 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
125 G4Exception("G4Torus::SetAllParameters()", << 114 G4cerr << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl
126 "GeomSolids0002", FatalExcepti << 115 << " Invalid values for radii !" << G4endl
>> 116 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
>> 117 G4Exception("G4Torus::SetAllParameters() - invalid radii");
127 } 118 }
128 119
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 120 // Check angles
136 // << 121
137 if ( pDPhi >= twopi ) { fDPhi = twopi ; } << 122 if ( pDPhi >= 2.0*M_PI ) fDPhi = 2*M_PI ;
138 else 123 else
139 { 124 {
140 if (pDPhi > 0) { fDPhi = pDPhi ; } << 125 if (pDPhi > 0) fDPhi = pDPhi ;
141 else 126 else
142 { 127 {
143 std::ostringstream message; << 128 G4cout << "ERROR - G4Torus::SetAllParameters(): " << GetName() << G4endl
144 message << "Invalid Z delta-Phi for Soli << 129 << " Negative delta-Phi ! - "
145 << " pDPhi = " << pDPhi; << 130 << pDPhi << G4endl;
146 G4Exception("G4Torus::SetAllParameters() << 131 G4cerr << "ERROR - G4Torus::SetAllParameters(): " << GetName() << G4endl
147 "GeomSolids0002", FatalExcep << 132 << " Negative Z delta-Phi ! - "
>> 133 << pDPhi << G4endl;
>> 134 G4Exception("G4Torus::SetAllParameters() - invalid dphi");
148 } 135 }
149 } 136 }
150 137
151 // Ensure psphi in 0-2PI or -2PI-0 range if 138 // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
152 // <<
153 fSPhi = pSPhi; <<
154 139
155 if (fSPhi < 0) { fSPhi = twopi-std::fmod(st << 140 fSPhi = pSPhi;
156 else { fSPhi = std::fmod(fSPhi,tw <<
157 141
158 if (fSPhi+fDPhi > twopi) { fSPhi-=twopi ; } << 142 if (fSPhi < 0) fSPhi = 2.0*M_PI-fmod(fabs(fSPhi),2.0*M_PI) ;
159 } << 143 else fSPhi = fmod(fSPhi,2.0*M_PI) ;
160 144
161 ////////////////////////////////////////////// << 145 if (fSPhi+fDPhi > 2.0*M_PI) fSPhi-=2.0*M_PI ;
162 // <<
163 // Fake default constructor - sets only member <<
164 // for usage restri <<
165 // <<
166 G4Torus::G4Torus( __void__& a ) <<
167 : G4CSGSolid(a) <<
168 { <<
169 } 146 }
170 147
171 ////////////////////////////////////////////// 148 //////////////////////////////////////////////////////////////////////
172 // 149 //
173 // Destructor 150 // Destructor
174 151
175 G4Torus::~G4Torus() = default; << 152 G4Torus::~G4Torus()
176 << 153 {}
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 154
209 ////////////////////////////////////////////// 155 //////////////////////////////////////////////////////////////////////
210 // 156 //
211 // Dispatch to parameterisation for replicatio 157 // Dispatch to parameterisation for replication mechanism dimension
212 // computation & modification. 158 // computation & modification.
213 159
214 void G4Torus::ComputeDimensions( G4VPVPa 160 void G4Torus::ComputeDimensions( G4VPVParameterisation* p,
215 const G4int n 161 const G4int n,
216 const G4VPhys 162 const G4VPhysicalVolume* pRep )
217 { 163 {
218 p->ComputeDimensions(*this,n,pRep); 164 p->ComputeDimensions(*this,n,pRep);
219 } 165 }
220 166
>> 167 ///////////////////////////////////////////////////////////////////////////
>> 168 //
>> 169 // Test function for study of intersections of a ray (starting from p along
>> 170 // v) with the torus
221 171
>> 172 G4int G4Torus::TorusRoots( G4double Ri,
>> 173 const G4ThreeVector& p,
>> 174 const G4ThreeVector& v ) const
>> 175 {
>> 176 // Define roots Si (generally real >=0) for intersection with
>> 177 // torus (Ri = fRmax or fRmin) of ray p +S*v . General equation is :
>> 178 // c[4]*S^4 + c[3]*S^3 +c[2]*S^2 + c[1]*S + c[0] = 0 .
>> 179
>> 180 G4double c[5],s[4] ;
>> 181 G4int num, i, j ;
>> 182 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
>> 183 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
>> 184 G4double Rtor2 = fRtor*fRtor, Ri2 = Ri*Ri ;
>> 185
>> 186 c[4] = 1.0 ;
>> 187 c[3] = 4*pDotV ;
>> 188 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - Ri2 + 2*Rtor2*v.z()*v.z()) ;
>> 189 c[1] = 4*(pDotV*(pRad2-Rtor2-Ri2) + 2*Rtor2*p.z()*v.z()) ;
>> 190 c[0] = pRad2*pRad2 - 2*pRad2*(Rtor2+Ri2)
>> 191 + 4*Rtor2*p.z()*p.z() + (Rtor2-Ri2)*(Rtor2-Ri2) ;
>> 192
>> 193 num = SolveBiQuadratic(c,s) ;
>> 194
>> 195 if(num)
>> 196 {
>> 197 for(i=0;i<num;i++) // leave only >=0 roots
>> 198 {
>> 199 if(s[i]<0)
>> 200 {
>> 201 for(j=i+1;j<num;j++) s[j-1] = s[j] ;
>> 202 i-- ;
>> 203 num-- ;
>> 204 }
>> 205 }
>> 206 if(num)
>> 207 {
>> 208 for(i=0;i<num;i++)
>> 209 {
>> 210 G4cout << i << " Root = " << s[i] << G4endl ;
>> 211 }
>> 212 }
>> 213 else G4cout << "All real roots are negative" << G4endl ;
>> 214 }
>> 215 else G4cout << "No real roots for intesection with torus" << G4endl;
222 216
223 ////////////////////////////////////////////// << 217 num = SolveBiQuadraticNew(c,s) ;
>> 218
>> 219 if(num)
>> 220 {
>> 221 for(i=0;i<num;i++) // leave only >=0 roots
>> 222 {
>> 223 if(s[i]<0)
>> 224 {
>> 225 for(j=i+1;j<num;j++) s[j-1] = s[j] ;
>> 226 i-- ;
>> 227 num-- ;
>> 228 }
>> 229 }
>> 230 if(num)
>> 231 {
>> 232 for(i=0;i<num;i++)
>> 233 {
>> 234 G4cout << i << " new Root = " << s[i] << G4endl ;
>> 235 }
>> 236 }
>> 237 else G4cout << "All real new roots are negative" << G4endl ;
>> 238 }
>> 239 else G4cout << "No real new roots for intesection with torus" << G4endl;
>> 240
>> 241 return num ;
>> 242 }
>> 243
>> 244 /////////////////////////////////////////////////////////////////////////
224 // 245 //
225 // Calculate the real roots to torus surface. << 246 // Auxiliary method for solving (in real numbers) biquadratic equation
226 // Returns negative solutions as well. << 247 // Algorithm based on : Graphics Gems I by Jochen Schwartz
227 248
228 void G4Torus::TorusRootsJT( const G4ThreeVecto << 249 G4int G4Torus::SolveBiQuadratic( G4double c[], G4double s[] ) const
229 const G4ThreeVecto <<
230 G4double r, <<
231 std::vector< <<
232 { 250 {
>> 251 G4double coeffs[ 4 ];
>> 252 G4double z, u, v, sub;
>> 253 G4double A, B, C, D;
>> 254 G4double A2, p, q, r;
>> 255 G4int i,j, num;
233 256
234 G4int i, num ; << 257 // normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0
235 G4double c[5], srd[4], si[4] ; <<
236 258
237 G4double Rtor2 = fRtor*fRtor, r2 = r*r ; << 259 A = c[ 3 ]; // c[ 4 ]; since always c[4]==1 !
>> 260 B = c[ 2 ]; // c[ 4 ];
>> 261 C = c[ 1 ]; // c[ 4 ];
>> 262 D = c[ 0 ]; // c[ 4 ];
238 263
239 G4double pDotV = p.x()*v.x() + p.y()*v.y() + << 264 // substitute x = y - A/4 to eliminate cubic term:
240 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + << 265 // y^4 + py^2 + qy + r = 0
241 266
242 G4double d=pRad2 - Rtor2; << 267 A2 = A*A;
243 c[0] = 1.0 ; << 268 p = - 0.375*A2 + B;
244 c[1] = 4*pDotV ; << 269 q = 0.125*A2*A - 0.5*A*B + C;
245 c[2] = 2*( (d + 2*pDotV*pDotV - r2) + 2*Rto << 270 r = - 3.0/256*A2*A2 + 1.0/16*A2*B - 0.25*A*C + D;
246 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 <<
248 271
249 G4JTPolynomialSolver torusEq; << 272 // y^4 + py^2 + r = 0 and z=y^2 so y = +-sqrt(z1) and y = +-sqrt(z2)
>> 273
>> 274 if(q==0)
>> 275 {
>> 276 coeffs[ 0 ] = r;
>> 277 coeffs[ 1 ] = p;
>> 278 coeffs[ 2 ] = 1;
>> 279 num = SolveQuadratic(coeffs, s) ;
250 280
251 num = torusEq.FindRoots( c, 4, srd, si ); << 281 if(num)
252 << 282 {
253 for ( i = 0; i < num; ++i ) << 283 if(num==2)
>> 284 {
>> 285 if(s[0]>=0)
>> 286 {
>> 287 if(s[0]==0) // Three roots and one of them == 0
>> 288 {
>> 289 s[2] = sqrt(s[1]) ;
>> 290 s[1] = s[0] ;
>> 291 s[0] = -s[2] ;
>> 292 num++ ;
>> 293 }
>> 294 else // Four roots
>> 295 {
>> 296 s[2] = sqrt(s[0]) ;
>> 297 s[3] = sqrt(s[1]) ;
>> 298 s[0] = -s[3] ;
>> 299 s[1] = -s[2] ;
>> 300 num +=2 ;
>> 301 }
>> 302 }
>> 303 else if(s[1]>=0)
>> 304 {
>> 305 if(s[1]==0) // One root == 0
>> 306 {
>> 307 s[0] = 0 ;
>> 308 num--;
>> 309 }
>> 310 else // Two roots
>> 311 {
>> 312 s[0] = -sqrt(s[1]) ;
>> 313 s[1] = -s[0] ;
>> 314 }
>> 315 }
>> 316 else return num = 0 ; // Both Quadratic roots are negative
>> 317 }
>> 318 else // num = 1 two equal roots from SolveQuadratic
>> 319 {
>> 320 if(s[0]>=0)
>> 321 {
>> 322 if(s[0]==0) ;
>> 323 else
>> 324 {
>> 325 s[1] = sqrt(s[0]) ;
>> 326 s[0] = -s[1] ;
>> 327 num +=1 ;
>> 328 }
>> 329 }
>> 330 else return num = 0 ;
>> 331 }
>> 332 }
>> 333 else return num ;
>> 334 }
>> 335 else if (r == 0) // no absolute term: y(y^3 + py + q) = 0
254 { 336 {
255 if( si[i] == 0. ) { roots.push_back(srd[i << 337 coeffs[ 0 ] = q ;
256 } << 338 coeffs[ 1 ] = p ;
>> 339 coeffs[ 2 ] = 0 ;
>> 340 coeffs[ 3 ] = 1 ;
>> 341 num = SolveCubic(coeffs, s) ;
>> 342
>> 343 s[ num++ ] = 0;
>> 344
>> 345 for(j=1;j<num;j++) // picksort of roots in ascending order
>> 346 {
>> 347 sub = s[j] ;
>> 348 i=j-1 ;
>> 349 while( i >= 0 && s[i] > sub )
>> 350 {
>> 351 i-- ;
>> 352 s[i+1] = s[i] ; // s[i--] ;
>> 353 }
>> 354 s[i+1] = sub ;
>> 355 }
>> 356 }
>> 357 else
>> 358 {
>> 359 // solve the resolvent cubic ...
>> 360
>> 361 coeffs[ 0 ] = 0.5*r*p - 0.125*q*q;
>> 362 coeffs[ 1 ] = - r;
>> 363 coeffs[ 2 ] = - 0.5*p;
>> 364 coeffs[ 3 ] = 1;
>> 365
>> 366 num = SolveCubic(coeffs, s);
>> 367
>> 368 // ... and take the one real solution ...
>> 369
>> 370 z = s[ 0 ];
>> 371
>> 372 // ... to Build two quadratic equations
>> 373
>> 374 u = z * z - r;
>> 375 v = 2 * z - p;
>> 376
>> 377 if (u==0) u = 0 ;
>> 378 else if (u > 0) u = sqrt(u) ;
>> 379 else return 0 ;
>> 380
>> 381 if (v==0) v = 0 ;
>> 382 else if (v > 0) v = sqrt(v);
>> 383 else return 0 ;
>> 384
>> 385 coeffs[ 0 ] = z - u;
>> 386 coeffs[ 1 ] = q < 0 ? -v : v;
>> 387 coeffs[ 2 ] = 1;
>> 388
>> 389 num = SolveQuadratic(coeffs, s);
>> 390
>> 391 coeffs[ 0 ]= z + u;
>> 392 coeffs[ 1 ] = q < 0 ? v : -v;
>> 393 coeffs[ 2 ] = 1;
>> 394
>> 395 num += SolveQuadratic(coeffs, s + num);
>> 396 }
>> 397
>> 398 // resubstitute
257 399
258 std::sort(roots.begin() , roots.end() ) ; / << 400 sub = 1.0/4 * A;
>> 401
>> 402 for (i = 0; i < num; ++i)
>> 403 s[ i ] -= sub;
>> 404
>> 405 return num;
259 } 406 }
260 407
261 ////////////////////////////////////////////// << 408 /////////////////////////////////////////////////////////////////////////////
262 // 409 //
263 // Interface for DistanceToIn and DistanceToOu << 410 // Auxiliary method for solving of cubic equation in real numbers
264 // Calls TorusRootsJT and returns the smalles << 411 // From Graphics Gems I bu Jochen Schwartz
265 // the surface. <<
266 // Attention: Difference in DistanceToIn/Out f <<
267 412
268 G4double G4Torus::SolveNumericJT( const G4Thre << 413 G4int G4Torus::SolveCubic( G4double c[], G4double s[] ) const
269 const G4Thre <<
270 G4doub <<
271 G4bool <<
272 { 414 {
273 G4double bigdist = 10*mm ; << 415 G4int i, num;
274 G4double tmin = kInfinity ; << 416 G4double sub;
275 G4double t, scal ; << 417 G4double A, B, C;
>> 418 G4double A2, p, q;
>> 419 G4double p3, D;
276 420
277 // calculate the distances to the intersecti << 421 // normal form: x^3 + Ax^2 + Bx + C = 0
278 // from a given point p and direction v. <<
279 // <<
280 std::vector<G4double> roots ; <<
281 std::vector<G4double> rootsrefined ; <<
282 TorusRootsJT(p,v,r,roots) ; <<
283 422
284 G4ThreeVector ptmp ; << 423 A = c[ 2 ]; // c[ 3 ]; since always c[3]==1 !
>> 424 B = c[ 1 ]; // c[ 3 ];
>> 425 C = c[ 0 ]; // c[ 3 ];
285 426
286 // determine the smallest non-negative solut << 427 // substitute x = y - A/3 to eliminate quadric term:
287 // << 428 // x^3 +px + q = 0
288 for ( std::size_t k = 0 ; k<roots.size() ; + << 429
289 { << 430 A2 = A*A;
290 t = roots[k] ; << 431 p = 1.0/3*(- 1.0/3*A2 + B);
>> 432 q = 1.0/2*(2.0/27*A*A2 - 1.0/3*A*B + C);
>> 433
>> 434 // use Cardano's formula
291 435
292 if ( t < -halfCarTolerance ) { continue ; << 436 p3 = p*p*p;
>> 437 D = q*q + p3;
293 438
294 if ( t > bigdist && t<kInfinity ) // pr << 439 if (D == 0)
>> 440 {
>> 441 if (q == 0) // one triple solution
295 { 442 {
296 ptmp = p + t*v ; << 443 s[ 0 ] = 0;
297 TorusRootsJT(ptmp,v,r,rootsrefined) ; << 444 num = 1;
298 if ( rootsrefined.size()==roots.size() ) << 445 }
299 { << 446 else // one single and one double solution
300 t = t + rootsrefined[k] ; << 447 {
301 } << 448 G4double u = pow(-q,1./3.);
>> 449 s[ 0 ] = 2 * u;
>> 450 s[ 1 ] = - u;
>> 451 num = 2;
302 } 452 }
>> 453 }
>> 454 else if (D < 0) // Casus irreducibilis: three real solutions
>> 455 {
>> 456 G4double phi = 1.0/3 * acos(-q / sqrt(-p3));
>> 457 G4double t = 2 * sqrt(-p);
>> 458
>> 459 s[ 0 ] = t * cos(phi);
>> 460 s[ 1 ] = - t * cos(phi + M_PI / 3);
>> 461 s[ 2 ] = - t * cos(phi - M_PI / 3);
>> 462 num = 3;
>> 463 }
>> 464 else // one real solution
>> 465 {
>> 466 G4double sqrt_D = sqrt(D);
>> 467 G4double u = pow(sqrt_D - q,1./3.);
>> 468 G4double v = - pow(sqrt_D + q,1./3.);
>> 469
>> 470 s[ 0 ] = u + v;
>> 471 num = 1;
>> 472 }
>> 473
>> 474 // resubstitute
>> 475
>> 476 sub = 1.0/3 * A;
>> 477
>> 478 for (i = 0; i < num; ++i)
>> 479 s[ i ] -= sub;
>> 480
>> 481 return num;
>> 482 }
>> 483
>> 484 // ---------------------------------------------------------------------
>> 485
>> 486 G4int G4Torus::SolveBiQuadraticNew( G4double c[], G4double s[] ) const
>> 487 {
>> 488 // From drte4 by McLareni; rewritten by O.Cremonesi
303 489
304 ptmp = p + t*v ; // calculate the positi << 490 G4double coeffs[ 4 ];
>> 491 G4double w1, w2, w3;
>> 492 G4double sub;
>> 493 G4double A, B, C, D;
>> 494 G4double A2, p, q, r ;
>> 495 G4int i,j, num;
305 496
306 G4double theta = std::atan2(ptmp.y(),ptmp. << 497 // normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0
>> 498
>> 499 A = c[ 3 ]; // c[ 4 ]; since always c[4]==1 !
>> 500 B = c[ 2 ]; // c[ 4 ];
>> 501 C = c[ 1 ]; // c[ 4 ];
>> 502 D = c[ 0 ]; // c[ 4 ];
>> 503
>> 504 if( B==0 && C==0 )
>> 505 {
>> 506 if( D==0 )
>> 507 {
>> 508 s[0] = -A;
>> 509 s[1] = s[2] = s[3] = 0;
>> 510 return 4;
>> 511 }
>> 512 }
>> 513 else if( A==0 )
>> 514 {
>> 515 if( D>0 ) return 0;
>> 516 else
>> 517 {
>> 518 s[0] = sqrt( sqrt( -D ) );
>> 519 s[1] = -s[0];
>> 520 return 2;
>> 521 }
>> 522 }
>> 523
>> 524 // substitute x = y - A/4 to eliminate cubic term:
>> 525 // y^4 + py^2 + qy + r = 0
>> 526
>> 527 A2 = A*A;
>> 528 p = B - 3.0*A2/8.0;
>> 529 q = C - 0.5*A*( B-A2/4.0 );
>> 530 r = D - (A*C-A2/4.0*(B-A2*3.0/16.0))/4.0;
>> 531 coeffs[ 0 ] = -q*q/64.;
>> 532 coeffs[ 1 ] = (p*p/4.0-r)/4.0;
>> 533 coeffs[ 2 ] = p/2.0;
>> 534 coeffs[ 3 ] = 1;
307 535
308 if ( fSPhi >= 0 ) << 536 G4double cubic_discr;
>> 537 num = SolveCubicNew(coeffs, s, cubic_discr);
>> 538
>> 539 sub = A/4.0;
>> 540 num = 0;
>> 541
>> 542 if( cubic_discr == 0 ) s[2] = s[1];
>> 543
>> 544 if( cubic_discr <= 0 )
>> 545 {
>> 546 num = 4;
>> 547 G4double v[3];
>> 548 G4double vm1 = -1.0e99, vm2 ;
>> 549 for( i=0; i<3; i++ )
309 { 550 {
310 if ( theta < - halfAngTolerance ) { the << 551 v[i] = fabs( s[i] ) ;
311 if ( (std::fabs(theta) < halfAngToleranc << 552 if( v[i] > vm1 ) vm1 = v[i] ;
312 && (std::fabs(fSPhi + fDPhi - twopi) < <<
313 { <<
314 theta += twopi ; // 0 <= theta < 2pi <<
315 } <<
316 } 553 }
317 if ((fSPhi <= -pi )&&(theta>halfAngToleran << 554 if( vm1 == v[0] )
318 << 555 {
319 // We have to verify if this root is insid << 556 i = 0;
320 // fSPhi and fSPhi + fDPhi << 557 if( v[1] > v[2] ) vm2 = v[1];
321 // << 558 else vm2 = v[2];
322 if ( (theta - fSPhi >= - halfAngTolerance) << 559 }
323 && (theta - (fSPhi + fDPhi) <= halfAngT << 560 else if( vm1 == v[1] )
324 { 561 {
325 // check if P is on the surface, and cal << 562 i = 1;
326 // DistanceToIn has to return 0.0 if par << 563 if( v[0] > v[2] ) vm2 = v[0];
>> 564 else vm2 = v[2];
>> 565 }
>> 566 else
>> 567 {
>> 568 i = 2;
>> 569 if( v[0] > v[1] ) vm2 = v[0];
>> 570 else vm2 = v[1];
>> 571 }
>> 572 if( vm2 == v[0] ) j = 0 ;
>> 573 else if( vm2 == v[1] ) j = 1 ;
>> 574 else j = 2 ;
327 575
328 if ( IsDistanceToIn ) << 576 w1 = sqrt( s[i] );
>> 577 w2 = sqrt( s[j] );
>> 578 }
>> 579 else
>> 580 {
>> 581 num = 2;
>> 582 w1 = w2 = sqrt( s[1] );
>> 583 }
>> 584 if( w1*w2 != 0. ) w3 = -q/( 8.0*w1*w2 ) ;
>> 585 else w3 = 0.0 ;
>> 586
>> 587 if( num == 4 )
>> 588 {
>> 589 s[0] = w1 + w2 + w3 - sub ;
>> 590 s[1] = -w1 - w2 + w3 - sub ;
>> 591 s[2] = -w1 + w2 - w3 - sub ;
>> 592 s[3] = w1 - w2 - w3 - sub ;
>> 593 }
>> 594 else if( num == 2 )
>> 595 {
>> 596 s[0] = w1 + w2 + w3 - sub ;
>> 597 s[1] = -w1 - w2 + w3 - sub ;
>> 598 }
>> 599 return num ;
>> 600 }
>> 601
>> 602 // -------------------------------------------------------------------------
>> 603
>> 604 G4int G4Torus::SolveCubicNew( G4double c[], G4double s[],
>> 605 G4double& cubic_discr ) const
>> 606 {
>> 607 // From drte3 by McLareni; rewritten by O.Cremonesi
>> 608
>> 609 const G4double eps = 1.e-6;
>> 610 const G4double delta = 1.e-15;
>> 611 G4int i, j;
>> 612 G4double sub;
>> 613 G4double y[3];
>> 614 G4double A, B, C;
>> 615 G4double A2, p, q;
>> 616 G4double h1,h2,h3;
>> 617 G4double u,v;
>> 618
>> 619 // normal form: x^3 + Ax^2 + Bx + C = 0
>> 620
>> 621 A = c[ 2 ]; // c[ 3 ]; since always c[3]==1 !
>> 622 B = c[ 1 ]; // c[ 3 ];
>> 623 C = c[ 0 ]; // c[ 3 ];
>> 624
>> 625 if( B==0 && C==0 )
>> 626 {
>> 627 s[0] = -A;
>> 628 s[1] = s[2] = 0.;
>> 629 cubic_discr = 0.;
>> 630 return 3;
>> 631 }
>> 632 A2 = A*A;
>> 633 p = B - A2/3.0;
>> 634 q = ( A2*2.0/27.-B/3.0 )*A + C;
>> 635 cubic_discr = q*q/4.0 + p*p*p/27.0;
>> 636 sub = A/3.0;
>> 637 h1 = q/2.0;
>> 638
>> 639 if( cubic_discr > delta )
>> 640 {
>> 641 h2 = sqrt( cubic_discr );
>> 642 u = -h1+h2;
>> 643 v = -h1-h2;
>> 644 if( u < 0 ) u = -pow(-u,1./3.);
>> 645 else u = pow(u,1./3.);
>> 646 if( v < 0 ) v = -pow(-v,1./3.);
>> 647 else v = pow(v,1./3.);
>> 648 s[0] = u+v-sub;
>> 649 s[1] = -(u+v)/2.0-sub;
>> 650 s[2] = fabs(u-v)*sqrt(3.0)/2.0;
>> 651 if( fabs(u) <= eps || fabs(v) <= eps )
>> 652 {
>> 653 y[0] = s[0] ;
>> 654 for( i=0; i<2; i++ )
329 { 655 {
330 if (std::fabs(t) < halfCarTolerance ) << 656 y[i+1] = y[i] - (((y[i]+A)*y[i]+B)*y[i]+C)/((3.*y[i]+2.*A)*y[i]+B);
331 { <<
332 // compute scalar product at positio <<
333 // ( n taken from SurfaceNormal, not <<
334 <<
335 scal = v* G4ThreeVector( p.x()*(1-fR <<
336 p.y()*(1-fR <<
337 p.z() ); <<
338 <<
339 // change sign in case of inner radi <<
340 // <<
341 if ( r == GetRmin() ) { scal = -sca <<
342 if ( scal < 0 ) { return 0.0 ; } <<
343 } <<
344 } 657 }
>> 658 s[0] = y[2];
>> 659 return 1;
>> 660 }
>> 661 }
>> 662 else if( fabs(cubic_discr) <= delta )
>> 663 {
>> 664 cubic_discr = 0.;
>> 665
>> 666 if( h1 < 0 ) u = pow(-h1,1./3.);
>> 667 else u = -pow(h1,1./3.);
345 668
346 // check if P is on the surface, and cal << 669 s[0] = u + u - sub ;
347 // DistanceToIn has to return 0.0 if par << 670 s[1] = -u - sub ;
>> 671 s[2] = s[1] ;
348 672
349 if ( !IsDistanceToIn ) << 673 if( fabs(h1) <= eps )
>> 674 {
>> 675 y[0] = s[0];
>> 676 for( i=0; i<2; i++ )
350 { 677 {
351 if (std::fabs(t) < halfCarTolerance ) << 678 h1 = (3.0*y[i]+2.*A)*y[i]+B;
352 { <<
353 // compute scalar product at positio <<
354 // <<
355 scal = v* G4ThreeVector( p.x()*(1-fR <<
356 p.y()*(1-fR <<
357 p.z() ); <<
358 679
359 // change sign in case of inner radi << 680 if( fabs(h1) > delta )
360 // << 681 y[i+1] = y[i]-(((y[i]+A)*y[i]+B)*y[i]+C)/h1;
361 if ( r == GetRmin() ) { scal = -sca << 682 else
362 if ( scal > 0 ) { return 0.0 ; } << 683 {
>> 684 s[0] = s[1] = s[2] = -A/3.;
>> 685 return 3;
363 } 686 }
364 } 687 }
>> 688 s[0] = y[2];
>> 689 s[1] = s[2] = -(A+s[0])/2.;
>> 690 return 3;
>> 691 }
>> 692 }
>> 693 else
>> 694 {
>> 695 h3 =fabs(p/3.);
>> 696 h3 = sqrt(h3*h3*h3);
>> 697 h2 = acos(-h1/h3)/3.;
>> 698 h1 = pow(h3,1./3.);
>> 699 u = h1*cos(h2);
>> 700 v = sqrt(3.)*h1*sin(h2);
>> 701 s[0] = u+u-sub;
>> 702 s[1] = -u-v-sub;
>> 703 s[2] = -u+v-sub;
365 704
366 // check if distance is larger than 1/2 << 705 if( h3 <= eps || s[0] <=eps || s[1] <= eps || s[2] <= eps )
367 // << 706 {
368 if( t > halfCarTolerance ) << 707 for( i=0; i<3; i++ )
369 { 708 {
370 tmin = t ; << 709 y[0] = s[i] ;
371 return tmin ; << 710
>> 711 for( j=0; j<2; j++ )
>> 712 {
>> 713 y[j+1] = y[j]-(((y[j]+A)*y[j]+B)*y[j]+C)/((3.*y[j]+2.*A)*y[j]+B);
>> 714 }
>> 715 s[i] = y[2] ;
372 } 716 }
373 } 717 }
374 } 718 }
375 << 719 return 3;
376 return tmin; <<
377 } 720 }
378 721
379 ////////////////////////////////////////////// << 722 ///////////////////////////////////////////////////////////////////////////
380 // 723 //
381 // Get bounding box << 724 // Auxiliary method for solving quadratic equations in real numbers
>> 725 // From Graphics Gems I by Jochen Schwartz
382 726
383 void G4Torus::BoundingLimits(G4ThreeVector& pM << 727 G4int G4Torus::SolveQuadratic( G4double c[], G4double s[] ) const
384 { 728 {
385 G4double rmax = GetRmax(); << 729 G4double p, q, D;
386 G4double rtor = GetRtor(); <<
387 G4double rint = rtor - rmax; <<
388 G4double rext = rtor + rmax; <<
389 G4double dz = rmax; <<
390 730
391 // Find bounding box << 731 // normal form: x^2 + px + q = 0
392 // << 732
393 if (GetDPhi() >= twopi) << 733 p = c[ 1 ]/2 ; // * c[ 2 ]); since always c[2]==1
>> 734 q = c[ 0 ] ; // c[ 2 ];
>> 735
>> 736 D = p * p - q;
>> 737
>> 738 if (D==0)
394 { 739 {
395 pMin.set(-rext,-rext,-dz); << 740 s[ 0 ] = - p; // Generally we have two equal roots ?!
396 pMax.set( rext, rext, dz); << 741 return 1; // But consider them as one for geometry
397 } 742 }
398 else << 743 else if (D > 0)
399 { 744 {
400 G4TwoVector vmin,vmax; << 745 G4double sqrt_D = sqrt(D);
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 746
409 // Check correctness of the bounding box << 747 s[ 0 ] = - p - sqrt_D ; // in ascending order !
410 // << 748 s[ 1 ] = - p + sqrt_D ;
411 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 749 return 2;
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 } 750 }
>> 751 return 0;
422 } 752 }
423 753
424 ////////////////////////////////////////////// 754 /////////////////////////////////////////////////////////////////////////////
425 // 755 //
426 // Calculate extent under transform and specif 756 // Calculate extent under transform and specified limit
427 757
428 G4bool G4Torus::CalculateExtent( const EAxis p 758 G4bool G4Torus::CalculateExtent( const EAxis pAxis,
429 const G4Voxel 759 const G4VoxelLimits& pVoxelLimit,
430 const G4Affin 760 const G4AffineTransform& pTransform,
431 G4doubl 761 G4double& pMin, G4double& pMax) const
432 { 762 {
433 G4ThreeVector bmin, bmax; << 763 if (!pTransform.IsRotated() && fDPhi==2.0*M_PI && fRmin==0)
434 G4bool exist; << 764 {
>> 765 // Special case handling for unrotated solid torus
>> 766 // Compute x/y/z mins and maxs for bounding box respecting limits,
>> 767 // with early returns if outside limits. Then switch() on pAxis,
>> 768 // and compute exact x and y limit for x/y case
>> 769
>> 770 G4double xoffset,xMin,xMax;
>> 771 G4double yoffset,yMin,yMax;
>> 772 G4double zoffset,zMin,zMax;
>> 773
>> 774 G4double diff1,diff2,maxDiff,newMin,newMax;
>> 775 G4double xoff1,xoff2,yoff1,yoff2;
>> 776
>> 777 xoffset = pTransform.NetTranslation().x();
>> 778 xMin = xoffset - fRmax - fRtor ;
>> 779 xMax = xoffset + fRmax + fRtor ;
>> 780
>> 781 if (pVoxelLimit.IsXLimited())
>> 782 {
>> 783 if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance)
>> 784 || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) )
>> 785 return false ;
>> 786 else
>> 787 {
>> 788 if (xMin < pVoxelLimit.GetMinXExtent())
>> 789 {
>> 790 xMin = pVoxelLimit.GetMinXExtent() ;
>> 791 }
>> 792 if (xMax > pVoxelLimit.GetMaxXExtent())
>> 793 {
>> 794 xMax = pVoxelLimit.GetMaxXExtent() ;
>> 795 }
>> 796 }
>> 797 }
>> 798 yoffset = pTransform.NetTranslation().y();
>> 799 yMin = yoffset - fRmax - fRtor ;
>> 800 yMax = yoffset + fRmax + fRtor ;
>> 801
>> 802 if (pVoxelLimit.IsYLimited())
>> 803 {
>> 804 if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance)
>> 805 || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) )
>> 806 return false ;
>> 807 else
>> 808 {
>> 809 if (yMin < pVoxelLimit.GetMinYExtent() )
>> 810 {
>> 811 yMin = pVoxelLimit.GetMinYExtent() ;
>> 812 }
>> 813 if (yMax > pVoxelLimit.GetMaxYExtent() )
>> 814 {
>> 815 yMax = pVoxelLimit.GetMaxYExtent() ;
>> 816 }
>> 817 }
>> 818 }
>> 819 zoffset = pTransform.NetTranslation().z() ;
>> 820 zMin = zoffset - fRmax ;
>> 821 zMax = zoffset + fRmax ;
>> 822
>> 823 if (pVoxelLimit.IsZLimited())
>> 824 {
>> 825 if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance)
>> 826 || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) )
>> 827 return false ;
>> 828 else
>> 829 {
>> 830 if (zMin < pVoxelLimit.GetMinZExtent() )
>> 831 {
>> 832 zMin = pVoxelLimit.GetMinZExtent() ;
>> 833 }
>> 834 if (zMax > pVoxelLimit.GetMaxZExtent() )
>> 835 {
>> 836 zMax = pVoxelLimit.GetMaxZExtent() ;
>> 837 }
>> 838 }
>> 839 }
>> 840
>> 841 // Known to cut cylinder
>> 842
>> 843 switch (pAxis)
>> 844 {
>> 845 case kXAxis:
>> 846 yoff1=yoffset-yMin;
>> 847 yoff2=yMax-yoffset;
>> 848 if ( yoff1 >= 0 && yoff2 >= 0 )
>> 849 {
>> 850 // Y limits cross max/min x => no change
>> 851 //
>> 852 pMin = xMin ;
>> 853 pMax = xMax ;
>> 854 }
>> 855 else
>> 856 {
>> 857 // Y limits don't cross max/min x => compute max delta x,
>> 858 // hence new mins/maxs
>> 859 //
>> 860 diff1 = sqrt(fRmax*fRmax - yoff1*yoff1) ;
>> 861 diff2 = sqrt(fRmax*fRmax - yoff2*yoff2) ;
>> 862 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
>> 863 newMin = xoffset - maxDiff ;
>> 864 newMax = xoffset + maxDiff ;
>> 865 pMin = (newMin < xMin) ? xMin : newMin ;
>> 866 pMax = (newMax > xMax) ? xMax : newMax ;
>> 867 }
>> 868 break;
>> 869
>> 870 case kYAxis:
>> 871 xoff1 = xoffset - xMin ;
>> 872 xoff2 = xMax - xoffset ;
>> 873 if (xoff1 >= 0 && xoff2 >= 0 )
>> 874 {
>> 875 // X limits cross max/min y => no change
>> 876 //
>> 877 pMin = yMin ;
>> 878 pMax = yMax ;
>> 879 }
>> 880 else
>> 881 {
>> 882 // X limits don't cross max/min y => compute max delta y,
>> 883 // hence new mins/maxs
>> 884 //
>> 885 diff1 = sqrt(fRmax*fRmax - xoff1*xoff1) ;
>> 886 diff2 = sqrt(fRmax*fRmax - xoff2*xoff2) ;
>> 887 maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
>> 888 newMin = yoffset - maxDiff ;
>> 889 newMax = yoffset + maxDiff ;
>> 890 pMin = (newMin < yMin) ? yMin : newMin ;
>> 891 pMax = (newMax > yMax) ? yMax : newMax ;
>> 892 }
>> 893 break;
435 894
436 // Get bounding box << 895 case kZAxis:
437 BoundingLimits(bmin,bmax); << 896 pMin=zMin;
>> 897 pMax=zMax;
>> 898 break;
438 899
439 // Check bounding box << 900 default:
440 G4BoundingEnvelope bbox(bmin,bmax); << 901 break;
441 #ifdef G4BBOX_EXTENT << 902 }
442 return bbox.CalculateExtent(pAxis,pVoxelLimi << 903 pMin -= kCarTolerance ;
443 #endif << 904 pMax += kCarTolerance ;
444 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox <<
445 { <<
446 return exist = pMin < pMax; <<
447 } <<
448 <<
449 // Get parameters of the solid <<
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 905
461 // Find bounding envelope and calculate exte << 906 return true;
462 // << 907 }
463 static const G4int NPHI = 24; // number of << 908 else
464 static const G4int NDISK = 16; // number of << 909 {
465 static const G4double sinHalfDisk = std::sin << 910 G4int i, noEntries, noBetweenSections4 ;
466 static const G4double cosHalfDisk = std::cos << 911 G4bool existsAfterClip = false ;
467 static const G4double sinStepDisk = 2.*sinHa << 912
468 static const G4double cosStepDisk = 1. - 2.* << 913 // Calculate rotated vertex coordinates
469 << 914
470 G4double astep = (360/NPHI)*deg; // max angl << 915 G4ThreeVectorList *vertices ;
471 G4int kphi = (dphi <= astep) ? 1 : (G4in << 916 G4int noPolygonVertices ; // will be 4
472 G4double ang = dphi/kphi; << 917 vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ;
473 << 918
474 G4double sinHalf = std::sin(0.5*ang); << 919 pMin = +kInfinity ;
475 G4double cosHalf = std::cos(0.5*ang); << 920 pMax = -kInfinity ;
476 G4double sinStep = 2.*sinHalf*cosHalf; << 921
477 G4double cosStep = 1. - 2.*sinHalf*sinHalf; << 922 noEntries = vertices->size() ;
478 << 923 noBetweenSections4 = noEntries - noPolygonVertices ;
479 // define vectors for bounding envelope << 924
480 G4ThreeVectorList pols[NDISK+1]; << 925 for (i=0;i<noEntries;i+=noPolygonVertices)
481 for (auto & pol : pols) pol.resize(4); << 926 {
482 << 927 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
483 std::vector<const G4ThreeVectorList *> polyg << 928 }
484 polygons.resize(NDISK+1); << 929 for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
485 for (G4int k=0; k<NDISK+1; ++k) polygons[k] << 930 {
486 << 931 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
487 // set internal and external reference circl << 932 }
488 G4TwoVector rzmin[NDISK]; << 933 if (pMin!=kInfinity||pMax!=-kInfinity)
489 G4TwoVector rzmax[NDISK]; << 934 {
490 << 935 existsAfterClip = true ; // Add 2*tolerance to avoid precision troubles
491 if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+ << 936 pMin -= kCarTolerance ;
492 rmax /= cosHalfDisk; << 937 pMax += kCarTolerance ;
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 } 938 }
526 else 939 else
527 { 940 {
528 sinCur1 = sinCur2; << 941 // Check for case where completely enveloping clipping volume
529 cosCur1 = cosCur2; << 942 // If point inside then we are confident that the solid completely
530 sinCur2 = (i == kphi) ? sinEnd : sinCur1 << 943 // envelopes the clipping volume. Hence set min/max extents according
531 cosCur2 = (i == kphi) ? cosEnd : cosCur1 << 944 // to clipping volume extents along the specified axis.
532 } << 945
533 for (G4int k=0; k<NDISK; ++k) << 946 G4ThreeVector clipCentre(
534 { << 947 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
535 G4double r1 = rzmin[k].x(), r2 = rzmax[k << 948 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
536 G4double z1 = rzmin[k].y(), z2 = rzmax[k << 949 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
537 pols[k][0].set(r1*cosCur1,r1*sinCur1,z1) << 950
538 pols[k][1].set(r2*cosCur1,r2*sinCur1,z2) << 951 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
539 pols[k][2].set(r2*cosCur2,r2*sinCur2,z2) << 952 {
540 pols[k][3].set(r1*cosCur2,r1*sinCur2,z1) << 953 existsAfterClip = true ;
541 } << 954 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
542 pols[NDISK] = pols[0]; << 955 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
543 << 956 }
544 // get bounding box of current slice << 957 }
545 G4TwoVector vmin,vmax; << 958 delete vertices;
546 G4GeomTools:: << 959 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 } 960 }
559 return (pMin < pMax); <<
560 } 961 }
561 962
562 ////////////////////////////////////////////// << 963 ////////////////////////////////////////////////////////////////////////////////
563 // 964 //
564 // Return whether point inside/outside/on surf 965 // Return whether point inside/outside/on surface
565 966
566 EInside G4Torus::Inside( const G4ThreeVector& 967 EInside G4Torus::Inside( const G4ThreeVector& p ) const
567 { 968 {
568 G4double r, pt2, pPhi, tolRMin, tolRMax ; << 969 G4double r2, pt2, pPhi, tolRMin, tolRMax ;
569 970
570 EInside in = kOutside ; 971 EInside in = kOutside ;
>> 972 // General precals
>> 973 r2 = p.x()*p.x() + p.y()*p.y() ;
>> 974 pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*sqrt(r2) ;
571 975
572 // General precals << 976 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 ; 977 else tolRMin = 0 ;
579 978
580 tolRMax = fRmax - fRmaxTolerance; << 979 tolRMax = fRmax - kRadTolerance*0.5;
581 980
582 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax 981 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
583 { 982 {
584 if ( fDPhi == twopi || pt2 == 0 ) // on t << 983 if ( fDPhi == 2*M_PI || pt2 == 0 ) // on torus swept axis
585 { 984 {
586 in = kInside ; 985 in = kInside ;
587 } 986 }
588 else 987 else
589 { 988 {
590 // Try inner tolerant phi boundaries (=> 989 // Try inner tolerant phi boundaries (=>inside)
591 // if not inside, try outer tolerant phi 990 // if not inside, try outer tolerant phi boundaries
592 991
593 pPhi = std::atan2(p.y(),p.x()) ; << 992 pPhi = atan2(p.y(),p.x()) ;
594 993
595 if ( pPhi < -halfAngTolerance ) { pPhi << 994 if ( pPhi < 0 ) pPhi += 2*M_PI ; // 0<=pPhi<2*M_PI
596 if ( fSPhi >= 0 ) 995 if ( fSPhi >= 0 )
597 { 996 {
598 if ( (std::fabs(pPhi) < halfAngToleran << 997 if ( (pPhi >= fSPhi+kAngTolerance*0.5)
599 && (std::fabs(fSPhi + fDPhi - twop << 998 && (pPhi <= fSPhi+fDPhi-kAngTolerance*0.5) )
600 { <<
601 pPhi += twopi ; // 0 <= pPhi < 2pi <<
602 } <<
603 if ( (pPhi >= fSPhi + halfAngTolerance <<
604 && (pPhi <= fSPhi + fDPhi - halfAn <<
605 { <<
606 in = kInside ; 999 in = kInside ;
607 } << 1000 else if ( (pPhi >= fSPhi-kAngTolerance*0.5)
608 else if ( (pPhi >= fSPhi - halfAngTo << 1001 && (pPhi <= fSPhi+fDPhi+kAngTolerance*0.5) )
609 && (pPhi <= fSPhi + fDPhi + h <<
610 { <<
611 in = kSurface ; 1002 in = kSurface ;
612 } <<
613 } 1003 }
614 else // fSPhi < 0 << 1004 else
615 { 1005 {
616 if ( (pPhi <= fSPhi + twopi - halfAn << 1006 if (pPhi < fSPhi+2*M_PI) pPhi += 2*M_PI ;
617 && (pPhi >= fSPhi + fDPhi + halfA << 1007
618 else << 1008 if ( (pPhi >= fSPhi+2*M_PI+kAngTolerance*0.5)
619 { << 1009 && (pPhi <= fSPhi+fDPhi+2*M_PI-kAngTolerance*0.5) )
620 in = kSurface ; << 1010 in = kInside ;
621 } << 1011 else if ( (pPhi >= fSPhi+2*M_PI-kAngTolerance*0.5)
>> 1012 && (pPhi <= fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5) )
>> 1013 in = kSurface ;
622 } 1014 }
623 } 1015 }
624 } 1016 }
625 else // Try generous boundaries 1017 else // Try generous boundaries
626 { 1018 {
627 tolRMin = fRmin - fRminTolerance ; << 1019 tolRMin = fRmin - kRadTolerance*0.5 ;
628 tolRMax = fRmax + fRmaxTolerance ; << 1020 tolRMax = fRmax + kRadTolerance*0.5 ;
629 1021
630 if (tolRMin < 0 ) { tolRMin = 0 ; } << 1022 if (tolRMin < 0 ) tolRMin = 0 ;
631 1023
632 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= t << 1024 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax)
633 { 1025 {
634 if ( (fDPhi == twopi) || (pt2 == 0) ) // << 1026 if (fDPhi == 2*M_PI || pt2 == 0 ) // Continuous in phi or on z-axis
635 { 1027 {
636 in = kSurface ; 1028 in = kSurface ;
637 } 1029 }
638 else // Try outer tolerant phi boundarie 1030 else // Try outer tolerant phi boundaries only
639 { 1031 {
640 pPhi = std::atan2(p.y(),p.x()) ; << 1032 pPhi = atan2(p.y(),p.x()) ;
641 1033
642 if ( pPhi < -halfAngTolerance ) { pPh << 1034 if ( pPhi < 0 ) pPhi += 2*M_PI ; // 0<=pPhi<2*M_PI
643 if ( fSPhi >= 0 ) << 1035
644 { << 1036 if (fSPhi >= 0 )
645 if ( (std::fabs(pPhi) < halfAngToler <<
646 && (std::fabs(fSPhi + fDPhi - twop <<
647 { <<
648 pPhi += twopi ; // 0 <= pPhi < 2pi <<
649 } <<
650 if ( (pPhi >= fSPhi - halfAngToleran <<
651 && (pPhi <= fSPhi + fDPhi + halfAn <<
652 { <<
653 in = kSurface; <<
654 } <<
655 } <<
656 else // fSPhi < 0 <<
657 { 1037 {
658 if ( (pPhi <= fSPhi + twopi - halfAn << 1038 if ( (pPhi >= fSPhi-kAngTolerance*0.5)
659 && (pPhi >= fSPhi + fDPhi + halfA << 1039 && (pPhi <= fSPhi+fDPhi+kAngTolerance*0.5) )
660 else <<
661 { <<
662 in = kSurface ; 1040 in = kSurface ;
663 } <<
664 } 1041 }
>> 1042 else
>> 1043 {
>> 1044 if (pPhi < fSPhi + 2*M_PI) pPhi += 2*M_PI ;
>> 1045
>> 1046 if ( (pPhi >= fSPhi+2*M_PI-kAngTolerance*0.5)
>> 1047 && (pPhi <= fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5) )
>> 1048 in = kSurface ;
>> 1049 }
665 } 1050 }
666 } 1051 }
667 } 1052 }
668 return in ; 1053 return in ;
669 } 1054 }
670 1055
671 ////////////////////////////////////////////// 1056 /////////////////////////////////////////////////////////////////////////////
672 // 1057 //
673 // Return unit normal of surface closest to p 1058 // Return unit normal of surface closest to p
674 // - note if point on z axis, ignore phi divid 1059 // - note if point on z axis, ignore phi divided sides
675 // - unsafe if point close to z axis a rmin=0 1060 // - unsafe if point close to z axis a rmin=0 - no explicit checks
676 1061
677 G4ThreeVector G4Torus::SurfaceNormal( const G4 1062 G4ThreeVector G4Torus::SurfaceNormal( const G4ThreeVector& p ) const
678 { 1063 {
679 G4int noSurfaces = 0; <<
680 G4double rho, pt, pPhi; <<
681 G4double distRMin = kInfinity; <<
682 G4double distSPhi = kInfinity, distEPhi = kI <<
683 <<
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; <<
691 G4ThreeVector norm, sumnorm(0.,0.,0.); <<
692 <<
693 rho = std::hypot(p.x(),p.y()); <<
694 pt = std::hypot(p.z(),rho-fRtor); <<
695 <<
696 G4double distRMax = std::fabs(pt - fRmax); <<
697 if(fRmin != 0.0) distRMin = std::fabs(pt - f <<
698 <<
699 if( rho > delta && pt != 0.0 ) <<
700 { <<
701 G4double redFactor= (rho-fRtor)/rho; <<
702 nR = G4ThreeVector( p.x()*redFactor, // p <<
703 p.y()*redFactor, // p <<
704 p.z() ); <<
705 nR *= 1.0/pt; <<
706 } <<
707 <<
708 if ( fDPhi < twopi ) // && rho ) // old limi <<
709 { <<
710 if ( rho != 0.0 ) <<
711 { <<
712 pPhi = std::atan2(p.y(),p.x()); <<
713 <<
714 if(pPhi < fSPhi-delta) { pPhi <<
715 else if(pPhi > fSPhi+fDPhi+delta) { pPhi <<
716 <<
717 distSPhi = std::fabs( pPhi - fSPhi ); <<
718 distEPhi = std::fabs(pPhi-fSPhi-fDPhi); <<
719 } <<
720 nPs = G4ThreeVector(std::sin(fSPhi),-std:: <<
721 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi) <<
722 } <<
723 if( distRMax <= delta ) <<
724 { <<
725 ++noSurfaces; <<
726 sumnorm += nR; <<
727 } <<
728 else if( (fRmin != 0.0) && (distRMin <= delt <<
729 { <<
730 ++noSurfaces; <<
731 sumnorm -= nR; <<
732 } <<
733 <<
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 { <<
739 if (distSPhi <= dAngle) <<
740 { <<
741 ++noSurfaces; <<
742 sumnorm += nPs; <<
743 } <<
744 if (distEPhi <= dAngle) <<
745 { <<
746 ++noSurfaces; <<
747 sumnorm += nPe; <<
748 } <<
749 } <<
750 if ( noSurfaces == 0 ) <<
751 { <<
752 #ifdef G4CSGDEBUG <<
753 G4ExceptionDescription ed; <<
754 ed.precision(16); <<
755 <<
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); <<
795 } <<
796 else if ( noSurfaces == 1 ) { norm = sumnor <<
797 else { norm = sumnor <<
798 <<
799 return norm ; <<
800 } <<
801 <<
802 ////////////////////////////////////////////// <<
803 // <<
804 // Algorithm for SurfaceNormal() following the <<
805 // for points not on the surface <<
806 <<
807 G4ThreeVector G4Torus::ApproxSurfaceNormal( co <<
808 { <<
809 ENorm side ; 1064 ENorm side ;
810 G4ThreeVector norm; 1065 G4ThreeVector norm;
811 G4double rho,pt,phi; << 1066 G4double rho2,rho,pt2,pt,phi;
812 G4double distRMin,distRMax,distSPhi,distEPhi 1067 G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
813 1068
814 rho = std::hypot(p.x(),p.y()); << 1069 rho2 = p.x()*p.x() + p.y()*p.y();
815 pt = std::hypot(p.z(),rho-fRtor); << 1070 rho = sqrt(rho2) ;
>> 1071 pt2 = fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho) ;
>> 1072 pt = sqrt(pt2) ;
816 1073
817 #ifdef G4CSGDEBUG << 1074 distRMax = fabs(pt - fRmax) ;
818 G4cout << " G4Torus::ApproximateSurfaceNorma <<
819 << G4endl; <<
820 #endif <<
821 <<
822 distRMax = std::fabs(pt - fRmax) ; <<
823 1075
824 if(fRmin != 0.0) // First minimum radius << 1076
>> 1077 if(fRmin) // First minimum radius
825 { 1078 {
826 distRMin = std::fabs(pt - fRmin) ; << 1079 distRMin = fabs(pt - fRmin) ;
827 1080
828 if (distRMin < distRMax) 1081 if (distRMin < distRMax)
829 { 1082 {
830 distMin = distRMin ; 1083 distMin = distRMin ;
831 side = kNRMin ; 1084 side = kNRMin ;
832 } 1085 }
833 else 1086 else
834 { 1087 {
835 distMin = distRMax ; 1088 distMin = distRMax ;
836 side = kNRMax ; 1089 side = kNRMax ;
837 } 1090 }
838 } 1091 }
839 else 1092 else
840 { 1093 {
841 distMin = distRMax ; 1094 distMin = distRMax ;
842 side = kNRMax ; 1095 side = kNRMax ;
843 } 1096 }
844 if ( (fDPhi < twopi) && (rho != 0.0) ) << 1097 if (fDPhi < 2.0*M_PI && rho )
845 { 1098 {
846 phi = std::atan2(p.y(),p.x()) ; // Protect << 1099 phi = atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho !=0)
847 1100
848 if (phi < 0) { phi += twopi ; } << 1101 if (phi < 0) phi += 2*M_PI ;
849 1102
850 if (fSPhi < 0 ) { distSPhi = std::fabs(ph << 1103 if (fSPhi < 0 ) distSPhi = fabs(phi-(fSPhi+2.0*M_PI))*rho ;
851 else { distSPhi = std::fabs(ph << 1104 else distSPhi = fabs(phi-fSPhi)*rho ;
852 1105
853 distEPhi = std::fabs(phi - fSPhi - fDPhi)* << 1106 distEPhi = fabs(phi - fSPhi - fDPhi)*rho ;
854 1107
855 if (distSPhi < distEPhi) // Find new minim 1108 if (distSPhi < distEPhi) // Find new minimum
856 { 1109 {
857 if (distSPhi<distMin) side = kNSPhi ; 1110 if (distSPhi<distMin) side = kNSPhi ;
858 } 1111 }
859 else 1112 else
860 { 1113 {
861 if (distEPhi < distMin) { side = kNEPhi << 1114 if (distEPhi < distMin) side = kNEPhi ;
862 } 1115 }
863 } 1116 }
864 switch (side) 1117 switch (side)
865 { 1118 {
866 case kNRMin: // Inner radius 1119 case kNRMin: // Inner radius
867 norm = G4ThreeVector( -p.x()*(1-fRtor/rh 1120 norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
868 -p.y()*(1-fRtor/rh 1121 -p.y()*(1-fRtor/rho)/pt,
869 -p.z()/pt 1122 -p.z()/pt ) ;
870 break ; 1123 break ;
871 case kNRMax: // Outer radius 1124 case kNRMax: // Outer radius
872 norm = G4ThreeVector( p.x()*(1-fRtor/rho 1125 norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
873 p.y()*(1-fRtor/rho 1126 p.y()*(1-fRtor/rho)/pt,
874 p.z()/pt 1127 p.z()/pt ) ;
875 break; 1128 break;
876 case kNSPhi: 1129 case kNSPhi:
877 norm = G4ThreeVector(std::sin(fSPhi),-st << 1130 norm = G4ThreeVector(sin(fSPhi),-cos(fSPhi),0) ;
878 break; 1131 break;
879 case kNEPhi: 1132 case kNEPhi:
880 norm = G4ThreeVector(-std::sin(fSPhi+fDP << 1133 norm = G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0) ;
881 break; 1134 break;
882 default: // Should never reach th << 1135 default:
883 DumpInfo(); 1136 DumpInfo();
884 G4Exception("G4Torus::ApproxSurfaceNorma << 1137 G4Exception("G4Torus::SurfaceNormal() - Logic error");
885 "GeomSolids1002", JustWarnin <<
886 "Undefined side for valid su <<
887 break ; 1138 break ;
888 } 1139 }
889 return norm ; 1140 return norm ;
890 } 1141 }
891 1142
892 ////////////////////////////////////////////// 1143 ///////////////////////////////////////////////////////////////////////
893 // 1144 //
894 // Calculate distance to shape from outside, a 1145 // Calculate distance to shape from outside, along normalised vector
895 // - return kInfinity if no intersection, or i 1146 // - return kInfinity if no intersection, or intersection distance <= tolerance
896 // 1147 //
897 // - Compute the intersection with the z plane 1148 // - Compute the intersection with the z planes
898 // - if at valid r, phi, return 1149 // - if at valid r, phi, return
899 // 1150 //
900 // -> If point is outer outer radius, compute 1151 // -> If point is outer outer radius, compute intersection with rmax
901 // - if at valid phi,z return 1152 // - if at valid phi,z return
902 // 1153 //
903 // -> Compute intersection with inner radius, 1154 // -> Compute intersection with inner radius, taking largest +ve root
904 // - if valid (phi), save intersction 1155 // - if valid (phi), save intersction
905 // 1156 //
906 // -> If phi segmented, compute intersectio 1157 // -> If phi segmented, compute intersections with phi half planes
907 // - return smallest of valid phi inter 1158 // - return smallest of valid phi intersections and
908 // inner radius intersection 1159 // inner radius intersection
909 // 1160 //
910 // NOTE: 1161 // NOTE:
911 // - Precalculations for phi trigonometry are 1162 // - Precalculations for phi trigonometry are Done `just in time'
912 // - `if valid' implies tolerant checking of i 1163 // - `if valid' implies tolerant checking of intersection points
913 1164
914 G4double G4Torus::DistanceToIn( const G4ThreeV 1165 G4double G4Torus::DistanceToIn( const G4ThreeVector& p,
915 const G4ThreeV 1166 const G4ThreeVector& v ) const
916 { 1167 {
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 1168
925 // Check if point is traveling away <<
926 // <<
927 G4double distX = std::abs(p.x()) - boxDx; <<
928 G4double distY = std::abs(p.y()) - boxDy; <<
929 G4double distZ = std::abs(p.z()) - boxDz; <<
930 if (distX >= -halfCarTolerance && p.x()*v.x( <<
931 if (distY >= -halfCarTolerance && p.y()*v.y( <<
932 if (distZ >= -halfCarTolerance && p.z()*v.z( <<
933 <<
934 // Calculate safety distance to bounding box <<
935 // If point is too far, move it closer and c <<
936 // <<
937 G4double Dmax = 32*boxMax; <<
938 G4double safe = std::max(std::max(distX,dist <<
939 if (safe > Dmax) <<
940 { <<
941 G4double dist = safe - 1.e-8*safe - boxMin <<
942 dist += DistanceToIn(p + dist*v, v); <<
943 return (dist >= kInfinity) ? kInfinity : d <<
944 } <<
945 <<
946 // Find intersection with torus <<
947 // <<
948 G4double snxt=kInfinity, sphi=kInfinity; // 1169 G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value
949 1170
950 G4double sd[4] ; << 1171 G4double c[5], s[4] ;
951 1172
952 // Precalculated trig for phi intersections 1173 // Precalculated trig for phi intersections - used by r,z intersections to
953 // 1174 // check validity
954 1175
955 G4bool seg; // true if segmented 1176 G4bool seg; // true if segmented
956 G4double hDPhi; // half dphi << 1177 G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0.;
>> 1178 // half dphi + outer tolerance
957 G4double cPhi,sinCPhi=0.,cosCPhi=0.; // cen 1179 G4double cPhi,sinCPhi=0.,cosCPhi=0.; // central phi
958 1180
959 G4double tolORMin2; // `generous' radii squ << 1181 G4double tolORMin2,tolIRMin2; // `generous' radii squared
960 G4double tolORMax2; << 1182 G4double tolORMax2,tolIRMax2 ;
>> 1183
>> 1184 G4double Dist,xi,yi,zi,rhoi2,it2,inum,cosPsi; // Intersection point variables
961 1185
962 G4double Dist,xi,yi,zi,rhoi,it2; // Intersec <<
963 1186
964 G4double Comp; 1187 G4double Comp;
965 G4double cosSPhi,sinSPhi; // Trig for 1188 G4double cosSPhi,sinSPhi; // Trig for phi start intersect
966 G4double ePhi,cosEPhi,sinEPhi; // for phi e 1189 G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect
967 1190
>> 1191 #if DEBUGTORUS
>> 1192 G4cout << "G4Torus::DistanceToIn " << p << ", " << v << G4endl;
>> 1193 #endif
>> 1194
968 // Set phi divided flag and precalcs 1195 // Set phi divided flag and precalcs
969 // << 1196
970 if ( fDPhi < twopi ) << 1197 if ( fDPhi < 2.0*M_PI )
971 { 1198 {
972 seg = true ; 1199 seg = true ;
973 hDPhi = 0.5*fDPhi ; // half delta 1200 hDPhi = 0.5*fDPhi ; // half delta phi
974 cPhi = fSPhi + hDPhi ; 1201 cPhi = fSPhi + hDPhi ;
975 sinCPhi = std::sin(cPhi) ; << 1202 hDPhiOT = hDPhi+0.5*kAngTolerance ; // outers tol' half delta phi
976 cosCPhi = std::cos(cPhi) ; << 1203 hDPhiIT = hDPhi - 0.5*kAngTolerance ;
977 } << 1204 sinCPhi = sin(cPhi) ;
978 else << 1205 cosCPhi = cos(cPhi) ;
979 { << 1206 cosHDPhiOT = cos(hDPhiOT) ;
980 seg = false ; << 1207 cosHDPhiIT = cos(hDPhiIT) ;
981 } 1208 }
>> 1209 else seg = false ;
982 1210
983 if (fRmin > fRminTolerance) // Calculate tol << 1211 if (fRmin > kRadTolerance) // Calculate tolerant rmin and rmax
984 { 1212 {
985 tolORMin2 = (fRmin - fRminTolerance)*(fRmi << 1213 tolORMin2 = (fRmin - 0.5*kRadTolerance)*(fRmin - 0.5*kRadTolerance) ;
>> 1214 tolIRMin2 = (fRmin + 0.5*kRadTolerance)*(fRmin + 0.5*kRadTolerance) ;
986 } 1215 }
987 else 1216 else
988 { 1217 {
989 tolORMin2 = 0 ; 1218 tolORMin2 = 0 ;
>> 1219 tolIRMin2 = 0 ;
990 } 1220 }
991 tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax << 1221 tolORMax2 = (fRmax + 0.5*kRadTolerance)*(fRmax + 0.5*kRadTolerance) ;
>> 1222 tolIRMax2 = (fRmax - kRadTolerance*0.5)*(fRmax - kRadTolerance*0.5) ;
992 1223
993 // Intersection with Rmax (possible return) 1224 // Intersection with Rmax (possible return) and Rmin (must also check phi)
994 1225
995 snxt = SolveNumericJT(p,v,fRmax,true); << 1226 G4int i, j, num ;
>> 1227 G4double Rtor2 = fRtor*fRtor, Rmax2 = fRmax*fRmax, Rmin2 = fRmin*fRmin ;
>> 1228 G4double rho2 = p.x()*p.x()+p.y()*p.y();
>> 1229 G4double rho = sqrt(rho2) ;
>> 1230 G4double pt2 = fabs(rho2+p.z()*p.z() +Rtor2 - 2*fRtor*rho) ;
>> 1231 // G4double pt = sqrt(pt2) ;
>> 1232 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
>> 1233 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
>> 1234 G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
>> 1235
>> 1236 // Inside outer radius :
>> 1237 // check not inside, and heading through tubs (-> 0 to in)
996 1238
997 if (fRmin != 0.0) // Possible Rmin intersec << 1239 if( pt2 <= tolORMax2 && pt2 >= tolIRMin2 && vDotNmax < 0 )
998 { 1240 {
999 sd[0] = SolveNumericJT(p,v,fRmin,true); << 1241 if (seg)
1000 if ( sd[0] < snxt ) { snxt = sd[0] ; } << 1242 {
>> 1243 inum = p.x()*cosCPhi + p.y()*sinCPhi ;
>> 1244 cosPsi = inum/rho ;
>> 1245
>> 1246 if (cosPsi>=cosHDPhiIT)
>> 1247 {
>> 1248 #if DEBUGTORUS
>> 1249 G4cout << "G4Torus::DistanceToIn (cosPsi>=cosHDPhiIT) "
>> 1250 << __LINE__ << G4endl << G4endl;
>> 1251 #endif
>> 1252 return snxt = 0 ;
>> 1253 }
>> 1254 }
>> 1255 else
>> 1256 {
>> 1257 #if DEBUGTORUS
>> 1258 G4cout << "G4Torus::DistanceToIn (seg) "
>> 1259 << __LINE__ << G4endl << G4endl;
>> 1260 #endif
>> 1261 return snxt = 0 ;
>> 1262 }
1001 } 1263 }
>> 1264 else // intersection with Rmax torus
>> 1265 {
>> 1266 c[4] = 1.0 ;
>> 1267 c[3] = 4*pDotV ;
>> 1268 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - Rmax2 + 2*Rtor2*v.z()*v.z()) ;
>> 1269
>> 1270 c[1] = 4*(pDotV*(pRad2 - Rtor2 - Rmax2) + 2*Rtor2*p.z()*v.z()) ;
>> 1271
>> 1272 c[0] = pRad2*pRad2 - 2*pRad2*(Rtor2+Rmax2)
>> 1273 + 4*Rtor2*p.z()*p.z() + (Rtor2-Rmax2)*(Rtor2-Rmax2) ;
>> 1274
>> 1275 // Uncomment the line below for activating the analytical method
>> 1276 //
>> 1277 // num = SolveBiQuadratic(c,s) ;
>> 1278
>> 1279 // Numerical root research
>> 1280 //
>> 1281 s[0] = SolveNumeric(p, v, true);
>> 1282 num = 1; // There is only one root: the correct one
>> 1283
>> 1284 #if DEBUGTORUS
>> 1285 G4cout << "G4Torus::DistanceToIn (" << __LINE__ << ") SolveNumeric : "
>> 1286 << s[0] << G4endl;
>> 1287 #endif
>> 1288
>> 1289 if(num)
>> 1290 {
>> 1291 for(i=0;i<num;i++) // leave only >=kRadTolerance/2 roots P?!
>> 1292 {
>> 1293 if(s[i]<kRadTolerance*0.5)
>> 1294 {
>> 1295 for(j=i+1;j<num;j++) s[j-1] = s[j] ;
>> 1296 i-- ;
>> 1297 num-- ;
>> 1298 }
>> 1299 }
>> 1300 if(num)
>> 1301 {
>> 1302 for(i=0;i<num;i++)
>> 1303 {
>> 1304 if (seg) // intersection point must have proper Phi
>> 1305 {
>> 1306 xi = p.x() + s[i]*v.x() ;
>> 1307 yi = p.y() + s[i]*v.y() ;
>> 1308 rhoi2 = xi*xi + yi*yi ;
>> 1309 inum = xi*cosCPhi + yi*sinCPhi ;
>> 1310 cosPsi = inum/sqrt(rhoi2) ;
>> 1311
>> 1312 if (cosPsi >= cosHDPhiIT)
>> 1313 {
>> 1314 snxt = s[i] ;
>> 1315 break ;
>> 1316 }
>> 1317 }
>> 1318 else
>> 1319 {
>> 1320 snxt = s[i] ;
>> 1321 break ;
>> 1322 }
>> 1323 }
>> 1324 }
>> 1325 }
>> 1326 }
>> 1327 if (fRmin) // Possible Rmin intersection
>> 1328 {
>> 1329 // Inside relative to inner radius :
>> 1330 // check not inside, and heading through tubs (-> 0 to in)
1002 1331
>> 1332 if( pt2 >= tolORMin2 && pt2 <= tolIRMax2 && vDotNmax > 0 )
>> 1333 {
>> 1334 if (seg)
>> 1335 {
>> 1336 inum = p.x()*cosCPhi + p.y()*sinCPhi;
>> 1337 cosPsi = inum/rho ;
>> 1338
>> 1339 if (cosPsi>=cosHDPhiIT)
>> 1340 {
>> 1341 #if DEBUGTORUS
>> 1342 G4cout << "G4Torus::DistanceToIn (cosPsi>=cosHDPhiIT) "
>> 1343 << __LINE__ << G4endl << G4endl;
>> 1344 #endif
>> 1345 return snxt = 0 ;
>> 1346 }
>> 1347 }
>> 1348 else
>> 1349 {
>> 1350 #if DEBUGTORUS
>> 1351 G4cout << "G4Torus::DistanceToIn (seg) "
>> 1352 << __LINE__ << G4endl << G4endl;
>> 1353 #endif
>> 1354 return snxt = 0 ;
>> 1355 }
>> 1356 }
>> 1357 else // intersection with Rmin torus
>> 1358 {
>> 1359 c[4] = 1.0 ;
>> 1360 c[3] = 4*pDotV ;
>> 1361 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - Rmin2 + 2*Rtor2*v.z()*v.z()) ;
>> 1362
>> 1363 c[1] = 4*(pDotV*(pRad2-Rtor2-Rmin2) + 2*Rtor2*p.z()*v.z()) ;
>> 1364
>> 1365 c[0] = pRad2*pRad2 - 2*pRad2*(Rtor2+Rmin2)
>> 1366 + 4*Rtor2*p.z()*p.z() + (Rtor2-Rmin2)*(Rtor2-Rmin2) ;
>> 1367
>> 1368 // Uncomment the line below to activate the analytical method
>> 1369 //
>> 1370 // num = SolveBiQuadratic(c,s) ;
>> 1371
>> 1372 // Numerical root research
>> 1373 // s[0] = s[0]; // We already take care of Rmin in SolveNumeric !
>> 1374 num = 1;
>> 1375
>> 1376 if(num)
>> 1377 {
>> 1378 for(i=0;i<num;i++) // leave only >=kRadTolerance/2 roots P?!
>> 1379 {
>> 1380 if(s[i] < kRadTolerance*0.5)
>> 1381 {
>> 1382 for(j=i+1;j<num;j++) s[j-1] = s[j] ;
>> 1383 i-- ;
>> 1384 num-- ;
>> 1385 }
>> 1386 }
>> 1387 if(num)
>> 1388 {
>> 1389 for(i = 0 ; i < num ; i++ )
>> 1390 {
>> 1391 if (seg) // intersection point must have proper Phi
>> 1392 {
>> 1393 xi = p.x() + s[i]*v.x() ;
>> 1394 yi = p.y() + s[i]*v.y() ;
>> 1395 rhoi2 = xi*xi + yi*yi ;
>> 1396 inum = xi*cosCPhi + yi*sinCPhi ;
>> 1397 cosPsi = inum/sqrt(rhoi2) ;
>> 1398
>> 1399 if ( cosPsi >= cosHDPhiIT && s[i] < snxt )
>> 1400 {
>> 1401 snxt = s[i] ;
>> 1402 break ;
>> 1403 }
>> 1404 }
>> 1405 else if(s[i] < snxt)
>> 1406 {
>> 1407 snxt = s[i] ;
>> 1408 break ;
>> 1409 }
>> 1410 }
>> 1411 }
>> 1412 }
>> 1413 }
>> 1414 } // if(Rmin)
>> 1415
1003 // 1416 //
1004 // Phi segment intersection 1417 // Phi segment intersection
1005 // 1418 //
1006 // o Tolerant of points inside phi planes b 1419 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1007 // 1420 //
1008 // o NOTE: Large duplication of code betwee 1421 // o NOTE: Large duplication of code between sphi & ephi checks
1009 // -> only diffs: sphi -> ephi, Com 1422 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1010 // intersection check <=0 -> >=0 1423 // intersection check <=0 -> >=0
1011 // -> use some form of loop Constru 1424 // -> use some form of loop Construct ?
1012 1425
1013 if (seg) 1426 if (seg)
1014 { << 1427 {
1015 sinSPhi = std::sin(fSPhi) ; // First phi << 1428 sinSPhi = sin(fSPhi) ; // First phi surface (`S'tarting phi)
1016 cosSPhi = std::cos(fSPhi) ; << 1429 cosSPhi = cos(fSPhi) ;
1017 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1430 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; // Component in outwards
1018 << 1431 // normal direction
1019 if (Comp < 0 ) 1432 if (Comp < 0 )
1020 { 1433 {
1021 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) 1434 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1022 1435
1023 if (Dist < halfCarTolerance) << 1436 if (Dist < kCarTolerance*0.5)
1024 { 1437 {
1025 sphi = Dist/Comp ; 1438 sphi = Dist/Comp ;
1026 if (sphi < snxt) 1439 if (sphi < snxt)
1027 { 1440 {
1028 if ( sphi < 0 ) { sphi = 0 ; } << 1441 if ( sphi < 0 ) sphi = 0 ;
1029 1442
1030 xi = p.x() + sphi*v.x() ; 1443 xi = p.x() + sphi*v.x() ;
1031 yi = p.y() + sphi*v.y() ; 1444 yi = p.y() + sphi*v.y() ;
1032 zi = p.z() + sphi*v.z() ; 1445 zi = p.z() + sphi*v.z() ;
1033 rhoi = std::hypot(xi,yi); << 1446 rhoi2 = xi*xi + yi*yi ;
1034 it2 = zi*zi + (rhoi-fRtor)*(rhoi-fR << 1447 it2 = fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*sqrt(rhoi2)) ;
1035 1448
1036 if ( it2 >= tolORMin2 && it2 <= tol 1449 if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
1037 { 1450 {
1038 // r intersection is good - check 1451 // r intersection is good - check intersecting
1039 // with correct half-plane 1452 // with correct half-plane
1040 // 1453 //
1041 if ((yi*cosCPhi-xi*sinCPhi)<=0) << 1454 if ((yi*cosCPhi-xi*sinCPhi)<=0) snxt=sphi;
1042 } << 1455 }
1043 } 1456 }
1044 } 1457 }
1045 } 1458 }
1046 ePhi=fSPhi+fDPhi; // Second phi surfac << 1459 ePhi=fSPhi+fDPhi; // Second phi surface (`E'nding phi)
1047 sinEPhi=std::sin(ePhi); << 1460 sinEPhi=sin(ePhi);
1048 cosEPhi=std::cos(ePhi); << 1461 cosEPhi=cos(ePhi);
1049 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi); 1462 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
1050 1463
1051 if ( Comp < 0 ) // Component in outward 1464 if ( Comp < 0 ) // Component in outwards normal dirn
1052 { 1465 {
1053 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) 1466 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1054 1467
1055 if (Dist < halfCarTolerance ) << 1468 if (Dist < kCarTolerance*0.5 )
1056 { 1469 {
1057 sphi = Dist/Comp ; 1470 sphi = Dist/Comp ;
1058 <<
1059 if (sphi < snxt ) 1471 if (sphi < snxt )
1060 { 1472 {
1061 if (sphi < 0 ) { sphi = 0 ; } << 1473 if (sphi < 0 ) sphi = 0 ;
1062 1474
1063 xi = p.x() + sphi*v.x() ; 1475 xi = p.x() + sphi*v.x() ;
1064 yi = p.y() + sphi*v.y() ; 1476 yi = p.y() + sphi*v.y() ;
1065 zi = p.z() + sphi*v.z() ; 1477 zi = p.z() + sphi*v.z() ;
1066 rhoi = std::hypot(xi,yi); << 1478 rhoi2 = xi*xi + yi*yi ;
1067 it2 = zi*zi + (rhoi-fRtor)*(rhoi-fR << 1479 it2 = fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*sqrt(rhoi2)) ;
1068 1480
1069 if (it2 >= tolORMin2 && it2 <= tolO 1481 if (it2 >= tolORMin2 && it2 <= tolORMax2)
1070 { 1482 {
1071 // z and r intersections good - c 1483 // z and r intersections good - check intersecting
1072 // with correct half-plane 1484 // with correct half-plane
1073 // 1485 //
1074 if ((yi*cosCPhi-xi*sinCPhi)>=0) << 1486 if ((yi*cosCPhi-xi*sinCPhi)>=0) snxt=sphi;
1075 } 1487 }
1076 } 1488 }
1077 } 1489 }
1078 } 1490 }
1079 } 1491 }
1080 if(snxt < halfCarTolerance) { snxt = 0.0 ; << 1492 if(snxt < 0.5*kCarTolerance) snxt = 0.0 ;
1081 1493
>> 1494 #if DEBUGTORUS
>> 1495 G4cout << "G4Torus::DistanceToIn Final Value is "
>> 1496 << snxt << G4endl << G4endl;
>> 1497 #endif
>> 1498
1082 return snxt ; 1499 return snxt ;
1083 } 1500 }
1084 1501
1085 ///////////////////////////////////////////// 1502 /////////////////////////////////////////////////////////////////////////////
1086 // 1503 //
1087 // Calculate distance (<= actual) to closest 1504 // Calculate distance (<= actual) to closest surface of shape from outside
1088 // - Calculate distance to z, radial planes 1505 // - Calculate distance to z, radial planes
1089 // - Only to phi planes if outside phi extent 1506 // - Only to phi planes if outside phi extent
1090 // - Return 0 if point inside 1507 // - Return 0 if point inside
1091 1508
1092 G4double G4Torus::DistanceToIn( const G4Three 1509 G4double G4Torus::DistanceToIn( const G4ThreeVector& p ) const
1093 { 1510 {
1094 G4double safe=0.0, safe1, safe2 ; << 1511 G4double safe, safe1, safe2 ;
1095 G4double phiC, cosPhiC, sinPhiC, safePhi, e 1512 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1096 G4double rho, pt ; << 1513 G4double rho2, rho, pt2, pt ;
1097 << 1514
1098 rho = std::hypot(p.x(),p.y()); << 1515 #if DEBUGTORUS
1099 pt = std::hypot(p.z(),rho-fRtor); << 1516 G4cout << G4endl ;
>> 1517 #endif
>> 1518
>> 1519 rho2 = p.x()*p.x() + p.y()*p.y() ;
>> 1520 rho = sqrt(rho2) ;
>> 1521 pt2 = fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
>> 1522 pt = sqrt(pt2) ;
>> 1523
1100 safe1 = fRmin - pt ; 1524 safe1 = fRmin - pt ;
1101 safe2 = pt - fRmax ; 1525 safe2 = pt - fRmax ;
1102 1526
1103 if (safe1 > safe2) { safe = safe1; } << 1527 if (safe1 > safe2) safe = safe1;
1104 else { safe = safe2; } << 1528 else safe = safe2;
1105 1529
1106 if ( fDPhi < twopi && (rho != 0.0) ) << 1530 if ( fDPhi < 2.0*M_PI && rho )
1107 { 1531 {
1108 phiC = fSPhi + fDPhi*0.5 ; 1532 phiC = fSPhi + fDPhi*0.5 ;
1109 cosPhiC = std::cos(phiC) ; << 1533 cosPhiC = cos(phiC) ;
1110 sinPhiC = std::sin(phiC) ; << 1534 sinPhiC = sin(phiC) ;
1111 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC) 1535 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1112 1536
1113 if (cosPsi < std::cos(fDPhi*0.5) ) // Psi << 1537 if (cosPsi < cos(fDPhi*0.5) ) // Psi=angle from central phi to point
1114 { // Poi << 1538 { // Point lies outside phi range
1115 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 1539 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1116 { 1540 {
1117 safePhi = std::fabs(p.x()*std::sin(fS << 1541 safePhi = fabs(p.x()*sin(fSPhi) - p.y()*cos(fSPhi)) ;
1118 } 1542 }
1119 else 1543 else
1120 { 1544 {
1121 ePhi = fSPhi + fDPhi ; 1545 ePhi = fSPhi + fDPhi ;
1122 safePhi = std::fabs(p.x()*std::sin(eP << 1546 safePhi = fabs(p.x()*sin(ePhi) - p.y()*cos(ePhi)) ;
1123 } 1547 }
1124 if (safePhi > safe) { safe = safePhi ; << 1548 if (safePhi > safe) safe = safePhi ;
1125 } 1549 }
1126 } 1550 }
1127 if (safe < 0 ) { safe = 0 ; } << 1551 if (safe < 0 ) safe = 0 ;
1128 return safe; 1552 return safe;
1129 } 1553 }
1130 1554
1131 ///////////////////////////////////////////// 1555 ///////////////////////////////////////////////////////////////////////////
1132 // 1556 //
1133 // Calculate distance to surface of shape fro 1557 // Calculate distance to surface of shape from `inside', allowing for tolerance
1134 // - Only Calc rmax intersection if no valid 1558 // - Only Calc rmax intersection if no valid rmin intersection
1135 // 1559 //
>> 1560 // Problem: if the ray exit the torus from the surface, the only solution
>> 1561 // is epsilon (~ 0).
>> 1562 // Then this solution is eliminated by the loop (>= kRadTolerance) we have
>> 1563 // nothing. This results in 'invalid enum'
>> 1564 // solution: apply DistanceToIn() instead DistanceToOut() ?
1136 1565
1137 G4double G4Torus::DistanceToOut( const G4Thre 1566 G4double G4Torus::DistanceToOut( const G4ThreeVector& p,
1138 const G4Thre 1567 const G4ThreeVector& v,
1139 const G4bool 1568 const G4bool calcNorm,
1140 G4bool << 1569 G4bool *validNorm,
1141 G4Thre << 1570 G4ThreeVector *n ) const
1142 { 1571 {
1143 ESide side = kNull, sidephi = kNull ; 1572 ESide side = kNull, sidephi = kNull ;
1144 G4double snxt = kInfinity, sphi, sd[4] ; << 1573 G4double snxt = kInfinity, sphi, c[5], s[4] ;
1145 1574
1146 // Vars for phi intersection 1575 // Vars for phi intersection
1147 // 1576 //
1148 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, c 1577 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1149 G4double cPhi, sinCPhi, cosCPhi ; 1578 G4double cPhi, sinCPhi, cosCPhi ;
1150 G4double pDistS, compS, pDistE, compE, sphi 1579 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1151 1580
1152 // Radial Intersections Defenitions & Gener 1581 // Radial Intersections Defenitions & General Precals
>> 1582
>> 1583 // Define roots Si (generally real >=0) for intersection with
>> 1584 // torus (Ri = fRmax or fRmin) of ray p +S*v . General equation is :
>> 1585 // c[4]*S^4 + c[3]*S^3 +c[2]*S^2 + c[1]*S + c[0] = 0 .
>> 1586
>> 1587 #if DEBUGTORUS
>> 1588 G4cout << G4endl ;
>> 1589 #endif
1153 1590
1154 //////////////////////// new calculation // << 1591 G4int i,j,num ;
1155 << 1592 G4double Rtor2 = fRtor*fRtor, Rmax2 = fRmax*fRmax, Rmin2 = fRmin*fRmin ;
1156 #if 1 << 1593 G4double rho2 = p.x()*p.x()+p.y()*p.y();
1157 << 1594 G4double rho = sqrt(rho2) ;
1158 // This is the version with the calculation << 1595 G4double pt2 = fabs(rho2 + p.z()*p.z() + Rtor2 - 2*fRtor*rho) ;
1159 // To be done: Check the precision of this << 1596 G4double pt = sqrt(pt2) ;
1160 // If you want return always validNorm = fa <<
1161 <<
1162 <<
1163 G4double rho = std::hypot(p.x(),p.y()); <<
1164 G4double pt = hypot(p.z(),rho-fRtor); <<
1165 <<
1166 G4double pDotV = p.x()*v.x() + p.y()*v.y() 1597 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
1167 << 1598 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
1168 G4double tolRMax = fRmax - fRmaxTolerance ; << 1599
>> 1600 G4double tolRMax = fRmax - kRadTolerance*0.5 ;
1169 1601
1170 G4double vDotNmax = pDotV - fRtor*(v.x()* 1602 G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
1171 G4double pDotxyNmax = (1 - fRtor/rho) ; 1603 G4double pDotxyNmax = (1 - fRtor/rho) ;
1172 1604
1173 if( (pt*pt > tolRMax*tolRMax) && (vDotNmax << 1605 #if DEBUGTORUS
>> 1606 G4cout << "G4Torus::DistanceToOut " << p << ", " << v << G4endl ;
>> 1607 #endif
>> 1608
>> 1609 if( pt2 > tolRMax*tolRMax && vDotNmax >= 0 )
1174 { 1610 {
1175 // On tolerant boundary & heading outward 1611 // On tolerant boundary & heading outwards (or perpendicular to) outer
1176 // radial surface -> leaving immediately 1612 // radial surface -> leaving immediately with *n for really convex part
1177 // only 1613 // only
1178 << 1614
1179 if ( calcNorm && (pDotxyNmax >= -2.*fRmax << 1615 if (calcNorm && pDotxyNmax >= -kRadTolerance)
1180 { 1616 {
1181 *n = G4ThreeVector( p.x()*(1 - fRtor/rh 1617 *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
1182 p.y()*(1 - fRtor/rh 1618 p.y()*(1 - fRtor/rho)/pt,
1183 p.z()/pt 1619 p.z()/pt ) ;
1184 *validNorm = true ; 1620 *validNorm = true ;
1185 } 1621 }
1186 << 1622 #if DEBUGTORUS
>> 1623 G4cout << "G4Torus::DistanceToOut Leaving by Rmax immediately"
>> 1624 << G4endl ;
>> 1625 #endif
1187 return snxt = 0 ; // Leaving by Rmax imme 1626 return snxt = 0 ; // Leaving by Rmax immediately
1188 } 1627 }
1189 << 1628 else // intersection with Rmax torus
1190 snxt = SolveNumericJT(p,v,fRmax,false); << 1629 {
1191 side = kRMax ; << 1630 c[4] = 1.0 ;
>> 1631 c[3] = 4*pDotV ;
>> 1632 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - Rmax2 + 2*Rtor2*v.z()*v.z()) ;
>> 1633 c[1] = 4*(pDotV*(pRad2-Rtor2-Rmax2) + 2*Rtor2*p.z()*v.z()) ;
>> 1634 c[0] = pRad2*pRad2 - 2*pRad2*(Rtor2+Rmax2)
>> 1635 + 4*Rtor2*p.z()*p.z() + (Rtor2-Rmax2)*(Rtor2-Rmax2) ;
1192 1636
1193 // rmin << 1637 // Uncomment the line below to activate the analytical method
>> 1638 //
>> 1639 // num = SolveBiQuadratic(c,s) ;
1194 1640
1195 if ( fRmin != 0.0 ) << 1641 // Numerical root research
1196 { << 1642 s[0] = SolveNumeric( p, v, false);
1197 G4double tolRMin = fRmin + fRminTolerance << 1643 num = 1; // There is only one root.
>> 1644
>> 1645 #if DEBUGTORUS
>> 1646 G4cout << "G4Torus::DistanceToOut (" << __LINE__
>> 1647 << ") SolveNumeric : " << s[0] << G4endl ;
>> 1648 #endif
1198 1649
1199 if ( (pt*pt < tolRMin*tolRMin) && (vDotNm << 1650 if(num)
1200 { <<
1201 if (calcNorm) { *validNorm = false ; } <<
1202 return snxt = 0 ; <<
1203 } <<
1204 <<
1205 sd[0] = SolveNumericJT(p,v,fRmin,false); <<
1206 if ( sd[0] < snxt ) <<
1207 { 1651 {
1208 snxt = sd[0] ; << 1652 for(i=0;i<num;i++) // leave only >=kRadTolerance/2 roots
1209 side = kRMin ; << 1653 {
>> 1654 if( s[i] < kRadTolerance*0.5 )
>> 1655 {
>> 1656 for(j=i+1;j<num;j++) s[j-1] = s[j] ;
>> 1657 i-- ;
>> 1658 num-- ;
>> 1659 }
>> 1660 }
>> 1661 if(num)
>> 1662 {
>> 1663 snxt = s[0] ;
>> 1664 side = kRMax ;
>> 1665 }
1210 } 1666 }
1211 } <<
1212 1667
1213 #else << 1668 if (fRmin) // Possible Rmin intersection
1214 <<
1215 // this is the "conservative" version which <<
1216 // NOTE: using this version the unit test t <<
1217 <<
1218 snxt = SolveNumericJT(p,v,fRmax,false); <<
1219 side = kRMax ; <<
1220 <<
1221 if ( fRmin ) <<
1222 { <<
1223 sd[0] = SolveNumericJT(p,v,fRmin,false); <<
1224 if ( sd[0] < snxt ) <<
1225 { 1669 {
1226 snxt = sd[0] ; << 1670 G4double tolRMin = fRmin + kRadTolerance*0.5 ;
1227 side = kRMin ; <<
1228 } <<
1229 } <<
1230 <<
1231 if ( calcNorm && (snxt == 0.0) ) <<
1232 { <<
1233 *validNorm = false ; // Leaving solid, <<
1234 return snxt ; <<
1235 } <<
1236 1671
1237 #endif << 1672 // Leaving via Rmin
1238 << 1673 // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1239 if (fDPhi < twopi) // Phi Intersections << 1674 // - avoid sqrt for efficiency
>> 1675 //
>> 1676 if (pt2 < tolRMin*tolRMin && vDotNmax < 0 )
>> 1677 {
>> 1678 if (calcNorm) *validNorm = false ; // Concave surface of the torus
>> 1679 #if DEBUGTORUS
>> 1680 G4cout << "G4Torus::DistanceToOut Leaving by Rmin immediately"
>> 1681 << G4endl ;
>> 1682 #endif
>> 1683 return snxt = 0 ; // Leaving by Rmin immediately
>> 1684 }
>> 1685 else // intersection with Rmin torus
>> 1686 {
>> 1687 c[4] = 1.0 ;
>> 1688 c[3] = 4*pDotV ;
>> 1689 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - Rmin2 + 2*Rtor2*v.z()*v.z()) ;
>> 1690 c[1] = 4*(pDotV*(pRad2-Rtor2-Rmin2) + 2*Rtor2*p.z()*v.z()) ;
>> 1691 c[0] = pRad2*pRad2 - 2*pRad2*(Rtor2+Rmin2)
>> 1692 + 4*Rtor2*p.z()*p.z() + (Rtor2-Rmin2)*(Rtor2-Rmin2) ;
>> 1693
>> 1694 // Uncomment the line below to activate the analytical method
>> 1695 //
>> 1696 // num = SolveBiQuadratic(c,s) ;
>> 1697 // s[0] = s[0]; // We already take care of Rmin in SolveNumeric
>> 1698 num = 1;
>> 1699
>> 1700 if(num)
>> 1701 {
>> 1702 for(i=0;i<num;i++) // leave only >=kRadTolerance/2 roots
>> 1703 {
>> 1704 if(s[i] < kRadTolerance*0.5)
>> 1705 {
>> 1706 for(j=i+1;j<num;j++) s[j-1] = s[j] ;
>> 1707 i-- ;
>> 1708 num-- ;
>> 1709 }
>> 1710 }
>> 1711 if(num && s[0]<snxt)
>> 1712 {
>> 1713 snxt = s[0] ;
>> 1714 side = kRMin ;
>> 1715 }
>> 1716 }
>> 1717 }
>> 1718 } // if(Rmin)
>> 1719 }
>> 1720 if (fDPhi < 2.0*M_PI) // Phi Intersections
1240 { 1721 {
1241 sinSPhi = std::sin(fSPhi) ; << 1722 sinSPhi = sin(fSPhi) ;
1242 cosSPhi = std::cos(fSPhi) ; << 1723 cosSPhi = cos(fSPhi) ;
1243 ePhi = fSPhi + fDPhi ; 1724 ePhi = fSPhi + fDPhi ;
1244 sinEPhi = std::sin(ePhi) ; << 1725 sinEPhi = sin(ePhi) ;
1245 cosEPhi = std::cos(ePhi) ; << 1726 cosEPhi = cos(ePhi) ;
1246 cPhi = fSPhi + fDPhi*0.5 ; 1727 cPhi = fSPhi + fDPhi*0.5 ;
1247 sinCPhi = std::sin(cPhi) ; << 1728 sinCPhi = sin(cPhi) ;
1248 cosCPhi = std::cos(cPhi) ; << 1729 cosCPhi = 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 1730
1258 if ( (p.x() != 0.0) || (p.y() != 0.0) ) / << 1731 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1259 { 1732 {
1260 pDistS = p.x()*sinSPhi - p.y()*cosSPhi 1733 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
1261 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi 1734 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1262 1735
1263 // Comp -ve when in direction of outwar 1736 // Comp -ve when in direction of outwards normal
1264 // 1737 //
1265 compS = -sinSPhi*v.x() + cosSPhi*v.y( 1738 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1266 compE = sinEPhi*v.x() - cosEPhi*v.y() 1739 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1267 sidephi = kNull ; 1740 sidephi = kNull ;
1268 << 1741
1269 if( ( (fDPhi <= pi) && ( (pDistS <= hal << 1742 if (pDistS <= 0 && pDistE <= 0 )
1270 && (pDistE <= hal <<
1271 || ( (fDPhi > pi) && ((pDistS <= hal <<
1272 || (pDistE <= ha <<
1273 { 1743 {
1274 // Inside both phi *full* planes 1744 // Inside both phi *full* planes
1275 1745
1276 if ( compS < 0 ) << 1746 if (compS<0)
>> 1747 {
>> 1748 sphi=pDistS/compS;
>> 1749 xi=p.x()+sphi*v.x();
>> 1750 yi=p.y()+sphi*v.y();
>> 1751
>> 1752 // Check intersecting with correct half-plane
>> 1753 // (if not -> no intersect)
>> 1754 //
>> 1755 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 1756 sphi=kInfinity;
>> 1757 else
>> 1758 {
>> 1759 sidephi=kSPhi;
>> 1760 if (pDistS>-kCarTolerance*0.5)
>> 1761 sphi=0;
>> 1762 // Leave by sphi immediately
>> 1763 }
>> 1764 }
>> 1765 else sphi=kInfinity;
>> 1766
>> 1767 if (compE<0)
1277 { 1768 {
1278 sphi = pDistS/compS ; << 1769 sphi2=pDistE/compE;
1279 << 1770
1280 if (sphi >= -halfCarTolerance) << 1771 // Only check further if < starting phi intersection
1281 { << 1772 //
1282 xi = p.x() + sphi*v.x() ; << 1773 if (sphi2<sphi)
1283 yi = p.y() + sphi*v.y() ; << 1774 {
1284 << 1775 xi=p.x()+sphi2*v.x();
>> 1776 yi=p.y()+sphi2*v.y();
>> 1777
1285 // Check intersecting with correc 1778 // Check intersecting with correct half-plane
1286 // (if not -> no intersect) << 1779 //
1287 // << 1780 if ((yi*cosCPhi-xi*sinCPhi)>=0)
1288 if ( (std::fabs(xi)<=kCarToleranc <<
1289 && (std::fabs(yi)<=kCarToleranc <<
1290 { 1781 {
1291 sidephi = kSPhi; << 1782 // Leaving via ending phi
1292 if ( ((fSPhi-halfAngTolerance)< << 1783 //
1293 && ((ePhi+halfAngTolerance)>= << 1784 sidephi=kEPhi;
>> 1785 if (pDistE<=-kCarTolerance*0.5)
>> 1786 {
>> 1787 sphi=sphi2;
>> 1788 }
>> 1789 else
1294 { 1790 {
1295 sphi = kInfinity; << 1791 sphi=0;
1296 } 1792 }
1297 } 1793 }
1298 else if ( yi*cosCPhi-xi*sinCPhi > << 1794 }
>> 1795 }
>> 1796 }
>> 1797 else if (pDistS>=0&&pDistE>=0)
>> 1798 {
>> 1799 // Outside both *full* phi planes
>> 1800
>> 1801 if (pDistS <= pDistE)
>> 1802 {
>> 1803 sidephi = kSPhi ;
>> 1804 }
>> 1805 else
>> 1806 {
>> 1807 sidephi = kEPhi ;
>> 1808 }
>> 1809 if (fDPhi>M_PI)
>> 1810 {
>> 1811 if (compS<0&&compE<0) sphi=0;
>> 1812 else sphi=kInfinity;
>> 1813 }
>> 1814 else
>> 1815 {
>> 1816 // if towards both >=0 then once inside (after error)
>> 1817 // will remain inside
>> 1818 //
>> 1819 if (compS>=0&&compE>=0)
>> 1820 {
>> 1821 sphi=kInfinity;
>> 1822 }
>> 1823 else
>> 1824 {
>> 1825 sphi=0;
>> 1826 }
>> 1827 }
>> 1828 }
>> 1829 else if (pDistS>0&&pDistE<0)
>> 1830 {
>> 1831 // Outside full starting plane, inside full ending plane
>> 1832
>> 1833 if (fDPhi>M_PI)
>> 1834 {
>> 1835 if (compE<0)
>> 1836 {
>> 1837 sphi=pDistE/compE;
>> 1838 xi=p.x()+sphi*v.x();
>> 1839 yi=p.y()+sphi*v.y();
>> 1840
>> 1841 // Check intersection in correct half-plane (if not -> not leaving phi extent)
>> 1842 //
>> 1843 if ((yi*cosCPhi-xi*sinCPhi)<=0)
1299 { 1844 {
1300 sphi = kInfinity ; << 1845 sphi=kInfinity;
1301 } 1846 }
1302 else 1847 else
1303 { 1848 {
1304 sidephi = kSPhi ; << 1849 // Leaving via Ending phi
1305 } << 1850 //
>> 1851 sidephi = kEPhi ;
>> 1852 if (pDistE>-kCarTolerance*0.5)
>> 1853 sphi=0;
>> 1854 }
1306 } 1855 }
1307 else 1856 else
1308 { 1857 {
1309 sphi = kInfinity ; << 1858 sphi=kInfinity;
1310 } 1859 }
1311 } 1860 }
1312 else 1861 else
1313 { 1862 {
1314 sphi = kInfinity ; << 1863 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 { 1864 {
1325 xi = p.x() + sphi2*v.x() ; << 1865 if (compE<0)
1326 yi = p.y() + sphi2*v.y() ; <<
1327 <<
1328 if ( (std::fabs(xi)<=kCarToleranc <<
1329 && (std::fabs(yi)<=kCarToleranc <<
1330 { 1866 {
1331 // Leaving via ending phi << 1867 sphi=pDistE/compE;
>> 1868 xi=p.x()+sphi*v.x();
>> 1869 yi=p.y()+sphi*v.y();
>> 1870
>> 1871 // Check intersection in correct half-plane
>> 1872 // (if not -> remain in extent)
1332 // 1873 //
1333 if( (fSPhi-halfAngTolerance > v << 1874 if ((yi*cosCPhi-xi*sinCPhi)<=0)
1334 || (ePhi+halfAngTolerance < <<
1335 { 1875 {
1336 sidephi = kEPhi ; << 1876 sphi=kInfinity;
1337 sphi = sphi2; <<
1338 } 1877 }
1339 } << 1878 else
1340 else // Check intersecting wit <<
1341 { <<
1342 if ( (yi*cosCPhi-xi*sinCPhi) >= <<
1343 { 1879 {
1344 // Leaving via ending phi << 1880 // otherwise leaving via Ending phi
1345 // 1881 //
1346 sidephi = kEPhi ; << 1882 sidephi=kEPhi;
1347 sphi = sphi2; <<
1348 <<
1349 } 1883 }
1350 } 1884 }
>> 1885 else sphi=kInfinity;
>> 1886 }
>> 1887 else
>> 1888 {
>> 1889 // leaving immediately by starting phi
>> 1890 //
>> 1891 sidephi=kSPhi;
>> 1892 sphi=0;
1351 } 1893 }
1352 } 1894 }
1353 } 1895 }
1354 else 1896 else
1355 { 1897 {
1356 sphi = kInfinity ; << 1898 // Must be pDistS<0&&pDistE>0
>> 1899 // Inside full starting plane, outside full ending plane
>> 1900
>> 1901 if (fDPhi>M_PI)
>> 1902 {
>> 1903 if (compS<0)
>> 1904 {
>> 1905 sphi=pDistS/compS;
>> 1906 xi=p.x()+sphi*v.x();
>> 1907 yi=p.y()+sphi*v.y();
>> 1908
>> 1909 // Check intersection in correct half-plane
>> 1910 // (if not -> not leaving phi extent)
>> 1911 //
>> 1912 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 1913 {
>> 1914 sphi=kInfinity;
>> 1915 }
>> 1916 else
>> 1917 {
>> 1918 // Leaving via Starting phi
>> 1919 //
>> 1920 sidephi = kSPhi ;
>> 1921 if (pDistS>-kCarTolerance*0.5)
>> 1922 sphi=0;
>> 1923 }
>> 1924 }
>> 1925 else
>> 1926 {
>> 1927 sphi=kInfinity;
>> 1928 }
>> 1929 }
>> 1930 else
>> 1931 {
>> 1932 if (compE>=0)
>> 1933 {
>> 1934 if (compS<0)
>> 1935 {
>> 1936 sphi=pDistS/compS;
>> 1937 xi=p.x()+sphi*v.x();
>> 1938 yi=p.y()+sphi*v.y();
>> 1939
>> 1940 // Check intersection in correct half-plane
>> 1941 // (if not -> remain in extent)
>> 1942 //
>> 1943 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 1944 {
>> 1945 sphi=kInfinity;
>> 1946 }
>> 1947 else
>> 1948 {
>> 1949 // otherwise leaving via Starting phi
>> 1950 //
>> 1951 sidephi=kSPhi;
>> 1952 }
>> 1953 }
>> 1954 else
>> 1955 {
>> 1956 sphi=kInfinity;
>> 1957 }
>> 1958 }
>> 1959 else
>> 1960 {
>> 1961 // leaving immediately by ending
>> 1962 //
>> 1963 sidephi=kEPhi;
>> 1964 sphi=0;
>> 1965 }
>> 1966 }
1357 } 1967 }
1358 } << 1968 }
1359 else 1969 else
1360 { 1970 {
1361 // On z axis + travel not || to z axis 1971 // On z axis + travel not || to z axis -> if phi of vector direction
1362 // within phi of shape, Step limited by 1972 // within phi of shape, Step limited by rmax, else Step =0
1363 1973
1364 vphi = std::atan2(v.y(),v.x()); << 1974 vphi=atan2(v.y(),v.x());
1365 << 1975 if (fSPhi<vphi&&vphi<fSPhi+fDPhi)
1366 if ( ( fSPhi-halfAngTolerance <= vphi ) <<
1367 ( vphi <= ( ePhi+halfAngTolerance <<
1368 { 1976 {
1369 sphi = kInfinity; << 1977 sphi=kInfinity;
1370 } 1978 }
1371 else 1979 else
1372 { 1980 {
1373 sidephi = kSPhi ; // arbitrary 1981 sidephi = kSPhi ; // arbitrary
1374 sphi=0; 1982 sphi=0;
1375 } 1983 }
1376 } 1984 }
1377 1985
1378 // Order intersections 1986 // Order intersections
1379 1987
1380 if (sphi<snxt) 1988 if (sphi<snxt)
1381 { 1989 {
1382 snxt=sphi; 1990 snxt=sphi;
1383 side=sidephi; 1991 side=sidephi;
1384 } << 1992 }
1385 } 1993 }
>> 1994 G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
1386 1995
1387 G4double rhoi,it,iDotxyNmax ; <<
1388 // Note: by numerical computation we know w 1996 // Note: by numerical computation we know where the ray hits the torus
1389 // So I propose to return the side where th 1997 // So I propose to return the side where the ray hits
1390 1998
1391 if (calcNorm) 1999 if (calcNorm)
1392 { 2000 {
1393 switch(side) 2001 switch(side)
1394 { 2002 {
1395 case kRMax: // n is 2003 case kRMax: // n is unit vector
>> 2004 #if DEBUGTORUS
>> 2005 G4cout << "G4Torus::DistanceToOut Side is RMax" << G4endl ;
>> 2006 #endif
1396 xi = p.x() + snxt*v.x() ; 2007 xi = p.x() + snxt*v.x() ;
1397 yi = p.y() + snxt*v.y() ; << 2008 yi =p.y() + snxt*v.y() ;
1398 zi = p.z() + snxt*v.z() ; 2009 zi = p.z() + snxt*v.z() ;
1399 rhoi = std::hypot(xi,yi); << 2010 rhoi2 = xi*xi + yi*yi ;
1400 it = hypot(zi,rhoi-fRtor); << 2011 rhoi = sqrt(rhoi2) ;
1401 << 2012 it2 = fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ;
>> 2013 it = sqrt(it2) ;
1402 iDotxyNmax = (1-fRtor/rhoi) ; 2014 iDotxyNmax = (1-fRtor/rhoi) ;
1403 if(iDotxyNmax >= -2.*fRmaxTolerance) << 2015 if(iDotxyNmax >= -kRadTolerance) // really convex part of Rmax
1404 { 2016 {
1405 *n = G4ThreeVector( xi*(1-fRtor/rho 2017 *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
1406 yi*(1-fRtor/rho 2018 yi*(1-fRtor/rhoi)/it,
1407 zi/it 2019 zi/it ) ;
1408 *validNorm = true ; 2020 *validNorm = true ;
1409 } 2021 }
1410 else << 2022 else *validNorm = false ; // concave-convex part of Rmax
1411 { <<
1412 *validNorm = false ; // concave-con <<
1413 } <<
1414 break ; 2023 break ;
1415 2024
1416 case kRMin: 2025 case kRMin:
>> 2026 #if DEBUGTORUS
>> 2027 G4cout << "G4Torus::DistanceToOut Side is RMin" << G4endl ;
>> 2028 #endif
1417 *validNorm = false ; // Rmin is conc 2029 *validNorm = false ; // Rmin is concave or concave-convex
1418 break; 2030 break;
1419 2031
1420 case kSPhi: 2032 case kSPhi:
1421 if (fDPhi <= pi ) << 2033 #if DEBUGTORUS
>> 2034 G4cout << "G4Torus::DistanceToOut Side is SPhi" << G4endl ;
>> 2035 #endif
>> 2036 if (fDPhi <= M_PI )
1422 { 2037 {
1423 *n=G4ThreeVector(std::sin(fSPhi),-s << 2038 *n=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0);
1424 *validNorm=true; 2039 *validNorm=true;
1425 } 2040 }
1426 else << 2041 else *validNorm = false ;
1427 { <<
1428 *validNorm = false ; <<
1429 } <<
1430 break ; 2042 break ;
1431 2043
1432 case kEPhi: 2044 case kEPhi:
1433 if (fDPhi <= pi) << 2045 #if DEBUGTORUS
>> 2046 G4cout << "G4Torus::DistanceToOut Side is EPhi" << G4endl ;
>> 2047 #endif
>> 2048 if (fDPhi <= M_PI)
1434 { 2049 {
1435 *n=G4ThreeVector(-std::sin(fSPhi+fD << 2050 *n=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0);
1436 *validNorm=true; 2051 *validNorm=true;
1437 } 2052 }
1438 else << 2053 else *validNorm = false ;
1439 { <<
1440 *validNorm = false ; <<
1441 } <<
1442 break; 2054 break;
1443 2055
1444 default: 2056 default:
1445 2057
1446 // It seems we go here from time to t 2058 // It seems we go here from time to time ...
>> 2059 //
>> 2060 // G4cout << "Side is " << side << G4endl ;
>> 2061 // G4cout << "Valid ESide are :" << kNull << " "
>> 2062 // << kRMin << " " << kRMax
>> 2063 // << " " << kSPhi << " " << kEPhi << G4endl;
1447 2064
>> 2065 G4cout.precision(16);
1448 G4cout << G4endl; 2066 G4cout << G4endl;
1449 DumpInfo(); 2067 DumpInfo();
1450 std::ostringstream message; << 2068 G4cout << "Position:" << G4endl << G4endl;
1451 G4long oldprc = message.precision(16) << 2069 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
1452 message << "Undefined side for valid << 2070 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
1453 << G4endl << 2071 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
1454 << "Position:" << G4endl << << 2072 G4cout << "Direction:" << G4endl << G4endl;
1455 << "p.x() = " << p.x()/mm < << 2073 G4cout << "v.x() = " << v.x() << G4endl;
1456 << "p.y() = " << p.y()/mm < << 2074 G4cout << "v.y() = " << v.y() << G4endl;
1457 << "p.z() = " << p.z()/mm < << 2075 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
1458 << "Direction:" << G4endl << << 2076 G4cout << "Proposed distance :" << G4endl << G4endl;
1459 << "v.x() = " << v.x() << G << 2077 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl;
1460 << "v.y() = " << v.y() << G << 2078 G4Exception("G4Torus::DistanceToOut() - Invalid enum");
1461 << "v.z() = " << v.z() << G <<
1462 << "Proposed distance :" << G <<
1463 << "snxt = " << snxt/mm << " <<
1464 message.precision(oldprc); <<
1465 G4Exception("G4Torus::DistanceToOut(p <<
1466 "GeomSolids1002",JustWarn <<
1467 break; 2079 break;
1468 } 2080 }
1469 } 2081 }
1470 if ( snxt<halfCarTolerance ) { snxt=0 ; } << 2082
>> 2083 #if DEBUGTORUS
>> 2084 G4cout << "G4Torus::DistanceToOut Final Value is "
>> 2085 << snxt << G4endl << G4endl;
>> 2086 #endif
1471 2087
1472 return snxt; 2088 return snxt;
1473 } 2089 }
1474 2090
1475 ///////////////////////////////////////////// 2091 /////////////////////////////////////////////////////////////////////////
1476 // 2092 //
1477 // Calculate distance (<=actual) to closest s 2093 // Calculate distance (<=actual) to closest surface of shape from inside
1478 2094
1479 G4double G4Torus::DistanceToOut( const G4Thre 2095 G4double G4Torus::DistanceToOut( const G4ThreeVector& p ) const
1480 { 2096 {
1481 G4double safe=0.0,safeR1,safeR2; << 2097 G4double safe,safeR1,safeR2;
1482 G4double rho,pt ; << 2098 G4double rho2,rho,pt2,pt ;
1483 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi; 2099 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1484 << 2100 rho2 = p.x()*p.x() + p.y()*p.y() ;
1485 rho = std::hypot(p.x(),p.y()); << 2101 rho = sqrt(rho2) ;
1486 pt = std::hypot(p.z(),rho-fRtor); << 2102 pt2 = fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
1487 << 2103 pt = sqrt(pt2) ;
>> 2104
1488 #ifdef G4CSGDEBUG 2105 #ifdef G4CSGDEBUG
1489 if( Inside(p) == kOutside ) 2106 if( Inside(p) == kOutside )
1490 { 2107 {
1491 G4long oldprc = G4cout.precision(16) ; << 2108 G4cout.precision(16) ;
1492 G4cout << G4endl ; 2109 G4cout << G4endl ;
1493 DumpInfo(); 2110 DumpInfo();
1494 G4cout << "Position:" << G4endl << G4en 2111 G4cout << "Position:" << G4endl << G4endl ;
1495 G4cout << "p.x() = " << p.x()/mm << " 2112 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1496 G4cout << "p.y() = " << p.y()/mm << " 2113 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1497 G4cout << "p.z() = " << p.z()/mm << " 2114 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1498 G4cout.precision(oldprc); << 2115 G4cout << "G4Torus::DistanceToOut(p) - point p is outside ?!" << G4endl ;
1499 G4Exception("G4Torus::DistanceToOut(p)", << 2116 G4cerr << "G4Torus::DistanceToOut(p) - point p is outside ?!" << G4endl ;
1500 JustWarning, "Point p is out <<
1501 } 2117 }
1502 #endif 2118 #endif
>> 2119 #if DEBUGTORUS
>> 2120 G4cout << G4endl ;
>> 2121 #endif
1503 2122
1504 if (fRmin != 0.0) << 2123 if (fRmin)
1505 { 2124 {
1506 safeR1 = pt - fRmin ; 2125 safeR1 = pt - fRmin ;
1507 safeR2 = fRmax - pt ; 2126 safeR2 = fRmax - pt ;
1508 2127
1509 if (safeR1 < safeR2) { safe = safeR1 ; } << 2128 if (safeR1 < safeR2) safe = safeR1 ;
1510 else { safe = safeR2 ; } << 2129 else safe = safeR2 ;
1511 } 2130 }
1512 else << 2131 else safe = fRmax - pt ;
1513 { <<
1514 safe = fRmax - pt ; <<
1515 } <<
1516 2132
1517 // Check if phi divided, Calc distances clo << 2133 // Check if phi divided, Calc distances closest phi plane
1518 // << 2134
1519 if (fDPhi < twopi) // Above/below central p << 2135 if (fDPhi<2.0*M_PI) // Above/below central phi of Torus?
1520 { 2136 {
1521 phiC = fSPhi + fDPhi*0.5 ; 2137 phiC = fSPhi + fDPhi*0.5 ;
1522 cosPhiC = std::cos(phiC) ; << 2138 cosPhiC = cos(phiC) ;
1523 sinPhiC = std::sin(phiC) ; << 2139 sinPhiC = sin(phiC) ;
1524 2140
1525 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) 2141 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1526 { 2142 {
1527 safePhi = -(p.x()*std::sin(fSPhi) - p.y << 2143 safePhi = -(p.x()*sin(fSPhi) - p.y()*cos(fSPhi)) ;
1528 } 2144 }
1529 else 2145 else
1530 { 2146 {
1531 ePhi = fSPhi + fDPhi ; 2147 ePhi = fSPhi + fDPhi ;
1532 safePhi = (p.x()*std::sin(ePhi) - p.y() << 2148 safePhi = (p.x()*sin(ePhi) - p.y()*cos(ePhi)) ;
1533 } 2149 }
1534 if (safePhi < safe) { safe = safePhi ; } << 2150 if (safePhi < safe) safe = safePhi ;
1535 } 2151 }
1536 if (safe < 0) { safe = 0 ; } << 2152 if (safe < 0) safe = 0 ;
1537 return safe ; 2153 return safe ;
1538 } 2154 }
1539 2155
1540 ///////////////////////////////////////////// << 2156 /////////////////////////////////////////////////////////////////////////////
1541 // 2157 //
1542 // Stream object contents to an output stream << 2158 // Create a List containing the transformed vertices
>> 2159 // Ordering [0-3] -fRtor cross section
>> 2160 // [4-7] +fRtor cross section such that [0] is below [4],
>> 2161 // [1] below [5] etc.
>> 2162 // Note:
>> 2163 // Caller has deletion resposibility
>> 2164 // Potential improvement: For last slice, use actual ending angle
>> 2165 // to avoid rounding error problems.
1543 2166
1544 G4GeometryType G4Torus::GetEntityType() const << 2167 G4ThreeVectorList*
>> 2168 G4Torus::CreateRotatedVertices( const G4AffineTransform& pTransform,
>> 2169 G4int& noPolygonVertices ) const
1545 { 2170 {
1546 return {"G4Torus"}; << 2171 G4ThreeVectorList *vertices;
>> 2172 G4ThreeVector vertex0,vertex1,vertex2,vertex3;
>> 2173 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
>> 2174 G4double rMaxX,rMaxY,rMinX,rMinY;
>> 2175 G4int crossSection,noCrossSections;
>> 2176
>> 2177 // Compute no of cross-sections necessary to mesh tube
>> 2178
>> 2179 noCrossSections = G4int (fDPhi/kMeshAngleDefault) + 1 ;
>> 2180
>> 2181 if (noCrossSections < kMinMeshSections)
>> 2182 {
>> 2183 noCrossSections = kMinMeshSections ;
>> 2184 }
>> 2185 else if (noCrossSections>kMaxMeshSections)
>> 2186 {
>> 2187 noCrossSections=kMaxMeshSections;
>> 2188 }
>> 2189 meshAngle = fDPhi/(noCrossSections - 1) ;
>> 2190 meshRMax = (fRtor + fRmax)/cos(meshAngle*0.5) ;
>> 2191
>> 2192 // If complete in phi, set start angle such that mesh will be at fRmax
>> 2193 // on the x axis. Will give better extent calculations when not rotated.
>> 2194
>> 2195 if ( fDPhi == M_PI*2.0 && fSPhi == 0 )
>> 2196 {
>> 2197 sAngle = -meshAngle*0.5 ;
>> 2198 }
>> 2199 else
>> 2200 {
>> 2201 sAngle = fSPhi ;
>> 2202 }
>> 2203 vertices = new G4ThreeVectorList();
>> 2204 vertices->reserve(noCrossSections*4) ;
>> 2205
>> 2206 if (vertices)
>> 2207 {
>> 2208 for (crossSection=0;crossSection<noCrossSections;crossSection++)
>> 2209 {
>> 2210 // Compute coordinates of cross section at section crossSection
>> 2211
>> 2212 crossAngle=sAngle+crossSection*meshAngle;
>> 2213 cosCrossAngle=cos(crossAngle);
>> 2214 sinCrossAngle=sin(crossAngle);
>> 2215
>> 2216 rMaxX=meshRMax*cosCrossAngle;
>> 2217 rMaxY=meshRMax*sinCrossAngle;
>> 2218 rMinX=(fRtor-fRmax)*cosCrossAngle;
>> 2219 rMinY=(fRtor-fRmax)*sinCrossAngle;
>> 2220 vertex0=G4ThreeVector(rMinX,rMinY,-fRmax);
>> 2221 vertex1=G4ThreeVector(rMaxX,rMaxY,-fRmax);
>> 2222 vertex2=G4ThreeVector(rMaxX,rMaxY,+fRmax);
>> 2223 vertex3=G4ThreeVector(rMinX,rMinY,+fRmax);
>> 2224
>> 2225 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 2226 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 2227 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 2228 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 2229 }
>> 2230 noPolygonVertices = 4 ;
>> 2231 }
>> 2232 else
>> 2233 {
>> 2234 DumpInfo();
>> 2235 G4Exception("G4Torus::CreateRotatedVertices() - Out of memory !");
>> 2236 }
>> 2237 return vertices;
1547 } 2238 }
1548 2239
1549 ///////////////////////////////////////////// 2240 //////////////////////////////////////////////////////////////////////////
1550 // 2241 //
1551 // Make a clone of the object << 2242 // Stream object contents to an output stream
1552 // << 2243
1553 G4VSolid* G4Torus::Clone() const << 2244 G4GeometryType G4Torus::GetEntityType() const
1554 { 2245 {
1555 return new G4Torus(*this); << 2246 return G4String("G4Torus");
1556 } 2247 }
1557 2248
1558 ///////////////////////////////////////////// 2249 //////////////////////////////////////////////////////////////////////////
1559 // 2250 //
1560 // Stream object contents to an output stream 2251 // Stream object contents to an output stream
1561 2252
1562 std::ostream& G4Torus::StreamInfo( std::ostre << 2253 G4std::ostream& G4Torus::StreamInfo( G4std::ostream& os ) const
1563 { 2254 {
1564 G4long oldprc = os.precision(16); <<
1565 os << "------------------------------------ 2255 os << "-----------------------------------------------------------\n"
1566 << " *** Dump for solid - " << GetNam 2256 << " *** Dump for solid - " << GetName() << " ***\n"
1567 << " ================================ 2257 << " ===================================================\n"
1568 << " Solid type: G4Torus\n" 2258 << " Solid type: G4Torus\n"
1569 << " Parameters: \n" 2259 << " Parameters: \n"
1570 << " inner radius: " << fRmin/mm << " 2260 << " inner radius: " << fRmin/mm << " mm \n"
1571 << " outer radius: " << fRmax/mm << " 2261 << " outer radius: " << fRmax/mm << " mm \n"
1572 << " swept radius: " << fRtor/mm << " 2262 << " swept radius: " << fRtor/mm << " mm \n"
1573 << " starting phi: " << fSPhi/degree 2263 << " starting phi: " << fSPhi/degree << " degrees \n"
1574 << " delta phi : " << fDPhi/degree 2264 << " delta phi : " << fDPhi/degree << " degrees \n"
1575 << "------------------------------------ 2265 << "-----------------------------------------------------------\n";
1576 os.precision(oldprc); <<
1577 2266
1578 return os; 2267 return os;
1579 } 2268 }
1580 2269
1581 ///////////////////////////////////////////// << 2270 ///////////////////////////////////////////////////////////////////////
1582 // 2271 //
1583 // GetPointOnSurface << 2272 // Visualisation Functions
1584 2273
1585 G4ThreeVector G4Torus::GetPointOnSurface() co << 2274 void G4Torus::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1586 { 2275 {
1587 G4double cosu, sinu,cosv, sinv, aOut, aIn, << 2276 scene.AddThis (*this);
1588 << 2277 }
1589 phi = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi << 2278
1590 theta = G4RandFlat::shoot(0.,twopi); << 2279 G4Polyhedron* G4Torus::CreatePolyhedron () const
>> 2280 {
>> 2281 return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi);
>> 2282 }
>> 2283
>> 2284 G4NURBS* G4Torus::CreateNURBS () const
>> 2285 {
>> 2286 G4NURBS* pNURBS;
>> 2287 if (fRmin != 0)
>> 2288 {
>> 2289 if (fDPhi >= 2.0 * M_PI)
>> 2290 {
>> 2291 pNURBS = new G4NURBStube(fRmin, fRmax, fRtor);
>> 2292 }
>> 2293 else
>> 2294 {
>> 2295 pNURBS = new G4NURBStubesector(fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi);
>> 2296 }
>> 2297 }
>> 2298 else
>> 2299 {
>> 2300 if (fDPhi >= 2.0 * M_PI)
>> 2301 {
>> 2302 pNURBS = new G4NURBScylinder (fRmax, fRtor);
>> 2303 }
>> 2304 else
>> 2305 {
>> 2306 const G4double epsilon = 1.e-4; // Cylinder sector not yet available!
>> 2307 pNURBS = new G4NURBStubesector (epsilon, fRmax, fRtor,
>> 2308 fSPhi, fSPhi + fDPhi);
>> 2309 }
>> 2310 }
>> 2311 return pNURBS;
>> 2312 }
>> 2313
>> 2314 // --------------------------------------------------------------------
>> 2315 // Numerical algorithms implementations
>> 2316 //
>> 2317 // Important : the precision could be tuned by TORUSPRECISION
>> 2318
>> 2319 #define TORUSPRECISION 1.0 // or whatever you want for precision
>> 2320 // (it is TorusEquation related)
>> 2321 #define HOLEBVM 0
>> 2322 #define NBPOINT 6
>> 2323 // --------------------------------------------------------------------
>> 2324
>> 2325 ///////////////////////////////////////////////////////////////////////
>> 2326 //
>> 2327 // Torus implementation with Newton Method and Bounding volume
>> 2328 //
>> 2329 // For speed issue, we lose time *only* when intersecting the BVM
>> 2330 // and SafeNewton when it is called.
>> 2331
>> 2332 G4double G4Torus::SolveNumeric( const G4ThreeVector& p,
>> 2333 const G4ThreeVector& v,
>> 2334 G4bool IsDistanceToIn ) const
>> 2335 {
>> 2336 // This method is a front-end to the numerical computation of roots
>> 2337 // In fact this computation takes care only of a perfect Torus
>> 2338 // So here we add Phi section/Tolerance/Rinterior
>> 2339
>> 2340 /*** SolveNumeric ***/
>> 2341
>> 2342 /** Conditions **/
>> 2343 // - if original point inside check for interior torus before
>> 2344 // - on surface it depends on the direction
>> 2345 // - the intersection point must be between fSPhi and fSPhi+fDPhi
>> 2346 // - If on the surface it depends on DistanceToOut or DistanceToIn :
>> 2347 // a ray from the surface to In called with DistanceToIn return 0.0
>> 2348 // and with DistanceToOut return the second intersection point
>> 2349
>> 2350 G4double lambda = 0;
>> 2351 G4double Value = TorusEquation(p.x(),p.y(),p.z(),GetRtor(),GetRmax());
>> 2352 EInside inside ;
>> 2353
>> 2354 #if DEBUGTORUS
>> 2355 G4cout << "G4Torus::SolveNumeric " << p << ", " << v << G4endl ;
>> 2356 G4cout << "G4Torus::SolveNumeric Value = " << Value << G4endl;
>> 2357 #endif
1591 2358
1592 cosu = std::cos(phi); sinu = std::sin( << 2359 if (Value < -TORUSPRECISION)
1593 cosv = std::cos(theta); sinv = std::sin( << 2360 {
>> 2361 inside = kInside ;
>> 2362 }
>> 2363 else
>> 2364 {
>> 2365 if (Value > TORUSPRECISION)
>> 2366 {
>> 2367 inside = kOutside;
>> 2368 }
>> 2369 else
>> 2370 {
>> 2371 inside = kSurface;
>> 2372 }
>> 2373 }
>> 2374
>> 2375 switch (inside)
>> 2376 {
>> 2377 case kInside:
1594 2378
1595 // compute the areas << 2379 #if DEBUGTORUS
>> 2380 G4cout << "G4Torus::SolveNumeric Point is Inside Rmax Torus "
>> 2381 << " Rtor = " << GetRtor()
>> 2382 << " Rmax = " << GetRmax() << G4endl ;
>> 2383 #endif
>> 2384
>> 2385 if (fabs(GetRmin()) > POLEPSILON)
>> 2386 {
>> 2387
>> 2388 #if DEBUGTORUS
>> 2389 G4cout << "G4Torus::SolveNumeric Testing interior torus .."
>> 2390 << G4endl ;
>> 2391 #endif
>> 2392 lambda = DistanceToTorus(p.x(),p.y(),p.z(),v.x(),v.y(),v.z(),
>> 2393 GetRtor(),GetRmin()); //Interior torus
>> 2394
>> 2395 #if DEBUGTORUS
>> 2396 G4cout << "G4Torus::SolveNumeric lambda to interior torus ="
>> 2397 << lambda << G4endl ;
>> 2398 G4cout << "G4Torus::SolveNumeric Tolerance is "
>> 2399 << kCarTolerance << G4endl ;
>> 2400 #endif
1596 2401
1597 aOut = (fDPhi)*twopi*fRtor*fRmax; << 2402 // Now check if on surface from interior torus
1598 aIn = (fDPhi)*twopi*fRtor*fRmin; << 2403 //
1599 aSide = pi*(fRmax*fRmax-fRmin*fRmin); << 2404 // PROBLEM: This may be a problem of precision
>> 2405 // if we are near kCarTolerance ...
>> 2406 //
>> 2407 if (fabs(lambda) < kCarTolerance)
>> 2408 {
>> 2409 G4double Lx,Ly,Lz;
>> 2410 G4double scal;
>> 2411 #if DEBUGTORUS
>> 2412 G4cout << "G4Torus::SolveNumeric"
>> 2413 << " In fact on the Surface of Rmin torus"
>> 2414 << G4endl ;
>> 2415 #endif
>> 2416
>> 2417 // Compute Surface point
>> 2418 //
>> 2419 Lx = p.x() + lambda*v.x();
>> 2420 Ly = p.y() + lambda*v.y();
>> 2421 Lz = p.z() + lambda*v.z();
>> 2422
>> 2423 // Scalar product
>> 2424 //
>> 2425 scal = v.x()*TorusDerivativeX(Lx,Ly,Lz,GetRtor(),GetRmin());
>> 2426 scal += v.y()*TorusDerivativeY(Lx,Ly,Lz,GetRtor(),GetRmin());
>> 2427 scal += v.z()*TorusDerivativeZ(Lx,Ly,Lz,GetRtor(),GetRmin());
>> 2428
>> 2429 // if entering and if it is DistToIn it is 0.0,
>> 2430 // but in fact it is the opposite because it is the interior torus
>> 2431 // beware that this could be DistanceToOut
>> 2432
>> 2433 if ( (IsDistanceToIn == true) && (scal > 0.0) )
>> 2434 {
>> 2435 #if DEBUGTORUS
>> 2436 G4cout << "G4Torus::SolveNumeric"
>> 2437 << " Entering Surface from Rmin Torus Gradient: "
>> 2438 << scal << G4endl ;
>> 2439 #endif
>> 2440
>> 2441 // DistanceToIn returns 0.0
>> 2442 //
>> 2443 lambda = 0.0;
>> 2444 }
>> 2445 else
>> 2446 {
>> 2447 #if DEBUGTORUS
>> 2448 G4cout << "G4Torus::SolveNumeric"
>> 2449 << " Exiting Surface (Recalculating) or DistanceToOut"
>> 2450 << " from surface"
>> 2451 << G4endl ;
>> 2452 G4cout << "G4Torus::SolveNumeric Recursive call lambda..."
>> 2453 << lambda << G4endl << G4endl;
>> 2454 #endif
>> 2455 // else it is not necessarily infinity !!
>> 2456 // (we could reach the opposite side..)
>> 2457 // To reach the opposite side we remark that from Surface
>> 2458 // the sphere of radius min((Rmax - Rmin)/2, Rmin) does not
>> 2459 // hit 2 surfaces of the torus so it is safe to do that way
>> 2460
>> 2461 if ((GetRmax() - GetRmin())/2.0 < GetRmin())
>> 2462 {
>> 2463 lambda =
>> 2464 SolveNumeric(p+((GetRmax()-GetRmin())/2.0)*v,v,IsDistanceToIn)
>> 2465 + (GetRmax() - GetRmin())/2.0;
>> 2466 }
>> 2467 else
>> 2468 {
>> 2469 lambda =
>> 2470 SolveNumeric(p+GetRmin()*v,v,IsDistanceToIn) + GetRmin();
>> 2471 }
>> 2472
>> 2473 #if DEBUGTORUS
>> 2474 G4cout << "G4Torus::SolveNumeric --> Recursive call: lambda = "
>> 2475 << lambda << G4endl;
>> 2476 #endif
>> 2477 }
>> 2478 }
>> 2479 else
>> 2480 {
>> 2481 // PROBLEM : could this be better done ?
1600 2482
1601 if ((fSPhi == 0) && (fDPhi == twopi)){ aSid << 2483 G4double lambdaToRmax = DistanceToTorus(p.x(),p.y(),p.z(),
1602 chose = G4RandFlat::shoot(0.,aOut + aIn + 2 << 2484 v.x(),v.y(),v.z(),
>> 2485 GetRtor(),GetRmax());
>> 2486 if (lambda >= lambdaToRmax)
>> 2487 {
>> 2488 #if DEBUGTORUS
>> 2489 G4cout << "G4Torus::SolveNumeric"
>> 2490 << " Point does not hit the Rmin torus from here"
>> 2491 << G4endl;
>> 2492 #endif
>> 2493 lambda = lambdaToRmax;
>> 2494 }
>> 2495 else
>> 2496 {
>> 2497 #if DEBUGTORUS
>> 2498 G4cout << "G4Torus::SolveNumeric We hit the Rmin torus with "
>> 2499 << lambda << G4endl;
>> 2500 G4cout << "G4Torus::SolveNumeric Note that this could be small "
>> 2501 << "and not in Tolerance resulting in wrong result "
>> 2502 << G4endl ;
>> 2503 #endif
>> 2504 }
>> 2505 }
>> 2506 }
>> 2507 else
>> 2508 {
>> 2509 // It is a whole torus
>> 2510 //
>> 2511 lambda = DistanceToTorus(p.x(),p.y(),p.z(),
>> 2512 v.x(),v.y(),v.z(),
>> 2513 GetRtor(),GetRmax());
>> 2514 }
>> 2515 break;
>> 2516 case kSurface:
>> 2517 {
>> 2518 G4double Lx,Ly,Lz;
>> 2519 G4double scal;
>> 2520
>> 2521 #if DEBUGTORUS
>> 2522 G4cout << "G4Torus::SolveNumeric Point is on the Rmax Surface"
>> 2523 << G4endl ;
>> 2524 #endif
>> 2525 // It is possible with Phi that this is not the correct point
>> 2526 //
>> 2527 lambda = DistanceToTorus(p.x(),p.y(),p.z(),
>> 2528 v.x(),v.y(),v.z(),
>> 2529 GetRtor(),GetRmax());
>> 2530 // Compute Surface point
>> 2531 //
>> 2532 Lx = p.x() + lambda*v.x();
>> 2533 Ly = p.y() + lambda*v.y();
>> 2534 Lz = p.z() + lambda*v.z();
>> 2535
>> 2536 // Scalar product
>> 2537 scal = v.x()*TorusDerivativeX(Lx,Ly,Lz,GetRtor(),GetRmax());
>> 2538 scal += v.y()*TorusDerivativeY(Lx,Ly,Lz,GetRtor(),GetRmax());
>> 2539 scal += v.z()*TorusDerivativeZ(Lx,Ly,Lz,GetRtor(),GetRmax());
>> 2540
>> 2541 // if entering it is < 0.0
>> 2542 //
>> 2543 if ( (IsDistanceToIn) && (scal < 0.0) )
>> 2544 {
>> 2545 #if DEBUGTORUS
>> 2546 G4cout << "G4Torus::SolveNumeric Point is Entering Surface "
>> 2547 << scal << G4endl ;
>> 2548 #endif
>> 2549 lambda = 0.0;
>> 2550 }
>> 2551 else
>> 2552 {
>> 2553 #if DEBUGTORUS
>> 2554 G4cout << "G4Torus::SolveNumeric Point is Exiting Surface "
>> 2555 << "or DistanceToOut " << scal << G4endl ;
>> 2556 G4cout << "Recursive call ..." << G4endl << G4endl ;
>> 2557 #endif
>> 2558 // To reach the opposite side we remark that from Surface the
>> 2559 // sphere of radius (Rmax - Rmin)/2 does not hit 2 surfaces of
>> 2560 // the torus so it is safe to do that way
>> 2561 // lambda = SolveNumeric(p+(lambda + kCarTolerance)*v,
>> 2562 // v,IsDistanceToIn);
>> 2563 //
>> 2564 lambda =
>> 2565 SolveNumeric(p+((GetRmax()-GetRmin())/2.0)*v, v, IsDistanceToIn)
>> 2566 + (GetRmax() - GetRmin())/2.0;
>> 2567 #if DEBUGTORUS
>> 2568 G4cout << "Recursive call ...END" << G4endl ;
>> 2569 #endif
>> 2570 }
>> 2571 }
>> 2572 break;
>> 2573 case kOutside:
>> 2574 #if DEBUGTORUS
>> 2575 G4cout << "G4Torus::SolveNumeric Point is Outside the Rmax torus"
>> 2576 << G4endl ;
>> 2577 #endif
>> 2578
>> 2579 lambda = DistanceToTorus(p.x(),p.y(),p.z(),
>> 2580 v.x(),v.y(),v.z(),
>> 2581 GetRtor(),GetRmax());
>> 2582 break;
>> 2583 }
>> 2584
>> 2585 if (lambda == kInfinity) return lambda;
>> 2586
>> 2587 #if DEBUGTORUS
>> 2588 G4cout << "G4Torus::SolveNumeric Intersection found. "
>> 2589 << "Now checking Phi angles" << G4endl ;
>> 2590 #endif
>> 2591
>> 2592 // Ok we have a lambda that is correct without Phi
>> 2593 // Now check Phi ..
>> 2594
>> 2595 // Eliminate the case of point (0,0,0)
>> 2596 //
>> 2597 if (((p.x()+ lambda*v.x())*(p.x()+ lambda*v.x()) +
>> 2598 (p.y()+ lambda*v.y())*(p.y()+ lambda*v.y()) +
>> 2599 (p.z()+ lambda*v.z())*(p.z()+ lambda*v.z())) > POLEPSILON)
>> 2600 {
>> 2601 G4double theta = atan2(p.y() + lambda*v.y(),p.x() + lambda*v.x());
>> 2602
>> 2603 #if DEBUGTORUS
>> 2604 G4cout << "G4Torus::SolveNumeric theta = " << theta << G4endl;
>> 2605 #endif
>> 2606
>> 2607 if (theta < 0) theta += 2*M_PI;
>> 2608
>> 2609 // We have to verify if this root is inside the region between
>> 2610 // fSPhi and fSPhi + fDPhi
>> 2611
>> 2612 #if DEBUGTORUS
>> 2613 G4cout << "G4Torus::SolveNumeric theta = " << theta
>> 2614 << " Phi = " << fSPhi
>> 2615 << " Phi + dPhi = " << fSPhi + fDPhi
>> 2616 << " kAngTolerance = " << kAngTolerance << G4endl ;
>> 2617 G4cout << " theta - Phi = " << theta - fSPhi << G4endl ;
>> 2618 #endif
>> 2619
>> 2620 if ( (theta - fSPhi >= - kAngTolerance*0.5)
>> 2621 && (theta - (fSPhi + fDPhi) <= kAngTolerance*0.5) )
>> 2622 {
>> 2623 // If this is the case we return this solution
>> 2624
>> 2625 #if DEBUGTORUS
>> 2626 G4cout << "G4Torus::SolveNumeric Correct Phi section" << G4endl ;
>> 2627 #endif
>> 2628
>> 2629 return lambda;
>> 2630 }
>> 2631 else
>> 2632 {
>> 2633 // Else we compute the intersection with the 2 half-plane [fSPhi]
>> 2634 // and [fSPhi + fDPhi]
>> 2635
>> 2636 G4double IntersectPlanar ;
>> 2637 IntersectPlanar = -(p.y()-p.x()*tan(fSPhi))/(v.y()-v.x()*tan(fSPhi));
>> 2638
>> 2639 #if DEBUGTORUS
>> 2640 G4cout << "G4Torus::SolveNumeric IntersectPlanar = "
>> 2641 << IntersectPlanar << G4endl ;
>> 2642 #endif
>> 2643
>> 2644 // If this is below lambda we check for the other plane
>> 2645 //
>> 2646 if (IntersectPlanar < lambda)
>> 2647 {
>> 2648 IntersectPlanar = - (p.y() - p.x()*tan(fSPhi + fDPhi))
>> 2649 / (v.y() - v.x()*tan(fSPhi + fDPhi)) ;
>> 2650 #if DEBUGTORUS
>> 2651 G4cout << "G4Torus::SolveNumeric IntersectPlanar (2) = "
>> 2652 << IntersectPlanar << G4endl ;
>> 2653 #endif
>> 2654 }
>> 2655
>> 2656 // If we does not hit the two plan then we does not hit the torus ..
>> 2657 //
>> 2658 if (IntersectPlanar < lambda)
>> 2659 {
>> 2660 #if DEBUGTORUS
>> 2661 G4cout << "G4Torus::SolveNumeric No intersection with planar Phi .."
>> 2662 << G4endl ;
>> 2663 #endif
>> 2664 return kInfinity;
>> 2665 }
>> 2666
>> 2667 #if DEBUGTORUS
>> 2668 G4cout << "G4Torus::SolveNumeric Incorrect Phi section" << G4endl ;
>> 2669 G4cout << "G4Torus::SolveNumeric point : " << p << " direction : "
>> 2670 << v << G4endl ;
>> 2671 G4cout << "G4Torus::SolveNumeric IntersectPlanar = "
>> 2672 << IntersectPlanar << G4endl ;
>> 2673 #endif
>> 2674
>> 2675 if ( (TorusEquation(p.x() + IntersectPlanar*v.x(),
>> 2676 p.y() + IntersectPlanar*v.y(),
>> 2677 p.z() + IntersectPlanar*v.z(),
>> 2678 GetRtor(),GetRmax()) < 0)
>> 2679 && (TorusEquation(p.x() + IntersectPlanar*v.x(),
>> 2680 p.y() + IntersectPlanar*v.y(),
>> 2681 p.z() + IntersectPlanar*v.z(),
>> 2682 GetRtor(),GetRmin()) > 0) )
>> 2683 {
>> 2684 // if this point is inside torus Rmax and outside torus Rmin
>> 2685 // then it is on the cut planar faces
>> 2686
>> 2687 #if DEBUGTORUS
>> 2688 G4cout << "G4Torus::SolveNumeric Hit planar section" << G4endl ;
>> 2689 #endif
>> 2690 return IntersectPlanar;
>> 2691 }
>> 2692 else
>> 2693 {
>> 2694 // else we continue from this new point (SolveNumeric)
>> 2695
>> 2696 #if DEBUGTORUS
>> 2697 G4cout << "G4Torus::SolveNumeric Recursive Phi call with "
>> 2698 << IntersectPlanar << " .." << G4endl << G4endl;
>> 2699 #endif
>> 2700
>> 2701 return IntersectPlanar + SolveNumeric(p+IntersectPlanar*v,
>> 2702 v,IsDistanceToIn);
>> 2703 }
>> 2704 }
>> 2705 }
>> 2706 else
>> 2707 {
>> 2708 #if DEBUGTORUS
>> 2709 G4cout << "G4Torus::SolveNumeric Phi not checked because point is "
>> 2710 << p + lambda*v << G4endl << G4endl;
>> 2711 #endif
>> 2712 }
>> 2713
>> 2714 return lambda;
>> 2715 }
>> 2716
>> 2717 ///////////////////////////////////////////////////////////////////////
>> 2718 //
>> 2719 // Utility function
>> 2720
>> 2721 void G4Torus::BVMIntersection( G4double x,G4double y,G4double z,
>> 2722 G4double dx,G4double dy,G4double dz,
>> 2723 G4double Rmax, G4double Rmin,
>> 2724 G4double *NewL,G4int *valid ) const
>> 2725 {
>> 2726
>> 2727 if (dz != 0)
>> 2728 {
>> 2729 G4double DistToZ ;
>> 2730
>> 2731 NewL[0] = (Rmin - z)/dz ; // z = + Rmin
>> 2732 NewL[1] = (-Rmin - z)/dz ; // z = - Rmin
>> 2733
>> 2734 // Test validity here (*** To be optimized ***)
>> 2735 //
>> 2736 if (NewL[0] < 0.0) valid[0] = 0;
>> 2737 if (NewL[1] < 0.0) valid[1] = 0;
>> 2738 DistToZ = (x+NewL[0]*dx)*(x+NewL[0]*dx) + (y+NewL[0]*dy)*(y+NewL[0]*dy);
>> 2739 if (DistToZ - (Rmax + Rmin)*(Rmax + Rmin) > 0)
>> 2740 valid[0] = 0;
>> 2741
>> 2742 #if HOLEBVM
>> 2743 if (DistToZ - (Rmax - Rmin)*(Rmax - Rmin) < 0)
>> 2744 valid[0] = 0;
>> 2745 #endif
>> 2746
>> 2747 DistToZ = (x+NewL[1]*dx)*(x+NewL[1]*dx) + (y+NewL[1]*dy)*(y+NewL[1]*dy);
>> 2748 if (DistToZ - (Rmax + Rmin)*(Rmax + Rmin) > 0)
>> 2749 valid[1] = 0;
>> 2750
>> 2751 #if HOLEBVM
>> 2752 if (DistToZ - (Rmax - Rmin)*(Rmax - Rmin) < 0)
>> 2753 valid[1] = 0;
>> 2754 #endif
>> 2755 }
>> 2756 else
>> 2757 {
>> 2758 // if dz == 0 we could know the exact solution
>> 2759 // Well, this is true but we have the expected precision
>> 2760 // issue from sqrt ...
>> 2761 NewL[0] = -1.0;
>> 2762 NewL[1] = -1.0;
>> 2763 valid[0] = 0;
>> 2764 valid[1] = 0;
>> 2765 }
>> 2766
>> 2767 // x² + y² = (Rmax + Rmin)²
>> 2768 //
>> 2769 if ((dx != 0) || (dy != 0))
>> 2770 {
>> 2771 G4double a,b,c,d;
>> 2772
>> 2773 a = dx*dx + dy*dy ;
>> 2774 b = 2*(x*dx + y*dy) ;
>> 2775 c = x*x + y*y - (Rmax + Rmin)*(Rmax + Rmin) ;
>> 2776 d = b*b - 4*a*c ;
>> 2777
>> 2778 if (d < 0)
>> 2779 {
>> 2780 valid[2] = 0;
>> 2781 valid[3] = 0;
>> 2782 NewL[2] = -1.0;
>> 2783 NewL[3] = -1.0;
>> 2784 }
>> 2785 else
>> 2786 {
>> 2787 d = sqrt(d) ;
>> 2788 NewL[2] = (d - b)/(2*a);
>> 2789 NewL[3] = (-d - b)/(2*a);
>> 2790 if (NewL[2] < 0.0) valid[2] = 0;
>> 2791 if (fabs(z + NewL[2]*dz) - Rmin > POLEPSILON) valid[2] = 0;
>> 2792 if (NewL[3] < 0.0) valid[3] = 0;
>> 2793 if (fabs(z + NewL[3]*dz) - Rmin > POLEPSILON) valid[3] = 0;
>> 2794 }
>> 2795 }
>> 2796 else
>> 2797 {
>> 2798 // only dz != 0 so we could know the exact solution
>> 2799 // this depends only for the distance to Z axis
>> 2800 // BUT big precision problem near the border..
>> 2801
>> 2802 NewL[2] = -1.0;
>> 2803 NewL[3] = -1.0;
>> 2804 valid[2] = 0;
>> 2805 valid[3] = 0;
>> 2806
>> 2807 /* SQRT Test ---------------------------------------------------
>> 2808
>> 2809 // Try This to see precision issue with sqrt(~ 0)
>> 2810 //
>> 2811 G4double DistToZ ;
>> 2812 G4double result;
>> 2813 G4double guess;
>> 2814
>> 2815 DistToZ = sqrt(x*x + y*y) ;
>> 2816
>> 2817 if ((DistToZ < (Rmax - Rmin)) || (DistToZ > (Rmax + Rmin)))
>> 2818 {
>> 2819 return -1.0 ;
>> 2820 }
>> 2821
>> 2822 result = sqrt((Rmin + Rmax - DistToZ)*(Rmin - Rmax + DistToZ));
>> 2823
>> 2824 if (dz < 0)
>> 2825 {
>> 2826 if (z > result)
>> 2827 {
>> 2828 return (result - z)/dz;
>> 2829 }
>> 2830 else
>> 2831 {
>> 2832 if (z > -result)
>> 2833 {
>> 2834 return (-result - z)/dz;
>> 2835 }
>> 2836 else
>> 2837 return -1.0;
>> 2838 }
>> 2839 }
>> 2840 else
>> 2841 {
>> 2842 if (z < -result)
>> 2843 {
>> 2844 return (z + result)/dz;
>> 2845 }
>> 2846 else
>> 2847 {
>> 2848 if (z < result)
>> 2849 {
>> 2850 return (z - result)/dz;
>> 2851 }
>> 2852 else
>> 2853 return -1.0;
>> 2854 }
>> 2855 }
>> 2856 --------------------------------------------------- END SQRT test */
>> 2857 }
>> 2858
>> 2859 // x² + y² = (Rmax - Rmin)²
>> 2860 //
>> 2861 #if HOLEBVM
>> 2862 if ((dx != 0) || (dy != 0))
>> 2863 {
>> 2864 G4double a,b,c,d;
>> 2865
>> 2866 a = dx*dx + dy*dy ;
>> 2867 b = 2*(x*dx + y*dy) ;
>> 2868 c = x*x + y*y - (Rmax - Rmin)*(Rmax - Rmin) ;
>> 2869 d = b*b - 4*a*c ;
>> 2870
>> 2871 if (d < 0)
>> 2872 {
>> 2873 valid[4] = 0;
>> 2874 valid[5] = 0;
>> 2875 NewL[4] = -1.0;
>> 2876 NewL[5] = -1.0;
>> 2877 }
>> 2878 else
>> 2879 {
>> 2880 d = sqrt(d) ;
>> 2881 NewL[4] = (d - b)/(2*a);
>> 2882 NewL[5] = (-d - b)/(2*a);
>> 2883 if (NewL[4] < 0.0) valid[4] = 0;
>> 2884 if (fabs(z + NewL[4]*dz) - Rmin > POLEPSILON) valid[4] = 0;
>> 2885 if (NewL[5] < 0.0) valid[5] = 0;
>> 2886 if (fabs(z + NewL[5]*dz) - Rmin > POLEPSILON) valid[5] = 0;
>> 2887 }
>> 2888 }
>> 2889 else
>> 2890 #endif
>> 2891 {
>> 2892 // only dz != 0 so we could know the exact solution
>> 2893 // OK but same as above ..
>> 2894 //
>> 2895 valid[4] = 0;
>> 2896 valid[5] = 0;
>> 2897 NewL[4] = -1.0;
>> 2898 NewL[5] = -1.0;
>> 2899 }
>> 2900 }
>> 2901
>> 2902 ///////////////////////////////////////////////////////////////////////
>> 2903 //
>> 2904 // Utility function
1603 2905
1604 if(chose < aOut) << 2906 void G4Torus::SortIntervals ( G4double *SortL, G4double *NewL,
>> 2907 G4int *valid, G4int *NbIntersection ) const
>> 2908 {
>> 2909 G4int i,j;
>> 2910 G4double swap;
>> 2911
>> 2912 (*NbIntersection) = 0;
>> 2913 SortL[0] = -kInfinity;
>> 2914
>> 2915 for (i=0;i<6;i++)
>> 2916 {
>> 2917 if (valid[i] != 0)
>> 2918 {
>> 2919 SortL[(*NbIntersection)] = NewL[i] ;
>> 2920 for (j=(*NbIntersection);j>0;j--)
>> 2921 {
>> 2922 if (SortL[j] < SortL[j-1])
>> 2923 {
>> 2924 swap = SortL[j-1] ;
>> 2925 SortL[j-1] = SortL[j];
>> 2926 SortL[j] = swap;
>> 2927 }
>> 2928 }
>> 2929 (*NbIntersection) ++;
>> 2930 }
>> 2931 }
>> 2932
>> 2933 // Delete double values
>> 2934 // When the ray hits a corner we have a double value
>> 2935 //
>> 2936 for (i=0;i<(*NbIntersection)-1;i++)
>> 2937 {
>> 2938 if (SortL[i+1] - SortL[i] < POLEPSILON)
>> 2939 {
>> 2940 if (((*NbIntersection) & (1)) == 1)
>> 2941 {
>> 2942 // If the NbIntersection is odd then we keep one value
>> 2943 //
>> 2944 for (j=i+1;j<(*NbIntersection);j++)
>> 2945 {
>> 2946 SortL[j-1] = SortL[j] ;
>> 2947 }
>> 2948 (*NbIntersection) --;
>> 2949 }
>> 2950 else
>> 2951 {
>> 2952 // If it is even we delete the 2 values
>> 2953 //
>> 2954 for (j=i+2;j<(*NbIntersection);j++)
>> 2955 {
>> 2956 SortL[j-2] = SortL[j] ;
>> 2957 }
>> 2958 (*NbIntersection) -= 2;
>> 2959 }
>> 2960 }
>> 2961 }
>> 2962 }
>> 2963
>> 2964 ///////////////////////////////////////////////////////////////////////
>> 2965 //
>> 2966 // Utility function
>> 2967
>> 2968 G4double G4Torus::DistanceToTorus ( G4double x, G4double y, G4double z,
>> 2969 G4double dx, G4double dy, G4double dz,
>> 2970 G4double Rmax, G4double Rmin ) const
>> 2971 {
>> 2972 G4double Lmin=0.;
>> 2973 G4double Lmax=0.;
>> 2974 G4double guess;
>> 2975 G4double SortL[4];
>> 2976
>> 2977 G4int NbIntersection = 0;
>> 2978
>> 2979 G4double NewL[NBPOINT];
>> 2980 G4int valid[] = {1,1,1,1,1,1} ;
>> 2981 G4int j;
>> 2982
>> 2983 j = 0;
>> 2984
>> 2985 // Compute Intervals from Bounding Volume
>> 2986 //
>> 2987 BVMIntersection(x,y,z,dx,dy,dz,Rmax,Rmin,NewL,valid);
>> 2988
>> 2989 // We could compute intervals value
>> 2990 // Sort all valid NewL to SortL.
>> 2991 // There must be 4 values at max and
>> 2992 // odd one if point is inside
>> 2993
>> 2994 SortIntervals(SortL,NewL,valid,&NbIntersection);
>> 2995
>> 2996 {
>> 2997 // Length check (Torus specific)
>> 2998 //
>> 2999 G4double LengthMin = 0.82842712*Rmin;
>> 3000
>> 3001 switch(NbIntersection)
>> 3002 {
>> 3003 case 1:
>> 3004 if (SortL[0] < POLEPSILON)
>> 3005 {
>> 3006 if (fabs(TorusEquation(x,y,z,Rmax,Rmin)) < TORUSPRECISION)
>> 3007 {
>> 3008 return 0.0;
>> 3009 }
>> 3010 else
>> 3011 {
>> 3012 return kInfinity;
>> 3013 }
>> 3014 }
>> 3015 break;
>> 3016 case 2:
>> 3017 if ((SortL[1] - SortL[0]) < LengthMin) NbIntersection = 0;
>> 3018 break;
>> 3019 case 3:
>> 3020 if (SortL[0] < POLEPSILON)
>> 3021 {
>> 3022 if (fabs(TorusEquation(x,y,z,Rmax,Rmin)) < TORUSPRECISION)
>> 3023 {
>> 3024 return 0.0;
>> 3025 }
>> 3026 else
>> 3027 {
>> 3028 NbIntersection --;
>> 3029 SortL[0] = SortL[1] ;
>> 3030 SortL[1] = SortL[2] ;
>> 3031 if ((SortL[1] - SortL[0]) < LengthMin) NbIntersection = 0;
>> 3032 }
>> 3033 }
>> 3034 else
>> 3035 {
>> 3036 if ((SortL[2] - SortL[1]) < LengthMin) NbIntersection -= 2;
>> 3037 }
>> 3038 break;
>> 3039 case 4:
>> 3040 if ((SortL[1] - SortL[0]) < LengthMin)
>> 3041 {
>> 3042 NbIntersection -= 2;
>> 3043 SortL[0] = SortL[2];
>> 3044 SortL[1] = SortL[3];
>> 3045 if ((SortL[1] - SortL[0]) < LengthMin) NbIntersection -= 2;
>> 3046 }
>> 3047 break;
>> 3048 }
>> 3049 }
>> 3050
>> 3051 #if DEBUGTORUS
1605 { 3052 {
1606 return { (fRtor+fRmax*cosv)*cosu, (fRtor+ << 3053 G4int i;
>> 3054 G4cout.precision(16);
>> 3055 G4cout << "G4Torus::DistanceToTorus INTERVALS" << G4endl ;
>> 3056 for (i=0;i<NbIntersection;i++)
>> 3057 {
>> 3058 G4cout << "G4Torus::DistanceToTorus " << SortL[i] << G4endl ;
>> 3059 }
1607 } 3060 }
1608 else if( (chose >= aOut) && (chose < aOut + << 3061 #endif
>> 3062
>> 3063 switch (NbIntersection)
1609 { 3064 {
1610 return { (fRtor+fRmin*cosv)*cosu, (fRtor+ << 3065 case 0:
>> 3066 return kInfinity ;
>> 3067 break;
>> 3068 case 1:
>> 3069 Lmin = 0.0 ;
>> 3070 Lmax = SortL[0] ;
>> 3071 break;
>> 3072 case 2:
>> 3073 Lmin = SortL[0] ;
>> 3074 Lmax = SortL[1] ;
>> 3075 break;
>> 3076 #if HOLEBVM
>> 3077 case 3:
>> 3078 Lmin = 0.0 ;
>> 3079 Lmax = SortL[0] ;
>> 3080
>> 3081 G4TorusEquation torus (Rmax,Rmin);
>> 3082 torus.setPosition(x,y,z);
>> 3083 torus.setDirection(dx,dy,dz);
>> 3084
>> 3085 G4PolynomialSolver<G4TorusEquation,
>> 3086 G4double(G4TorusEquation::*)(G4double)>
>> 3087 PolySolver(&torus, &G4TorusEquation::Function,
>> 3088 &G4TorusEquation::Derivative, TORUSPRECISION) ;
>> 3089 guess = PolySolver.solve(Lmin,Lmax);
>> 3090 if ((guess >= (Lmin - POLEPSILON)) && (guess <= (Lmax + POLEPSILON)))
>> 3091 {
>> 3092 return guess ;
>> 3093 }
>> 3094 else
>> 3095 {
>> 3096 Lmin = SortL[1] ;
>> 3097 Lmax = SortL[2] ;
>> 3098 }
>> 3099 break;
>> 3100 case 4:
>> 3101 Lmin = SortL[0] ;
>> 3102 Lmax = SortL[1] ;
>> 3103
>> 3104 G4TorusEquation torus (Rmax,Rmin);
>> 3105 torus.setPosition(x,y,z);
>> 3106 torus.setDirection(dx,dy,dz);
>> 3107
>> 3108 G4PolynomialSolver<G4TorusEquation,
>> 3109 G4double(G4TorusEquation::*)(G4double)>
>> 3110 PolySolver(&torus, &G4TorusEquation::Function,
>> 3111 &G4TorusEquation::Derivative, TORUSPRECISION) ;
>> 3112 guess = PolySolver.solve(Lmin,Lmax);
>> 3113 if ((guess >= (Lmin - POLEPSILON)) && (guess <= (Lmax + POLEPSILON)))
>> 3114 {
>> 3115 return guess ;
>> 3116 }
>> 3117 else
>> 3118 {
>> 3119 Lmin = SortL[2] ;
>> 3120 Lmax = SortL[3] ;
>> 3121 }
>> 3122 break;
>> 3123 #endif
>> 3124
>> 3125 default:
>> 3126 G4cerr << "G4Torus::DistanceToTorus NbIntersection = "
>> 3127 << NbIntersection << G4endl;
>> 3128 break;
1611 } 3129 }
1612 else if( (chose >= aOut + aIn) && (chose < << 3130
>> 3131 G4TorusEquation torus (Rmax,Rmin);
>> 3132 torus.setPosition(x,y,z);
>> 3133 torus.setDirection(dx,dy,dz);
>> 3134
>> 3135 G4PolynomialSolver<G4TorusEquation,
>> 3136 G4double(G4TorusEquation::*)(G4double)>
>> 3137 PolySolver(&torus, &G4TorusEquation::Function,
>> 3138 &G4TorusEquation::Derivative, TORUSPRECISION) ;
>> 3139 guess = PolySolver.solve(Lmin,Lmax);
>> 3140 if ((guess >= (Lmin - POLEPSILON)) && (guess <= (Lmax + POLEPSILON)))
1613 { 3141 {
1614 rRand = GetRadiusInRing(fRmin,fRmax); << 3142 #if DEBUGTORUS
1615 return { (fRtor+rRand*cosv)*std::cos(fSPh << 3143 G4cout << "G4Torus::DistanceToTorus distance = " << guess << G4endl ;
1616 (fRtor+rRand*cosv)*std::sin(fSPh << 3144 #endif
>> 3145 return guess ;
1617 } 3146 }
1618 else 3147 else
1619 { << 3148 {
1620 rRand = GetRadiusInRing(fRmin,fRmax); << 3149 #if DEBUGTORUS
1621 return { (fRtor+rRand*cosv)*std::cos(fSPh << 3150 G4cout << "G4Torus::DistanceToTorus : kInfinity" << G4endl ;
1622 (fRtor+rRand*cosv)*std::sin(fSPh << 3151 #endif
1623 } << 3152 return kInfinity;
>> 3153 }
1624 } 3154 }
1625 3155
1626 ///////////////////////////////////////////// 3156 ///////////////////////////////////////////////////////////////////////
1627 // 3157 //
1628 // Visualisation Functions << 3158 // G4TorusEquation definition
1629 3159
1630 void G4Torus::DescribeYourselfTo ( G4VGraphic << 3160 G4TorusEquation::G4TorusEquation()
1631 { 3161 {
1632 scene.AddSolid (*this); <<
1633 } 3162 }
1634 3163
1635 G4Polyhedron* G4Torus::CreatePolyhedron () co << 3164 G4TorusEquation::G4TorusEquation(G4double Rmax, G4double Rmin)
1636 { 3165 {
1637 return new G4PolyhedronTorus (fRmin, fRmax, << 3166 R0 = Rmax;
>> 3167 R1 = Rmin;
1638 } 3168 }
1639 3169
1640 #endif // !defined(G4GEOM_USE_TORUS) || !defi << 3170 G4TorusEquation::~G4TorusEquation()
>> 3171 {
>> 3172 }
1641 3173