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