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