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