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