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