| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
>> 23 //
>> 24 // $Id: G4Sphere.cc,v 1.19 2003/06/16 16:53:38 gunter Exp $
>> 25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
>> 26 //
>> 27 // class G4Sphere
>> 28 //
26 // Implementation for G4Sphere class 29 // Implementation for G4Sphere class
27 // 30 //
>> 31 // History:
28 // 28.03.94 P.Kent: old C++ code converted to 32 // 28.03.94 P.Kent: old C++ code converted to tolerant geometry
29 // 17.09.96 V.Grichine: final modifications to 33 // 17.09.96 V.Grichine: final modifications to commit
30 // 30.10.03 J.Apostolakis: new algorithm in In << 34 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
31 // 03.05.05 V.Grichine: SurfaceNormal(p) accor << 35 // 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections
32 // 22.07.05 O.Link: Added check for intersecti << 36 // 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections
33 // 26.03.09 G.Cosmo: optimisations and uniform << 37 // 18.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
34 // 26.10.16 E.Tcherniaev: re-implemented Calcu << 38 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
35 // G4BoundingEnvelope, << 39 // 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451
36 // ------------------------------------------- << 40 // 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance
>> 41 //
>> 42 // ********************************************************************
37 43
38 #include "G4Sphere.hh" << 44 #include <assert.h>
39 45
40 #if !defined(G4GEOM_USE_USPHERE) << 46 #include "G4Sphere.hh"
41 47
42 #include "G4GeomTools.hh" <<
43 #include "G4VoxelLimits.hh" 48 #include "G4VoxelLimits.hh"
44 #include "G4AffineTransform.hh" 49 #include "G4AffineTransform.hh"
45 #include "G4GeometryTolerance.hh" <<
46 #include "G4BoundingEnvelope.hh" <<
47 50
48 #include "G4VPVParameterisation.hh" 51 #include "G4VPVParameterisation.hh"
49 52
50 #include "G4QuickRand.hh" <<
51 <<
52 #include "meshdefs.hh" 53 #include "meshdefs.hh"
53 54
54 #include "G4VGraphicsScene.hh" 55 #include "G4VGraphicsScene.hh"
55 #include "G4VisExtent.hh" << 56 #include "G4Polyhedron.hh"
56 << 57 #include "G4NURBS.hh"
57 using namespace CLHEP; << 58 #include "G4NURBSbox.hh"
58 59
59 // Private enum: Not for external use - used b 60 // Private enum: Not for external use - used by distanceToOut
60 61
61 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTh 62 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTheta,kETheta};
62 63
63 // used by normal 64 // used by normal
64 65
65 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSThe 66 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta};
66 67
>> 68 /////////////////////////////////////////////////////////////////////
>> 69 //
>> 70 // Destructor
>> 71
>> 72 G4Sphere::~G4Sphere()
>> 73 {
>> 74 }
>> 75
67 ////////////////////////////////////////////// 76 ////////////////////////////////////////////////////////////////////////
68 // 77 //
69 // constructor - check parameters, convert ang 78 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
70 // - note if pDPhi>2PI then reset 79 // - note if pDPhi>2PI then reset to 2PI
71 80
72 G4Sphere::G4Sphere( const G4String& pName, 81 G4Sphere::G4Sphere( const G4String& pName,
73 G4double pRmin, G4do 82 G4double pRmin, G4double pRmax,
74 G4double pSPhi, G4do 83 G4double pSPhi, G4double pDPhi,
75 G4double pSTheta, G4 84 G4double pSTheta, G4double pDTheta )
76 : G4CSGSolid(pName), fSPhi(0.0), fFullPhiSph << 85 : G4CSGSolid(pName)
77 { 86 {
78 kAngTolerance = G4GeometryTolerance::GetInst <<
79 kRadTolerance = G4GeometryTolerance::GetInst <<
80 <<
81 halfCarTolerance = 0.5*kCarTolerance; <<
82 halfAngTolerance = 0.5*kAngTolerance; <<
83 87
84 // Check radii and set radial tolerances << 88 // Check radii
85 89
86 if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTo << 90 if (pRmin<pRmax&&pRmin>=0)
>> 91 {
>> 92 fRmin=pRmin; fRmax=pRmax;
>> 93 }
>> 94 else
87 { 95 {
88 std::ostringstream message; << 96 G4cout << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
89 message << "Invalid radii for Solid: " << << 97 << " Invalid values for radii !" << G4endl
90 << " pRmin = " << pRmin << << 98 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
91 G4Exception("G4Sphere::G4Sphere()", "GeomS << 99 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
92 FatalException, message); << 100 << " Invalide values for radii ! - "
>> 101 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
>> 102 G4Exception("G4Sphere::G4Sphere() - invalid radii");
93 } 103 }
94 fRmin=pRmin; fRmax=pRmax; <<
95 fRminTolerance = (fRmin) != 0.0 ? std::max( <<
96 fRmaxTolerance = std::max( kRadTolerance, fE <<
97 104
98 // Check angles 105 // Check angles
99 106
100 CheckPhiAngles(pSPhi, pDPhi); << 107 if (pDPhi>=2.0*M_PI)
101 CheckThetaAngles(pSTheta, pDTheta); << 108 {
102 } << 109 fDPhi=2*M_PI;
103 << 110 }
104 ////////////////////////////////////////////// << 111 else if (pDPhi>0)
105 // << 112 {
106 // Fake default constructor - sets only member << 113 fDPhi=pDPhi;
107 // for usage restri << 114 }
108 // << 115 else
109 G4Sphere::G4Sphere( __void__& a ) << 116 {
110 : G4CSGSolid(a) << 117 G4cout << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
111 { << 118 << " Negative delta-Phi ! - "
112 } << 119 << pDPhi << G4endl;
113 << 120 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
114 ////////////////////////////////////////////// << 121 << " Negative Z delta-Phi ! - "
115 // << 122 << pDPhi << G4endl;
116 // Destructor << 123 G4Exception("G4Sphere::G4Sphere() - invalid DPhi");
>> 124 }
117 125
118 G4Sphere::~G4Sphere() = default; << 126 // Convert fSPhi to 0-2PI
119 127
120 ////////////////////////////////////////////// << 128 if (pSPhi<0)
121 // << 129 {
122 // Copy constructor << 130 fSPhi=2.0*M_PI-fmod(fabs(pSPhi),2.0*M_PI);
>> 131 }
>> 132 else
>> 133 {
>> 134 fSPhi=fmod(pSPhi,2.0*M_PI);
>> 135 }
123 136
124 G4Sphere::G4Sphere(const G4Sphere&) = default; << 137 // Sphere is placed such that fSPhi+fDPhi>2.0*M_PI !
>> 138 // fSPhi could be < 0 !!?
>> 139 //
>> 140 if (fSPhi+fDPhi>2.0*M_PI) fSPhi-=2.0*M_PI;
125 141
126 ////////////////////////////////////////////// << 142 // Check theta angles
127 // <<
128 // Assignment operator <<
129 143
130 G4Sphere& G4Sphere::operator = (const G4Sphere << 144 if (pSTheta<0 || pSTheta>M_PI)
131 { << 145 {
132 // Check assignment to self << 146 G4cout << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl;
133 // << 147 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl;
134 if (this == &rhs) { return *this; } << 148 G4Exception("G4Sphere::G4Sphere() - stheta outside 0-PI range");
135 << 149 }
136 // Copy base class data << 150 else
137 // << 151 {
138 G4CSGSolid::operator=(rhs); << 152 fSTheta=pSTheta;
139 << 153 }
140 // Copy data <<
141 // <<
142 fRminTolerance = rhs.fRminTolerance; fRmaxT <<
143 kAngTolerance = rhs.kAngTolerance; kRadTole <<
144 fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; <<
145 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSThe <<
146 fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPh <<
147 cosHDPhi = rhs.cosHDPhi; <<
148 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r <<
149 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh <<
150 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh <<
151 hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = <<
152 sinSTheta = rhs.sinSTheta; cosSTheta = rhs. <<
153 sinETheta = rhs.sinETheta; cosETheta = rhs. <<
154 tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs <<
155 tanETheta = rhs.tanETheta; tanETheta2 = rhs <<
156 eTheta = rhs.eTheta; fFullPhiSphere = rhs.f <<
157 fFullThetaSphere = rhs.fFullThetaSphere; fF <<
158 halfCarTolerance = rhs.halfCarTolerance; <<
159 halfAngTolerance = rhs.halfAngTolerance; <<
160 154
161 return *this; << 155 if (pDTheta+pSTheta>=M_PI)
>> 156 {
>> 157 fDTheta=M_PI-pSTheta;
>> 158 }
>> 159 else if (pDTheta>0)
>> 160 {
>> 161 fDTheta=pDTheta;
>> 162 }
>> 163 else
>> 164 {
>> 165 G4cout << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
>> 166 << " Negative delta-Theta ! - "
>> 167 << pDTheta << G4endl;
>> 168 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
>> 169 << " Negative delta-Theta ! - "
>> 170 << pDTheta << G4endl;
>> 171 G4Exception("G4Sphere::G4Sphere() - invalid pDTheta");
>> 172 }
162 } 173 }
163 174
164 ////////////////////////////////////////////// 175 //////////////////////////////////////////////////////////////////////////
165 // 176 //
166 // Dispatch to parameterisation for replicatio 177 // Dispatch to parameterisation for replication mechanism dimension
167 // computation & modification. 178 // computation & modification.
168 179
169 void G4Sphere::ComputeDimensions( G4VPVP 180 void G4Sphere::ComputeDimensions( G4VPVParameterisation* p,
170 const G4int 181 const G4int n,
171 const G4VPhy 182 const G4VPhysicalVolume* pRep)
172 { 183 {
173 p->ComputeDimensions(*this,n,pRep); 184 p->ComputeDimensions(*this,n,pRep);
174 } 185 }
175 186
176 ////////////////////////////////////////////// <<
177 // <<
178 // Get bounding box <<
179 <<
180 void G4Sphere::BoundingLimits(G4ThreeVector& p <<
181 { <<
182 G4double rmin = GetInnerRadius(); <<
183 G4double rmax = GetOuterRadius(); <<
184 <<
185 // Find bounding box <<
186 // <<
187 if (GetDeltaThetaAngle() >= pi && GetDeltaPh <<
188 { <<
189 pMin.set(-rmax,-rmax,-rmax); <<
190 pMax.set( rmax, rmax, rmax); <<
191 } <<
192 else <<
193 { <<
194 G4double sinStart = GetSinStartTheta(); <<
195 G4double cosStart = GetCosStartTheta(); <<
196 G4double sinEnd = GetSinEndTheta(); <<
197 G4double cosEnd = GetCosEndTheta(); <<
198 <<
199 G4double stheta = GetStartThetaAngle(); <<
200 G4double etheta = stheta + GetDeltaThetaAn <<
201 G4double rhomin = rmin*std::min(sinStart,s <<
202 G4double rhomax = rmax; <<
203 if (stheta > halfpi) rhomax = rmax*sinStar <<
204 if (etheta < halfpi) rhomax = rmax*sinEnd; <<
205 <<
206 G4TwoVector xymin,xymax; <<
207 G4GeomTools::DiskExtent(rhomin,rhomax, <<
208 GetSinStartPhi(),G <<
209 GetSinEndPhi(),Get <<
210 xymin,xymax); <<
211 <<
212 G4double zmin = std::min(rmin*cosEnd,rmax* <<
213 G4double zmax = std::max(rmin*cosStart,rma <<
214 pMin.set(xymin.x(),xymin.y(),zmin); <<
215 pMax.set(xymax.x(),xymax.y(),zmax); <<
216 } <<
217 <<
218 // Check correctness of the bounding box <<
219 // <<
220 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
221 { <<
222 std::ostringstream message; <<
223 message << "Bad bounding box (min >= max) <<
224 << GetName() << " !" <<
225 << "\npMin = " << pMin <<
226 << "\npMax = " << pMax; <<
227 G4Exception("G4Sphere::BoundingLimits()", <<
228 JustWarning, message); <<
229 DumpInfo(); <<
230 } <<
231 } <<
232 <<
233 ////////////////////////////////////////////// 187 ////////////////////////////////////////////////////////////////////////////
234 // 188 //
235 // Calculate extent under transform and specif 189 // Calculate extent under transform and specified limit
236 190
237 G4bool G4Sphere::CalculateExtent( const EAxis 191 G4bool G4Sphere::CalculateExtent( const EAxis pAxis,
238 const G4Voxe 192 const G4VoxelLimits& pVoxelLimit,
239 const G4Affi 193 const G4AffineTransform& pTransform,
240 G4doub 194 G4double& pMin, G4double& pMax ) const
241 { 195 {
242 G4ThreeVector bmin, bmax; << 196 if ( fDPhi==2.0*M_PI && fDTheta==M_PI) // !pTransform.IsRotated() &&
243 <<
244 // Get bounding box <<
245 BoundingLimits(bmin,bmax); <<
246 <<
247 // Find extent <<
248 G4BoundingEnvelope bbox(bmin,bmax); <<
249 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
250 } <<
251 <<
252 ////////////////////////////////////////////// <<
253 // <<
254 // Return whether point inside/outside/on surf <<
255 // Split into radius, phi, theta checks <<
256 // Each check modifies 'in', or returns as app <<
257 <<
258 EInside G4Sphere::Inside( const G4ThreeVector& <<
259 { <<
260 G4double rho,rho2,rad2,tolRMin,tolRMax; <<
261 G4double pPhi,pTheta; <<
262 EInside in = kOutside; <<
263 <<
264 const G4double halfRmaxTolerance = fRmaxTole <<
265 const G4double halfRminTolerance = fRminTole <<
266 const G4double Rmax_minus = fRmax - halfRmax <<
267 const G4double Rmin_plus = (fRmin > 0) ? fR <<
268 <<
269 rho2 = p.x()*p.x() + p.y()*p.y() ; <<
270 rad2 = rho2 + p.z()*p.z() ; <<
271 <<
272 // Check radial surfaces. Sets 'in' <<
273 <<
274 tolRMin = Rmin_plus; <<
275 tolRMax = Rmax_minus; <<
276 <<
277 if(rad2 == 0.0) <<
278 { 197 {
279 if (fRmin > 0.0) << 198 // Special case handling for solid spheres-shells
280 { << 199 // (rotation doesn't influence).
281 return in = kOutside; << 200 // Compute x/y/z mins and maxs for bounding box respecting limits,
282 } << 201 // with early returns if outside limits. Then switch() on pAxis,
283 if ( (!fFullPhiSphere) || (!fFullThetaSphe << 202 // and compute exact x and y limit for x/y case
>> 203
>> 204 G4double xoffset,xMin,xMax;
>> 205 G4double yoffset,yMin,yMax;
>> 206 G4double zoffset,zMin,zMax;
>> 207
>> 208 G4double diff1,diff2,maxDiff,newMin,newMax;
>> 209 G4double xoff1,xoff2,yoff1,yoff2;
>> 210
>> 211 xoffset=pTransform.NetTranslation().x();
>> 212 xMin=xoffset-fRmax;
>> 213 xMax=xoffset+fRmax;
>> 214 if (pVoxelLimit.IsXLimited())
284 { 215 {
285 return in = kSurface; << 216 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
>> 217 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
>> 218 {
>> 219 return false;
>> 220 }
>> 221 else
>> 222 {
>> 223 if (xMin<pVoxelLimit.GetMinXExtent())
>> 224 {
>> 225 xMin=pVoxelLimit.GetMinXExtent();
>> 226 }
>> 227 if (xMax>pVoxelLimit.GetMaxXExtent())
>> 228 {
>> 229 xMax=pVoxelLimit.GetMaxXExtent();
>> 230 }
>> 231 }
286 } 232 }
287 else << 233
>> 234 yoffset=pTransform.NetTranslation().y();
>> 235 yMin=yoffset-fRmax;
>> 236 yMax=yoffset+fRmax;
>> 237 if (pVoxelLimit.IsYLimited())
288 { 238 {
289 return in = kInside; << 239 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
>> 240 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
>> 241 {
>> 242 return false;
>> 243 }
>> 244 else
>> 245 {
>> 246 if (yMin<pVoxelLimit.GetMinYExtent())
>> 247 {
>> 248 yMin=pVoxelLimit.GetMinYExtent();
>> 249 }
>> 250 if (yMax>pVoxelLimit.GetMaxYExtent())
>> 251 {
>> 252 yMax=pVoxelLimit.GetMaxYExtent();
>> 253 }
>> 254 }
290 } 255 }
291 } <<
292 256
293 if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad << 257 zoffset=pTransform.NetTranslation().z();
294 { << 258 zMin=zoffset-fRmax;
295 in = kInside; << 259 zMax=zoffset+fRmax;
296 } << 260 if (pVoxelLimit.IsZLimited())
297 else <<
298 { <<
299 tolRMax = fRmax + halfRmaxTolerance; <<
300 tolRMin = std::max(fRmin-halfRminTolerance <<
301 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= <<
302 { 261 {
303 in = kSurface; << 262 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
>> 263 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
>> 264 {
>> 265 return false;
>> 266 }
>> 267 else
>> 268 {
>> 269 if (zMin<pVoxelLimit.GetMinZExtent())
>> 270 {
>> 271 zMin=pVoxelLimit.GetMinZExtent();
>> 272 }
>> 273 if (zMax>pVoxelLimit.GetMaxZExtent())
>> 274 {
>> 275 zMax=pVoxelLimit.GetMaxZExtent();
>> 276 }
>> 277 }
304 } 278 }
305 else << 279
>> 280 // Known to cut sphere
>> 281
>> 282 switch (pAxis)
306 { 283 {
307 return in = kOutside; << 284 case kXAxis:
>> 285 yoff1=yoffset-yMin;
>> 286 yoff2=yMax-yoffset;
>> 287 if (yoff1>=0&&yoff2>=0)
>> 288 {
>> 289 // Y limits cross max/min x => no change
>> 290 //
>> 291 pMin=xMin;
>> 292 pMax=xMax;
>> 293 }
>> 294 else
>> 295 {
>> 296 // Y limits don't cross max/min x => compute max delta x,
>> 297 // hence new mins/maxs
>> 298 //
>> 299 diff1=sqrt(fRmax*fRmax-yoff1*yoff1);
>> 300 diff2=sqrt(fRmax*fRmax-yoff2*yoff2);
>> 301 maxDiff=(diff1>diff2) ? diff1:diff2;
>> 302 newMin=xoffset-maxDiff;
>> 303 newMax=xoffset+maxDiff;
>> 304 pMin=(newMin<xMin) ? xMin : newMin;
>> 305 pMax=(newMax>xMax) ? xMax : newMax;
>> 306 }
>> 307 break;
>> 308 case kYAxis:
>> 309 xoff1=xoffset-xMin;
>> 310 xoff2=xMax-xoffset;
>> 311 if (xoff1>=0&&xoff2>=0)
>> 312 {
>> 313 // X limits cross max/min y => no change
>> 314 //
>> 315 pMin=yMin;
>> 316 pMax=yMax;
>> 317 }
>> 318 else
>> 319 {
>> 320 // X limits don't cross max/min y => compute max delta y,
>> 321 // hence new mins/maxs
>> 322 //
>> 323 diff1=sqrt(fRmax*fRmax-xoff1*xoff1);
>> 324 diff2=sqrt(fRmax*fRmax-xoff2*xoff2);
>> 325 maxDiff=(diff1>diff2) ? diff1:diff2;
>> 326 newMin=yoffset-maxDiff;
>> 327 newMax=yoffset+maxDiff;
>> 328 pMin=(newMin<yMin) ? yMin : newMin;
>> 329 pMax=(newMax>yMax) ? yMax : newMax;
>> 330 }
>> 331 break;
>> 332 case kZAxis:
>> 333 pMin=zMin;
>> 334 pMax=zMax;
>> 335 break;
>> 336 default:
>> 337 break;
308 } 338 }
>> 339 pMin-=kCarTolerance;
>> 340 pMax+=kCarTolerance;
>> 341
>> 342 return true;
309 } 343 }
>> 344 else // Transformed cutted sphere
>> 345 {
>> 346 G4int i,j,noEntries,noBetweenSections;
>> 347 G4bool existsAfterClip=false;
310 348
311 // Phi boundaries : Do not check if it has << 349 // Calculate rotated vertex coordinates
312 350
313 if ( !fFullPhiSphere && (rho2 != 0.0) ) // << 351 G4ThreeVectorList* vertices;
314 { << 352 G4int noPolygonVertices ;
315 pPhi = std::atan2(p.y(),p.x()) ; << 353 vertices=CreateRotatedVertices(pTransform,noPolygonVertices);
316 354
317 if ( pPhi < fSPhi - halfAngTolerance << 355 pMin=+kInfinity;
318 else if ( pPhi > ePhi + halfAngTolerance ) << 356 pMax=-kInfinity;
319 357
320 if ( (pPhi < fSPhi - halfAngTolerance) << 358 noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections
321 || (pPhi > ePhi + halfAngTolerance) ) << 359 noBetweenSections=noEntries-noPolygonVertices;
322 360
323 else if (in == kInside) // else it's kSur << 361 G4ThreeVectorList ThetaPolygon ;
>> 362 for (i=0;i<noEntries;i+=noPolygonVertices)
324 { 363 {
325 if ( (pPhi < fSPhi + halfAngTolerance) << 364 for(j=0;j<(noPolygonVertices/2)-1;j++)
326 || (pPhi > ePhi - halfAngTolerance) ) << 365 {
>> 366 ThetaPolygon.push_back((*vertices)[i+j]) ;
>> 367 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
>> 368 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;
>> 369 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;
>> 370 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
>> 371 ThetaPolygon.clear() ;
>> 372 }
327 } 373 }
328 } << 374 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
329 <<
330 // Theta bondaries <<
331 <<
332 if ( ((rho2 != 0.0) || (p.z() != 0.0)) && (! <<
333 { <<
334 rho = std::sqrt(rho2); <<
335 pTheta = std::atan2(rho,p.z()); <<
336 <<
337 if ( in == kInside ) <<
338 { 375 {
339 if ( ((fSTheta > 0.0) && (pTheta < fSThe << 376 for(j=0;j<noPolygonVertices-1;j++)
340 || ((eTheta < pi) && (pTheta > eTheta <<
341 { 377 {
342 if ( (( (fSTheta>0.0)&&(pTheta>=fSThet << 378 ThetaPolygon.push_back((*vertices)[i+j]) ;
343 || (fSTheta == 0.0) ) << 379 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
344 && ((eTheta==pi)||(pTheta <= eTheta << 380 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;
345 { << 381 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;
346 in = kSurface; << 382 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
347 } << 383 ThetaPolygon.clear() ;
348 else <<
349 { <<
350 in = kOutside; <<
351 } <<
352 } 384 }
>> 385 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;
>> 386 ThetaPolygon.push_back((*vertices)[i]) ;
>> 387 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;
>> 388 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;
>> 389 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
>> 390 ThetaPolygon.clear() ;
>> 391 }
>> 392
>> 393 if (pMin!=kInfinity || pMax!=-kInfinity)
>> 394 {
>> 395 existsAfterClip=true;
>> 396
>> 397 // Add 2*tolerance to avoid precision troubles
>> 398 //
>> 399 pMin-=kCarTolerance;
>> 400 pMax+=kCarTolerance;
353 } 401 }
354 else 402 else
355 { 403 {
356 if ( ((fSTheta > 0.0)&&(pTheta < fSThe << 404 // Check for case where completely enveloping clipping volume
357 ||((eTheta < pi )&&(pTheta > eThet << 405 // If point inside then we are confident that the solid completely
358 { << 406 // envelopes the clipping volume. Hence set min/max extents according
359 in = kOutside; << 407 // to clipping volume extents along the specified axis.
>> 408
>> 409 G4ThreeVector clipCentre(
>> 410 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 411 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 412 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
>> 413
>> 414 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
>> 415 {
>> 416 existsAfterClip=true;
>> 417 pMin=pVoxelLimit.GetMinExtent(pAxis);
>> 418 pMax=pVoxelLimit.GetMaxExtent(pAxis);
360 } 419 }
361 } 420 }
>> 421 delete vertices;
>> 422 return existsAfterClip;
362 } 423 }
363 return in; <<
364 } 424 }
365 425
366 ////////////////////////////////////////////// << 426 ///////////////////////////////////////////////////////////////////////////
367 // 427 //
368 // Return unit normal of surface closest to p << 428 // Return whether point inside/outside/on surface
369 // - note if point on z axis, ignore phi divid << 429 // Split into radius, phi, theta checks
370 // - unsafe if point close to z axis a rmin=0 << 430 // Each check modifies `in', or returns as approprate
371 431
372 G4ThreeVector G4Sphere::SurfaceNormal( const G << 432 EInside G4Sphere::Inside( const G4ThreeVector& p ) const
373 { 433 {
374 G4int noSurfaces = 0; << 434 G4double rho,rho2,rad2,tolRMin,tolRMax;
375 G4double rho, rho2, radius, pTheta, pPhi=0.; << 435 G4double pPhi,pTheta;
376 G4double distRMin = kInfinity; << 436 EInside in;
377 G4double distSPhi = kInfinity, distEPhi = kI <<
378 G4double distSTheta = kInfinity, distETheta <<
379 G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0., <<
380 G4ThreeVector norm, sumnorm(0.,0.,0.); <<
381 <<
382 rho2 = p.x()*p.x()+p.y()*p.y(); <<
383 radius = std::sqrt(rho2+p.z()*p.z()); <<
384 rho = std::sqrt(rho2); <<
385 437
386 G4double distRMax = std::fabs(radius-fRma << 438 rho2 = p.x()*p.x() + p.y()*p.y() ;
387 if (fRmin != 0.0) distRMin = std::fabs(radi << 439 rad2 = rho2 + p.z()*p.z() ;
388 440
389 if ( (rho != 0.0) && !fFullSphere ) << 441 // G4double rad = sqrt(rad2);
>> 442 // Check radial surfaces
>> 443 // sets `in'
>> 444
>> 445 if ( fRmin ) tolRMin = fRmin + kRadTolerance*0.5;
>> 446 else tolRMin = 0 ;
>> 447
>> 448 tolRMax = fRmax - kRadTolerance*0.5 ;
>> 449
>> 450 if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin ) in = kInside ;
>> 451 else
390 { 452 {
391 pPhi = std::atan2(p.y(),p.x()); << 453 tolRMax = fRmax + kRadTolerance*0.5 ;
>> 454 tolRMin = fRmin - kRadTolerance*0.5 ;
392 455
393 if (pPhi < fSPhi-halfAngTolerance) { p << 456 if ( tolRMin < 0.0 ) tolRMin = 0.0 ;
394 else if (pPhi > ePhi+halfAngTolerance) { p << 457
>> 458 if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin ) in = kSurface ;
>> 459 else return in = kOutside ;
395 } 460 }
396 if ( !fFullPhiSphere ) <<
397 { <<
398 if ( rho != 0.0 ) <<
399 { <<
400 distSPhi = std::fabs( pPhi-fSPhi ); <<
401 distEPhi = std::fabs( pPhi-ePhi ); <<
402 } <<
403 else if( fRmin == 0.0 ) <<
404 { <<
405 distSPhi = 0.; <<
406 distEPhi = 0.; <<
407 } <<
408 nPs = G4ThreeVector(sinSPhi,-cosSPhi,0); <<
409 nPe = G4ThreeVector(-sinEPhi,cosEPhi,0); <<
410 } <<
411 if ( !fFullThetaSphere ) <<
412 { <<
413 if ( rho != 0.0 ) <<
414 { <<
415 pTheta = std::atan2(rho,p.z()); <<
416 distSTheta = std::fabs(pTheta-fSTheta); <<
417 distETheta = std::fabs(pTheta-eTheta); <<
418 461
419 nTs = G4ThreeVector(-cosSTheta*p.x()/rho << 462 // Phi boundaries : Do not check if it has no phi boundary!
420 -cosSTheta*p.y()/rho << 463 // (in != kOutside)
421 sinSTheta <<
422 464
423 nTe = G4ThreeVector( cosETheta*p.x()/rho << 465 if ( ( fDPhi < 2*M_PI - kAngTolerance ) &&
424 cosETheta*p.y()/rho << 466 ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
425 -sinETheta << 467 {
426 } << 468 pPhi = atan2(p.y(),p.x()) ;
427 else if( fRmin == 0.0 ) << 469 if ( pPhi < 0.0) pPhi += 2*M_PI ; // 0<=pPhi<2pi
>> 470 if ( fSPhi >= 0.0 )
428 { 471 {
429 if ( fSTheta != 0.0 ) << 472 if (in == kInside)
430 { 473 {
431 distSTheta = 0.; << 474 if ( (pPhi < fSPhi + kAngTolerance*0.5)
432 nTs = G4ThreeVector(0.,0.,-1.); << 475 || (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) )
>> 476 {
>> 477 // Not `inside' tolerant bounds
>> 478 //
>> 479 if ( (pPhi >= fSPhi - kAngTolerance*0.5)
>> 480 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
>> 481 {
>> 482 in = kSurface ;
>> 483 }
>> 484 else
>> 485 {
>> 486 return in = kOutside ;
>> 487 }
>> 488 }
433 } 489 }
434 if ( eTheta < pi ) << 490 else // in==kSurface
435 { 491 {
436 distETheta = 0.; << 492 if ( (pPhi < fSPhi - kAngTolerance*0.5)
437 nTe = G4ThreeVector(0.,0.,1.); << 493 || (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) ) // && bug216
>> 494 {
>> 495 return in = kOutside ;
>> 496 }
438 } 497 }
439 } 498 }
440 } << 499 else
441 if( radius != 0.0 ) { nR = G4ThreeVector(p. <<
442 <<
443 if( distRMax <= halfCarTolerance ) <<
444 { <<
445 ++noSurfaces; <<
446 sumnorm += nR; <<
447 } <<
448 if( (fRmin != 0.0) && (distRMin <= halfCarTo <<
449 { <<
450 ++noSurfaces; <<
451 sumnorm -= nR; <<
452 } <<
453 if( !fFullPhiSphere ) <<
454 { <<
455 if (distSPhi <= halfAngTolerance) <<
456 { <<
457 ++noSurfaces; <<
458 sumnorm += nPs; <<
459 } <<
460 if (distEPhi <= halfAngTolerance) <<
461 { 500 {
462 ++noSurfaces; << 501 if ( pPhi < fSPhi + 2*M_PI ) pPhi += 2*M_PI ;
463 sumnorm += nPe; << 502
>> 503 if ( (pPhi < fSPhi + 2*M_PI + kAngTolerance*0.5)
>> 504 || (pPhi > fSPhi + fDPhi + 2*M_PI - kAngTolerance*0.5) )
>> 505 {
>> 506 // Not `inside' tolerant bounds
>> 507 //
>> 508 if ( (pPhi >= fSPhi + 2*M_PI - kAngTolerance*0.5)
>> 509 && (pPhi <= fSPhi + fDPhi + 2*M_PI + kAngTolerance*0.5) )
>> 510 {
>> 511 in = kSurface ;
>> 512 }
>> 513 else
>> 514 {
>> 515 return in = kOutside ;
>> 516 }
>> 517 }
464 } 518 }
465 } 519 }
466 if ( !fFullThetaSphere ) << 520 // Theta bondaries
>> 521 // (in!=kOutside)
>> 522
>> 523 if ( (rho2 || p.z()) && fDTheta < M_PI - kAngTolerance*0.5 )
467 { 524 {
468 if ((distSTheta <= halfAngTolerance) && (f << 525 rho = sqrt(rho2);
>> 526 pTheta = atan2(rho,p.z());
>> 527
>> 528 if ( in == kInside )
469 { 529 {
470 ++noSurfaces; << 530 if ( (pTheta < fSTheta + kAngTolerance*0.5)
471 if ((radius <= halfCarTolerance) && fFul << 531 || (pTheta > fSTheta + fDTheta - kAngTolerance*0.5) )
472 else << 532 {
>> 533 if ( (pTheta >= fSTheta - kAngTolerance*0.5)
>> 534 && (pTheta <= fSTheta + fDTheta + kAngTolerance*0.5) )
>> 535 {
>> 536 in = kSurface ;
>> 537 }
>> 538 else
>> 539 {
>> 540 in = kOutside ;
>> 541 }
>> 542 }
473 } 543 }
474 if ((distETheta <= halfAngTolerance) && (e << 544 else
475 { 545 {
476 ++noSurfaces; << 546 if ( (pTheta < fSTheta - kAngTolerance*0.5)
477 if ((radius <= halfCarTolerance) && fFul << 547 || (pTheta > fSTheta + fDTheta + kAngTolerance*0.5) )
478 else << 548 {
479 if(sumnorm.z() == 0.) { sumnorm += nZ; << 549 in = kOutside ;
>> 550 }
480 } 551 }
481 } 552 }
482 if ( noSurfaces == 0 ) << 553 return in;
483 { <<
484 #ifdef G4CSGDEBUG <<
485 G4Exception("G4Sphere::SurfaceNormal(p)", <<
486 JustWarning, "Point p is not o <<
487 #endif <<
488 norm = ApproxSurfaceNormal(p); <<
489 } <<
490 else if ( noSurfaces == 1 ) { norm = sumnor <<
491 else { norm = sumnor <<
492 return norm; <<
493 } 554 }
494 555
495 <<
496 ////////////////////////////////////////////// 556 /////////////////////////////////////////////////////////////////////
497 // 557 //
498 // Algorithm for SurfaceNormal() following the << 558 // Return unit normal of surface closest to p
499 // for points not on the surface << 559 // - note if point on z axis, ignore phi divided sides
>> 560 // - unsafe if point close to z axis a rmin=0 - no explicit checks
500 561
501 G4ThreeVector G4Sphere::ApproxSurfaceNormal( c << 562 G4ThreeVector G4Sphere::SurfaceNormal( const G4ThreeVector& p ) const
502 { 563 {
503 ENorm side; 564 ENorm side;
504 G4ThreeVector norm; 565 G4ThreeVector norm;
505 G4double rho,rho2,radius,pPhi,pTheta; << 566 G4double rho,rho2,rad,pPhi,pTheta;
506 G4double distRMin,distRMax,distSPhi,distEPhi 567 G4double distRMin,distRMax,distSPhi,distEPhi,
507 distSTheta,distETheta,distMin; 568 distSTheta,distETheta,distMin;
508 569
509 rho2=p.x()*p.x()+p.y()*p.y(); 570 rho2=p.x()*p.x()+p.y()*p.y();
510 radius=std::sqrt(rho2+p.z()*p.z()); << 571 rad=sqrt(rho2+p.z()*p.z());
511 rho=std::sqrt(rho2); << 572 rho=sqrt(rho2);
512 573
513 // 574 //
514 // Distance to r shells 575 // Distance to r shells
515 // 576 //
516 577
517 distRMax=std::fabs(radius-fRmax); << 578 distRMax=fabs(rad-fRmax);
518 if (fRmin != 0.0) << 579 if (fRmin)
519 { 580 {
520 distRMin=std::fabs(radius-fRmin); << 581 distRMin=fabs(rad-fRmin);
521 << 582
522 if (distRMin<distRMax) 583 if (distRMin<distRMax)
523 { 584 {
524 distMin=distRMin; 585 distMin=distRMin;
525 side=kNRMin; 586 side=kNRMin;
526 } 587 }
527 else 588 else
528 { << 589 {
529 distMin=distRMax; 590 distMin=distRMax;
530 side=kNRMax; 591 side=kNRMax;
531 } 592 }
532 } 593 }
533 else 594 else
534 { 595 {
535 distMin=distRMax; 596 distMin=distRMax;
536 side=kNRMax; 597 side=kNRMax;
537 } 598 }
538 599
539 // 600 //
540 // Distance to phi planes 601 // Distance to phi planes
541 // 602 //
542 // Protected against (0,0,z) << 603 if (fDPhi<2.0*M_PI&&rho)
543 <<
544 pPhi = std::atan2(p.y(),p.x()); <<
545 if (pPhi<0) { pPhi += twopi; } <<
546 <<
547 if (!fFullPhiSphere && (rho != 0.0)) <<
548 { 604 {
>> 605 // Protected against (0,0,z) (above)
>> 606 //
>> 607 pPhi=atan2(p.y(),p.x());
>> 608 if (pPhi<0) pPhi+=2*M_PI;
549 if (fSPhi<0) 609 if (fSPhi<0)
550 { 610 {
551 distSPhi=std::fabs(pPhi-(fSPhi+twopi))*r << 611 distSPhi=fabs(pPhi-(fSPhi+2.0*M_PI))*rho;
552 } 612 }
553 else 613 else
554 { 614 {
555 distSPhi=std::fabs(pPhi-fSPhi)*rho; << 615 distSPhi=fabs(pPhi-fSPhi)*rho;
556 } 616 }
557 617
558 distEPhi=std::fabs(pPhi-fSPhi-fDPhi)*rho; << 618 distEPhi=fabs(pPhi-fSPhi-fDPhi)*rho;
559 619
560 // Find new minimum 620 // Find new minimum
561 // 621 //
562 if (distSPhi<distEPhi) 622 if (distSPhi<distEPhi)
563 { 623 {
564 if (distSPhi<distMin) 624 if (distSPhi<distMin)
565 { 625 {
566 distMin = distSPhi; << 626 distMin=distSPhi;
567 side = kNSPhi; << 627 side=kNSPhi;
568 } 628 }
569 } 629 }
570 else 630 else
571 { 631 {
572 if (distEPhi<distMin) 632 if (distEPhi<distMin)
573 { 633 {
574 distMin = distEPhi; << 634 distMin=distEPhi;
575 side = kNEPhi; << 635 side=kNEPhi;
576 } 636 }
577 } 637 }
578 } 638 }
579 639
580 // 640 //
581 // Distance to theta planes 641 // Distance to theta planes
582 // 642 //
583 643
584 if (!fFullThetaSphere && (radius != 0.0)) << 644 if (fDTheta<M_PI&&rad)
585 { 645 {
586 pTheta=std::atan2(rho,p.z()); << 646 pTheta=atan2(rho,p.z());
587 distSTheta=std::fabs(pTheta-fSTheta)*radiu << 647 distSTheta=fabs(pTheta-fSTheta)*rad;
588 distETheta=std::fabs(pTheta-fSTheta-fDThet << 648 distETheta=fabs(pTheta-fSTheta-fDTheta)*rad;
589 649
590 // Find new minimum 650 // Find new minimum
591 // 651 //
592 if (distSTheta<distETheta) 652 if (distSTheta<distETheta)
593 { 653 {
594 if (distSTheta<distMin) 654 if (distSTheta<distMin)
595 { 655 {
596 distMin = distSTheta ; 656 distMin = distSTheta ;
597 side = kNSTheta ; 657 side = kNSTheta ;
598 } 658 }
599 } 659 }
600 else 660 else
601 { 661 {
602 if (distETheta<distMin) 662 if (distETheta<distMin)
603 { 663 {
604 distMin = distETheta ; 664 distMin = distETheta ;
605 side = kNETheta ; 665 side = kNETheta ;
606 } 666 }
607 } 667 }
608 } 668 }
609 669
610 switch (side) 670 switch (side)
611 { 671 {
612 case kNRMin: // Inner radius 672 case kNRMin: // Inner radius
613 norm=G4ThreeVector(-p.x()/radius,-p.y()/ << 673 norm=G4ThreeVector(-p.x()/rad,-p.y()/rad,-p.z()/rad);
614 break; 674 break;
615 case kNRMax: // Outer radius 675 case kNRMax: // Outer radius
616 norm=G4ThreeVector(p.x()/radius,p.y()/ra << 676 norm=G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);
617 break; 677 break;
618 case kNSPhi: 678 case kNSPhi:
619 norm=G4ThreeVector(sinSPhi,-cosSPhi,0); << 679 norm=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0);
620 break; 680 break;
621 case kNEPhi: 681 case kNEPhi:
622 norm=G4ThreeVector(-sinEPhi,cosEPhi,0); << 682 norm=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0);
623 break; 683 break;
>> 684
624 case kNSTheta: 685 case kNSTheta:
625 norm=G4ThreeVector(-cosSTheta*std::cos(p << 686 norm=G4ThreeVector(-cos(fSTheta)*sin(fSPhi),
626 -cosSTheta*std::sin(p << 687 -cos(fSTheta)*cos(fSPhi),sin(fSTheta));
627 sinSTheta <<
628 break; 688 break;
629 case kNETheta: 689 case kNETheta:
630 norm=G4ThreeVector( cosETheta*std::cos(p << 690 norm=G4ThreeVector(-cos(fSTheta+fDTheta)*cos(fSPhi+fDPhi),
631 cosETheta*std::sin(p << 691 -cos(fSTheta+fDTheta)*sin(fSPhi+fDPhi),
632 -sinETheta << 692 -sin(fSTheta+fDTheta));
633 break; 693 break;
634 default: // Should never reach th << 694 default:
635 DumpInfo(); 695 DumpInfo();
636 G4Exception("G4Sphere::ApproxSurfaceNorm << 696 G4Exception("G4Sphere::SurfaceNormal() - Logic error");
637 "GeomSolids1002", JustWarnin << 697 break;
638 "Undefined side for valid su << 698 } // end case
639 break; <<
640 } <<
641 699
642 return norm; 700 return norm;
643 } 701 }
644 702
645 ////////////////////////////////////////////// 703 ///////////////////////////////////////////////////////////////////////////////
646 // 704 //
647 // Calculate distance to shape from outside, a 705 // Calculate distance to shape from outside, along normalised vector
648 // - return kInfinity if no intersection, or i 706 // - return kInfinity if no intersection, or intersection distance <= tolerance
649 // 707 //
650 // -> If point is outside outer radius, comput 708 // -> If point is outside outer radius, compute intersection with rmax
651 // - if no intersection return 709 // - if no intersection return
652 // - if valid phi,theta return interse 710 // - if valid phi,theta return intersection Dist
653 // 711 //
654 // -> If shell, compute intersection with inne 712 // -> If shell, compute intersection with inner radius, taking largest +ve root
655 // - if valid phi,theta, save intersect 713 // - if valid phi,theta, save intersection
656 // 714 //
657 // -> If phi segmented, compute intersection w 715 // -> If phi segmented, compute intersection with phi half planes
658 // - if valid intersection(r,theta), re 716 // - if valid intersection(r,theta), return smallest intersection of
659 // inner shell & phi intersection 717 // inner shell & phi intersection
660 // 718 //
661 // -> If theta segmented, compute intersection 719 // -> If theta segmented, compute intersection with theta cones
662 // - if valid intersection(r,phi), retu 720 // - if valid intersection(r,phi), return smallest intersection of
663 // inner shell & theta intersection 721 // inner shell & theta intersection
664 // 722 //
665 // 723 //
666 // NOTE: 724 // NOTE:
667 // - `if valid' (above) implies tolerant check 725 // - `if valid' (above) implies tolerant checking of intersection points
668 // 726 //
669 // OPT: 727 // OPT:
670 // Move tolIO/ORmin/RMax2 precalcs to where th 728 // Move tolIO/ORmin/RMax2 precalcs to where they are needed -
671 // not required for most cases. 729 // not required for most cases.
672 // Avoid atan2 for non theta cut G4Sphere. 730 // Avoid atan2 for non theta cut G4Sphere.
673 731
674 G4double G4Sphere::DistanceToIn( const G4Three 732 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p,
675 const G4Three 733 const G4ThreeVector& v ) const
676 { 734 {
677 G4double snxt = kInfinity ; // snxt = d 735 G4double snxt = kInfinity ; // snxt = default return value
>> 736
678 G4double rho2, rad2, pDotV2d, pDotV3d, pThet 737 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
679 G4double tolSTheta=0., tolETheta=0. ; <<
680 const G4double dRmax = 100.*fRmax; <<
681 738
682 const G4double halfRmaxTolerance = fRmaxTole << 739 G4double tolIRMin2, tolORMin2, tolORMax2, tolIRMax2 ;
683 const G4double halfRminTolerance = fRminTole << 740 G4double tolSTheta=0., tolETheta=0. ;
684 const G4double tolORMin2 = (fRmin>halfRminTo <<
685 ? (fRmin-halfRminTolerance)*(fR <<
686 const G4double tolIRMin2 = <<
687 (fRmin+halfRminTolerance)*(fRmi <<
688 const G4double tolORMax2 = <<
689 (fRmax+halfRmaxTolerance)*(fRma <<
690 const G4double tolIRMax2 = <<
691 (fRmax-halfRmaxTolerance)*(fRma <<
692 741
693 // Intersection point 742 // Intersection point
694 // << 743
695 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTh 744 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
696 745
697 // Phi intersection 746 // Phi intersection
698 // <<
699 G4double Comp ; <<
700 747
701 // Phi precalcs << 748 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi , Comp ;
702 // << 749
>> 750 // Phi flag and precalcs
>> 751
>> 752 G4bool segPhi ;
>> 753 G4double hDPhi, hDPhiOT, hDPhiIT, cPhi, sinCPhi=0., cosCPhi=0. ;
>> 754 G4double cosHDPhiOT=0., cosHDPhiIT=0. ;
703 G4double Dist, cosPsi ; 755 G4double Dist, cosPsi ;
704 756
705 // Theta precalcs << 757 G4bool segTheta ; // Theta flag and precals
706 // << 758 G4double tanSTheta, tanETheta ;
>> 759 G4double tanSTheta2, tanETheta2 ;
707 G4double dist2STheta, dist2ETheta ; 760 G4double dist2STheta, dist2ETheta ;
708 G4double t1, t2, b, c, d2, d, sd = kInfinity << 761 G4double t1, t2, b, c, d2, d, s = kInfinity ;
709 762
710 // General Precalcs 763 // General Precalcs
711 // << 764
712 rho2 = p.x()*p.x() + p.y()*p.y() ; 765 rho2 = p.x()*p.x() + p.y()*p.y() ;
713 rad2 = rho2 + p.z()*p.z() ; 766 rad2 = rho2 + p.z()*p.z() ;
714 pTheta = std::atan2(std::sqrt(rho2),p.z()) ; << 767 pTheta = atan2(sqrt(rho2),p.z()) ;
715 768
716 pDotV2d = p.x()*v.x() + p.y()*v.y() ; 769 pDotV2d = p.x()*v.x() + p.y()*v.y() ;
717 pDotV3d = pDotV2d + p.z()*v.z() ; 770 pDotV3d = pDotV2d + p.z()*v.z() ;
718 771
719 // Theta precalcs << 772 // Radial Precalcs
720 // << 773
721 if (!fFullThetaSphere) << 774 if (fRmin > kRadTolerance*0.5)
>> 775 {
>> 776 tolORMin2=(fRmin-kRadTolerance*0.5)*(fRmin-kRadTolerance*0.5);
>> 777 }
>> 778 else
722 { 779 {
723 tolSTheta = fSTheta - halfAngTolerance ; << 780 tolORMin2 = 0 ;
724 tolETheta = eTheta + halfAngTolerance ; << 781 }
>> 782 tolIRMin2 = (fRmin+kRadTolerance*0.5)*(fRmin+kRadTolerance*0.5) ;
>> 783 tolORMax2 = (fRmax+kRadTolerance*0.5)*(fRmax+kRadTolerance*0.5) ;
>> 784 tolIRMax2 = (fRmax-kRadTolerance*0.5)*(fRmax-kRadTolerance*0.5) ;
725 785
726 // Special case rad2 = 0 comparing with di << 786 // Set phi divided flag and precalcs
727 // << 787
728 if ((rad2!=0.0) || (fRmin!=0.0)) << 788 if (fDPhi < 2.0*M_PI)
729 { << 789 {
730 // Keep going for computation of distanc << 790 segPhi = true ;
731 } << 791 hDPhi = 0.5*fDPhi ; // half delta phi
732 else // Positioned on the sphere's origin << 792 cPhi = fSPhi + hDPhi ;
733 { << 793
734 G4double vTheta = std::atan2(std::sqrt(v << 794 hDPhiOT = hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi
735 if ( (vTheta < tolSTheta) || (vTheta > t << 795 hDPhiIT = hDPhi-0.5*kAngTolerance;
736 { << 796
737 return snxt ; // kInfinity << 797 sinCPhi = sin(cPhi) ;
738 } << 798 cosCPhi = cos(cPhi) ;
739 return snxt = 0.0 ; << 799 cosHDPhiOT = cos(hDPhiOT) ;
740 } << 800 cosHDPhiIT = cos(hDPhiIT) ;
>> 801 }
>> 802 else
>> 803 {
>> 804 segPhi = false ;
>> 805 }
>> 806
>> 807 // Theta precalcs
>> 808
>> 809 if (fDTheta < M_PI )
>> 810 {
>> 811 segTheta = true ;
>> 812 tolSTheta = fSTheta - kAngTolerance*0.5 ;
>> 813 tolETheta = fSTheta + fDTheta + kAngTolerance*0.5 ;
>> 814 }
>> 815 else
>> 816 {
>> 817 segTheta = false ;
741 } 818 }
742 819
743 // Outer spherical shell intersection 820 // Outer spherical shell intersection
744 // - Only if outside tolerant fRmax 821 // - Only if outside tolerant fRmax
745 // - Check for if inside and outer G4Sphere 822 // - Check for if inside and outer G4Sphere heading through solid (-> 0)
746 // - No intersect -> no intersection with G4 823 // - No intersect -> no intersection with G4Sphere
747 // 824 //
748 // Shell eqn: x^2+y^2+z^2=RSPH^2 825 // Shell eqn: x^2+y^2+z^2=RSPH^2
749 // 826 //
750 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 827 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
751 // 828 //
752 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+ << 829 // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
753 // => rad2 +2sd(pDotV3d) + << 830 // => rad2 +2s(pDotV3d) +s^2 =R^2
754 // 831 //
755 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2 << 832 // => s=-pDotV3d+-sqrt(pDotV3d^2-(rad2-R^2))
756 833
757 c = rad2 - fRmax*fRmax ; 834 c = rad2 - fRmax*fRmax ;
758 835
759 if (c > fRmaxTolerance*fRmax) << 836 if (c > kRadTolerance*fRmax)
760 { 837 {
761 // If outside tolerant boundary of outer G << 838 // If outside toleranct boundary of outer G4Sphere
762 // [should be std::sqrt(rad2)-fRmax > half << 839 // [should be sqrt(rad2)-fRmax > kRadTolerance*0.5]
763 840
764 d2 = pDotV3d*pDotV3d - c ; 841 d2 = pDotV3d*pDotV3d - c ;
765 842
766 if ( d2 >= 0 ) 843 if ( d2 >= 0 )
767 { 844 {
768 sd = -pDotV3d - std::sqrt(d2) ; << 845 s = -pDotV3d - sqrt(d2) ;
769 846
770 if (sd >= 0 ) << 847 if (s >= 0 )
771 { 848 {
772 if ( sd>dRmax ) // Avoid rounding erro << 849 xi = p.x() + s*v.x() ;
773 { // 64 bits systems. Sp << 850 yi = p.y() + s*v.y() ;
774 G4double fTerm = sd-std::fmod(sd,dRm << 851 rhoi = sqrt(xi*xi + yi*yi) ;
775 sd = fTerm + DistanceToIn(p+fTerm*v, <<
776 } <<
777 xi = p.x() + sd*v.x() ; <<
778 yi = p.y() + sd*v.y() ; <<
779 rhoi = std::sqrt(xi*xi + yi*yi) ; <<
780 852
781 if (!fFullPhiSphere && (rhoi != 0.0)) << 853 if (segPhi && rhoi) // Check phi intersection
782 { 854 {
783 cosPsi = (xi*cosCPhi + yi*sinCPhi)/r 855 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
784 856
785 if (cosPsi >= cosHDPhiOT) 857 if (cosPsi >= cosHDPhiOT)
786 { 858 {
787 if (!fFullThetaSphere) // Check << 859 if (segTheta) // Check theta intersection
788 { 860 {
789 zi = p.z() + sd*v.z() ; << 861 zi = p.z() + s*v.z() ;
790 862
791 // rhoi & zi can never both be 0 863 // rhoi & zi can never both be 0
792 // (=>intersect at origin =>fRma 864 // (=>intersect at origin =>fRmax=0)
793 // 865 //
794 iTheta = std::atan2(rhoi,zi) ; << 866 iTheta = atan2(rhoi,zi) ;
795 if ( (iTheta >= tolSTheta) && (i 867 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
796 { 868 {
797 return snxt = sd ; << 869 return snxt = s ;
798 } 870 }
799 } 871 }
800 else 872 else
801 { 873 {
802 return snxt=sd; << 874 return snxt=s;
803 } 875 }
804 } 876 }
805 } 877 }
806 else 878 else
807 { 879 {
808 if (!fFullThetaSphere) // Check t << 880 if (segTheta) // Check theta intersection
809 { 881 {
810 zi = p.z() + sd*v.z() ; << 882 zi = p.z() + s*v.z() ;
811 883
812 // rhoi & zi can never both be 0 884 // rhoi & zi can never both be 0
813 // (=>intersect at origin => fRmax 885 // (=>intersect at origin => fRmax=0 !)
814 // 886 //
815 iTheta = std::atan2(rhoi,zi) ; << 887 iTheta = atan2(rhoi,zi) ;
816 if ( (iTheta >= tolSTheta) && (iTh 888 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
817 { 889 {
818 return snxt=sd; << 890 return snxt=s;
819 } 891 }
820 } 892 }
821 else 893 else
822 { 894 {
823 return snxt = sd; << 895 return snxt = s ;
824 } 896 }
825 } << 897 }
826 } 898 }
827 } 899 }
828 else // No intersection with G4Sphere 900 else // No intersection with G4Sphere
829 { 901 {
830 return snxt=kInfinity; 902 return snxt=kInfinity;
831 } 903 }
832 } 904 }
833 else 905 else
834 { 906 {
835 // Inside outer radius 907 // Inside outer radius
836 // check not inside, and heading through G 908 // check not inside, and heading through G4Sphere (-> 0 to in)
837 909
838 d2 = pDotV3d*pDotV3d - c ; << 910 if (rad2 > tolIRMin2 && pDotV3d < 0 )
839 <<
840 if ( (rad2 > tolIRMax2) <<
841 && ( (d2 >= fRmaxTolerance*fRmax) && (pD <<
842 { 911 {
843 if (!fFullPhiSphere) << 912 if (segPhi)
844 { 913 {
845 // Use inner phi tolerant boundary -> 914 // Use inner phi tolerant boundary -> if on tolerant
846 // phi boundaries, phi intersect code 915 // phi boundaries, phi intersect code handles leaving/entering checks
847 916
848 cosPsi = (p.x()*cosCPhi + p.y()*sinCPh << 917 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/sqrt(rho2) ;
849 918
850 if (cosPsi>=cosHDPhiIT) 919 if (cosPsi>=cosHDPhiIT)
851 { << 920 {
852 // inside radii, delta r -ve, inside 921 // inside radii, delta r -ve, inside phi
853 922
854 if ( !fFullThetaSphere ) << 923 if (segTheta)
855 { 924 {
856 if ( (pTheta >= tolSTheta + kAngTo 925 if ( (pTheta >= tolSTheta + kAngTolerance)
857 && (pTheta <= tolETheta - kAngTo 926 && (pTheta <= tolETheta - kAngTolerance) )
858 { 927 {
859 return snxt=0; 928 return snxt=0;
860 } 929 }
861 } 930 }
862 else // strictly inside Theta in 931 else // strictly inside Theta in both cases
863 { 932 {
864 return snxt=0; 933 return snxt=0;
865 } 934 }
866 } 935 }
867 } 936 }
868 else 937 else
869 { 938 {
870 if ( !fFullThetaSphere ) << 939 if ( segTheta )
871 { 940 {
872 if ( (pTheta >= tolSTheta + kAngTole 941 if ( (pTheta >= tolSTheta + kAngTolerance)
873 && (pTheta <= tolETheta - kAngTole 942 && (pTheta <= tolETheta - kAngTolerance) )
874 { 943 {
875 return snxt=0; 944 return snxt=0;
876 } 945 }
877 } 946 }
878 else // strictly inside Theta in bot 947 else // strictly inside Theta in both cases
879 { 948 {
880 return snxt=0; 949 return snxt=0;
881 } 950 }
882 } 951 }
883 } 952 }
884 } 953 }
885 954
886 // Inner spherical shell intersection 955 // Inner spherical shell intersection
887 // - Always farthest root, because would hav 956 // - Always farthest root, because would have passed through outer
888 // surface first. 957 // surface first.
889 // - Tolerant check if travelling through so << 958 // - Tolerant check for if travelling through solid
890 959
891 if (fRmin != 0.0) << 960 if (fRmin)
892 { 961 {
893 c = rad2 - fRmin*fRmin ; << 962 c = rad2 - fRmin*fRmin ;
894 d2 = pDotV3d*pDotV3d - c ; <<
895 963
896 // Within tolerance inner radius of inner 964 // Within tolerance inner radius of inner G4Sphere
897 // Check for immediate entry/already insid 965 // Check for immediate entry/already inside and travelling outwards
898 966
899 if ( (c > -halfRminTolerance) && (rad2 < t << 967 if (c >- kRadTolerance*0.5 && pDotV3d >= 0 && rad2 < tolIRMax2 )
900 && ( (d2 < fRmin*kCarTolerance) || (pDot <<
901 { 968 {
902 if ( !fFullPhiSphere ) << 969 if (segPhi)
903 { 970 {
904 // Use inner phi tolerant boundary -> 971 // Use inner phi tolerant boundary -> if on tolerant
905 // phi boundaries, phi intersect code 972 // phi boundaries, phi intersect code handles leaving/entering checks
906 973
907 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi) << 974 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/sqrt(rho2) ;
908 if (cosPsi >= cosHDPhiIT) 975 if (cosPsi >= cosHDPhiIT)
909 { << 976 {
910 // inside radii, delta r -ve, inside 977 // inside radii, delta r -ve, inside phi
911 // 978 //
912 if ( !fFullThetaSphere ) << 979 if (segTheta)
913 { 980 {
914 if ( (pTheta >= tolSTheta + kAngTo 981 if ( (pTheta >= tolSTheta + kAngTolerance)
915 && (pTheta <= tolETheta - kAngTo 982 && (pTheta <= tolETheta - kAngTolerance) )
916 { 983 {
917 return snxt=0; 984 return snxt=0;
918 } 985 }
919 } 986 }
920 else 987 else
921 { 988 {
922 return snxt = 0 ; 989 return snxt = 0 ;
923 } 990 }
924 } 991 }
925 } 992 }
926 else 993 else
927 { 994 {
928 if ( !fFullThetaSphere ) << 995 if (segTheta)
929 { 996 {
930 if ( (pTheta >= tolSTheta + kAngTole 997 if ( (pTheta >= tolSTheta + kAngTolerance)
931 && (pTheta <= tolETheta - kAngTole 998 && (pTheta <= tolETheta - kAngTolerance) )
932 { 999 {
933 return snxt = 0 ; 1000 return snxt = 0 ;
934 } 1001 }
935 } 1002 }
936 else 1003 else
937 { 1004 {
938 return snxt=0; 1005 return snxt=0;
939 } 1006 }
940 } 1007 }
941 } 1008 }
942 else // Not special tolerant case 1009 else // Not special tolerant case
943 { 1010 {
>> 1011 d2 = pDotV3d*pDotV3d - c ;
944 if (d2 >= 0) 1012 if (d2 >= 0)
945 { 1013 {
946 sd = -pDotV3d + std::sqrt(d2) ; << 1014 s = -pDotV3d + sqrt(d2) ;
947 if ( sd >= halfRminTolerance ) // It << 1015 if ( s >= kRadTolerance*0.5 ) // It was >= 0 ??
948 { 1016 {
949 xi = p.x() + sd*v.x() ; << 1017 xi = p.x() + s*v.x() ;
950 yi = p.y() + sd*v.y() ; << 1018 yi = p.y() + s*v.y() ;
951 rhoi = std::sqrt(xi*xi+yi*yi) ; << 1019 rhoi = sqrt(xi*xi+yi*yi) ;
952 1020
953 if ( !fFullPhiSphere && (rhoi != 0.0 << 1021 if ( segPhi && rhoi ) // Check phi intersection
954 { 1022 {
955 cosPsi = (xi*cosCPhi + yi*sinCPhi) 1023 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
956 1024
957 if (cosPsi >= cosHDPhiOT) 1025 if (cosPsi >= cosHDPhiOT)
958 { 1026 {
959 if ( !fFullThetaSphere ) // Che << 1027 if (segTheta) // Check theta intersection
960 { 1028 {
961 zi = p.z() + sd*v.z() ; << 1029 zi = p.z() + s*v.z() ;
962 1030
963 // rhoi & zi can never both be 1031 // rhoi & zi can never both be 0
964 // (=>intersect at origin =>fR 1032 // (=>intersect at origin =>fRmax=0)
965 // 1033 //
966 iTheta = std::atan2(rhoi,zi) ; << 1034 iTheta = atan2(rhoi,zi) ;
967 if ( (iTheta >= tolSTheta) && 1035 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
968 { 1036 {
969 snxt = sd; << 1037 snxt = s ;
970 } 1038 }
971 } 1039 }
972 else 1040 else
973 { 1041 {
974 snxt=sd; << 1042 snxt=s;
975 } 1043 }
976 } 1044 }
977 } 1045 }
978 else 1046 else
979 { 1047 {
980 if ( !fFullThetaSphere ) // Chec << 1048 if (segTheta) // Check theta intersection
981 { 1049 {
982 zi = p.z() + sd*v.z() ; << 1050 zi = p.z() + s*v.z() ;
983 1051
984 // rhoi & zi can never both be 0 1052 // rhoi & zi can never both be 0
985 // (=>intersect at origin => fRm 1053 // (=>intersect at origin => fRmax=0 !)
986 // 1054 //
987 iTheta = std::atan2(rhoi,zi) ; << 1055 iTheta = atan2(rhoi,zi) ;
988 if ( (iTheta >= tolSTheta) && (i 1056 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
989 { 1057 {
990 snxt = sd; << 1058 snxt = s ;
991 } 1059 }
992 } 1060 }
993 else 1061 else
994 { 1062 {
995 snxt = sd; << 1063 snxt=s;
996 } 1064 }
997 } 1065 }
998 } 1066 }
999 } 1067 }
1000 } 1068 }
1001 } 1069 }
1002 1070
1003 // Phi segment intersection 1071 // Phi segment intersection
1004 // 1072 //
1005 // o Tolerant of points inside phi planes b 1073 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1006 // 1074 //
1007 // o NOTE: Large duplication of code betwee 1075 // o NOTE: Large duplication of code between sphi & ephi checks
1008 // -> only diffs: sphi -> ephi, Com 1076 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1009 // intersection check <=0 -> >=0 1077 // intersection check <=0 -> >=0
1010 // -> Should use some form of loop 1078 // -> Should use some form of loop Construct
1011 // 1079 //
1012 if ( !fFullPhiSphere ) << 1080 if ( segPhi )
1013 { 1081 {
1014 // First phi surface ('S'tarting phi) << 1082 // First phi surface (`S'tarting phi)
>> 1083
>> 1084 sinSPhi = sin(fSPhi) ;
>> 1085 cosSPhi = cos(fSPhi) ;
>> 1086
1015 // Comp = Component in outwards normal di 1087 // Comp = Component in outwards normal dirn
1016 // 1088 //
1017 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; << 1089 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1018 << 1090
1019 if ( Comp < 0 ) 1091 if ( Comp < 0 )
1020 { 1092 {
1021 Dist = p.y()*cosSPhi - p.x()*sinSPhi ; 1093 Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
1022 1094
1023 if (Dist < halfCarTolerance) << 1095 if (Dist < kCarTolerance*0.5)
1024 { 1096 {
1025 sd = Dist/Comp ; << 1097 s = Dist/Comp ;
1026 1098
1027 if (sd < snxt) << 1099 if (s < snxt)
1028 { 1100 {
1029 if ( sd > 0 ) << 1101 if ( s > 0 )
1030 { 1102 {
1031 xi = p.x() + sd*v.x() ; << 1103 xi = p.x() + s*v.x() ;
1032 yi = p.y() + sd*v.y() ; << 1104 yi = p.y() + s*v.y() ;
1033 zi = p.z() + sd*v.z() ; << 1105 zi = p.z() + s*v.z() ;
1034 rhoi2 = xi*xi + yi*yi ; 1106 rhoi2 = xi*xi + yi*yi ;
1035 radi2 = rhoi2 + zi*zi ; 1107 radi2 = rhoi2 + zi*zi ;
1036 } 1108 }
1037 else 1109 else
1038 { 1110 {
1039 sd = 0 ; << 1111 s = 0 ;
1040 xi = p.x() ; 1112 xi = p.x() ;
1041 yi = p.y() ; 1113 yi = p.y() ;
1042 zi = p.z() ; 1114 zi = p.z() ;
1043 rhoi2 = rho2 ; 1115 rhoi2 = rho2 ;
1044 radi2 = rad2 ; 1116 radi2 = rad2 ;
1045 } 1117 }
1046 if ( (radi2 <= tolORMax2) 1118 if ( (radi2 <= tolORMax2)
1047 && (radi2 >= tolORMin2) 1119 && (radi2 >= tolORMin2)
1048 && ((yi*cosCPhi-xi*sinCPhi) <= 0) 1120 && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
1049 { 1121 {
1050 // Check theta intersection 1122 // Check theta intersection
1051 // rhoi & zi can never both be 0 1123 // rhoi & zi can never both be 0
1052 // (=>intersect at origin =>fRmax 1124 // (=>intersect at origin =>fRmax=0)
1053 // 1125 //
1054 if ( !fFullThetaSphere ) << 1126 if ( segTheta )
1055 { 1127 {
1056 iTheta = std::atan2(std::sqrt(r << 1128 iTheta = atan2(sqrt(rhoi2),zi) ;
1057 if ( (iTheta >= tolSTheta) && ( 1129 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1058 { 1130 {
1059 // r and theta intersections 1131 // r and theta intersections good
1060 // - check intersecting with 1132 // - check intersecting with correct half-plane
1061 1133
1062 if ((yi*cosCPhi-xi*sinCPhi) < 1134 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
1063 { 1135 {
1064 snxt = sd; << 1136 snxt = s ;
1065 } 1137 }
1066 } 1138 }
1067 } 1139 }
1068 else 1140 else
1069 { 1141 {
1070 snxt = sd; << 1142 snxt = s ;
1071 } 1143 }
1072 } 1144 }
1073 } 1145 }
1074 } 1146 }
1075 } 1147 }
1076 1148
1077 // Second phi surface ('E'nding phi) << 1149 // Second phi surface (`E'nding phi)
1078 // Component in outwards normal dirn <<
1079 1150
1080 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ; << 1151 ePhi = fSPhi + fDPhi ;
>> 1152 sinEPhi = sin(ePhi) ;
>> 1153 cosEPhi = cos(ePhi) ;
1081 1154
>> 1155 // Compnent in outwards normal dirn
>> 1156
>> 1157 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
>> 1158
1082 if (Comp < 0) 1159 if (Comp < 0)
1083 { 1160 {
1084 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ; 1161 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
1085 if ( Dist < halfCarTolerance ) << 1162 if ( Dist < kCarTolerance*0.5 )
1086 { 1163 {
1087 sd = Dist/Comp ; << 1164 s = Dist/Comp ;
1088 1165
1089 if ( sd < snxt ) << 1166 if ( s < snxt )
1090 { 1167 {
1091 if (sd > 0) << 1168 if (s > 0)
1092 { 1169 {
1093 xi = p.x() + sd*v.x() ; << 1170 xi = p.x() + s*v.x() ;
1094 yi = p.y() + sd*v.y() ; << 1171 yi = p.y() + s*v.y() ;
1095 zi = p.z() + sd*v.z() ; << 1172 zi = p.z() + s*v.z() ;
1096 rhoi2 = xi*xi + yi*yi ; 1173 rhoi2 = xi*xi + yi*yi ;
1097 radi2 = rhoi2 + zi*zi ; 1174 radi2 = rhoi2 + zi*zi ;
1098 } 1175 }
1099 else 1176 else
1100 { 1177 {
1101 sd = 0 ; << 1178 s = 0 ;
1102 xi = p.x() ; 1179 xi = p.x() ;
1103 yi = p.y() ; 1180 yi = p.y() ;
1104 zi = p.z() ; 1181 zi = p.z() ;
1105 rhoi2 = rho2 ; 1182 rhoi2 = rho2 ;
1106 radi2 = rad2 ; 1183 radi2 = rad2 ;
1107 } 1184 }
1108 if ( (radi2 <= tolORMax2) 1185 if ( (radi2 <= tolORMax2)
1109 && (radi2 >= tolORMin2) 1186 && (radi2 >= tolORMin2)
1110 && ((yi*cosCPhi-xi*sinCPhi) >= 0) 1187 && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
1111 { 1188 {
1112 // Check theta intersection 1189 // Check theta intersection
1113 // rhoi & zi can never both be 0 1190 // rhoi & zi can never both be 0
1114 // (=>intersect at origin =>fRmax 1191 // (=>intersect at origin =>fRmax=0)
1115 // 1192 //
1116 if ( !fFullThetaSphere ) << 1193 if ( segTheta )
1117 { 1194 {
1118 iTheta = std::atan2(std::sqrt(r << 1195 iTheta = atan2(sqrt(rhoi2),zi) ;
1119 if ( (iTheta >= tolSTheta) && ( 1196 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1120 { 1197 {
1121 // r and theta intersections 1198 // r and theta intersections good
1122 // - check intersecting with 1199 // - check intersecting with correct half-plane
1123 1200
1124 if ((yi*cosCPhi-xi*sinCPhi) > 1201 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
1125 { 1202 {
1126 snxt = sd; << 1203 snxt = s ;
1127 } 1204 }
1128 } 1205 }
1129 } 1206 }
1130 else 1207 else
1131 { 1208 {
1132 snxt = sd; << 1209 snxt = s ;
1133 } 1210 }
1134 } 1211 }
1135 } 1212 }
1136 } 1213 }
1137 } 1214 }
1138 } 1215 }
1139 1216
1140 // Theta segment intersection 1217 // Theta segment intersection
1141 1218
1142 if ( !fFullThetaSphere ) << 1219 if ( segTheta )
1143 { 1220 {
1144 1221
1145 // Intersection with theta surfaces 1222 // Intersection with theta surfaces
1146 // Known failure cases: 1223 // Known failure cases:
1147 // o Inside tolerance of stheta surface, 1224 // o Inside tolerance of stheta surface, skim
1148 // ~parallel to cone and Hit & enter e 1225 // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1149 // 1226 //
1150 // To solve: Check 2nd root of etheta 1227 // To solve: Check 2nd root of etheta surface in addition to stheta
1151 // 1228 //
1152 // o start/end theta is exactly pi/2 << 1229 // o start/end theta is exactly pi/2
1153 // Intersections with cones 1230 // Intersections with cones
1154 // 1231 //
1155 // Cone equation: x^2+y^2=z^2tan^2(t) 1232 // Cone equation: x^2+y^2=z^2tan^2(t)
1156 // 1233 //
1157 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan 1234 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1158 // 1235 //
1159 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+p << 1236 // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t))
1160 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = << 1237 // + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1161 // 1238 //
1162 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d << 1239 // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1163 // + (rho2-pz^2tan^2(t)) = 0 <<
1164 1240
1165 if (fSTheta != 0.0) << 1241 tanSTheta = tan(fSTheta) ;
>> 1242 tanSTheta2 = tanSTheta*tanSTheta ;
>> 1243 tanETheta = tan(fSTheta+fDTheta) ;
>> 1244 tanETheta2 = tanETheta*tanETheta ;
>> 1245
>> 1246 if (fSTheta)
1166 { 1247 {
1167 dist2STheta = rho2 - p.z()*p.z()*tanSTh 1248 dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
1168 } 1249 }
1169 else 1250 else
1170 { 1251 {
1171 dist2STheta = kInfinity ; 1252 dist2STheta = kInfinity ;
1172 } 1253 }
1173 if ( eTheta < pi ) << 1254 if ( fSTheta + fDTheta < M_PI )
1174 { 1255 {
1175 dist2ETheta=rho2-p.z()*p.z()*tanETheta2 1256 dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
1176 } 1257 }
1177 else << 1258 else
1178 { 1259 {
1179 dist2ETheta=kInfinity; 1260 dist2ETheta=kInfinity;
1180 } << 1261 }
1181 if ( pTheta < tolSTheta ) << 1262 if ( pTheta < tolSTheta) // dist2STheta<-kRadTolerance*0.5 && dist2ETheta>0)
1182 { 1263 {
1183 // Inside (theta<stheta-tol) stheta con << 1264 // Inside (theta<stheta-tol) s theta cone
1184 // First root of stheta cone, second if 1265 // First root of stheta cone, second if first root -ve
1185 1266
1186 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1267 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1187 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1268 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1188 if (t1 != 0.0) << 1269
>> 1270 b = t2/t1 ;
>> 1271 c = dist2STheta/t1 ;
>> 1272 d2 = b*b - c ;
>> 1273
>> 1274 if ( d2 >= 0 )
1189 { 1275 {
>> 1276 d = sqrt(d2) ;
>> 1277 s = -b - d ; // First root
>> 1278
>> 1279 if ( s < 0 )
>> 1280 {
>> 1281 s=-b+d; // Second root
>> 1282 }
>> 1283 if (s >= 0 && s < snxt)
>> 1284 {
>> 1285 xi = p.x() + s*v.x() ;
>> 1286 yi = p.y() + s*v.y() ;
>> 1287 zi = p.z() + s*v.z() ;
>> 1288 rhoi2 = xi*xi + yi*yi ;
>> 1289 radi2 = rhoi2 + zi*zi ;
>> 1290 if ( (radi2 <= tolORMax2)
>> 1291 && (radi2 >= tolORMin2)
>> 1292 && (zi*(fSTheta - 0.5*pi) <= 0) )
>> 1293 {
>> 1294 if ( segPhi && rhoi2 ) // Check phi intersection
>> 1295 {
>> 1296 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
>> 1297 if (cosPsi >= cosHDPhiOT)
>> 1298 {
>> 1299 snxt = s ;
>> 1300 }
>> 1301 }
>> 1302 else
>> 1303 {
>> 1304 snxt = s ;
>> 1305 }
>> 1306 }
>> 1307 }
>> 1308 }
>> 1309
>> 1310 // Possible intersection with ETheta cone.
>> 1311 // Second >= 0 root should be considered
>> 1312
>> 1313 if ( fSTheta + fDTheta < M_PI )
>> 1314 {
>> 1315 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
>> 1316 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
>> 1317
1190 b = t2/t1 ; 1318 b = t2/t1 ;
1191 c = dist2STheta/t1 ; << 1319 c = dist2ETheta/t1 ;
1192 d2 = b*b - c ; 1320 d2 = b*b - c ;
1193 1321
1194 if ( d2 >= 0 ) << 1322 if (d2 >= 0)
1195 { 1323 {
1196 d = std::sqrt(d2) ; << 1324 d = sqrt(d2) ;
1197 sd = -b - d ; // First root << 1325 s = -b + d ; // Second root
1198 zi = p.z() + sd*v.z(); <<
1199 1326
1200 if ( (sd < 0) || (zi*(fSTheta - hal << 1327 if (s >= 0 && s < snxt)
1201 { 1328 {
1202 sd = -b+d; // Second root << 1329 xi = p.x() + s*v.x() ;
1203 } << 1330 yi = p.y() + s*v.y() ;
1204 if ((sd >= 0) && (sd < snxt)) << 1331 zi = p.z() + s*v.z() ;
1205 { << 1332 rhoi2 = xi*xi + yi*yi ;
1206 xi = p.x() + sd*v.x(); << 1333 radi2 = rhoi2 + zi*zi ;
1207 yi = p.y() + sd*v.y(); << 1334
1208 zi = p.z() + sd*v.z(); <<
1209 rhoi2 = xi*xi + yi*yi; <<
1210 radi2 = rhoi2 + zi*zi; <<
1211 if ( (radi2 <= tolORMax2) 1335 if ( (radi2 <= tolORMax2)
1212 && (radi2 >= tolORMin2) 1336 && (radi2 >= tolORMin2)
1213 && (zi*(fSTheta - halfpi) <= 0) << 1337 && (zi*(fSTheta + fDTheta - 0.5*pi) <= 0) )
1214 { 1338 {
1215 if ( !fFullPhiSphere && (rhoi2 << 1339 if (segPhi && rhoi2) // Check phi intersection
1216 { 1340 {
1217 cosPsi = (xi*cosCPhi + yi*sin << 1341 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
1218 if (cosPsi >= cosHDPhiOT) 1342 if (cosPsi >= cosHDPhiOT)
1219 { 1343 {
1220 snxt = sd; << 1344 snxt = s ;
1221 } 1345 }
1222 } 1346 }
1223 else 1347 else
1224 { 1348 {
1225 snxt = sd; << 1349 snxt = s ;
1226 } 1350 }
1227 } 1351 }
1228 } 1352 }
1229 } 1353 }
1230 } 1354 }
>> 1355 }
>> 1356 else if (pTheta > tolETheta)
>> 1357 { // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
>> 1358 // Inside (theta>etheta+tol) e theta cone
>> 1359 // First root of etheta cone, second if first root `imaginary'
1231 1360
1232 // Possible intersection with ETheta co << 1361 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1233 // Second >= 0 root should be considere << 1362 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
>> 1363
>> 1364 b = t2/t1 ;
>> 1365 c = dist2ETheta/t1 ;
>> 1366 d2 = b*b - c ;
1234 1367
1235 if ( eTheta < pi ) << 1368 if (d2 >= 0)
1236 { 1369 {
1237 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) << 1370 d = sqrt(d2) ;
1238 t2 = pDotV2d - p.z()*v.z()*tanETheta2 << 1371 s = -b - d ; // First root
1239 if (t1 != 0.0) << 1372 if (s < 0)
1240 { << 1373 {
1241 b = t2/t1 ; << 1374 s = -b + d ; // second root
1242 c = dist2ETheta/t1 ; << 1375 }
1243 d2 = b*b - c ; << 1376 if (s >= 0 && s < snxt)
>> 1377 {
>> 1378 xi = p.x() + s*v.x() ;
>> 1379 yi = p.y() + s*v.y() ;
>> 1380 zi = p.z() + s*v.z() ;
>> 1381 rhoi2 = xi*xi + yi*yi ;
>> 1382 radi2 = rhoi2 + zi*zi ;
1244 1383
1245 if (d2 >= 0) << 1384 if ( (radi2 <= tolORMax2)
>> 1385 && (radi2 >= tolORMin2)
>> 1386 && (zi*(fSTheta + fDTheta - 0.5*pi) <= 0) )
1246 { 1387 {
1247 d = std::sqrt(d2) ; << 1388 if (segPhi && rhoi2) // Check phi intersection
1248 sd = -b + d ; // Second root <<
1249 <<
1250 if ( (sd >= 0) && (sd < snxt) ) <<
1251 { 1389 {
1252 xi = p.x() + sd*v.x() ; << 1390 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
1253 yi = p.y() + sd*v.y() ; << 1391 if (cosPsi >= cosHDPhiOT)
1254 zi = p.z() + sd*v.z() ; <<
1255 rhoi2 = xi*xi + yi*yi ; <<
1256 radi2 = rhoi2 + zi*zi ; <<
1257 <<
1258 if ( (radi2 <= tolORMax2) <<
1259 && (radi2 >= tolORMin2) <<
1260 && (zi*(eTheta - halfpi) <= 0 <<
1261 { 1392 {
1262 if (!fFullPhiSphere && (rhoi2 << 1393 snxt = s ;
1263 { <<
1264 cosPsi = (xi*cosCPhi + yi*s <<
1265 if (cosPsi >= cosHDPhiOT) <<
1266 { <<
1267 snxt = sd; <<
1268 } <<
1269 } <<
1270 else <<
1271 { <<
1272 snxt = sd; <<
1273 } <<
1274 } 1394 }
1275 } 1395 }
>> 1396 else
>> 1397 {
>> 1398 snxt = s ;
>> 1399 }
1276 } 1400 }
1277 } 1401 }
1278 } 1402 }
1279 } <<
1280 else if ( pTheta > tolETheta ) <<
1281 { <<
1282 // dist2ETheta<-kRadTolerance*0.5 && di <<
1283 // Inside (theta > etheta+tol) e-theta <<
1284 // First root of etheta cone, second if <<
1285 1403
1286 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; << 1404 // Possible intersection with STheta cone.
1287 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; << 1405 // Second >= 0 root should be considered
1288 if (t1 != 0.0) << 1406
>> 1407 if ( fSTheta )
1289 { 1408 {
>> 1409 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
>> 1410 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
>> 1411
1290 b = t2/t1 ; 1412 b = t2/t1 ;
1291 c = dist2ETheta/t1 ; << 1413 c = dist2STheta/t1 ;
1292 d2 = b*b - c ; 1414 d2 = b*b - c ;
1293 1415
1294 if (d2 >= 0) 1416 if (d2 >= 0)
1295 { 1417 {
1296 d = std::sqrt(d2) ; << 1418 d = sqrt(d2) ;
1297 sd = -b - d ; // First root << 1419 s = -b + d ; // Second root
1298 zi = p.z() + sd*v.z(); <<
1299 1420
1300 if ( (sd < 0) || (zi*(eTheta - half << 1421 if ( (s >= 0) && (s < snxt) )
1301 { 1422 {
1302 sd = -b + d ; // second << 1423 xi = p.x() + s*v.x() ;
1303 } << 1424 yi = p.y() + s*v.y() ;
1304 if ( (sd >= 0) && (sd < snxt) ) << 1425 zi = p.z() + s*v.z() ;
1305 { <<
1306 xi = p.x() + sd*v.x() ; <<
1307 yi = p.y() + sd*v.y() ; <<
1308 zi = p.z() + sd*v.z() ; <<
1309 rhoi2 = xi*xi + yi*yi ; 1426 rhoi2 = xi*xi + yi*yi ;
1310 radi2 = rhoi2 + zi*zi ; 1427 radi2 = rhoi2 + zi*zi ;
1311 1428
1312 if ( (radi2 <= tolORMax2) 1429 if ( (radi2 <= tolORMax2)
1313 && (radi2 >= tolORMin2) 1430 && (radi2 >= tolORMin2)
1314 && (zi*(eTheta - halfpi) <= 0) << 1431 && (zi*(fSTheta - 0.5*pi) <= 0) )
1315 { 1432 {
1316 if (!fFullPhiSphere && (rhoi2 ! << 1433 if (segPhi && rhoi2) // Check phi intersection
1317 { 1434 {
1318 cosPsi = (xi*cosCPhi + yi*sin << 1435 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
1319 if (cosPsi >= cosHDPhiOT) 1436 if (cosPsi >= cosHDPhiOT)
1320 { 1437 {
1321 snxt = sd; << 1438 snxt = s ;
1322 } 1439 }
1323 } 1440 }
1324 else 1441 else
1325 { 1442 {
1326 snxt = sd; << 1443 snxt = s ;
1327 } 1444 }
1328 } 1445 }
1329 } 1446 }
1330 } 1447 }
1331 } << 1448 }
1332 << 1449 }
1333 // Possible intersection with STheta co << 1450 else if ( (pTheta <tolSTheta + kAngTolerance)
1334 // Second >= 0 root should be considere << 1451 && (fSTheta > kAngTolerance) )
1335 <<
1336 if ( fSTheta != 0.0 ) <<
1337 { <<
1338 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) <<
1339 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 <<
1340 if (t1 != 0.0) <<
1341 { <<
1342 b = t2/t1 ; <<
1343 c = dist2STheta/t1 ; <<
1344 d2 = b*b - c ; <<
1345 <<
1346 if (d2 >= 0) <<
1347 { <<
1348 d = std::sqrt(d2) ; <<
1349 sd = -b + d ; // Second root <<
1350 <<
1351 if ( (sd >= 0) && (sd < snxt) ) <<
1352 { <<
1353 xi = p.x() + sd*v.x() ; <<
1354 yi = p.y() + sd*v.y() ; <<
1355 zi = p.z() + sd*v.z() ; <<
1356 rhoi2 = xi*xi + yi*yi ; <<
1357 radi2 = rhoi2 + zi*zi ; <<
1358 <<
1359 if ( (radi2 <= tolORMax2) <<
1360 && (radi2 >= tolORMin2) <<
1361 && (zi*(fSTheta - halfpi) <= <<
1362 { <<
1363 if (!fFullPhiSphere && (rhoi2 <<
1364 { <<
1365 cosPsi = (xi*cosCPhi + yi*s <<
1366 if (cosPsi >= cosHDPhiOT) <<
1367 { <<
1368 snxt = sd; <<
1369 } <<
1370 } <<
1371 else <<
1372 { <<
1373 snxt = sd; <<
1374 } <<
1375 } <<
1376 } <<
1377 } <<
1378 } <<
1379 } <<
1380 } <<
1381 else if ( (pTheta < tolSTheta + kAngToler <<
1382 && (fSTheta > halfAngTolerance) ) <<
1383 { 1452 {
1384 // In tolerance of stheta 1453 // In tolerance of stheta
1385 // If entering through solid [r,phi] => 1454 // If entering through solid [r,phi] => 0 to in
1386 // else try 2nd root 1455 // else try 2nd root
1387 1456
1388 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1457 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1389 if ( (t2>=0 && tolIRMin2<rad2 && rad2<t << 1458 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<M_PI*.5)
1390 || (t2<0 && tolIRMin2<rad2 && rad2<t << 1459 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>M_PI*.5)
1391 || (v.z()<0 && tolIRMin2<rad2 && rad2 << 1460 || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==M_PI*.5) )
1392 { 1461 {
1393 if (!fFullPhiSphere && (rho2 != 0.0)) << 1462 if (segPhi && rho2) // Check phi intersection
1394 { 1463 {
1395 cosPsi = (p.x()*cosCPhi + p.y()*sin << 1464 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/sqrt(rho2) ;
1396 if (cosPsi >= cosHDPhiIT) 1465 if (cosPsi >= cosHDPhiIT)
1397 { 1466 {
1398 return 0 ; 1467 return 0 ;
1399 } 1468 }
1400 } 1469 }
1401 else 1470 else
1402 { 1471 {
1403 return 0 ; 1472 return 0 ;
1404 } 1473 }
1405 } 1474 }
1406 1475
1407 // Not entering immediately/travelling 1476 // Not entering immediately/travelling through
1408 1477
1409 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1478 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1410 if (t1 != 0.0) << 1479 b = t2/t1 ;
1411 { << 1480 c = dist2STheta/t1 ;
1412 b = t2/t1 ; << 1481 d2 = b*b - c ;
1413 c = dist2STheta/t1 ; <<
1414 d2 = b*b - c ; <<
1415 1482
1416 if (d2 >= 0) << 1483 if (d2 >= 0)
1417 { << 1484 {
1418 d = std::sqrt(d2) ; << 1485 d = sqrt(d2) ;
1419 sd = -b + d ; << 1486 s = -b + d ;
1420 if ( (sd >= halfCarTolerance) && (s << 1487 if ( (s >= kCarTolerance*0.5) && (s < snxt) && (fSTheta < M_PI*0.5) )
1421 { // ^^^^^^^^^^^^^^^^^^^^^ shoul << 1488 {
1422 xi = p.x() + sd*v.x() ; << 1489 xi = p.x() + s*v.x() ;
1423 yi = p.y() + sd*v.y() ; << 1490 yi = p.y() + s*v.y() ;
1424 zi = p.z() + sd*v.z() ; << 1491 zi = p.z() + s*v.z() ;
1425 rhoi2 = xi*xi + yi*yi ; << 1492 rhoi2 = xi*xi + yi*yi ;
1426 radi2 = rhoi2 + zi*zi ; << 1493 radi2 = rhoi2 + zi*zi ;
1427 1494
1428 if ( (radi2 <= tolORMax2) << 1495 if ( (radi2 <= tolORMax2)
1429 && (radi2 >= tolORMin2) << 1496 && (radi2 >= tolORMin2)
1430 && (zi*(fSTheta - halfpi) <= 0) << 1497 && (zi*(fSTheta - 0.5*pi) <= 0) )
>> 1498 {
>> 1499 if ( segPhi && rhoi2 ) // Check phi intersection
1431 { 1500 {
1432 if ( !fFullPhiSphere && (rhoi2 << 1501 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
1433 { << 1502 if ( cosPsi >= cosHDPhiOT )
1434 cosPsi = (xi*cosCPhi + yi*sin <<
1435 if ( cosPsi >= cosHDPhiOT ) <<
1436 { <<
1437 snxt = sd; <<
1438 } <<
1439 } <<
1440 else <<
1441 { 1503 {
1442 snxt = sd; << 1504 snxt = s ;
1443 } 1505 }
1444 } 1506 }
>> 1507 else
>> 1508 {
>> 1509 snxt = s ;
>> 1510 }
1445 } 1511 }
1446 } 1512 }
1447 } 1513 }
1448 } << 1514 }
1449 else if ((pTheta > tolETheta-kAngToleranc << 1515 else if ( (pTheta > tolETheta - kAngTolerance)
>> 1516 && ((fSTheta + fDTheta) < M_PI-kAngTolerance) )
1450 { 1517 {
1451 1518
1452 // In tolerance of etheta 1519 // In tolerance of etheta
1453 // If entering through solid [r,phi] => 1520 // If entering through solid [r,phi] => 0 to in
1454 // else try 2nd root 1521 // else try 2nd root
1455 1522
1456 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; 1523 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1457 1524
1458 if ( ((t2<0) && (eTheta < halfpi) << 1525 if (
1459 && (tolIRMin2 < rad2) && (rad2 < to << 1526 (t2<0 && (fSTheta+fDTheta) <M_PI*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
1460 || ((t2>=0) && (eTheta > halfpi) << 1527 || (t2>=0 && (fSTheta+fDTheta) >M_PI*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
1461 && (tolIRMin2 < rad2) && (rad2 < to << 1528 || (v.z()>0 && (fSTheta+fDTheta)==M_PI*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
1462 || ((v.z()>0) && (eTheta == halfpi) << 1529 )
1463 && (tolIRMin2 < rad2) && (rad2 < to <<
1464 { 1530 {
1465 if (!fFullPhiSphere && (rho2 != 0.0)) << 1531 if (segPhi && rho2) // Check phi intersection
1466 { 1532 {
1467 cosPsi = (p.x()*cosCPhi + p.y()*sin << 1533 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/sqrt(rho2) ;
1468 if (cosPsi >= cosHDPhiIT) 1534 if (cosPsi >= cosHDPhiIT)
1469 { 1535 {
1470 return 0 ; 1536 return 0 ;
1471 } 1537 }
1472 } 1538 }
1473 else 1539 else
1474 { 1540 {
1475 return 0 ; 1541 return 0 ;
1476 } 1542 }
1477 } 1543 }
1478 1544
1479 // Not entering immediately/travelling 1545 // Not entering immediately/travelling through
1480 1546
1481 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1547 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1482 if (t1 != 0.0) << 1548 b = t2/t1 ;
1483 { << 1549 c = dist2ETheta/t1 ;
1484 b = t2/t1 ; << 1550 d2 = b*b - c ;
1485 c = dist2ETheta/t1 ; <<
1486 d2 = b*b - c ; <<
1487 1551
1488 if (d2 >= 0) << 1552 if (d2 >= 0)
1489 { << 1553 {
1490 d = std::sqrt(d2) ; << 1554 d = sqrt(d2) ;
1491 sd = -b + d ; << 1555 s = -b + d ;
>> 1556
>> 1557 if ( (s >= kCarTolerance*0.5)
>> 1558 && (s < snxt) && ((fSTheta + fDTheta) > M_PI*0.5) )
>> 1559 {
>> 1560 xi = p.x() + s*v.x() ;
>> 1561 yi = p.y() + s*v.y() ;
>> 1562 zi = p.z() + s*v.z() ;
>> 1563 rhoi2 = xi*xi + yi*yi ;
>> 1564 radi2 = rhoi2 + zi*zi ;
1492 1565
1493 if ( (sd >= halfCarTolerance) << 1566 if ( (radi2 <= tolORMax2)
1494 && (sd < snxt) && (eTheta > halfp << 1567 && (radi2 >= tolORMin2)
>> 1568 && (zi*(fSTheta + fDTheta - 0.5*pi) <= 0) )
1495 { 1569 {
1496 xi = p.x() + sd*v.x() ; << 1570 if (segPhi && rhoi2) // Check phi intersection
1497 yi = p.y() + sd*v.y() ; <<
1498 zi = p.z() + sd*v.z() ; <<
1499 rhoi2 = xi*xi + yi*yi ; <<
1500 radi2 = rhoi2 + zi*zi ; <<
1501 <<
1502 if ( (radi2 <= tolORMax2) <<
1503 && (radi2 >= tolORMin2) <<
1504 && (zi*(eTheta - halfpi) <= 0) <<
1505 { 1571 {
1506 if (!fFullPhiSphere && (rhoi2 ! << 1572 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
1507 { << 1573 if (cosPsi>=cosHDPhiOT)
1508 cosPsi = (xi*cosCPhi + yi*sin <<
1509 if (cosPsi >= cosHDPhiOT) <<
1510 { <<
1511 snxt = sd; <<
1512 } <<
1513 } <<
1514 else <<
1515 { 1574 {
1516 snxt = sd; << 1575 snxt = s ;
1517 } 1576 }
1518 } 1577 }
>> 1578 else
>> 1579 {
>> 1580 snxt = s ;
>> 1581 }
1519 } 1582 }
1520 } 1583 }
1521 } << 1584 }
1522 } << 1585 }
1523 else 1586 else
1524 { 1587 {
1525 // stheta+tol<theta<etheta-tol 1588 // stheta+tol<theta<etheta-tol
1526 // For BOTH stheta & etheta check 2nd r 1589 // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1527 1590
1528 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1591 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1529 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1592 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1530 if (t1 != 0.0) << 1593
>> 1594 b = t2/t1;
>> 1595 c = dist2STheta/t1 ;
>> 1596 d2 = b*b - c ;
>> 1597
>> 1598 if (d2 >= 0)
1531 { 1599 {
1532 b = t2/t1; << 1600 d = sqrt(d2) ;
1533 c = dist2STheta/t1 ; << 1601 s = -b + d ; // second root
1534 d2 = b*b - c ; <<
1535 1602
1536 if (d2 >= 0) << 1603 if (s >= 0 && s < snxt)
1537 { 1604 {
1538 d = std::sqrt(d2) ; << 1605 xi = p.x() + s*v.x() ;
1539 sd = -b + d ; // second root << 1606 yi = p.y() + s*v.y() ;
>> 1607 zi = p.z() + s*v.z() ;
>> 1608 rhoi2 = xi*xi + yi*yi ;
>> 1609 radi2 = rhoi2 + zi*zi ;
1540 1610
1541 if ((sd >= 0) && (sd < snxt)) << 1611 if ( (radi2 <= tolORMax2)
>> 1612 && (radi2 >= tolORMin2)
>> 1613 && (zi*(fSTheta - 0.5*pi) <= 0) )
1542 { 1614 {
1543 xi = p.x() + sd*v.x() ; << 1615 if (segPhi && rhoi2) // Check phi intersection
1544 yi = p.y() + sd*v.y() ; <<
1545 zi = p.z() + sd*v.z() ; <<
1546 rhoi2 = xi*xi + yi*yi ; <<
1547 radi2 = rhoi2 + zi*zi ; <<
1548 <<
1549 if ( (radi2 <= tolORMax2) <<
1550 && (radi2 >= tolORMin2) <<
1551 && (zi*(fSTheta - halfpi) <= 0) <<
1552 { 1616 {
1553 if (!fFullPhiSphere && (rhoi2 ! << 1617 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
1554 { << 1618 if (cosPsi >= cosHDPhiOT)
1555 cosPsi = (xi*cosCPhi + yi*sin <<
1556 if (cosPsi >= cosHDPhiOT) <<
1557 { <<
1558 snxt = sd; <<
1559 } <<
1560 } <<
1561 else <<
1562 { 1619 {
1563 snxt = sd; << 1620 snxt = s ;
1564 } 1621 }
1565 } 1622 }
>> 1623 else
>> 1624 {
>> 1625 snxt = s ;
>> 1626 }
1566 } 1627 }
1567 } 1628 }
1568 } << 1629 }
1569 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1630 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1570 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; 1631 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1571 if (t1 != 0.0) << 1632
>> 1633 b = t2/t1 ;
>> 1634 c = dist2ETheta/t1 ;
>> 1635 d2 = b*b - c ;
>> 1636
>> 1637 if (d2 >= 0)
1572 { 1638 {
1573 b = t2/t1 ; << 1639 d = sqrt(d2) ;
1574 c = dist2ETheta/t1 ; << 1640 s = -b + d; // second root
1575 d2 = b*b - c ; <<
1576 1641
1577 if (d2 >= 0) << 1642 if (s >= 0 && s < snxt)
1578 { 1643 {
1579 d = std::sqrt(d2) ; << 1644 xi = p.x() + s*v.x() ;
1580 sd = -b + d; // second root << 1645 yi = p.y() + s*v.y() ;
>> 1646 zi = p.z() + s*v.z() ;
>> 1647 rhoi2 = xi*xi + yi*yi ;
>> 1648 radi2 = rhoi2 + zi*zi ;
1581 1649
1582 if ((sd >= 0) && (sd < snxt)) << 1650 if ( (radi2 <= tolORMax2)
>> 1651 && (radi2 >= tolORMin2)
>> 1652 && (zi*(fSTheta + fDTheta - 0.5*pi) <= 0) )
1583 { 1653 {
1584 xi = p.x() + sd*v.x() ; << 1654 if (segPhi && rhoi2) // Check phi intersection
1585 yi = p.y() + sd*v.y() ; <<
1586 zi = p.z() + sd*v.z() ; <<
1587 rhoi2 = xi*xi + yi*yi ; <<
1588 radi2 = rhoi2 + zi*zi ; <<
1589 <<
1590 if ( (radi2 <= tolORMax2) <<
1591 && (radi2 >= tolORMin2) <<
1592 && (zi*(eTheta - halfpi) <= 0) <<
1593 { 1655 {
1594 if (!fFullPhiSphere && (rhoi2 ! << 1656 cosPsi = (xi*cosCPhi + yi*sinCPhi)/sqrt(rhoi2) ;
>> 1657 if ( cosPsi >= cosHDPhiOT )
1595 { 1658 {
1596 cosPsi = (xi*cosCPhi + yi*sin << 1659 snxt=s;
1597 if ( cosPsi >= cosHDPhiOT ) <<
1598 { <<
1599 snxt = sd; <<
1600 } <<
1601 } <<
1602 else <<
1603 { <<
1604 snxt = sd; <<
1605 } 1660 }
1606 } 1661 }
>> 1662 else
>> 1663 {
>> 1664 snxt = s ;
>> 1665 }
1607 } 1666 }
1608 } 1667 }
1609 } 1668 }
1610 } << 1669 }
1611 } 1670 }
1612 return snxt; 1671 return snxt;
1613 } 1672 }
1614 1673
1615 ///////////////////////////////////////////// 1674 //////////////////////////////////////////////////////////////////////
1616 // 1675 //
1617 // Calculate distance (<= actual) to closest 1676 // Calculate distance (<= actual) to closest surface of shape from outside
1618 // - Calculate distance to radial planes 1677 // - Calculate distance to radial planes
1619 // - Only to phi planes if outside phi extent 1678 // - Only to phi planes if outside phi extent
1620 // - Only to theta planes if outside theta ex 1679 // - Only to theta planes if outside theta extent
1621 // - Return 0 if point inside 1680 // - Return 0 if point inside
1622 1681
1623 G4double G4Sphere::DistanceToIn( const G4Thre 1682 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p ) const
1624 { 1683 {
1625 G4double safe=0.0,safeRMin,safeRMax,safePhi << 1684 G4double safe,safeRMin,safeRMax,safePhi,safeTheta;
1626 G4double rho2,rds,rho; << 1685 G4double rho2,rad,rho;
1627 G4double cosPsi; << 1686 G4double phiC,cosPhiC,sinPhiC,cosPsi,ePhi;
1628 G4double pTheta,dTheta1,dTheta2; 1687 G4double pTheta,dTheta1,dTheta2;
1629 rho2=p.x()*p.x()+p.y()*p.y(); 1688 rho2=p.x()*p.x()+p.y()*p.y();
1630 rds=std::sqrt(rho2+p.z()*p.z()); << 1689 rad=sqrt(rho2+p.z()*p.z());
1631 rho=std::sqrt(rho2); << 1690 rho=sqrt(rho2);
1632 1691
1633 // 1692 //
1634 // Distance to r shells 1693 // Distance to r shells
1635 // << 1694 //
1636 if (fRmin != 0.0) << 1695 if (fRmin)
1637 { 1696 {
1638 safeRMin=fRmin-rds; << 1697 safeRMin=fRmin-rad;
1639 safeRMax=rds-fRmax; << 1698 safeRMax=rad-fRmax;
1640 if (safeRMin>safeRMax) 1699 if (safeRMin>safeRMax)
1641 { 1700 {
1642 safe=safeRMin; 1701 safe=safeRMin;
1643 } 1702 }
1644 else 1703 else
1645 { 1704 {
1646 safe=safeRMax; 1705 safe=safeRMax;
1647 } 1706 }
1648 } 1707 }
1649 else 1708 else
1650 { 1709 {
1651 safe=rds-fRmax; << 1710 safe=rad-fRmax;
1652 } 1711 }
1653 1712
1654 // 1713 //
1655 // Distance to phi extent 1714 // Distance to phi extent
1656 // 1715 //
1657 if (!fFullPhiSphere && (rho != 0.0)) << 1716 if (fDPhi<2.0*M_PI&&rho)
1658 { 1717 {
>> 1718 phiC=fSPhi+fDPhi*0.5;
>> 1719 cosPhiC=cos(phiC);
>> 1720 sinPhiC=sin(phiC);
>> 1721
1659 // Psi=angle from central phi to point 1722 // Psi=angle from central phi to point
1660 // 1723 //
1661 cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho; << 1724 cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho;
1662 if (cosPsi<cosHDPhi) << 1725 if (cosPsi<cos(fDPhi*0.5))
1663 { 1726 {
1664 // Point lies outside phi range 1727 // Point lies outside phi range
1665 // 1728 //
1666 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0) << 1729 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1667 { 1730 {
1668 safePhi=std::fabs(p.x()*sinSPhi-p.y() << 1731 safePhi=fabs(p.x()*sin(fSPhi)-p.y()*cos(fSPhi));
1669 } 1732 }
1670 else 1733 else
1671 { 1734 {
1672 safePhi=std::fabs(p.x()*sinEPhi-p.y() << 1735 ePhi=fSPhi+fDPhi;
>> 1736 safePhi=fabs(p.x()*sin(ePhi)-p.y()*cos(ePhi));
1673 } 1737 }
1674 if (safePhi>safe) { safe=safePhi; } << 1738 if (safePhi>safe) safe=safePhi;
1675 } 1739 }
1676 } 1740 }
1677 // 1741 //
1678 // Distance to Theta extent 1742 // Distance to Theta extent
1679 // << 1743 //
1680 if ((rds!=0.0) && (!fFullThetaSphere)) << 1744 if ((rad!=0.0) && (fDTheta<M_PI))
1681 { 1745 {
1682 pTheta=std::acos(p.z()/rds); << 1746 pTheta=acos(p.z()/rad);
1683 if (pTheta<0) { pTheta+=pi; } << 1747 if (pTheta<0) pTheta+=M_PI;
1684 dTheta1=fSTheta-pTheta; 1748 dTheta1=fSTheta-pTheta;
1685 dTheta2=pTheta-eTheta; << 1749 dTheta2=pTheta-(fSTheta+fDTheta);
1686 if (dTheta1>dTheta2) 1750 if (dTheta1>dTheta2)
1687 { 1751 {
1688 if (dTheta1>=0) // WHY ???? 1752 if (dTheta1>=0) // WHY ???????????
1689 { 1753 {
1690 safeTheta=rds*std::sin(dTheta1); << 1754 safeTheta=rad*sin(dTheta1);
1691 if (safe<=safeTheta) 1755 if (safe<=safeTheta)
1692 { 1756 {
1693 safe=safeTheta; 1757 safe=safeTheta;
1694 } 1758 }
1695 } 1759 }
1696 } 1760 }
1697 else 1761 else
1698 { 1762 {
1699 if (dTheta2>=0) 1763 if (dTheta2>=0)
1700 { 1764 {
1701 safeTheta=rds*std::sin(dTheta2); << 1765 safeTheta=rad*sin(dTheta2);
1702 if (safe<=safeTheta) 1766 if (safe<=safeTheta)
1703 { 1767 {
1704 safe=safeTheta; 1768 safe=safeTheta;
1705 } 1769 }
1706 } 1770 }
1707 } 1771 }
1708 } 1772 }
1709 1773
1710 if (safe<0) { safe=0; } << 1774 if (safe<0) safe=0;
1711 return safe; 1775 return safe;
1712 } 1776 }
1713 1777
1714 ///////////////////////////////////////////// 1778 /////////////////////////////////////////////////////////////////////
1715 // 1779 //
1716 // Calculate distance to surface of shape fro << 1780 // Calculate distance to surface of shape from `inside', allowing for tolerance
1717 // - Only Calc rmax intersection if no valid 1781 // - Only Calc rmax intersection if no valid rmin intersection
1718 1782
1719 G4double G4Sphere::DistanceToOut( const G4Thr 1783 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p,
1720 const G4Thr 1784 const G4ThreeVector& v,
1721 const G4boo 1785 const G4bool calcNorm,
1722 G4boo << 1786 G4bool *validNorm,
1723 G4Thr << 1787 G4ThreeVector *n ) const
1724 { 1788 {
1725 G4double snxt = kInfinity; // snxt is d 1789 G4double snxt = kInfinity; // snxt is default return value
1726 G4double sphi= kInfinity,stheta= kInfinity; 1790 G4double sphi= kInfinity,stheta= kInfinity;
1727 ESide side=kNull,sidephi=kNull,sidetheta=kN << 1791 ESide side=kNull,sidephi=kNull,sidetheta=kNull;
1728 1792
1729 const G4double halfRmaxTolerance = fRmaxTol <<
1730 const G4double halfRminTolerance = fRminTol <<
1731 const G4double Rmax_plus = fRmax + halfRma <<
1732 const G4double Rmin_minus = (fRmin) != 0.0 <<
1733 G4double t1,t2; 1793 G4double t1,t2;
1734 G4double b,c,d; 1794 G4double b,c,d;
1735 1795
1736 // Variables for phi intersection: 1796 // Variables for phi intersection:
1737 1797
>> 1798 G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi;
>> 1799 G4double cPhi,sinCPhi,cosCPhi;
1738 G4double pDistS,compS,pDistE,compE,sphi2,vp 1800 G4double pDistS,compS,pDistE,compE,sphi2,vphi;
>> 1801
>> 1802 G4double rho2,rad2,pDotV2d,pDotV3d,pTheta;
1739 1803
1740 G4double rho2,rad2,pDotV2d,pDotV3d; << 1804 G4double tolSTheta=0.,tolETheta=0.;
1741 <<
1742 G4double xi,yi,zi; // Intersection poi 1805 G4double xi,yi,zi; // Intersection point
1743 1806
1744 // Theta precals << 1807 // G4double Comp; // Phi intersection
1745 // << 1808
1746 G4double rhoSecTheta; << 1809 G4bool segPhi; // Phi flag and precalcs
1747 G4double dist2STheta, dist2ETheta, distThet << 1810 G4double hDPhi,hDPhiOT,hDPhiIT;
1748 G4double d2,sd; << 1811 G4double cosHDPhiOT,cosHDPhiIT;
>> 1812
>> 1813 G4bool segTheta; // Theta flag and precals
>> 1814 G4double tanSTheta=0.,tanETheta, rhoSecTheta;
>> 1815 G4double tanSTheta2=0.,tanETheta2=0.;
>> 1816 G4double dist2STheta,dist2ETheta;
>> 1817 G4double d2,s;
1749 1818
1750 // General Precalcs 1819 // General Precalcs
1751 // <<
1752 rho2 = p.x()*p.x()+p.y()*p.y(); <<
1753 rad2 = rho2+p.z()*p.z(); <<
1754 1820
1755 pDotV2d = p.x()*v.x()+p.y()*v.y(); << 1821 rho2=p.x()*p.x()+p.y()*p.y();
1756 pDotV3d = pDotV2d+p.z()*v.z(); << 1822 rad2=rho2+p.z()*p.z();
>> 1823 pTheta=atan2(sqrt(rho2),p.z());
>> 1824
>> 1825 pDotV2d=p.x()*v.x()+p.y()*v.y();
>> 1826 pDotV3d=pDotV2d+p.z()*v.z();
>> 1827
>> 1828 // Set phi divided flag and precalcs
>> 1829
>> 1830 if(fDPhi<2.0*M_PI)
>> 1831 {
>> 1832 segPhi=true;
>> 1833 hDPhi=0.5*fDPhi; // half delta phi
>> 1834 cPhi=fSPhi+hDPhi;;
>> 1835 hDPhiOT=hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi
>> 1836 hDPhiIT=hDPhi-0.5*kAngTolerance;
>> 1837 sinCPhi=sin(cPhi);
>> 1838 cosCPhi=cos(cPhi);
>> 1839 cosHDPhiOT=cos(hDPhiOT);
>> 1840 cosHDPhiIT=cos(hDPhiIT);
>> 1841 }
>> 1842 else
>> 1843 {
>> 1844 segPhi=false;
>> 1845 }
1757 1846
>> 1847 // Theta precalcs
>> 1848
>> 1849 if (fDTheta < M_PI)
>> 1850 {
>> 1851 segTheta=true;
>> 1852 tolSTheta=fSTheta-kAngTolerance*0.5;
>> 1853 tolETheta=fSTheta+fDTheta+kAngTolerance*0.5;
>> 1854 }
>> 1855 else
>> 1856 {
>> 1857 segTheta=false;
>> 1858 }
>> 1859
1758 // Radial Intersections from G4Sphere::Dist 1860 // Radial Intersections from G4Sphere::DistanceToIn
1759 // 1861 //
1760 // Outer spherical shell intersection 1862 // Outer spherical shell intersection
1761 // - Only if outside tolerant fRmax 1863 // - Only if outside tolerant fRmax
1762 // - Check for if inside and outer G4Sphere 1864 // - Check for if inside and outer G4Sphere heading through solid (-> 0)
1763 // - No intersect -> no intersection with G 1865 // - No intersect -> no intersection with G4Sphere
1764 // 1866 //
1765 // Shell eqn: x^2+y^2+z^2=RSPH^2 1867 // Shell eqn: x^2+y^2+z^2=RSPH^2
1766 // 1868 //
1767 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 1869 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
1768 // 1870 //
1769 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz) << 1871 // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
1770 // => rad2 +2sd(pDotV3d) << 1872 // => rad2 +2s(pDotV3d) +s^2 =R^2
>> 1873 //
>> 1874 // => s=-pDotV3d+-sqrt(pDotV3d^2-(rad2-R^2))
1771 // 1875 //
1772 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad << 1876 const G4double fractionTolerance = 1.0e-14;
>> 1877 const G4double flexRadMaxTolerance = std::max(kRadTolerance,
>> 1878 fractionTolerance * fRmax);
>> 1879 const G4double Rmax_plus = fRmax + flexRadMaxTolerance*0.5;
>> 1880 const G4double flexRadMinTolerance = std::max(kRadTolerance,
>> 1881 fractionTolerance * fRmin);
>> 1882 const G4double Rmin_minus= (fRmin>0) ? fRmin-flexRadMinTolerance*0.5 : 0 ;
1773 1883
1774 if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 << 1884 if(rad2 <= Rmax_plus*Rmax_plus && rad2 >= Rmin_minus*Rmin_minus)
1775 { 1885 {
1776 c = rad2 - fRmax*fRmax; << 1886 c=rad2-fRmax*fRmax ;
1777 1887
1778 if (c < fRmaxTolerance*fRmax) << 1888 if (c < flexRadMaxTolerance*fRmax)
1779 { 1889 {
1780 // Within tolerant Outer radius << 1890 // Within tolerant Outer radius
1781 // << 1891 //
1782 // The test is 1892 // The test is
1783 // rad - fRmax < 0.5*kRadTolerance 1893 // rad - fRmax < 0.5*kRadTolerance
1784 // => rad < fRmax + 0.5*kRadTol 1894 // => rad < fRmax + 0.5*kRadTol
1785 // => rad2 < (fRmax + 0.5*kRadTol)^2 1895 // => rad2 < (fRmax + 0.5*kRadTol)^2
1786 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kR 1896 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
1787 // => rad2 - fRmax^2 <~ fRmax*kR << 1897 // => rad2 - fRmax^2 <~ fRmax*kRadTol
1788 1898
1789 d2 = pDotV3d*pDotV3d - c; << 1899 d2=pDotV3d*pDotV3d-c ;
1790 1900
1791 if( (c >- fRmaxTolerance*fRmax) / << 1901 if( (c >- flexRadMaxTolerance*fRmax) // on tolerant surface
1792 && ((pDotV3d >=0) || (d2 < 0)) ) / << 1902 && ((pDotV3d >=0) || (d2 < 0)) ) // leaving outside from Rmax
1793 / << 1903 // not re-entering
1794 { 1904 {
1795 if(calcNorm) 1905 if(calcNorm)
1796 { 1906 {
1797 *validNorm = true ; 1907 *validNorm = true ;
1798 *n = G4ThreeVector(p.x()/fR 1908 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
1799 } 1909 }
1800 return snxt = 0; 1910 return snxt = 0;
1801 } 1911 }
1802 else << 1912 else
1803 { 1913 {
1804 snxt = -pDotV3d+std::sqrt(d2); // << 1914 snxt=-pDotV3d+sqrt(d2); // second root since inside Rmax
1805 side = kRMax ; << 1915 side = kRMax ;
1806 } 1916 }
1807 } 1917 }
1808 1918
1809 // Inner spherical shell intersection: 1919 // Inner spherical shell intersection:
1810 // Always first >=0 root, because would h 1920 // Always first >=0 root, because would have passed
1811 // from outside of Rmin surface . 1921 // from outside of Rmin surface .
1812 1922
1813 if (fRmin != 0.0) << 1923 if (fRmin)
1814 { 1924 {
1815 c = rad2 - fRmin*fRmin; << 1925 c=rad2-fRmin*fRmin ;
1816 d2 = pDotV3d*pDotV3d - c; <<
1817 1926
1818 if (c >- fRminTolerance*fRmin) // 2.0 * << 1927 if (c >- flexRadMinTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
1819 { 1928 {
1820 if ( (c < fRminTolerance*fRmin) << 1929 if( c < flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin
1821 && (d2 >= fRminTolerance*fRmin) && <<
1822 { 1930 {
1823 if(calcNorm) { *validNorm = false; << 1931 if(calcNorm)
>> 1932 {
>> 1933 *validNorm = false ; // Rmin surface is concave
>> 1934 }
1824 return snxt = 0 ; 1935 return snxt = 0 ;
1825 } 1936 }
1826 else 1937 else
1827 { << 1938 {
1828 if ( d2 >= 0. ) << 1939 d2=pDotV3d*pDotV3d-c ;
>> 1940 if (d2 >= 0)
1829 { 1941 {
1830 sd = -pDotV3d-std::sqrt(d2); << 1942 s = -pDotV3d-sqrt(d2) ;
1831 << 1943 if (s>=0) // Always intersect Rmin first
1832 if ( sd >= 0. ) // Always int <<
1833 { 1944 {
1834 snxt = sd ; << 1945 snxt = s ;
1835 side = kRMin ; 1946 side = kRMin ;
1836 } 1947 }
1837 } 1948 }
1838 } 1949 }
1839 } 1950 }
1840 } 1951 }
1841 } 1952 }
1842 1953
1843 // Theta segment intersection 1954 // Theta segment intersection
1844 1955
1845 if ( !fFullThetaSphere ) << 1956 if (segTheta)
1846 { 1957 {
1847 // Intersection with theta surfaces 1958 // Intersection with theta surfaces
1848 // 1959 //
1849 // Known failure cases: 1960 // Known failure cases:
1850 // o Inside tolerance of stheta surface, 1961 // o Inside tolerance of stheta surface, skim
1851 // ~parallel to cone and Hit & enter e 1962 // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1852 // 1963 //
1853 // To solve: Check 2nd root of etheta 1964 // To solve: Check 2nd root of etheta surface in addition to stheta
1854 // 1965 //
1855 // o start/end theta is exactly pi/2 << 1966 // o start/end theta is exactly pi/2
1856 // 1967 //
1857 // Intersections with cones 1968 // Intersections with cones
1858 // 1969 //
1859 // Cone equation: x^2+y^2=z^2tan^2(t) 1970 // Cone equation: x^2+y^2=z^2tan^2(t)
1860 // 1971 //
1861 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan 1972 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1862 // 1973 //
1863 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+p << 1974 // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t))
1864 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = << 1975 // + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1865 // 1976 //
1866 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d << 1977 // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1867 // + (rho2-pz^2tan^2(t)) = 0 <<
1868 // 1978 //
1869 << 1979 tanSTheta=tan(fSTheta);
1870 if(fSTheta != 0.0) // intersection with f << 1980 tanSTheta2=tanSTheta*tanSTheta;
>> 1981 tanETheta=tan(fSTheta+fDTheta);
>> 1982 tanETheta2=tanETheta*tanETheta;
>> 1983
>> 1984 if (fSTheta)
>> 1985 {
>> 1986 dist2STheta=rho2-p.z()*p.z()*tanSTheta2;
>> 1987 }
>> 1988 else
>> 1989 {
>> 1990 dist2STheta = kInfinity;
>> 1991 }
>> 1992 if (fSTheta + fDTheta < M_PI)
>> 1993 {
>> 1994 dist2ETheta = rho2-p.z()*p.z()*tanETheta2;
>> 1995 }
>> 1996 else
1871 { 1997 {
1872 if( std::fabs(tanSTheta) > 5./kAngToler << 1998 dist2ETheta = kInfinity ;
>> 1999 }
>> 2000 if (pTheta > tolSTheta && pTheta < tolETheta) // Inside theta
>> 2001 {
>> 2002 // In tolerance of STheta and possible leaving out to small thetas N-
>> 2003
>> 2004 if(pTheta < tolSTheta + kAngTolerance && fSTheta > kAngTolerance)
1873 { 2005 {
1874 if( v.z() > 0. ) << 2006 t2=pDotV2d-p.z()*v.z()*tanSTheta2 ; // =(VdotN+)*rhoSecSTheta
>> 2007
>> 2008 if( fSTheta < M_PI*0.5 && t2 < 0)
>> 2009 {
>> 2010 if(calcNorm) *validNorm = false ;
>> 2011 return snxt = 0 ;
>> 2012 }
>> 2013 else if(fSTheta > M_PI*0.5 && t2 >= 0)
1875 { 2014 {
1876 if ( std::fabs( p.z() ) <= halfRmax << 2015 if(calcNorm)
1877 { 2016 {
1878 if(calcNorm) << 2017 rhoSecTheta = sqrt(rho2*(1+tanSTheta2)) ;
1879 { << 2018 *validNorm = true ;
1880 *validNorm = true; << 2019 *n = G4ThreeVector(-p.x()/rhoSecTheta, // N-
1881 *n = G4ThreeVector(0.,0.,1.); << 2020 -p.y()/rhoSecTheta,
1882 } << 2021 tanSTheta/sqrt(1+tanSTheta2) ) ;
1883 return snxt = 0 ; <<
1884 } 2022 }
1885 stheta = -p.z()/v.z(); << 2023 return snxt = 0 ;
1886 sidetheta = kSTheta; << 2024 }
>> 2025 else if( fSTheta == M_PI*0.5 && v.z() > 0)
>> 2026 {
>> 2027 if(calcNorm)
>> 2028 {
>> 2029 *validNorm = true ;
>> 2030 *n = G4ThreeVector(0,0,1) ;
>> 2031 }
>> 2032 return snxt = 0 ;
1887 } 2033 }
1888 } 2034 }
1889 else // kons is not plane <<
1890 { <<
1891 t1 = 1-v.z()*v.z()*(1+tanSTh <<
1892 t2 = pDotV2d-p.z()*v.z()*tan <<
1893 dist2STheta = rho2-p.z()*p.z()*tanSTh <<
1894 2035
1895 distTheta = std::sqrt(rho2)-p.z()*tan << 2036 // In tolerance of ETheta and possible leaving out to larger thetas N+
1896 2037
1897 if( std::fabs(t1) < halfAngTolerance << 2038 if ( (pTheta > tolETheta - kAngTolerance)
1898 { << 2039 && (( fSTheta + fDTheta) < M_PI - kAngTolerance) )
1899 if( v.z() > 0. ) << 2040 {
1900 { << 2041 t2=pDotV2d-p.z()*v.z()*tanETheta2 ;
1901 if(std::fabs(distTheta) < halfRma << 2042 if((fSTheta+fDTheta)>M_PI*0.5 && t2<0)
1902 { <<
1903 if( (fSTheta < halfpi) && (p.z( <<
1904 { <<
1905 if( calcNorm ) { *validNorm <<
1906 return snxt = 0.; <<
1907 } <<
1908 else if( (fSTheta > halfpi) && <<
1909 { <<
1910 if( calcNorm ) <<
1911 { <<
1912 *validNorm = true; <<
1913 if (rho2 != 0.0) <<
1914 { <<
1915 rhoSecTheta = std::sqrt(r <<
1916 <<
1917 *n = G4ThreeVector( p.x() <<
1918 p.y() <<
1919 std:: <<
1920 } <<
1921 else *n = G4ThreeVector(0., <<
1922 } <<
1923 return snxt = 0.; <<
1924 } <<
1925 } <<
1926 stheta = -0.5*dist2STheta/t2; <<
1927 sidetheta = kSTheta; <<
1928 } <<
1929 } // 2nd order equation, 1st roo <<
1930 else // 2nd if 1st root -ve <<
1931 { 2043 {
1932 if( std::fabs(distTheta) < halfRmax << 2044 if(calcNorm) *validNorm = false ;
>> 2045 return snxt = 0 ;
>> 2046 }
>> 2047 else if( (fSTheta+fDTheta) < M_PI*0.5 && t2 >= 0 )
>> 2048 {
>> 2049 if(calcNorm)
1933 { 2050 {
1934 if( (fSTheta > halfpi) && (t2 >= << 2051 rhoSecTheta = sqrt(rho2*(1+tanETheta2)) ;
1935 { << 2052 *validNorm = true ;
1936 if( calcNorm ) << 2053 *n = G4ThreeVector( p.x()/rhoSecTheta, // N+
1937 { << 2054 p.y()/rhoSecTheta,
1938 *validNorm = true; << 2055 -tanETheta/sqrt(1+tanETheta2) ) ;
1939 if (rho2 != 0.0) <<
1940 { <<
1941 rhoSecTheta = std::sqrt(rho <<
1942 <<
1943 *n = G4ThreeVector( p.x()/r <<
1944 p.y()/r <<
1945 std::si <<
1946 } <<
1947 else { *n = G4ThreeVector(0. <<
1948 } <<
1949 return snxt = 0.; <<
1950 } <<
1951 else if( (fSTheta < halfpi) && (t <<
1952 { <<
1953 if( calcNorm ) { *validNorm = <<
1954 return snxt = 0.; <<
1955 } <<
1956 } 2056 }
1957 b = t2/t1; << 2057 return snxt = 0 ;
1958 c = dist2STheta/t1; << 2058 }
1959 d2 = b*b - c ; << 2059 else if( ( fSTheta+fDTheta) == M_PI*0.5 && v.z() < 0 )
1960 << 2060 {
1961 if ( d2 >= 0. ) << 2061 if(calcNorm)
1962 { 2062 {
1963 d = std::sqrt(d2); << 2063 *validNorm = true ;
1964 << 2064 *n = G4ThreeVector(0,0,-1) ;
1965 if( fSTheta > halfpi ) <<
1966 { <<
1967 sd = -b - d; // First r <<
1968 <<
1969 if ( ((std::fabs(s) < halfRmaxT <<
1970 || (sd < 0.) || ( (sd > 0.) <<
1971 { <<
1972 sd = -b + d ; // 2nd root <<
1973 } <<
1974 if( (sd > halfRmaxTolerance) && <<
1975 { <<
1976 stheta = sd; <<
1977 sidetheta = kSTheta; <<
1978 } <<
1979 } <<
1980 else // sTheta < pi/2, concave su <<
1981 { <<
1982 sd = -b - d; // First r <<
1983 <<
1984 if ( ( (std::fabs(sd) < halfRma <<
1985 || (sd < 0.) || ( (sd > 0.) & <<
1986 { <<
1987 sd = -b + d ; // 2nd root <<
1988 } <<
1989 if( (sd > halfRmaxTolerance) && <<
1990 { <<
1991 stheta = sd; <<
1992 sidetheta = kSTheta; <<
1993 } <<
1994 } <<
1995 } 2065 }
>> 2066 return snxt = 0 ;
1996 } 2067 }
1997 } 2068 }
1998 } << 2069 if( fSTheta > 0 )
1999 if (eTheta < pi) // intersection with sec << 2070 {
2000 { << 2071 // First root of fSTheta cone, second if first root -ve
2001 if( std::fabs(tanETheta) > 5./kAngToler << 2072
2002 { << 2073 t1 = 1-v.z()*v.z()*(1+tanSTheta2);
2003 if( v.z() < 0. ) << 2074 t2 = pDotV2d-p.z()*v.z()*tanSTheta2;
>> 2075
>> 2076 b = t2/t1;
>> 2077 c = dist2STheta/t1;
>> 2078 d2 = b*b - c ;
>> 2079
>> 2080 if ( d2 >= 0 )
2004 { 2081 {
2005 if ( std::fabs( p.z() ) <= halfRmax << 2082 d = sqrt(d2) ;
>> 2083 s = -b - d ; // First root
>> 2084
>> 2085 if ( s < 0 )
2006 { 2086 {
2007 if(calcNorm) << 2087 s = -b + d ; // Second root
2008 { <<
2009 *validNorm = true; <<
2010 *n = G4ThreeVector(0.,0.,-1.); <<
2011 } <<
2012 return snxt = 0 ; <<
2013 } 2088 }
2014 sd = -p.z()/v.z(); << 2089 if (s > flexRadMaxTolerance*0.5 ) // && s<sr)
2015 <<
2016 if( sd < stheta ) <<
2017 { 2090 {
2018 stheta = sd; << 2091 stheta = s ;
2019 sidetheta = kETheta; << 2092 sidetheta = kSTheta ;
2020 } 2093 }
2021 } 2094 }
2022 } 2095 }
2023 else // kons is not plane <<
2024 { <<
2025 t1 = 1-v.z()*v.z()*(1+tanETh <<
2026 t2 = pDotV2d-p.z()*v.z()*tan <<
2027 dist2ETheta = rho2-p.z()*p.z()*tanETh <<
2028 2096
2029 distTheta = std::sqrt(rho2)-p.z()*tan << 2097 // Possible intersection with ETheta cone
2030 << 2098
2031 if( std::fabs(t1) < halfAngTolerance << 2099 if (fSTheta + fDTheta < M_PI)
2032 { << 2100 {
2033 if( v.z() < 0. ) << 2101 t1 = 1-v.z()*v.z()*(1+tanETheta2);
2034 { << 2102 t2 = pDotV2d-p.z()*v.z()*tanETheta2;
2035 if(std::fabs(distTheta) < halfRma << 2103 b = t2/t1;
2036 { << 2104 c = dist2ETheta/t1;
2037 if( (eTheta > halfpi) && (p.z() << 2105 d2 = b*b-c ;
2038 { <<
2039 if( calcNorm ) { *validNorm <<
2040 return snxt = 0.; <<
2041 } <<
2042 else if ( (eTheta < halfpi) && <<
2043 { <<
2044 if( calcNorm ) <<
2045 { <<
2046 *validNorm = true; <<
2047 if (rho2 != 0.0) <<
2048 { <<
2049 rhoSecTheta = std::sqrt(r <<
2050 *n = G4ThreeVector( p.x() <<
2051 p.y() <<
2052 -sinE <<
2053 } <<
2054 else { *n = G4ThreeVector( <<
2055 } <<
2056 return snxt = 0.; <<
2057 } <<
2058 } <<
2059 sd = -0.5*dist2ETheta/t2; <<
2060 2106
2061 if( sd < stheta ) << 2107 if ( d2 >= 0 )
2062 { <<
2063 stheta = sd; <<
2064 sidetheta = kETheta; <<
2065 } <<
2066 } <<
2067 } // 2nd order equation, 1st roo <<
2068 else // 2nd if 1st root -ve <<
2069 { 2108 {
2070 if ( std::fabs(distTheta) < halfRma << 2109 d = sqrt(d2);
>> 2110 s = -b - d ; // First root
>> 2111
>> 2112 if ( s < 0 )
2071 { 2113 {
2072 if( (eTheta < halfpi) && (t2 >= 0 << 2114 s=-b+d; // Second root
2073 { <<
2074 if( calcNorm ) <<
2075 { <<
2076 *validNorm = true; <<
2077 if (rho2 != 0.0) <<
2078 { <<
2079 rhoSecTheta = std::sqrt(r <<
2080 *n = G4ThreeVector( p.x() <<
2081 p.y() <<
2082 -sinE <<
2083 } <<
2084 else *n = G4ThreeVector(0.,0. <<
2085 } <<
2086 return snxt = 0.; <<
2087 } <<
2088 else if ( (eTheta > halfpi) <<
2089 && (t2 < 0.) && (p.z() <=0 <<
2090 { <<
2091 if( calcNorm ) { *validNorm = <<
2092 return snxt = 0.; <<
2093 } <<
2094 } 2115 }
2095 b = t2/t1; << 2116 if (s > flexRadMaxTolerance*0.5 && s < stheta )
2096 c = dist2ETheta/t1; <<
2097 d2 = b*b - c ; <<
2098 if ( (d2 <halfRmaxTolerance) && (d2 <<
2099 { 2117 {
2100 d2 = 0.; << 2118 stheta = s ;
2101 } << 2119 sidetheta = kETheta ;
2102 if ( d2 >= 0. ) <<
2103 { <<
2104 d = std::sqrt(d2); <<
2105 <<
2106 if( eTheta < halfpi ) <<
2107 { <<
2108 sd = -b - d; // First r <<
2109 <<
2110 if( ((std::fabs(sd) < halfRmaxT <<
2111 || (sd < 0.) ) <<
2112 { <<
2113 sd = -b + d ; // 2nd root <<
2114 } <<
2115 if( sd > halfRmaxTolerance ) <<
2116 { <<
2117 if( sd < stheta ) <<
2118 { <<
2119 stheta = sd; <<
2120 sidetheta = kETheta; <<
2121 } <<
2122 } <<
2123 } <<
2124 else // sTheta+fDTheta > pi/2, co <<
2125 { <<
2126 sd = -b - d; // First r <<
2127 <<
2128 if ( ((std::fabs(sd) < halfRmax <<
2129 || (sd < 0.) <<
2130 || ( (sd > 0.) && (p.z() + sd <<
2131 { <<
2132 sd = -b + d ; // 2nd root <<
2133 } <<
2134 if ( ( sd>halfRmaxTolerance ) <<
2135 && ( p.z()+sd*v.z() <= halfRm <<
2136 { <<
2137 if( sd < stheta ) <<
2138 { <<
2139 stheta = sd; <<
2140 sidetheta = kETheta; <<
2141 } <<
2142 } <<
2143 } <<
2144 } 2120 }
2145 } 2121 }
2146 } 2122 }
2147 } << 2123 }
2148 << 2124 }
2149 } // end theta intersections <<
2150 2125
2151 // Phi Intersection 2126 // Phi Intersection
2152 << 2127
2153 if ( !fFullPhiSphere ) << 2128 if ( fDPhi < 2.0*M_PI)
2154 { 2129 {
2155 if ( (p.x() != 0.0) || (p.y() != 0.0) ) / << 2130 sinSPhi=sin(fSPhi);
>> 2131 cosSPhi=cos(fSPhi);
>> 2132 ePhi=fSPhi+fDPhi;
>> 2133 sinEPhi=sin(ePhi);
>> 2134 cosEPhi=cos(ePhi);
>> 2135 cPhi=fSPhi+fDPhi*0.5;
>> 2136 sinCPhi=sin(cPhi);
>> 2137 cosCPhi=cos(cPhi);
>> 2138
>> 2139 if ( p.x()||p.y() ) // Check if on z axis (rho not needed later)
2156 { 2140 {
2157 // pDist -ve when inside 2141 // pDist -ve when inside
2158 2142
2159 pDistS=p.x()*sinSPhi-p.y()*cosSPhi; 2143 pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
2160 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi; 2144 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
2161 2145
2162 // Comp -ve when in direction of outwar 2146 // Comp -ve when in direction of outwards normal
2163 2147
2164 compS = -sinSPhi*v.x()+cosSPhi*v.y() 2148 compS = -sinSPhi*v.x()+cosSPhi*v.y() ;
2165 compE = sinEPhi*v.x()-cosEPhi*v.y() 2149 compE = sinEPhi*v.x()-cosEPhi*v.y() ;
2166 sidephi = kNull ; 2150 sidephi = kNull ;
2167 2151
2168 if ( (pDistS <= 0) && (pDistE <= 0) ) << 2152 if ( pDistS <= 0 && pDistE <= 0 )
2169 { 2153 {
2170 // Inside both phi *full* planes 2154 // Inside both phi *full* planes
2171 2155
2172 if ( compS < 0 ) 2156 if ( compS < 0 )
2173 { 2157 {
2174 sphi = pDistS/compS ; 2158 sphi = pDistS/compS ;
2175 xi = p.x()+sphi*v.x() ; 2159 xi = p.x()+sphi*v.x() ;
2176 yi = p.y()+sphi*v.y() ; 2160 yi = p.y()+sphi*v.y() ;
2177 2161
2178 // Check intersection with correct << 2162 // Check intersecting with correct half-plane
2179 // << 2163 // (if not -> no intersect)
2180 if( (std::fabs(xi)<=kCarTolerance) << 2164
2181 { << 2165 if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2182 vphi = std::atan2(v.y(),v.x()); <<
2183 sidephi = kSPhi; <<
2184 if ( ( (fSPhi-halfAngTolerance) < <<
2185 && ( (ePhi+halfAngTolerance) > <<
2186 { <<
2187 sphi = kInfinity; <<
2188 } <<
2189 } <<
2190 else if ( ( yi*cosCPhi - xi*sinCPhi <<
2191 { 2166 {
2192 sphi=kInfinity; 2167 sphi=kInfinity;
2193 } 2168 }
2194 else 2169 else
2195 { 2170 {
2196 sidephi = kSPhi ; 2171 sidephi = kSPhi ;
2197 if ( pDistS > -halfCarTolerance) << 2172 if ( pDistS > -0.5*kCarTolerance) sphi =0 ; // Leave by sphi
2198 } 2173 }
2199 } 2174 }
2200 else { sphi = kInfinity; } << 2175 else sphi = kInfinity ;
2201 2176
2202 if ( compE < 0 ) 2177 if ( compE < 0 )
2203 { 2178 {
2204 sphi2=pDistE/compE ; 2179 sphi2=pDistE/compE ;
2205 if (sphi2 < sphi) // Only check fur 2180 if (sphi2 < sphi) // Only check further if < starting phi intersection
2206 { 2181 {
2207 xi = p.x()+sphi2*v.x() ; 2182 xi = p.x()+sphi2*v.x() ;
2208 yi = p.y()+sphi2*v.y() ; 2183 yi = p.y()+sphi2*v.y() ;
2209 2184
2210 // Check intersection with correc << 2185 // Check intersecting with correct half-plane
2211 // << 2186
2212 if ( (std::fabs(xi)<=kCarToleranc << 2187 if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
2213 && (std::fabs(yi)<=kCarToleranc <<
2214 { <<
2215 // Leaving via ending phi <<
2216 // <<
2217 vphi = std::atan2(v.y(),v.x()) <<
2218 <<
2219 if( (fSPhi-halfAngTolerance > v <<
2220 ||(fSPhi+fDPhi+halfAngToler <<
2221 { <<
2222 sidephi = kEPhi; <<
2223 if ( pDistE <= -halfCarTolera <<
2224 else <<
2225 } <<
2226 } <<
2227 else if ((yi*cosCPhi-xi*sinCPhi)> <<
2228 { 2188 {
2229 sidephi = kEPhi ; 2189 sidephi = kEPhi ;
2230 if ( pDistE <= -halfCarToleranc << 2190 if ( pDistE <= -0.5*kCarTolerance )
2231 { 2191 {
2232 sphi=sphi2; 2192 sphi=sphi2;
2233 } 2193 }
2234 else << 2194 else
2235 { 2195 {
2236 sphi = 0 ; 2196 sphi = 0 ;
2237 } 2197 }
2238 } 2198 }
2239 } 2199 }
2240 } << 2200 }
2241 } 2201 }
2242 else if ((pDistS >= 0) && (pDistE >= 0) << 2202 else if ( pDistS >= 0 && pDistE >= 0 ) // Outside both *full* phi planes
2243 { 2203 {
2244 if ( pDistS <= pDistE ) 2204 if ( pDistS <= pDistE )
2245 { 2205 {
2246 sidephi = kSPhi ; 2206 sidephi = kSPhi ;
2247 } 2207 }
2248 else 2208 else
2249 { 2209 {
2250 sidephi = kEPhi ; 2210 sidephi = kEPhi ;
2251 } 2211 }
2252 if ( fDPhi > pi ) << 2212 if ( fDPhi > M_PI )
2253 { 2213 {
2254 if ( (compS < 0) && (compE < 0) ) << 2214 if ( compS < 0 && compE < 0 ) sphi = 0 ;
2255 else << 2215 else sphi = kInfinity ;
2256 } 2216 }
2257 else 2217 else
2258 { 2218 {
2259 // if towards both >=0 then once in 2219 // if towards both >=0 then once inside (after error)
2260 // will remain inside 2220 // will remain inside
2261 2221
2262 if ( (compS >= 0) && (compE >= 0) ) << 2222 if ( compS >= 0 && compE >= 0 )
2263 else << 2223 {
2264 } << 2224 sphi=kInfinity;
>> 2225 }
>> 2226 else
>> 2227 {
>> 2228 sphi=0;
>> 2229 }
>> 2230 }
2265 } 2231 }
2266 else if ( (pDistS > 0) && (pDistE < 0) << 2232 else if ( pDistS > 0 && pDistE < 0 )
2267 { 2233 {
2268 // Outside full starting plane, insid 2234 // Outside full starting plane, inside full ending plane
2269 2235
2270 if ( fDPhi > pi ) << 2236 if ( fDPhi > M_PI )
2271 { 2237 {
2272 if ( compE < 0 ) 2238 if ( compE < 0 )
2273 { 2239 {
2274 sphi = pDistE/compE ; 2240 sphi = pDistE/compE ;
2275 xi = p.x() + sphi*v.x() ; 2241 xi = p.x() + sphi*v.x() ;
2276 yi = p.y() + sphi*v.y() ; 2242 yi = p.y() + sphi*v.y() ;
2277 2243
2278 // Check intersection in correct 2244 // Check intersection in correct half-plane
2279 // (if not -> not leaving phi ext 2245 // (if not -> not leaving phi extent)
2280 // 2246 //
2281 if( (std::fabs(xi)<=kCarTolerance << 2247 if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
2282 { <<
2283 vphi = std::atan2(v.y(),v.x()); <<
2284 sidephi = kSPhi; <<
2285 if ( ( (fSPhi-halfAngTolerance) <<
2286 && ( (ePhi+halfAngTolerance) <<
2287 { <<
2288 sphi = kInfinity; <<
2289 } <<
2290 } <<
2291 else if ( ( yi*cosCPhi - xi*sinCP <<
2292 { 2248 {
2293 sphi = kInfinity ; 2249 sphi = kInfinity ;
2294 } 2250 }
2295 else // Leaving via Ending phi 2251 else // Leaving via Ending phi
2296 { 2252 {
2297 sidephi = kEPhi ; 2253 sidephi = kEPhi ;
2298 if ( pDistE > -halfCarTolerance << 2254 if ( pDistE > -0.5*kCarTolerance ) sphi = 0. ;
2299 } 2255 }
2300 } 2256 }
2301 else 2257 else
2302 { 2258 {
2303 sphi = kInfinity ; 2259 sphi = kInfinity ;
2304 } 2260 }
2305 } 2261 }
2306 else 2262 else
2307 { 2263 {
2308 if ( compS >= 0 ) 2264 if ( compS >= 0 )
2309 { 2265 {
2310 if ( compE < 0 ) 2266 if ( compE < 0 )
2311 { << 2267 {
2312 sphi = pDistE/compE ; 2268 sphi = pDistE/compE ;
2313 xi = p.x() + sphi*v.x() ; 2269 xi = p.x() + sphi*v.x() ;
2314 yi = p.y() + sphi*v.y() ; 2270 yi = p.y() + sphi*v.y() ;
2315 2271
2316 // Check intersection in correc 2272 // Check intersection in correct half-plane
2317 // (if not -> remain in extent) 2273 // (if not -> remain in extent)
2318 // 2274 //
2319 if( (std::fabs(xi)<=kCarToleran << 2275 if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
2320 && (std::fabs(yi)<=kCarToleran <<
2321 { <<
2322 vphi = std::atan2(v.y(),v.x() <<
2323 sidephi = kSPhi; <<
2324 if ( ( (fSPhi-halfAngToleranc <<
2325 && ( (ePhi+halfAngTolerance <<
2326 { <<
2327 sphi = kInfinity; <<
2328 } <<
2329 } <<
2330 else if ( ( yi*cosCPhi - xi*sin <<
2331 { 2276 {
2332 sphi=kInfinity; 2277 sphi=kInfinity;
2333 } 2278 }
2334 else // otherwise leaving via E 2279 else // otherwise leaving via Ending phi
2335 { 2280 {
2336 sidephi = kEPhi ; 2281 sidephi = kEPhi ;
2337 } 2282 }
2338 } 2283 }
2339 else sphi=kInfinity; 2284 else sphi=kInfinity;
2340 } 2285 }
2341 else // leaving immediately by star 2286 else // leaving immediately by starting phi
2342 { 2287 {
2343 sidephi = kSPhi ; 2288 sidephi = kSPhi ;
2344 sphi = 0 ; 2289 sphi = 0 ;
2345 } 2290 }
2346 } 2291 }
2347 } 2292 }
2348 else 2293 else
2349 { 2294 {
2350 // Must be pDistS < 0 && pDistE > 0 2295 // Must be pDistS < 0 && pDistE > 0
2351 // Inside full starting plane, outsid 2296 // Inside full starting plane, outside full ending plane
2352 2297
2353 if ( fDPhi > pi ) << 2298 if ( fDPhi > M_PI )
2354 { 2299 {
2355 if ( compS < 0 ) 2300 if ( compS < 0 )
2356 { 2301 {
2357 sphi=pDistS/compS; 2302 sphi=pDistS/compS;
2358 xi=p.x()+sphi*v.x(); 2303 xi=p.x()+sphi*v.x();
2359 yi=p.y()+sphi*v.y(); 2304 yi=p.y()+sphi*v.y();
2360 2305
2361 // Check intersection in correct 2306 // Check intersection in correct half-plane
2362 // (if not -> not leaving phi ext 2307 // (if not -> not leaving phi extent)
2363 // 2308 //
2364 if( (std::fabs(xi)<=kCarTolerance << 2309 if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2365 { <<
2366 vphi = std::atan2(v.y(),v.x()) <<
2367 sidephi = kSPhi; <<
2368 if ( ( (fSPhi-halfAngTolerance) <<
2369 && ( (ePhi+halfAngTolerance) <<
2370 { <<
2371 sphi = kInfinity; <<
2372 } <<
2373 } <<
2374 else if ( ( yi*cosCPhi - xi*sinCP <<
2375 { 2310 {
2376 sphi = kInfinity ; 2311 sphi = kInfinity ;
2377 } 2312 }
2378 else // Leaving via Starting phi << 2313 else // Leaving via Starting phi
2379 { << 2314 {
2380 sidephi = kSPhi ; 2315 sidephi = kSPhi ;
2381 if ( pDistS > -halfCarTolerance << 2316 if ( pDistS > -0.5*kCarTolerance ) sphi = 0 ;
2382 } 2317 }
2383 } 2318 }
2384 else 2319 else
2385 { 2320 {
2386 sphi = kInfinity ; 2321 sphi = kInfinity ;
2387 } 2322 }
2388 } 2323 }
2389 else 2324 else
2390 { 2325 {
2391 if ( compE >= 0 ) 2326 if ( compE >= 0 )
2392 { 2327 {
2393 if ( compS < 0 ) 2328 if ( compS < 0 )
2394 { 2329 {
2395 sphi = pDistS/compS ; 2330 sphi = pDistS/compS ;
2396 xi = p.x()+sphi*v.x() ; 2331 xi = p.x()+sphi*v.x() ;
2397 yi = p.y()+sphi*v.y() ; 2332 yi = p.y()+sphi*v.y() ;
2398 2333
2399 // Check intersection in correc 2334 // Check intersection in correct half-plane
2400 // (if not -> remain in extent) 2335 // (if not -> remain in extent)
2401 // 2336 //
2402 if( (std::fabs(xi)<=kCarToleran << 2337 if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2403 && (std::fabs(yi)<=kCarToleran <<
2404 { <<
2405 vphi = std::atan2(v.y(),v.x() <<
2406 sidephi = kSPhi; <<
2407 if ( ( (fSPhi-halfAngToleranc <<
2408 && ( (ePhi+halfAngTolerance <<
2409 { <<
2410 sphi = kInfinity; <<
2411 } <<
2412 } <<
2413 else if ( ( yi*cosCPhi - xi*sin <<
2414 { 2338 {
2415 sphi = kInfinity ; 2339 sphi = kInfinity ;
2416 } 2340 }
2417 else // otherwise leaving via S 2341 else // otherwise leaving via Starting phi
2418 { 2342 {
2419 sidephi = kSPhi ; 2343 sidephi = kSPhi ;
2420 } 2344 }
2421 } 2345 }
2422 else 2346 else
2423 { 2347 {
2424 sphi = kInfinity ; 2348 sphi = kInfinity ;
2425 } 2349 }
2426 } 2350 }
2427 else // leaving immediately by endi 2351 else // leaving immediately by ending
2428 { 2352 {
2429 sidephi = kEPhi ; 2353 sidephi = kEPhi ;
2430 sphi = 0 ; 2354 sphi = 0 ;
2431 } 2355 }
2432 } 2356 }
2433 } << 2357 }
2434 } 2358 }
2435 else 2359 else
2436 { 2360 {
2437 // On z axis + travel not || to z axis 2361 // On z axis + travel not || to z axis -> if phi of vector direction
2438 // within phi of shape, Step limited by 2362 // within phi of shape, Step limited by rmax, else Step =0
2439 2363
2440 if ( (v.x() != 0.0) || (v.y() != 0.0) ) << 2364 if ( v.x() || v.y() )
2441 { 2365 {
2442 vphi = std::atan2(v.y(),v.x()) ; << 2366 vphi = atan2(v.y(),v.x()) ;
2443 if ((fSPhi-halfAngTolerance < vphi) & << 2367 if ( fSPhi < vphi && vphi < fSPhi + fDPhi )
2444 { 2368 {
2445 sphi = kInfinity; << 2369 sphi=kInfinity;
2446 } 2370 }
2447 else 2371 else
2448 { 2372 {
2449 sidephi = kSPhi ; // arbitrary << 2373 sidephi = kSPhi ; // arbitrary
2450 sphi = 0 ; 2374 sphi = 0 ;
2451 } 2375 }
2452 } 2376 }
2453 else // travel along z - no phi inters << 2377 else // travel along z - no phi intersaction
2454 { 2378 {
2455 sphi = kInfinity ; 2379 sphi = kInfinity ;
2456 } 2380 }
2457 } 2381 }
2458 if ( sphi < snxt ) // Order intersecttio 2382 if ( sphi < snxt ) // Order intersecttions
2459 { 2383 {
2460 snxt = sphi ; 2384 snxt = sphi ;
2461 side = sidephi ; 2385 side = sidephi ;
2462 } 2386 }
2463 } 2387 }
2464 if (stheta < snxt ) // Order intersections 2388 if (stheta < snxt ) // Order intersections
2465 { 2389 {
2466 snxt = stheta ; 2390 snxt = stheta ;
2467 side = sidetheta ; 2391 side = sidetheta ;
2468 } 2392 }
2469 2393
2470 if (calcNorm) // Output switch operator 2394 if (calcNorm) // Output switch operator
2471 { 2395 {
2472 switch( side ) 2396 switch( side )
2473 { 2397 {
2474 case kRMax: 2398 case kRMax:
2475 xi=p.x()+snxt*v.x(); 2399 xi=p.x()+snxt*v.x();
2476 yi=p.y()+snxt*v.y(); 2400 yi=p.y()+snxt*v.y();
2477 zi=p.z()+snxt*v.z(); 2401 zi=p.z()+snxt*v.z();
2478 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi 2402 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
2479 *validNorm=true; 2403 *validNorm=true;
2480 break; 2404 break;
2481 <<
2482 case kRMin: 2405 case kRMin:
2483 *validNorm=false; // Rmin is concave 2406 *validNorm=false; // Rmin is concave
2484 break; 2407 break;
2485 <<
2486 case kSPhi: 2408 case kSPhi:
2487 if ( fDPhi <= pi ) // Normal to P << 2409 if (fDPhi<=M_PI) // Normal to Phi-
2488 { 2410 {
2489 *n=G4ThreeVector(sinSPhi,-cosSPhi,0 << 2411 *n=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0);
2490 *validNorm=true; 2412 *validNorm=true;
2491 } 2413 }
2492 else { *validNorm=false; } << 2414 else *validNorm=false;
2493 break ; 2415 break ;
2494 <<
2495 case kEPhi: 2416 case kEPhi:
2496 if ( fDPhi <= pi ) // Normal to << 2417 if (fDPhi<=M_PI) // Normal to Phi+
2497 { 2418 {
2498 *n=G4ThreeVector(-sinEPhi,cosEPhi,0 << 2419 *n=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0);
2499 *validNorm=true; 2420 *validNorm=true;
2500 } 2421 }
2501 else { *validNorm=false; } << 2422 else *validNorm=false;
2502 break; 2423 break;
2503 <<
2504 case kSTheta: 2424 case kSTheta:
2505 if( fSTheta == halfpi ) << 2425 if( fSTheta == M_PI*0.5 )
2506 { 2426 {
2507 *n=G4ThreeVector(0.,0.,1.); << 2427 *n=G4ThreeVector(0,0,1);
2508 *validNorm=true; 2428 *validNorm=true;
2509 } 2429 }
2510 else if ( fSTheta > halfpi ) << 2430 else if ( fSTheta > M_PI )
2511 { 2431 {
2512 xi = p.x() + snxt*v.x(); << 2432 xi=p.x()+snxt*v.x();
2513 yi = p.y() + snxt*v.y(); << 2433 yi=p.y()+snxt*v.y();
2514 rho2=xi*xi+yi*yi; << 2434 rhoSecTheta = sqrt((xi*xi+yi*yi)*(1+tanSTheta2)) ;
2515 if (rho2 != 0.0) << 2435 *n = G4ThreeVector(-xi/rhoSecTheta, // N-
2516 { << 2436 -yi/rhoSecTheta,
2517 rhoSecTheta = std::sqrt(rho2*(1+t << 2437 tanSTheta/sqrt(1+tanSTheta2)) ;
2518 *n = G4ThreeVector( xi/rhoSecThet <<
2519 -tanSTheta/std <<
2520 } <<
2521 else <<
2522 { <<
2523 *n = G4ThreeVector(0.,0.,1.); <<
2524 } <<
2525 *validNorm=true; 2438 *validNorm=true;
2526 } 2439 }
2527 else { *validNorm=false; } // Conca << 2440 else *validNorm=false; // Concave STheta cone
2528 break; 2441 break;
2529 <<
2530 case kETheta: 2442 case kETheta:
2531 if( eTheta == halfpi ) << 2443 if( ( fSTheta + fDTheta ) == M_PI*0.5 )
2532 { 2444 {
2533 *n = G4ThreeVector(0.,0.,-1 << 2445 *n = G4ThreeVector(0,0,-1);
2534 *validNorm = true; << 2446 *validNorm = true ;
2535 } 2447 }
2536 else if ( eTheta < halfpi ) << 2448 else if ( ( fSTheta + fDTheta ) < M_PI )
2537 { 2449 {
2538 xi=p.x()+snxt*v.x(); 2450 xi=p.x()+snxt*v.x();
2539 yi=p.y()+snxt*v.y(); 2451 yi=p.y()+snxt*v.y();
2540 rho2=xi*xi+yi*yi; << 2452 rhoSecTheta = sqrt((xi*xi+yi*yi)*(1+tanETheta2)) ;
2541 if (rho2 != 0.0) << 2453 *n = G4ThreeVector( xi/rhoSecTheta, // N+
2542 { << 2454 yi/rhoSecTheta,
2543 rhoSecTheta = std::sqrt(rho2*(1+t << 2455 -tanSTheta/sqrt(1+tanSTheta2) ) ;
2544 *n = G4ThreeVector( xi/rhoSecThet <<
2545 -tanETheta/std <<
2546 } <<
2547 else <<
2548 { <<
2549 *n = G4ThreeVector(0.,0.,-1.); <<
2550 } <<
2551 *validNorm=true; 2456 *validNorm=true;
2552 } 2457 }
2553 else { *validNorm=false; } // Conc << 2458 else *validNorm=false; // Concave ETheta cone
2554 break; 2459 break;
2555 <<
2556 default: 2460 default:
>> 2461 G4cout.precision(16);
2557 G4cout << G4endl; 2462 G4cout << G4endl;
2558 DumpInfo(); 2463 DumpInfo();
2559 std::ostringstream message; << 2464 G4cout << "Position:" << G4endl << G4endl;
2560 G4long oldprc = message.precision(16) << 2465 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
2561 message << "Undefined side for valid << 2466 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
2562 << G4endl << 2467 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
2563 << "Position:" << G4endl << << 2468 G4cout << "Direction:" << G4endl << G4endl;
2564 << "p.x() = " << p.x()/mm < << 2469 G4cout << "v.x() = " << v.x() << G4endl;
2565 << "p.y() = " << p.y()/mm < << 2470 G4cout << "v.y() = " << v.y() << G4endl;
2566 << "p.z() = " << p.z()/mm < << 2471 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
2567 << "Direction:" << G4endl << << 2472 G4cout << "Proposed distance :" << G4endl << G4endl;
2568 << "v.x() = " << v.x() << G << 2473 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl;
2569 << "v.y() = " << v.y() << G << 2474 G4Exception("G4Sphere::DistanceToOut() - Invalid enum");
2570 << "v.z() = " << v.z() << G <<
2571 << "Proposed distance :" << G <<
2572 << "snxt = " << snxt/mm << <<
2573 message.precision(oldprc); <<
2574 G4Exception("G4Sphere::DistanceToOut( <<
2575 "GeomSolids1002", JustWar <<
2576 break; 2475 break;
2577 } 2476 }
2578 } 2477 }
2579 if (snxt == kInfinity) <<
2580 { <<
2581 G4cout << G4endl; <<
2582 DumpInfo(); <<
2583 std::ostringstream message; <<
2584 G4long oldprc = message.precision(16); <<
2585 message << "Logic error: snxt = kInfinity <<
2586 << "Position:" << G4endl << G4en <<
2587 << "p.x() = " << p.x()/mm << " <<
2588 << "p.y() = " << p.y()/mm << " <<
2589 << "p.z() = " << p.z()/mm << " <<
2590 << "Rp = "<< std::sqrt( p.x()*p.x <<
2591 << " mm" << G4endl << G4endl <<
2592 << "Direction:" << G4endl << G4en <<
2593 << "v.x() = " << v.x() << G4end <<
2594 << "v.y() = " << v.y() << G4end <<
2595 << "v.z() = " << v.z() << G4end <<
2596 << "Proposed distance :" << G4end <<
2597 << "snxt = " << snxt/mm << " m <<
2598 message.precision(oldprc); <<
2599 G4Exception("G4Sphere::DistanceToOut(p,v, <<
2600 "GeomSolids1002", JustWarning <<
2601 } <<
2602 <<
2603 return snxt; 2478 return snxt;
2604 } 2479 }
2605 2480
2606 ///////////////////////////////////////////// 2481 /////////////////////////////////////////////////////////////////////////
2607 // 2482 //
2608 // Calculate distance (<=actual) to closest s << 2483 // Calcluate distance (<=actual) to closest surface of shape from inside
2609 2484
2610 G4double G4Sphere::DistanceToOut( const G4Thr 2485 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
2611 { 2486 {
2612 G4double safe=0.0,safeRMin,safeRMax,safePhi << 2487 G4double safe,safeRMin,safeRMax,safePhi,safeTheta;
2613 G4double rho2,rds,rho; << 2488 G4double rho2,rad,rho;
2614 G4double pTheta,dTheta1 = kInfinity,dTheta2 << 2489 G4double phiC,cosPhiC,sinPhiC,ePhi;
>> 2490 G4double pTheta,dTheta1,dTheta2;
2615 rho2=p.x()*p.x()+p.y()*p.y(); 2491 rho2=p.x()*p.x()+p.y()*p.y();
2616 rds=std::sqrt(rho2+p.z()*p.z()); << 2492 rad=sqrt(rho2+p.z()*p.z());
2617 rho=std::sqrt(rho2); << 2493 rho=sqrt(rho2);
2618 2494
2619 #ifdef G4CSGDEBUG 2495 #ifdef G4CSGDEBUG
2620 if( Inside(p) == kOutside ) 2496 if( Inside(p) == kOutside )
2621 { 2497 {
2622 G4long old_prc = G4cout.precision(16); << 2498 G4cout.precision(16) ;
2623 G4cout << G4endl; << 2499 G4cout << G4endl ;
2624 DumpInfo(); 2500 DumpInfo();
2625 G4cout << "Position:" << G4endl << G4en 2501 G4cout << "Position:" << G4endl << G4endl ;
2626 G4cout << "p.x() = " << p.x()/mm << " 2502 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2627 G4cout << "p.y() = " << p.y()/mm << " 2503 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2628 G4cout << "p.z() = " << p.z()/mm << " 2504 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2629 G4cout.precision(old_prc) ; << 2505 G4cout << "G4Sphere::DistanceToOut(p) - point p is outside ?!" << G4endl ;
2630 G4Exception("G4Sphere::DistanceToOut(p)" << 2506 G4cerr << "G4Sphere::DistanceToOut(p) - point p is outside ?!" << G4endl ;
2631 "GeomSolids1002", JustWarnin <<
2632 } 2507 }
2633 #endif 2508 #endif
2634 2509
2635 // Distance to r shells <<
2636 // 2510 //
2637 safeRMax = fRmax-rds; << 2511 // Distance to r shells
2638 safe = safeRMax; << 2512 //
2639 if (fRmin != 0.0) << 2513 if (fRmin)
>> 2514 {
>> 2515 safeRMin=rad-fRmin;
>> 2516 safeRMax=fRmax-rad;
>> 2517 if (safeRMin<safeRMax)
>> 2518 {
>> 2519 safe=safeRMin;
>> 2520 }
>> 2521 else
>> 2522 {
>> 2523 safe=safeRMax;
>> 2524 }
>> 2525 }
>> 2526 else
2640 { 2527 {
2641 safeRMin = rds-fRmin; << 2528 safe=fRmax-rad;
2642 safe = std::min( safeRMin, safeRMax ); <<
2643 } 2529 }
2644 2530
>> 2531 //
2645 // Distance to phi extent 2532 // Distance to phi extent
2646 // 2533 //
2647 if ( !fFullPhiSphere ) << 2534 if (fDPhi<2.0*M_PI && rho)
2648 { 2535 {
2649 if (rho>0.0) << 2536 phiC=fSPhi+fDPhi*0.5;
2650 { << 2537 cosPhiC=cos(phiC);
2651 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0) << 2538 sinPhiC=sin(phiC);
2652 { << 2539 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
2653 safePhi=-(p.x()*sinSPhi-p.y()*cosS << 2540 {
2654 } << 2541 safePhi=-(p.x()*sin(fSPhi)-p.y()*cos(fSPhi));
2655 else << 2542 }
2656 { << 2543 else
2657 safePhi=(p.x()*sinEPhi-p.y()*cosEP << 2544 {
2658 } << 2545 ePhi=fSPhi+fDPhi;
2659 } << 2546 safePhi=(p.x()*sin(ePhi)-p.y()*cos(ePhi));
2660 else << 2547 }
2661 { << 2548 if (safePhi<safe) safe=safePhi;
2662 safePhi = 0.0; // Distance to both P <<
2663 } <<
2664 // Both cases above can be improved - in <<
2665 // although it may be costlier (good fo <<
2666 <<
2667 safe= std::min(safe, safePhi); <<
2668 } 2549 }
2669 2550
2670 // Distance to Theta extent <<
2671 // 2551 //
2672 if ( !fFullThetaSphere ) << 2552 // Distance to Theta extent
>> 2553 //
>> 2554 if (rad)
2673 { 2555 {
2674 if( rds > 0.0 ) << 2556 pTheta=acos(p.z()/rad);
>> 2557 if (pTheta<0) pTheta+=M_PI;
>> 2558 dTheta1=pTheta-fSTheta;
>> 2559 dTheta2=(fSTheta+fDTheta)-pTheta;
>> 2560 if (dTheta1<dTheta2)
2675 { 2561 {
2676 pTheta=std::acos(p.z()/rds); << 2562 safeTheta=rad*sin(dTheta1);
2677 if (pTheta<0) { pTheta+=pi; } << 2563 if (safe>safeTheta)
2678 if(fSTheta>0.) << 2564 {
2679 { dTheta1=pTheta-fSTheta;} << 2565 safe=safeTheta;
2680 if(eTheta<pi) << 2566 }
2681 { dTheta2=eTheta-pTheta;} <<
2682 <<
2683 safeTheta=rds*std::sin(std::min(dTheta <<
2684 } 2567 }
2685 else 2568 else
2686 { 2569 {
2687 safeTheta= 0.0; << 2570 safeTheta=rad*sin(dTheta2);
2688 // An improvement will be to return << 2571 if (safe>safeTheta)
>> 2572 {
>> 2573 safe=safeTheta;
>> 2574 }
2689 } 2575 }
2690 safe = std::min( safe, safeTheta ); <<
2691 } 2576 }
2692 2577
2693 if (safe<0.0) { safe=0; } << 2578 if (safe<0) safe=0;
2694 // An improvement to return negative answ << 2579 return safe;
2695 <<
2696 return safe; <<
2697 } 2580 }
2698 2581
2699 ///////////////////////////////////////////// 2582 //////////////////////////////////////////////////////////////////////////
2700 // 2583 //
2701 // G4EntityType << 2584 // Create a List containing the transformed vertices
>> 2585 // Ordering [0-3] -fDz cross section
>> 2586 // [4-7] +fDz cross section such that [0] is below [4],
>> 2587 // [1] below [5] etc.
>> 2588 // Note:
>> 2589 // Caller has deletion resposibility
>> 2590 // Potential improvement: For last slice, use actual ending angle
>> 2591 // to avoid rounding error problems.
>> 2592
>> 2593 G4ThreeVectorList*
>> 2594 G4Sphere::CreateRotatedVertices( const G4AffineTransform& pTransform,
>> 2595 G4int& noPolygonVertices ) const
>> 2596 {
>> 2597 G4ThreeVectorList *vertices;
>> 2598 G4ThreeVector vertex;
>> 2599 G4double meshAnglePhi,meshRMax,crossAnglePhi,
>> 2600 coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
>> 2601 G4double meshTheta,crossTheta,startTheta;
>> 2602 G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
>> 2603 G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;
>> 2604
>> 2605 // Phi cross sections
>> 2606
>> 2607 noPhiCrossSections=G4int (fDPhi/kMeshAngleDefault)+1;
>> 2608
>> 2609 if (noPhiCrossSections<kMinMeshSections)
>> 2610 {
>> 2611 noPhiCrossSections=kMinMeshSections;
>> 2612 }
>> 2613 else if (noPhiCrossSections>kMaxMeshSections)
>> 2614 {
>> 2615 noPhiCrossSections=kMaxMeshSections;
>> 2616 }
>> 2617 meshAnglePhi=fDPhi/(noPhiCrossSections-1);
>> 2618
>> 2619 // If complete in phi, set start angle such that mesh will be at fRMax
>> 2620 // on the x axis. Will give better extent calculations when not rotated.
>> 2621
>> 2622 if (fDPhi==M_PI*2.0 && fSPhi==0)
>> 2623 {
>> 2624 sAnglePhi = -meshAnglePhi*0.5;
>> 2625 }
>> 2626 else
>> 2627 {
>> 2628 sAnglePhi=fSPhi;
>> 2629 }
2702 2630
2703 G4GeometryType G4Sphere::GetEntityType() cons << 2631 // Theta cross sections
2704 { << 2632
2705 return {"G4Sphere"}; << 2633 noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1;
>> 2634
>> 2635 if (noThetaSections<kMinMeshSections)
>> 2636 {
>> 2637 noThetaSections=kMinMeshSections;
>> 2638 }
>> 2639 else if (noThetaSections>kMaxMeshSections)
>> 2640 {
>> 2641 noThetaSections=kMaxMeshSections;
>> 2642 }
>> 2643 meshTheta=fDTheta/(noThetaSections-1);
>> 2644
>> 2645 // If complete in Theta, set start angle such that mesh will be at fRMax
>> 2646 // on the z axis. Will give better extent calculations when not rotated.
>> 2647
>> 2648 if (fDTheta==M_PI && fSTheta==0)
>> 2649 {
>> 2650 startTheta = -meshTheta*0.5;
>> 2651 }
>> 2652 else
>> 2653 {
>> 2654 startTheta=fSTheta;
>> 2655 }
>> 2656
>> 2657 meshRMax = (meshAnglePhi >= meshTheta) ?
>> 2658 fRmax/cos(meshAnglePhi*0.5) : fRmax/cos(meshTheta*0.5);
>> 2659 G4double* cosCrossTheta = new G4double[noThetaSections];
>> 2660 G4double* sinCrossTheta = new G4double[noThetaSections];
>> 2661 vertices=new G4ThreeVectorList();
>> 2662 vertices->reserve(noPhiCrossSections*(noThetaSections*2));
>> 2663 if (vertices && cosCrossTheta && sinCrossTheta)
>> 2664 {
>> 2665 for (crossSectionPhi=0;
>> 2666 crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
>> 2667 {
>> 2668 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
>> 2669 coscrossAnglePhi=cos(crossAnglePhi);
>> 2670 sincrossAnglePhi=sin(crossAnglePhi);
>> 2671 for (crossSectionTheta=0;
>> 2672 crossSectionTheta<noThetaSections;crossSectionTheta++)
>> 2673 {
>> 2674 // Compute coordinates of cross section at section crossSectionPhi
>> 2675 //
>> 2676 crossTheta=startTheta+crossSectionTheta*meshTheta;
>> 2677 cosCrossTheta[crossSectionTheta]=cos(crossTheta);
>> 2678 sinCrossTheta[crossSectionTheta]=sin(crossTheta);
>> 2679
>> 2680 rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
>> 2681 rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
>> 2682 rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
>> 2683
>> 2684 vertex=G4ThreeVector(rMinX,rMinY,rMinZ);
>> 2685 vertices->push_back(pTransform.TransformPoint(vertex));
>> 2686
>> 2687 } // Theta forward
>> 2688
>> 2689 for (crossSectionTheta=noThetaSections-1;
>> 2690 crossSectionTheta>=0; crossSectionTheta--)
>> 2691 {
>> 2692 rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
>> 2693 rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
>> 2694 rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
>> 2695
>> 2696 vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ);
>> 2697 vertices->push_back(pTransform.TransformPoint(vertex));
>> 2698
>> 2699 } // Theta back
>> 2700 } // Phi
>> 2701 noPolygonVertices = noThetaSections*2 ;
>> 2702 }
>> 2703 else
>> 2704 {
>> 2705 DumpInfo();
>> 2706 G4Exception("G4Sphere::CreateRotatedVertices() - Out of memory !");
>> 2707 }
>> 2708
>> 2709 delete[] cosCrossTheta;
>> 2710 delete[] sinCrossTheta;
>> 2711
>> 2712 return vertices;
2706 } 2713 }
2707 2714
2708 ///////////////////////////////////////////// 2715 //////////////////////////////////////////////////////////////////////////
2709 // 2716 //
2710 // Make a clone of the object << 2717 // G4EntityType
2711 // << 2718
2712 G4VSolid* G4Sphere::Clone() const << 2719 G4GeometryType G4Sphere::GetEntityType() const
2713 { 2720 {
2714 return new G4Sphere(*this); << 2721 return G4String("G4Sphere");
2715 } 2722 }
2716 2723
2717 ///////////////////////////////////////////// 2724 //////////////////////////////////////////////////////////////////////////
2718 // 2725 //
2719 // Stream object contents to an output stream 2726 // Stream object contents to an output stream
2720 2727
2721 std::ostream& G4Sphere::StreamInfo( std::ostr 2728 std::ostream& G4Sphere::StreamInfo( std::ostream& os ) const
2722 { 2729 {
2723 G4long oldprc = os.precision(16); <<
2724 os << "------------------------------------ 2730 os << "-----------------------------------------------------------\n"
2725 << " *** Dump for solid - " << GetNam 2731 << " *** Dump for solid - " << GetName() << " ***\n"
2726 << " ================================ 2732 << " ===================================================\n"
2727 << " Solid type: G4Sphere\n" 2733 << " Solid type: G4Sphere\n"
2728 << " Parameters: \n" 2734 << " Parameters: \n"
2729 << " inner radius: " << fRmin/mm << " 2735 << " inner radius: " << fRmin/mm << " mm \n"
2730 << " outer radius: " << fRmax/mm << " 2736 << " outer radius: " << fRmax/mm << " mm \n"
2731 << " starting phi of segment : " << 2737 << " starting phi of segment : " << fSPhi/degree << " degrees \n"
2732 << " delta phi of segment : " << 2738 << " delta phi of segment : " << fDPhi/degree << " degrees \n"
2733 << " starting theta of segment: " << 2739 << " starting theta of segment: " << fSTheta/degree << " degrees \n"
2734 << " delta theta of segment : " << 2740 << " delta theta of segment : " << fDTheta/degree << " degrees \n"
2735 << "------------------------------------ 2741 << "-----------------------------------------------------------\n";
2736 os.precision(oldprc); <<
2737 2742
2738 return os; 2743 return os;
2739 } 2744 }
2740 2745
2741 ///////////////////////////////////////////// <<
2742 // <<
2743 // Get volume <<
2744 <<
2745 G4double G4Sphere::GetCubicVolume() <<
2746 { <<
2747 if (fCubicVolume == 0.) <<
2748 { <<
2749 G4double RRR = fRmax*fRmax*fRmax; <<
2750 G4double rrr = fRmin*fRmin*fRmin; <<
2751 fCubicVolume = fDPhi*(cosSTheta - cosEThe <<
2752 } <<
2753 return fCubicVolume; <<
2754 } <<
2755 <<
2756 ///////////////////////////////////////////// <<
2757 // <<
2758 // Get surface area <<
2759 <<
2760 G4double G4Sphere::GetSurfaceArea() <<
2761 { <<
2762 if (fSurfaceArea == 0.) <<
2763 { <<
2764 G4double RR = fRmax*fRmax; <<
2765 G4double rr = fRmin*fRmin; <<
2766 fSurfaceArea = fDPhi*(RR + rr)*(cosSTheta <<
2767 if (!fFullPhiSphere) fSurfaceArea += f <<
2768 if (fSTheta > 0) fSurfaceArea += 0 <<
2769 if (eTheta < CLHEP::pi) fSurfaceArea += 0 <<
2770 } <<
2771 return fSurfaceArea; <<
2772 } <<
2773 <<
2774 ///////////////////////////////////////////// <<
2775 // <<
2776 // Return a point randomly and uniformly sele <<
2777 <<
2778 G4ThreeVector G4Sphere::GetPointOnSurface() c <<
2779 { <<
2780 G4double RR = fRmax*fRmax; <<
2781 G4double rr = fRmin*fRmin; <<
2782 <<
2783 // Find surface areas <<
2784 // <<
2785 G4double aInner = fDPhi*rr*(cosSTheta - c <<
2786 G4double aOuter = fDPhi*RR*(cosSTheta - c <<
2787 G4double aPhi = (!fFullPhiSphere) ? fDT <<
2788 G4double aSTheta = (fSTheta > 0) ? 0.5*fDP <<
2789 G4double aETheta = (eTheta < pi) ? 0.5*fDP <<
2790 G4double aTotal = aInner + aOuter + aPhi <<
2791 <<
2792 // Select surface and generate a point <<
2793 // <<
2794 G4double select = aTotal*G4QuickRand(); <<
2795 G4double u = G4QuickRand(); <<
2796 G4double v = G4QuickRand(); <<
2797 if (select < aInner + aOuter) // <<
2798 { <<
2799 G4double r = (select < aInner) ? fRmin <<
2800 G4double z = cosSTheta + (cosETheta - c <<
2801 G4double rho = std::sqrt(1. - z*z); <<
2802 G4double phi = fDPhi*v + fSPhi; <<
2803 return { r*rho*std::cos(phi), r*rho*std:: <<
2804 } <<
2805 else if (select < aInner + aOuter + aPhi) / <<
2806 { <<
2807 G4double phi = (select < aInner + aOute <<
2808 G4double r = std::sqrt((RR - rr)*u + <<
2809 G4double theta = fDTheta*v + fSTheta; <<
2810 G4double z = std::cos(theta); <<
2811 G4double rho = std::sin(theta); <<
2812 return { r*rho*std::cos(phi), r*rho*std:: <<
2813 } <<
2814 else // <<
2815 { <<
2816 G4double theta = (select < aTotal - aEThe <<
2817 G4double r = std::sqrt((RR - rr)*u + <<
2818 G4double phi = fDPhi*v + fSPhi; <<
2819 G4double z = std::cos(theta); <<
2820 G4double rho = std::sin(theta); <<
2821 return { r*rho*std::cos(phi), r*rho*std:: <<
2822 } <<
2823 } <<
2824 <<
2825 ///////////////////////////////////////////// 2746 /////////////////////////////////////////////////////////////////////////////
2826 // 2747 //
2827 // Methods for visualisation 2748 // Methods for visualisation
2828 2749
2829 G4VisExtent G4Sphere::GetExtent() const <<
2830 { <<
2831 return { -fRmax, fRmax,-fRmax, fRmax,-fRmax <<
2832 } <<
2833 <<
2834 <<
2835 void G4Sphere::DescribeYourselfTo ( G4VGraphi 2750 void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
2836 { 2751 {
2837 scene.AddSolid (*this); << 2752 scene.AddThis (*this);
2838 } 2753 }
2839 2754
2840 G4Polyhedron* G4Sphere::CreatePolyhedron () c 2755 G4Polyhedron* G4Sphere::CreatePolyhedron () const
2841 { 2756 {
2842 return new G4PolyhedronSphere (fRmin, fRmax 2757 return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
2843 } 2758 }
2844 2759
2845 #endif << 2760 G4NURBS* G4Sphere::CreateNURBS () const
>> 2761 {
>> 2762 return new G4NURBSbox (fRmax, fRmax, fRmax); // Box for now!!!
>> 2763 }
>> 2764
>> 2765
>> 2766 // ****************************** End of G4Sphere.cc ****************************************
2846 2767