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