| Geant4 Cross Reference |
>> 1 //
1 // ******************************************* 2 // ********************************************************************
2 // * License and Disclaimer 3 // * License and Disclaimer *
3 // * 4 // * *
4 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
9 // * 10 // * *
10 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
16 // * 17 // * *
17 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
23 // ******************************************* 24 // ********************************************************************
24 // 25 //
>> 26 // $Id: G4Orb.cc 69788 2013-05-15 12:06:57Z gcosmo $
>> 27 //
>> 28 // class G4Orb
>> 29 //
25 // Implementation for G4Orb class 30 // Implementation for G4Orb class
26 // 31 //
>> 32 // History:
>> 33 //
>> 34 // 05.04.12 M.Kelsey - GetPointOnSurface() throw flat in cos(theta)
>> 35 // 30.06.04 V.Grichine - bug fixed in DistanceToIn(p,v) on Rmax surface
27 // 20.08.03 V.Grichine - created 36 // 20.08.03 V.Grichine - created
28 // 08.08.17 E.Tcherniaev - complete revision, << 37 //
29 // ------------------------------------------- << 38 //////////////////////////////////////////////////////////////
30 39
31 #include "G4Orb.hh" << 40 #include <assert.h>
32 41
33 #if !defined(G4GEOM_USE_UORB) << 42 #include "G4Orb.hh"
34 43
35 #include "G4TwoVector.hh" <<
36 #include "G4VoxelLimits.hh" 44 #include "G4VoxelLimits.hh"
37 #include "G4AffineTransform.hh" 45 #include "G4AffineTransform.hh"
38 #include "G4GeometryTolerance.hh" 46 #include "G4GeometryTolerance.hh"
39 #include "G4BoundingEnvelope.hh" <<
40 47
41 #include "G4VPVParameterisation.hh" 48 #include "G4VPVParameterisation.hh"
42 49
43 #include "G4RandomDirection.hh" <<
44 #include "Randomize.hh" 50 #include "Randomize.hh"
45 51
>> 52 #include "meshdefs.hh"
>> 53
46 #include "G4VGraphicsScene.hh" 54 #include "G4VGraphicsScene.hh"
47 #include "G4VisExtent.hh" << 55 #include "G4Polyhedron.hh"
>> 56 #include "G4NURBS.hh"
>> 57 #include "G4NURBSbox.hh"
48 58
49 using namespace CLHEP; 59 using namespace CLHEP;
50 60
51 ////////////////////////////////////////////// << 61 // Private enum: Not for external use - used by distanceToOut
>> 62
>> 63 enum ESide {kNull,kRMax};
>> 64
>> 65 // used by normal
>> 66
>> 67 enum ENorm {kNRMax};
>> 68
>> 69
>> 70 ////////////////////////////////////////////////////////////////////////
52 // 71 //
53 // Constructor << 72 // constructor - check positive radius
>> 73 //
54 74
55 G4Orb::G4Orb( const G4String& pName, G4double 75 G4Orb::G4Orb( const G4String& pName, G4double pRmax )
56 : G4CSGSolid(pName), fRmax(pRmax) << 76 : G4CSGSolid(pName), fRmax(pRmax)
57 { 77 {
58 Initialize(); << 78
>> 79 const G4double fEpsilon = 2.e-11; // relative tolerance of fRmax
>> 80
>> 81 G4double kRadTolerance
>> 82 = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
>> 83
>> 84 // Check radius
>> 85 //
>> 86 if ( pRmax < 10*kCarTolerance )
>> 87 {
>> 88 G4Exception("G4Orb::G4Orb()", "GeomSolids0002", FatalException,
>> 89 "Invalid radius > 10*kCarTolerance.");
>> 90 }
>> 91 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax);
>> 92
59 } 93 }
60 94
61 ////////////////////////////////////////////// << 95 ///////////////////////////////////////////////////////////////////////
62 // 96 //
63 // Fake default constructor - sets only member 97 // Fake default constructor - sets only member data and allocates memory
64 // for usage restri << 98 // for usage restricted to object persistency.
65 << 99 //
66 G4Orb::G4Orb( __void__& a ) 100 G4Orb::G4Orb( __void__& a )
67 : G4CSGSolid(a) << 101 : G4CSGSolid(a), fRmax(0.), fRmaxTolerance(0.)
68 { 102 {
69 } 103 }
70 104
71 ////////////////////////////////////////////// << 105 /////////////////////////////////////////////////////////////////////
72 // 106 //
73 // Destructor 107 // Destructor
74 108
75 G4Orb::~G4Orb() = default; << 109 G4Orb::~G4Orb()
>> 110 {
>> 111 }
76 112
77 ////////////////////////////////////////////// 113 //////////////////////////////////////////////////////////////////////////
78 // 114 //
79 // Copy constructor 115 // Copy constructor
80 116
81 G4Orb::G4Orb(const G4Orb&) = default; << 117 G4Orb::G4Orb(const G4Orb& rhs)
>> 118 : G4CSGSolid(rhs), fRmax(rhs.fRmax), fRmaxTolerance(rhs.fRmaxTolerance)
>> 119 {
>> 120 }
82 121
83 ////////////////////////////////////////////// 122 //////////////////////////////////////////////////////////////////////////
84 // 123 //
85 // Assignment operator 124 // Assignment operator
86 125
87 G4Orb& G4Orb::operator = (const G4Orb& rhs) << 126 G4Orb& G4Orb::operator = (const G4Orb& rhs)
88 { 127 {
89 // Check assignment to self 128 // Check assignment to self
90 // 129 //
91 if (this == &rhs) { return *this; } 130 if (this == &rhs) { return *this; }
92 131
93 // Copy base class data 132 // Copy base class data
94 // 133 //
95 G4CSGSolid::operator=(rhs); 134 G4CSGSolid::operator=(rhs);
96 135
97 // Copy data 136 // Copy data
98 // 137 //
99 fRmax = rhs.fRmax; 138 fRmax = rhs.fRmax;
100 halfRmaxTol = rhs.halfRmaxTol; << 139 fRmaxTolerance = rhs.fRmaxTolerance;
101 sqrRmaxPlusTol = rhs.sqrRmaxPlusTol; <<
102 sqrRmaxMinusTol = rhs.sqrRmaxMinusTol; <<
103 140
104 return *this; 141 return *this;
105 } 142 }
106 143
107 ////////////////////////////////////////////// 144 //////////////////////////////////////////////////////////////////////////
108 // 145 //
109 // Check radius and initialize dada members <<
110 <<
111 void G4Orb::Initialize() <<
112 { <<
113 const G4double fEpsilon = 2.e-11; // relati <<
114 <<
115 // Check radius <<
116 // <<
117 if ( fRmax < 10*kCarTolerance ) <<
118 { <<
119 G4Exception("G4Orb::Initialize()", "GeomSo <<
120 "Invalid radius < 10*kCarToler <<
121 } <<
122 halfRmaxTol = 0.5 * std::max(kCarTolerance, <<
123 G4double rmaxPlusTol = fRmax + halfRmaxTol; <<
124 G4double rmaxMinusTol = fRmax - halfRmaxTol; <<
125 sqrRmaxPlusTol = rmaxPlusTol*rmaxPlusTol; <<
126 sqrRmaxMinusTol = rmaxMinusTol*rmaxMinusTol; <<
127 } <<
128 <<
129 ////////////////////////////////////////////// <<
130 // <<
131 // Dispatch to parameterisation for replicatio 146 // Dispatch to parameterisation for replication mechanism dimension
132 // computation & modification << 147 // computation & modification.
133 148
134 void G4Orb::ComputeDimensions( G4VPVPara 149 void G4Orb::ComputeDimensions( G4VPVParameterisation* p,
135 const G4int n, 150 const G4int n,
136 const G4VPhysic 151 const G4VPhysicalVolume* pRep )
137 { 152 {
138 p->ComputeDimensions(*this,n,pRep); 153 p->ComputeDimensions(*this,n,pRep);
139 } 154 }
140 155
141 ////////////////////////////////////////////// << 156 ////////////////////////////////////////////////////////////////////////////
142 // 157 //
143 // Get bounding box << 158 // Calculate extent under transform and specified limit
144 159
145 void G4Orb::BoundingLimits(G4ThreeVector& pMin << 160 G4bool G4Orb::CalculateExtent( const EAxis pAxis,
146 { << 161 const G4VoxelLimits& pVoxelLimit,
147 G4double radius = GetRadius(); << 162 const G4AffineTransform& pTransform,
148 pMin.set(-radius,-radius,-radius); << 163 G4double& pMin, G4double& pMax ) const
149 pMax.set( radius, radius, radius); << 164 {
>> 165 // Compute x/y/z mins and maxs for bounding box respecting limits,
>> 166 // with early returns if outside limits. Then switch() on pAxis,
>> 167 // and compute exact x and y limit for x/y case
>> 168
>> 169 G4double xoffset,xMin,xMax;
>> 170 G4double yoffset,yMin,yMax;
>> 171 G4double zoffset,zMin,zMax;
>> 172
>> 173 G4double diff1,diff2,delta,maxDiff,newMin,newMax;
>> 174 G4double xoff1,xoff2,yoff1,yoff2;
>> 175
>> 176 xoffset=pTransform.NetTranslation().x();
>> 177 xMin=xoffset-fRmax;
>> 178 xMax=xoffset+fRmax;
150 179
151 // Check correctness of the bounding box << 180 if (pVoxelLimit.IsXLimited())
152 // << 181 {
153 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 182 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
154 { << 183 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
155 std::ostringstream message; << 184 {
156 message << "Bad bounding box (min >= max) << 185 return false;
157 << GetName() << " !" << 186 }
158 << "\npMin = " << pMin << 187 else
159 << "\npMax = " << pMax; << 188 {
160 G4Exception("G4Orb::BoundingLimits()", "Ge << 189 if (xMin<pVoxelLimit.GetMinXExtent())
161 JustWarning, message); << 190 {
162 DumpInfo(); << 191 xMin=pVoxelLimit.GetMinXExtent();
163 } << 192 }
164 } << 193 if (xMax>pVoxelLimit.GetMaxXExtent())
>> 194 {
>> 195 xMax=pVoxelLimit.GetMaxXExtent();
>> 196 }
>> 197 }
>> 198 }
>> 199 yoffset=pTransform.NetTranslation().y();
>> 200 yMin=yoffset-fRmax;
>> 201 yMax=yoffset+fRmax;
165 202
166 ////////////////////////////////////////////// << 203 if (pVoxelLimit.IsYLimited())
167 // << 204 {
168 // Calculate extent under transform and specif << 205 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
>> 206 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
>> 207 {
>> 208 return false;
>> 209 }
>> 210 else
>> 211 {
>> 212 if (yMin<pVoxelLimit.GetMinYExtent())
>> 213 {
>> 214 yMin=pVoxelLimit.GetMinYExtent();
>> 215 }
>> 216 if (yMax>pVoxelLimit.GetMaxYExtent())
>> 217 {
>> 218 yMax=pVoxelLimit.GetMaxYExtent();
>> 219 }
>> 220 }
>> 221 }
>> 222 zoffset=pTransform.NetTranslation().z();
>> 223 zMin=zoffset-fRmax;
>> 224 zMax=zoffset+fRmax;
169 225
170 G4bool G4Orb::CalculateExtent(const EAxis pAxi << 226 if (pVoxelLimit.IsZLimited())
171 const G4VoxelLim << 227 {
172 const G4AffineTr << 228 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
173 G4doub << 229 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
174 { << 230 {
175 G4ThreeVector bmin, bmax; << 231 return false;
176 G4bool exist; << 232 }
177 << 233 else
178 // Get bounding box << 234 {
179 BoundingLimits(bmin,bmax); << 235 if (zMin<pVoxelLimit.GetMinZExtent())
180 << 236 {
181 // Check bounding box << 237 zMin=pVoxelLimit.GetMinZExtent();
182 G4BoundingEnvelope bbox(bmin,bmax); << 238 }
183 #ifdef G4BBOX_EXTENT << 239 if (zMax>pVoxelLimit.GetMaxZExtent())
184 return bbox.CalculateExtent(pAxis,pVoxelLimi << 240 {
185 #endif << 241 zMax=pVoxelLimit.GetMaxZExtent();
186 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 242 }
187 { << 243 }
188 return exist = pMin < pMax; << 244 }
189 } <<
190 245
191 // Find bounding envelope and calculate exte << 246 // Known to cut sphere
192 // << 247
193 static const G4int NTHETA = 8; // number of << 248 switch (pAxis)
194 static const G4int NPHI = 16; // number of <<
195 static const G4double sinHalfTheta = std::si <<
196 static const G4double cosHalfTheta = std::co <<
197 static const G4double sinHalfPhi = std::si <<
198 static const G4double cosHalfPhi = std::co <<
199 static const G4double sinStepTheta = 2.*sinH <<
200 static const G4double cosStepTheta = 1. - 2. <<
201 static const G4double sinStepPhi = 2.*sinH <<
202 static const G4double cosStepPhi = 1. - 2. <<
203 <<
204 G4double radius = GetRadius(); <<
205 G4double rtheta = radius/cosHalfTheta; <<
206 G4double rphi = rtheta/cosHalfPhi; <<
207 <<
208 // set reference circle <<
209 G4TwoVector xy[NPHI]; <<
210 G4double sinCurPhi = sinHalfPhi; <<
211 G4double cosCurPhi = cosHalfPhi; <<
212 for (auto & k : xy) <<
213 { <<
214 k.set(cosCurPhi,sinCurPhi); <<
215 G4double sinTmpPhi = sinCurPhi; <<
216 sinCurPhi = sinCurPhi*cosStepPhi + cosCurP <<
217 cosCurPhi = cosCurPhi*cosStepPhi - sinTmpP <<
218 } <<
219 <<
220 // set bounding circles <<
221 G4ThreeVectorList circles[NTHETA]; <<
222 for (auto & circle : circles) { circle.resiz <<
223 <<
224 G4double sinCurTheta = sinHalfTheta; <<
225 G4double cosCurTheta = cosHalfTheta; <<
226 for (auto & circle : circles) <<
227 { <<
228 G4double z = rtheta*cosCurTheta; <<
229 G4double rho = rphi*sinCurTheta; <<
230 for (G4int k=0; k<NPHI; ++k) <<
231 { 249 {
232 circle[k].set(rho*xy[k].x(),rho*xy[k].y( << 250 case kXAxis:
>> 251 yoff1=yoffset-yMin;
>> 252 yoff2=yMax-yoffset;
>> 253
>> 254 if ( yoff1 >= 0 && yoff2 >= 0 )
>> 255 {
>> 256 // Y limits cross max/min x => no change
>> 257 //
>> 258 pMin=xMin;
>> 259 pMax=xMax;
>> 260 }
>> 261 else
>> 262 {
>> 263 // Y limits don't cross max/min x => compute max delta x,
>> 264 // hence new mins/maxs
>> 265 //
>> 266 delta=fRmax*fRmax-yoff1*yoff1;
>> 267 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
>> 268 delta=fRmax*fRmax-yoff2*yoff2;
>> 269 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
>> 270 maxDiff=(diff1>diff2) ? diff1:diff2;
>> 271 newMin=xoffset-maxDiff;
>> 272 newMax=xoffset+maxDiff;
>> 273 pMin=(newMin<xMin) ? xMin : newMin;
>> 274 pMax=(newMax>xMax) ? xMax : newMax;
>> 275 }
>> 276 break;
>> 277 case kYAxis:
>> 278 xoff1=xoffset-xMin;
>> 279 xoff2=xMax-xoffset;
>> 280 if (xoff1>=0&&xoff2>=0)
>> 281 {
>> 282 // X limits cross max/min y => no change
>> 283 //
>> 284 pMin=yMin;
>> 285 pMax=yMax;
>> 286 }
>> 287 else
>> 288 {
>> 289 // X limits don't cross max/min y => compute max delta y,
>> 290 // hence new mins/maxs
>> 291 //
>> 292 delta=fRmax*fRmax-xoff1*xoff1;
>> 293 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
>> 294 delta=fRmax*fRmax-xoff2*xoff2;
>> 295 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
>> 296 maxDiff=(diff1>diff2) ? diff1:diff2;
>> 297 newMin=yoffset-maxDiff;
>> 298 newMax=yoffset+maxDiff;
>> 299 pMin=(newMin<yMin) ? yMin : newMin;
>> 300 pMax=(newMax>yMax) ? yMax : newMax;
>> 301 }
>> 302 break;
>> 303 case kZAxis:
>> 304 pMin=zMin;
>> 305 pMax=zMax;
>> 306 break;
>> 307 default:
>> 308 break;
233 } 309 }
234 G4double sinTmpTheta = sinCurTheta; << 310 pMin -= fRmaxTolerance;
235 sinCurTheta = sinCurTheta*cosStepTheta + c << 311 pMax += fRmaxTolerance;
236 cosCurTheta = cosCurTheta*cosStepTheta - s <<
237 } <<
238 312
239 // set envelope and calculate extent << 313 return true;
240 std::vector<const G4ThreeVectorList *> polyg << 314
241 polygons.resize(NTHETA); <<
242 for (G4int i=0; i<NTHETA; ++i) { polygons[i] <<
243 <<
244 G4BoundingEnvelope benv(bmin,bmax,polygons); <<
245 exist = benv.CalculateExtent(pAxis,pVoxelLim <<
246 return exist; <<
247 } 315 }
248 316
249 ////////////////////////////////////////////// << 317 ///////////////////////////////////////////////////////////////////////////
250 // 318 //
251 // Return whether point is inside/outside/on s << 319 // Return whether point inside/outside/on surface
>> 320 // Split into radius checks
>> 321 //
252 322
253 EInside G4Orb::Inside( const G4ThreeVector& p 323 EInside G4Orb::Inside( const G4ThreeVector& p ) const
254 { 324 {
255 G4double rr = p.mag2(); << 325 G4double rad2,tolRMax;
256 if (rr > sqrRmaxPlusTol) return kOutside; << 326 EInside in;
257 return (rr > sqrRmaxMinusTol) ? kSurface : k << 327
>> 328
>> 329 rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z();
>> 330
>> 331 G4double radius = std::sqrt(rad2);
>> 332
>> 333 // G4double radius = std::sqrt(rad2);
>> 334 // Check radial surface
>> 335 // sets `in'
>> 336
>> 337 tolRMax = fRmax - fRmaxTolerance*0.5;
>> 338
>> 339 if ( radius <= tolRMax ) { in = kInside; }
>> 340 else
>> 341 {
>> 342 tolRMax = fRmax + fRmaxTolerance*0.5;
>> 343 if ( radius <= tolRMax ) { in = kSurface; }
>> 344 else { in = kOutside; }
>> 345 }
>> 346 return in;
258 } 347 }
259 348
260 ////////////////////////////////////////////// << 349 /////////////////////////////////////////////////////////////////////
261 // 350 //
262 // Return unit normal of surface closest to p 351 // Return unit normal of surface closest to p
>> 352 // - note if point on z axis, ignore phi divided sides
>> 353 // - unsafe if point close to z axis a rmin=0 - no explicit checks
263 354
264 G4ThreeVector G4Orb::SurfaceNormal( const G4Th 355 G4ThreeVector G4Orb::SurfaceNormal( const G4ThreeVector& p ) const
265 { 356 {
266 return (1/p.mag())*p; << 357 ENorm side = kNRMax;
267 } << 358 G4ThreeVector norm;
>> 359 G4double radius = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
268 360
269 ////////////////////////////////////////////// << 361 switch (side)
270 // << 362 {
271 // Calculate distance to the surface of the or << 363 case kNRMax:
272 // - return kInfinity if no intersection or << 364 norm = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius);
273 // intersection distance <= tolerance << 365 break;
>> 366 default: // Should never reach this case ...
>> 367 DumpInfo();
>> 368 G4Exception("G4Orb::SurfaceNormal()", "GeomSolids1002", JustWarning,
>> 369 "Undefined side for valid surface normal to solid.");
>> 370 break;
>> 371 }
>> 372
>> 373 return norm;
>> 374 }
>> 375
>> 376 ///////////////////////////////////////////////////////////////////////////////
>> 377 //
>> 378 // Calculate distance to shape from outside, along normalised vector
>> 379 // - return kInfinity if no intersection, or intersection distance <= tolerance
>> 380 //
>> 381 // -> If point is outside outer radius, compute intersection with rmax
>> 382 // - if no intersection return
>> 383 // - if valid phi,theta return intersection Dist
274 384
275 G4double G4Orb::DistanceToIn( const G4ThreeVec 385 G4double G4Orb::DistanceToIn( const G4ThreeVector& p,
276 const G4ThreeVec 386 const G4ThreeVector& v ) const
277 { 387 {
278 // Check if point is on the surface and trav << 388 G4double snxt = kInfinity; // snxt = default return value
279 // << 389
280 G4double rr = p.mag2(); << 390 G4double radius, pDotV3d; // , tolORMax2, tolIRMax2;
281 G4double pv = p.dot(v); << 391 G4double c, d2, sd = kInfinity;
282 if (rr >= sqrRmaxMinusTol && pv >= 0) return << 392
>> 393 const G4double dRmax = 100.*fRmax;
>> 394
>> 395 // General Precalcs
>> 396
>> 397 radius = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
>> 398 pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z();
283 399
284 // Find intersection << 400 // Radial Precalcs
>> 401
>> 402 // tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5);
>> 403 // tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5);
>> 404
>> 405 // Outer spherical shell intersection
>> 406 // - Only if outside tolerant fRmax
>> 407 // - Check for if inside and outer G4Orb heading through solid (-> 0)
>> 408 // - No intersect -> no intersection with G4Orb
>> 409 //
>> 410 // Shell eqn: x^2+y^2+z^2 = RSPH^2
285 // 411 //
286 // Sphere eqn: x^2 + y^2 + z^2 = R^2 << 412 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
287 // 413 //
288 // => (px + t*vx)^2 + (py + t*vy)^2 + (pz << 414 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
289 // => r^2 + 2t(p.v) + t^2 = R^2 << 415 // => rad2 +2sd(pDotV3d) +sd^2 =R^2
290 // => tmin = -(p.v) - Sqrt((p.v)^2 - (r^2 <<
291 // 416 //
292 G4double D = pv*pv - rr + fRmax*fRmax; << 417 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
293 if (D < 0) return kInfinity; // no intersect <<
294 418
295 G4double sqrtD = std::sqrt(D); << 419 c = (radius - fRmax)*(radius + fRmax);
296 G4double dist = -pv - sqrtD; <<
297 420
298 // Avoid rounding errors due to precision is << 421 if( radius > fRmax-fRmaxTolerance*0.5 ) // not inside in terms of Inside(p)
299 // Split long distances and recompute << 422 {
300 // << 423 if ( c > fRmaxTolerance*fRmax )
301 G4double Dmax = 32*fRmax; << 424 {
302 if (dist > Dmax) << 425 // If outside tolerant boundary of outer G4Orb in terms of c
>> 426 // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ]
>> 427
>> 428 d2 = pDotV3d*pDotV3d - c;
>> 429
>> 430 if ( d2 >= 0 )
>> 431 {
>> 432 sd = -pDotV3d - std::sqrt(d2);
>> 433 if ( sd >= 0 )
>> 434 {
>> 435 if ( sd > dRmax ) // Avoid rounding errors due to precision issues seen on
>> 436 { // 64 bits systems. Split long distances and recompute
>> 437 G4double fTerm = sd - std::fmod(sd,dRmax);
>> 438 sd = fTerm + DistanceToIn(p+fTerm*v,v);
>> 439 }
>> 440 return snxt = sd;
>> 441 }
>> 442 }
>> 443 else // No intersection with G4Orb
>> 444 {
>> 445 return snxt = kInfinity;
>> 446 }
>> 447 }
>> 448 else // not outside in terms of c
>> 449 {
>> 450 if ( c > -fRmaxTolerance*fRmax ) // on surface
>> 451 {
>> 452 d2 = pDotV3d*pDotV3d - c;
>> 453 if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
>> 454 {
>> 455 return snxt = kInfinity;
>> 456 }
>> 457 else
>> 458 {
>> 459 return snxt = 0.;
>> 460 }
>> 461 }
>> 462 }
>> 463 }
>> 464 #ifdef G4CSGDEBUG
>> 465 else // inside ???
303 { 466 {
304 dist = dist - 1.e-8*dist - fRmax; // to s << 467 G4Exception("G4Orb::DistanceToIn(p,v)", "GeomSolids1002",
305 dist += DistanceToIn(p + dist*v, v); << 468 JustWarning, "Point p is inside !?");
306 return (dist >= kInfinity) ? kInfinity : d <<
307 } 469 }
>> 470 #endif
308 471
309 if (sqrtD*2 <= halfRmaxTol) return kInfinity << 472 return snxt;
310 return (dist < halfRmaxTol) ? 0. : dist; <<
311 } 473 }
312 474
313 ////////////////////////////////////////////// << 475 //////////////////////////////////////////////////////////////////////
314 // 476 //
315 // Calculate shortest distance to the boundary << 477 // Calculate distance (<= actual) to closest surface of shape from outside
316 // - Return 0 if point is inside << 478 // - Calculate distance to radial plane
>> 479 // - Return 0 if point inside
317 480
318 G4double G4Orb::DistanceToIn( const G4ThreeVec 481 G4double G4Orb::DistanceToIn( const G4ThreeVector& p ) const
319 { 482 {
320 G4double dist = p.mag() - fRmax; << 483 G4double safe = 0.0,
321 return (dist > 0) ? dist : 0.; << 484 radius = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
>> 485 safe = radius - fRmax;
>> 486 if( safe < 0 ) { safe = 0.; }
>> 487 return safe;
322 } 488 }
323 489
324 ////////////////////////////////////////////// << 490 /////////////////////////////////////////////////////////////////////
325 // 491 //
326 // Calculate distance to the surface of the or << 492 // Calculate distance to surface of shape from `inside', allowing for tolerance
327 // find normal at exit point, if required << 493 //
328 // - when leaving the surface, return 0 <<
329 494
330 G4double G4Orb::DistanceToOut( const G4ThreeVe 495 G4double G4Orb::DistanceToOut( const G4ThreeVector& p,
331 const G4ThreeVe 496 const G4ThreeVector& v,
332 const G4bool ca 497 const G4bool calcNorm,
333 G4bool* v << 498 G4bool *validNorm,
334 G4ThreeVe << 499 G4ThreeVector *n ) const
335 { 500 {
336 // Check if point is on the surface and trav << 501 G4double snxt = kInfinity; // ??? snxt is default return value
>> 502 ESide side = kNull;
>> 503
>> 504 G4double rad2,pDotV3d;
>> 505 G4double xi,yi,zi; // Intersection point
>> 506 G4double c,d2;
>> 507
>> 508 rad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z();
>> 509 pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z();
>> 510
>> 511 // Radial Intersection from G4Orb::DistanceToIn
>> 512 //
>> 513 // Outer spherical shell intersection
>> 514 // - Only if outside tolerant fRmax
>> 515 // - Check for if inside and outer G4Orb heading through solid (-> 0)
>> 516 // - No intersect -> no intersection with G4Orb
>> 517 //
>> 518 // Shell eqn: x^2+y^2+z^2=RSPH^2
337 // 519 //
338 G4double rr = p.mag2(); << 520 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
339 G4double pv = p.dot(v); << 521 //
340 if (rr >= sqrRmaxMinusTol && pv > 0) << 522 // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
>> 523 // => rad2 +2s(pDotV3d) +s^2 =R^2
>> 524 //
>> 525 // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
>> 526
>> 527 const G4double Rmax_plus = fRmax + fRmaxTolerance*0.5;
>> 528 G4double radius = std::sqrt(rad2);
>> 529
>> 530 if ( radius <= Rmax_plus )
341 { 531 {
342 if (calcNorm) << 532 c = (radius - fRmax)*(radius + fRmax);
>> 533
>> 534 if ( c < fRmaxTolerance*fRmax )
343 { 535 {
344 *validNorm = true; << 536 // Within tolerant Outer radius
345 *n = p*(1./std::sqrt(rr)); << 537 //
>> 538 // The test is
>> 539 // radius - fRmax < 0.5*fRmaxTolerance
>> 540 // => radius < fRmax + 0.5*kRadTol
>> 541 // => rad2 < (fRmax + 0.5*kRadTol)^2
>> 542 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
>> 543 // => rad2 - fRmax^2 <~ fRmax*kRadTol
>> 544
>> 545 d2 = pDotV3d*pDotV3d - c;
>> 546
>> 547 if( ( c > -fRmaxTolerance*fRmax) && // on tolerant surface
>> 548 ( ( pDotV3d >= 0 ) || ( d2 < 0 )) ) // leaving outside from Rmax
>> 549 // not re-entering
>> 550 {
>> 551 if(calcNorm)
>> 552 {
>> 553 *validNorm = true;
>> 554 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax);
>> 555 }
>> 556 return snxt = 0;
>> 557 }
>> 558 else
>> 559 {
>> 560 snxt = -pDotV3d + std::sqrt(d2); // second root since inside Rmax
>> 561 side = kRMax;
>> 562 }
346 } 563 }
347 return 0.; <<
348 } 564 }
349 << 565 else // p is outside ???
350 // Find intersection << 566 {
351 // << 567 G4cout << G4endl;
352 // Sphere eqn: x^2 + y^2 + z^2 = R^2 << 568 DumpInfo();
353 // << 569 std::ostringstream message;
354 // => (px + t*vx)^2 + (py + t*vy)^2 + (pz << 570 G4int oldprc = message.precision(16);
355 // => r^2 + 2t(p.v) + t^2 = R^2 << 571 message << "Logic error: snxt = kInfinity ???" << G4endl
356 // => tmax = -(p.v) + Sqrt((p.v)^2 - (r^2 << 572 << "Position:" << G4endl << G4endl
357 // << 573 << "p.x() = " << p.x()/mm << " mm" << G4endl
358 G4double D = pv*pv - rr + fRmax*fRmax; << 574 << "p.y() = " << p.y()/mm << " mm" << G4endl
359 G4double tmax = (D <= 0) ? 0. : std::sqrt(D) << 575 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
360 if (tmax < halfRmaxTol) tmax = 0.; << 576 << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
361 if (calcNorm) << 577 << " mm" << G4endl << G4endl
>> 578 << "Direction:" << G4endl << G4endl
>> 579 << "v.x() = " << v.x() << G4endl
>> 580 << "v.y() = " << v.y() << G4endl
>> 581 << "v.z() = " << v.z() << G4endl << G4endl
>> 582 << "Proposed distance :" << G4endl << G4endl
>> 583 << "snxt = " << snxt/mm << " mm" << G4endl;
>> 584 message.precision(oldprc);
>> 585 G4Exception("G4Orb::DistanceToOut(p,v,..)", "GeomSolids1002",
>> 586 JustWarning, message);
>> 587 }
>> 588 if (calcNorm) // Output switch operator
362 { 589 {
363 *validNorm = true; << 590 switch( side )
364 G4ThreeVector ptmax = p + tmax*v; << 591 {
365 *n = ptmax*(1./ptmax.mag()); << 592 case kRMax:
>> 593 xi=p.x()+snxt*v.x();
>> 594 yi=p.y()+snxt*v.y();
>> 595 zi=p.z()+snxt*v.z();
>> 596 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
>> 597 *validNorm=true;
>> 598 break;
>> 599 default:
>> 600 G4cout << G4endl;
>> 601 DumpInfo();
>> 602 std::ostringstream message;
>> 603 G4int oldprc = message.precision(16);
>> 604 message << "Undefined side for valid surface normal to solid."
>> 605 << G4endl
>> 606 << "Position:" << G4endl << G4endl
>> 607 << "p.x() = " << p.x()/mm << " mm" << G4endl
>> 608 << "p.y() = " << p.y()/mm << " mm" << G4endl
>> 609 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
>> 610 << "Direction:" << G4endl << G4endl
>> 611 << "v.x() = " << v.x() << G4endl
>> 612 << "v.y() = " << v.y() << G4endl
>> 613 << "v.z() = " << v.z() << G4endl << G4endl
>> 614 << "Proposed distance :" << G4endl << G4endl
>> 615 << "snxt = " << snxt/mm << " mm" << G4endl;
>> 616 message.precision(oldprc);
>> 617 G4Exception("G4Orb::DistanceToOut(p,v,..)","GeomSolids1002",
>> 618 JustWarning, message);
>> 619 break;
>> 620 }
366 } 621 }
367 return tmax; << 622 return snxt;
368 } 623 }
369 624
370 ////////////////////////////////////////////// << 625 /////////////////////////////////////////////////////////////////////////
371 // 626 //
372 // Calculate distance (<=actual) to closest su 627 // Calculate distance (<=actual) to closest surface of shape from inside
373 628
374 G4double G4Orb::DistanceToOut( const G4ThreeVe 629 G4double G4Orb::DistanceToOut( const G4ThreeVector& p ) const
375 { 630 {
>> 631 G4double safe=0.0,radius = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
>> 632
376 #ifdef G4CSGDEBUG 633 #ifdef G4CSGDEBUG
377 if( Inside(p) == kOutside ) 634 if( Inside(p) == kOutside )
378 { 635 {
379 std::ostringstream message; << 636 G4int oldprc = G4cout.precision(16);
380 G4int oldprc = message.precision(16); << 637 G4cout << G4endl;
381 message << "Point p is outside (!?) of sol << 638 DumpInfo();
382 message << "Position:\n"; << 639 G4cout << "Position:" << G4endl << G4endl;
383 message << " p.x() = " << p.x()/mm << " << 640 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
384 message << " p.y() = " << p.y()/mm << " << 641 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
385 message << " p.z() = " << p.z()/mm << " << 642 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
386 G4cout.precision(oldprc); << 643 G4cout.precision(oldprc);
387 G4Exception("G4Trap::DistanceToOut(p)", "G << 644 G4Exception("G4Orb::DistanceToOut(p)", "GeomSolids1002",
388 JustWarning, message ); << 645 JustWarning, "Point p is outside !?" );
389 DumpInfo(); <<
390 } 646 }
391 #endif 647 #endif
392 G4double dist = fRmax - p.mag(); << 648
393 return (dist > 0) ? dist : 0.; << 649 safe = fRmax - radius;
>> 650 if ( safe < 0. ) safe = 0.;
>> 651 return safe;
394 } 652 }
395 653
396 ////////////////////////////////////////////// 654 //////////////////////////////////////////////////////////////////////////
397 // 655 //
398 // G4EntityType 656 // G4EntityType
399 657
400 G4GeometryType G4Orb::GetEntityType() const 658 G4GeometryType G4Orb::GetEntityType() const
401 { 659 {
402 return {"G4Orb"}; << 660 return G4String("G4Orb");
403 } 661 }
404 662
405 ////////////////////////////////////////////// 663 //////////////////////////////////////////////////////////////////////////
406 // 664 //
407 // Make a clone of the object 665 // Make a clone of the object
408 << 666 //
409 G4VSolid* G4Orb::Clone() const 667 G4VSolid* G4Orb::Clone() const
410 { 668 {
411 return new G4Orb(*this); 669 return new G4Orb(*this);
412 } 670 }
413 671
414 ////////////////////////////////////////////// 672 //////////////////////////////////////////////////////////////////////////
415 // 673 //
416 // Stream object contents to an output stream 674 // Stream object contents to an output stream
417 675
418 std::ostream& G4Orb::StreamInfo( std::ostream& 676 std::ostream& G4Orb::StreamInfo( std::ostream& os ) const
419 { 677 {
420 G4long oldprc = os.precision(16); << 678 G4int oldprc = os.precision(16);
421 os << "------------------------------------- 679 os << "-----------------------------------------------------------\n"
422 << " *** Dump for solid - " << GetName 680 << " *** Dump for solid - " << GetName() << " ***\n"
423 << " ================================= 681 << " ===================================================\n"
424 << " Solid type: G4Orb\n" 682 << " Solid type: G4Orb\n"
425 << " Parameters: \n" 683 << " Parameters: \n"
>> 684
426 << " outer radius: " << fRmax/mm << " 685 << " outer radius: " << fRmax/mm << " mm \n"
427 << "------------------------------------- 686 << "-----------------------------------------------------------\n";
428 os.precision(oldprc); 687 os.precision(oldprc);
>> 688
429 return os; 689 return os;
430 } 690 }
431 691
432 ////////////////////////////////////////////// << 692 /////////////////////////////////////////////////////////////////////////
433 // 693 //
434 // GetPointOnSurface 694 // GetPointOnSurface
435 695
436 G4ThreeVector G4Orb::GetPointOnSurface() const 696 G4ThreeVector G4Orb::GetPointOnSurface() const
437 { 697 {
438 return fRmax * G4RandomDirection(); << 698 // generate a random number from zero to 2pi...
>> 699 //
>> 700 G4double phi = RandFlat::shoot(0.,2.*pi);
>> 701 G4double cosphi = std::cos(phi);
>> 702 G4double sinphi = std::sin(phi);
>> 703
>> 704 // generate a random point uniform in area
>> 705 G4double costheta = RandFlat::shoot(-1.,1.);
>> 706 G4double sintheta = std::sqrt(1.-sqr(costheta));
>> 707
>> 708 return G4ThreeVector (fRmax*sintheta*cosphi,
>> 709 fRmax*sintheta*sinphi, fRmax*costheta);
439 } 710 }
440 711
441 ////////////////////////////////////////////// << 712 ////////////////////////////////////////////////////////////////////////
442 // 713 //
443 // Methods for visualisation 714 // Methods for visualisation
444 715
445 void G4Orb::DescribeYourselfTo ( G4VGraphicsSc 716 void G4Orb::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
446 { 717 {
447 scene.AddSolid (*this); 718 scene.AddSolid (*this);
448 } 719 }
449 720
450 G4VisExtent G4Orb::GetExtent() const <<
451 { <<
452 return {-fRmax, fRmax, -fRmax, fRmax, -fRmax <<
453 } <<
454 <<
455 G4Polyhedron* G4Orb::CreatePolyhedron () const 721 G4Polyhedron* G4Orb::CreatePolyhedron () const
456 { 722 {
457 return new G4PolyhedronSphere (0., fRmax, 0. 723 return new G4PolyhedronSphere (0., fRmax, 0., 2*pi, 0., pi);
458 } 724 }
459 725
460 #endif << 726 G4NURBS* G4Orb::CreateNURBS () const
>> 727 {
>> 728 return new G4NURBSbox (fRmax, fRmax, fRmax); // Box for now!!!
>> 729 }
461 730