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