| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 //
>> 27 // $Id: G4Box.cc 99469 2016-09-22 15:04:36Z gcosmo $
>> 28 //
>> 29 //
>> 30 //
26 // Implementation for G4Box class 31 // Implementation for G4Box class
27 // 32 //
28 // 30.06.95 - P.Kent: First version << 33 // 24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
29 // 20.09.98 - V.Grichine: new algorithm of Di 34 // 20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
30 // 18.04.17 - E.Tcherniaev: complete revision << 35 // 07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
>> 36 // 09.06.00 - V.Grichine: safety in DistanceToIn(p) against Inside(p)=kOutside
>> 37 // and information before exception in DistanceToOut(p,v,...)
>> 38 // 15.11.00 - D.Williams, V.Grichine: bug fixed in CalculateExtent - change
>> 39 // algorithm for rotated vertices
>> 40 // 23.08.16 - E.Tcherniaev: use G4BoundingEnvelope for CalculateExtent()
>> 41 // 20.09.16 - E.Tcherniaev: added Extent(pmin,pmax)
31 // ------------------------------------------- 42 // --------------------------------------------------------------------
32 43
33 #include "G4Box.hh" 44 #include "G4Box.hh"
34 45
35 #if !defined(G4GEOM_USE_UBOX) 46 #if !defined(G4GEOM_USE_UBOX)
36 47
37 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh"
38 #include "G4VoxelLimits.hh" 49 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 50 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh" 51 #include "G4BoundingEnvelope.hh"
41 #include "G4QuickRand.hh" << 52 #include "Randomize.hh"
42 53
43 #include "G4VPVParameterisation.hh" 54 #include "G4VPVParameterisation.hh"
44 55
45 #include "G4VGraphicsScene.hh" 56 #include "G4VGraphicsScene.hh"
46 #include "G4VisExtent.hh" 57 #include "G4VisExtent.hh"
47 58
48 ////////////////////////////////////////////// 59 ////////////////////////////////////////////////////////////////////////
49 // 60 //
50 // Constructor - check & set half widths 61 // Constructor - check & set half widths
51 62
52 G4Box::G4Box(const G4String& pName, 63 G4Box::G4Box(const G4String& pName,
53 G4double pX, 64 G4double pX,
54 G4double pY, 65 G4double pY,
55 G4double pZ) 66 G4double pZ)
56 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(p 67 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
57 { 68 {
58 delta = 0.5*kCarTolerance; 69 delta = 0.5*kCarTolerance;
59 if (pX < 2*kCarTolerance || << 70 if ( (pX < 2*kCarTolerance)
60 pY < 2*kCarTolerance || << 71 || (pY < 2*kCarTolerance)
61 pZ < 2*kCarTolerance) // limit to thick << 72 || (pZ < 2*kCarTolerance) ) // limit to thickness of surfaces
62 { 73 {
63 std::ostringstream message; 74 std::ostringstream message;
64 message << "Dimensions too small for Solid 75 message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
65 << " hX, hY, hZ = " << pX << " 76 << " hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
66 G4Exception("G4Box::G4Box()", "GeomSolids0 77 G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
67 } 78 }
68 } 79 }
69 80
70 ////////////////////////////////////////////// 81 //////////////////////////////////////////////////////////////////////////
71 // 82 //
72 // Fake default constructor - sets only member 83 // Fake default constructor - sets only member data and allocates memory
73 // for usage restri << 84 // for usage restricted to object persistency.
74 85
75 G4Box::G4Box( __void__& a ) 86 G4Box::G4Box( __void__& a )
76 : G4CSGSolid(a) << 87 : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.), delta(0.)
77 { 88 {
78 } 89 }
79 90
80 ////////////////////////////////////////////// 91 //////////////////////////////////////////////////////////////////////////
81 // 92 //
82 // Destructor 93 // Destructor
83 94
84 G4Box::~G4Box() = default; << 95 G4Box::~G4Box()
>> 96 {
>> 97 }
85 98
86 ////////////////////////////////////////////// 99 //////////////////////////////////////////////////////////////////////////
87 // 100 //
88 // Copy constructor 101 // Copy constructor
89 102
90 G4Box::G4Box(const G4Box&) = default; << 103 G4Box::G4Box(const G4Box& rhs)
>> 104 : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), delta(rhs.delta)
>> 105 {
>> 106 }
91 107
92 ////////////////////////////////////////////// 108 //////////////////////////////////////////////////////////////////////////
93 // 109 //
94 // Assignment operator 110 // Assignment operator
95 111
96 G4Box& G4Box::operator = (const G4Box& rhs) << 112 G4Box& G4Box::operator = (const G4Box& rhs)
97 { 113 {
98 // Check assignment to self 114 // Check assignment to self
99 // 115 //
100 if (this == &rhs) { return *this; } 116 if (this == &rhs) { return *this; }
101 117
102 // Copy base class data 118 // Copy base class data
103 // 119 //
104 G4CSGSolid::operator=(rhs); 120 G4CSGSolid::operator=(rhs);
105 121
106 // Copy data 122 // Copy data
107 // 123 //
108 fDx = rhs.fDx; 124 fDx = rhs.fDx;
109 fDy = rhs.fDy; 125 fDy = rhs.fDy;
110 fDz = rhs.fDz; 126 fDz = rhs.fDz;
111 delta = rhs.delta; 127 delta = rhs.delta;
112 128
113 return *this; 129 return *this;
114 } 130 }
115 131
116 ////////////////////////////////////////////// << 132 //////////////////////////////////////////////////////////////////////////////
117 // <<
118 // Set X dimension <<
119 133
120 void G4Box::SetXHalfLength(G4double dx) 134 void G4Box::SetXHalfLength(G4double dx)
121 { 135 {
122 if(dx > 2*kCarTolerance) // limit to thickn 136 if(dx > 2*kCarTolerance) // limit to thickness of surfaces
123 { 137 {
124 fDx = dx; 138 fDx = dx;
125 } 139 }
126 else 140 else
127 { 141 {
128 std::ostringstream message; 142 std::ostringstream message;
129 message << "Dimension X too small for soli 143 message << "Dimension X too small for solid: " << GetName() << "!"
130 << G4endl 144 << G4endl
131 << " hX = " << dx; 145 << " hX = " << dx;
132 G4Exception("G4Box::SetXHalfLength()", "Ge 146 G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
133 FatalException, message); 147 FatalException, message);
134 } 148 }
135 fCubicVolume = 0.; << 149 fCubicVolume= 0.;
136 fSurfaceArea = 0.; << 150 fSurfaceArea= 0.;
137 fRebuildPolyhedron = true; 151 fRebuildPolyhedron = true;
138 } << 152 }
139 <<
140 ////////////////////////////////////////////// <<
141 // <<
142 // Set Y dimension <<
143 153
144 void G4Box::SetYHalfLength(G4double dy) << 154 void G4Box::SetYHalfLength(G4double dy)
145 { 155 {
146 if(dy > 2*kCarTolerance) // limit to thickn 156 if(dy > 2*kCarTolerance) // limit to thickness of surfaces
147 { 157 {
148 fDy = dy; 158 fDy = dy;
149 } 159 }
150 else 160 else
151 { 161 {
152 std::ostringstream message; 162 std::ostringstream message;
153 message << "Dimension Y too small for soli << 163 message << "Dimension Y too small for solid: " << GetName() << "!"
>> 164 << G4endl
154 << " hY = " << dy; 165 << " hY = " << dy;
155 G4Exception("G4Box::SetYHalfLength()", "Ge 166 G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
156 FatalException, message); 167 FatalException, message);
157 } 168 }
158 fCubicVolume = 0.; << 169 fCubicVolume= 0.;
159 fSurfaceArea = 0.; << 170 fSurfaceArea= 0.;
160 fRebuildPolyhedron = true; 171 fRebuildPolyhedron = true;
161 } << 172 }
162 <<
163 ////////////////////////////////////////////// <<
164 // <<
165 // Set Z dimension <<
166 173
167 void G4Box::SetZHalfLength(G4double dz) << 174 void G4Box::SetZHalfLength(G4double dz)
168 { 175 {
169 if(dz > 2*kCarTolerance) // limit to thickn 176 if(dz > 2*kCarTolerance) // limit to thickness of surfaces
170 { 177 {
171 fDz = dz; 178 fDz = dz;
172 } 179 }
173 else 180 else
174 { 181 {
175 std::ostringstream message; 182 std::ostringstream message;
176 message << "Dimension Z too small for soli << 183 message << "Dimension Z too small for solid: " << GetName() << "!"
>> 184 << G4endl
177 << " hZ = " << dz; 185 << " hZ = " << dz;
178 G4Exception("G4Box::SetZHalfLength()", "Ge 186 G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
179 FatalException, message); 187 FatalException, message);
180 } 188 }
181 fCubicVolume = 0.; << 189 fCubicVolume= 0.;
182 fSurfaceArea = 0.; << 190 fSurfaceArea= 0.;
183 fRebuildPolyhedron = true; 191 fRebuildPolyhedron = true;
184 } << 192 }
185 193
186 ////////////////////////////////////////////// << 194 ////////////////////////////////////////////////////////////////////////
187 // 195 //
188 // Dispatch to parameterisation for replicatio 196 // Dispatch to parameterisation for replication mechanism dimension
189 // computation & modification. 197 // computation & modification.
190 198
191 void G4Box::ComputeDimensions(G4VPVParameteris 199 void G4Box::ComputeDimensions(G4VPVParameterisation* p,
192 const G4int n, 200 const G4int n,
193 const G4VPhysica 201 const G4VPhysicalVolume* pRep)
194 { 202 {
195 p->ComputeDimensions(*this,n,pRep); 203 p->ComputeDimensions(*this,n,pRep);
196 } 204 }
197 205
198 ////////////////////////////////////////////// 206 //////////////////////////////////////////////////////////////////////////
199 // 207 //
200 // Get bounding box 208 // Get bounding box
201 209
202 void G4Box::BoundingLimits(G4ThreeVector& pMin << 210 void G4Box::Extent(G4ThreeVector& pMin, G4ThreeVector& pMax) const
203 { 211 {
204 pMin.set(-fDx,-fDy,-fDz); 212 pMin.set(-fDx,-fDy,-fDz);
205 pMax.set( fDx, fDy, fDz); 213 pMax.set( fDx, fDy, fDz);
206 214
207 // Check correctness of the bounding box 215 // Check correctness of the bounding box
208 // 216 //
209 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 217 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
210 { 218 {
211 std::ostringstream message; 219 std::ostringstream message;
212 message << "Bad bounding box (min >= max) << 220 message << "Bad bounding box (min >= max) for solid: "
213 << GetName() << " !" 221 << GetName() << " !"
214 << "\npMin = " << pMin 222 << "\npMin = " << pMin
215 << "\npMax = " << pMax; 223 << "\npMax = " << pMax;
216 G4Exception("G4Box::BoundingLimits()", "Ge << 224 G4Exception("G4Box::Extent()", "GeomMgt0001", JustWarning, message);
217 DumpInfo(); 225 DumpInfo();
218 } 226 }
219 } << 227 }
220 228
221 ////////////////////////////////////////////// 229 //////////////////////////////////////////////////////////////////////////
222 // 230 //
223 // Calculate extent under transform and specif 231 // Calculate extent under transform and specified limit
224 232
225 G4bool G4Box::CalculateExtent(const EAxis pAxi 233 G4bool G4Box::CalculateExtent(const EAxis pAxis,
226 const G4VoxelLim 234 const G4VoxelLimits& pVoxelLimit,
227 const G4AffineTr 235 const G4AffineTransform& pTransform,
228 G4double& 236 G4double& pMin, G4double& pMax) const
229 { 237 {
230 G4ThreeVector bmin, bmax; 238 G4ThreeVector bmin, bmax;
231 239
232 // Get bounding box 240 // Get bounding box
233 BoundingLimits(bmin,bmax); << 241 Extent(bmin,bmax);
234 242
235 // Find extent 243 // Find extent
236 G4BoundingEnvelope bbox(bmin,bmax); 244 G4BoundingEnvelope bbox(bmin,bmax);
237 return bbox.CalculateExtent(pAxis,pVoxelLimi 245 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
238 } << 246 }
239 247
240 ////////////////////////////////////////////// << 248 /////////////////////////////////////////////////////////////////////////
241 // 249 //
242 // Return whether point inside/outside/on surf 250 // Return whether point inside/outside/on surface, using tolerance
243 251
244 EInside G4Box::Inside(const G4ThreeVector& p) 252 EInside G4Box::Inside(const G4ThreeVector& p) const
245 { 253 {
246 G4double dist = std::max(std::max( << 254 EInside in = kOutside ;
247 std::abs(p.x())-fDx, << 255 G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
248 std::abs(p.y())-fDy), << 256
249 std::abs(p.z())-fDz); << 257 if ( q.x() <= (fDx - delta) )
250 return (dist > delta) ? kOutside : << 258 {
251 ((dist > -delta) ? kSurface : kInside); << 259 if (q.y() <= (fDy - delta) )
>> 260 {
>> 261 if ( q.z() <= (fDz - delta) ) { in = kInside ; }
>> 262 else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
>> 263 }
>> 264 else if ( q.y() <= (fDy + delta) )
>> 265 {
>> 266 if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
>> 267 }
>> 268 }
>> 269 else if ( q.x() <= (fDx + delta) )
>> 270 {
>> 271 if ( q.y() <= (fDy + delta) )
>> 272 {
>> 273 if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
>> 274 }
>> 275 }
>> 276 return in ;
252 } 277 }
253 278
254 ////////////////////////////////////////////// << 279 ///////////////////////////////////////////////////////////////////////
255 // 280 //
256 // Detect the side(s) and return corresponding << 281 // Calculate side nearest to p, and return normal
>> 282 // If two sides are equidistant, normal of first side (x/y/z)
>> 283 // encountered returned
257 284
258 G4ThreeVector G4Box::SurfaceNormal( const G4Th 285 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
259 { 286 {
260 G4ThreeVector norm(0,0,0); << 287 G4double distx, disty, distz ;
261 G4double px = p.x(); << 288 G4ThreeVector norm(0.,0.,0.);
262 if (std::abs(std::abs(px) - fDx) <= delta) n << 289
263 G4double py = p.y(); << 290 // Calculate distances as if in 1st octant
264 if (std::abs(std::abs(py) - fDy) <= delta) n << 291
265 G4double pz = p.z(); << 292 distx = std::fabs(std::fabs(p.x()) - fDx) ;
266 if (std::abs(std::abs(pz) - fDz) <= delta) n << 293 disty = std::fabs(std::fabs(p.y()) - fDy) ;
267 << 294 distz = std::fabs(std::fabs(p.z()) - fDz) ;
268 G4double nside = norm.mag2(); // number of s << 295
269 if (nside == 1) << 296 // New code for particle on surface including edges and corners with specific
270 return norm; << 297 // normals
271 else if (nside > 1) << 298
272 return norm.unit(); // edge or corner << 299 const G4ThreeVector nX = G4ThreeVector( 1.0, 0,0 );
>> 300 const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0 );
>> 301 const G4ThreeVector nY = G4ThreeVector( 0, 1.0,0 );
>> 302 const G4ThreeVector nmY = G4ThreeVector( 0,-1.0,0 );
>> 303 const G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
>> 304 const G4ThreeVector nmZ = G4ThreeVector( 0, 0,- 1.0);
>> 305
>> 306 G4ThreeVector normX(0.,0.,0.), normY(0.,0.,0.), normZ(0.,0.,0.);
>> 307 G4ThreeVector sumnorm(0., 0., 0.);
>> 308 G4int noSurfaces=0;
>> 309
>> 310 if (distx <= delta) // on X/mX surface and around
>> 311 {
>> 312 noSurfaces ++;
>> 313 if ( p.x() >= 0. ) { normX= nX ; } // on +X surface : (1,0,0)
>> 314 else { normX= nmX; } // (-1,0,0)
>> 315 sumnorm= normX;
>> 316 }
>> 317
>> 318 if (disty <= delta) // on one of the +Y or -Y surfaces
>> 319 {
>> 320 noSurfaces ++;
>> 321 if ( p.y() >= 0. ) { normY= nY; } // on +Y surface
>> 322 else { normY= nmY; }
>> 323 sumnorm += normY;
>> 324 }
>> 325
>> 326 if (distz <= delta) // on one of the +Z or -Z surfaces
>> 327 {
>> 328 noSurfaces ++;
>> 329 if ( p.z() >= 0. ) { normZ= nZ; } // on +Z surface
>> 330 else { normZ= nmZ; }
>> 331 sumnorm += normZ;
>> 332 }
>> 333
>> 334 static const G4double invSqrt2 = 1.0 / std::sqrt(2.0);
>> 335 static const G4double invSqrt3 = 1.0 / std::sqrt(3.0);
>> 336
>> 337 if( noSurfaces > 0 )
>> 338 {
>> 339 if( noSurfaces == 1 )
>> 340 {
>> 341 norm= sumnorm;
>> 342 }
>> 343 else
>> 344 {
>> 345 // norm = sumnorm . unit();
>> 346 if( noSurfaces == 2 )
>> 347 {
>> 348 // 2 surfaces -> on edge
>> 349 norm = invSqrt2 * sumnorm;
>> 350 }
>> 351 else
>> 352 {
>> 353 // 3 surfaces (on corner)
>> 354 norm = invSqrt3 * sumnorm;
>> 355 }
>> 356 }
>> 357 }
273 else 358 else
274 { 359 {
275 // Point is not on the surface <<
276 // <<
277 #ifdef G4CSGDEBUG 360 #ifdef G4CSGDEBUG
278 std::ostringstream message; << 361 G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning,
279 G4int oldprc = message.precision(16); << 362 "Point p is not on surface !?" );
280 message << "Point p is not on surface (!?) << 363 #endif
281 << GetName() << G4endl; << 364 norm = ApproxSurfaceNormal(p);
282 message << "Position:\n"; <<
283 message << " p.x() = " << p.x()/mm << " <<
284 message << " p.y() = " << p.y()/mm << " <<
285 message << " p.z() = " << p.z()/mm << " <<
286 G4cout.precision(oldprc); <<
287 G4Exception("G4Box::SurfaceNormal(p)", "Ge <<
288 JustWarning, message ); <<
289 DumpInfo(); <<
290 #endif <<
291 return ApproxSurfaceNormal(p); <<
292 } 365 }
>> 366
>> 367 return norm;
293 } 368 }
294 369
295 ////////////////////////////////////////////// 370 //////////////////////////////////////////////////////////////////////////
296 // 371 //
297 // Algorithm for SurfaceNormal() following the 372 // Algorithm for SurfaceNormal() following the original specification
298 // for points not on the surface 373 // for points not on the surface
299 374
300 G4ThreeVector G4Box::ApproxSurfaceNormal(const << 375 G4ThreeVector G4Box::ApproxSurfaceNormal( const G4ThreeVector& p ) const
301 { 376 {
302 G4double distx = std::abs(p.x()) - fDx; << 377 G4double distx, disty, distz ;
303 G4double disty = std::abs(p.y()) - fDy; << 378 G4ThreeVector norm(0.,0.,0.);
304 G4double distz = std::abs(p.z()) - fDz; << 379
305 << 380 // Calculate distances as if in 1st octant
306 if (distx >= disty && distx >= distz) << 381
307 return {std::copysign(1.,p.x()), 0., 0.}; << 382 distx = std::fabs(std::fabs(p.x()) - fDx) ;
308 if (disty >= distx && disty >= distz) << 383 disty = std::fabs(std::fabs(p.y()) - fDy) ;
309 return {0., std::copysign(1.,p.y()), 0.}; << 384 distz = std::fabs(std::fabs(p.z()) - fDz) ;
>> 385
>> 386 if ( distx <= disty )
>> 387 {
>> 388 if ( distx <= distz ) // Closest to X
>> 389 {
>> 390 if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
>> 391 else { norm = G4ThreeVector( 1.0,0,0) ; }
>> 392 }
>> 393 else // Closest to Z
>> 394 {
>> 395 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
>> 396 else { norm = G4ThreeVector(0,0, 1.0) ; }
>> 397 }
>> 398 }
310 else 399 else
311 return {0., 0., std::copysign(1.,p.z())}; << 400 {
>> 401 if ( disty <= distz ) // Closest to Y
>> 402 {
>> 403 if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
>> 404 else { norm = G4ThreeVector(0, 1.0,0) ; }
>> 405 }
>> 406 else // Closest to Z
>> 407 {
>> 408 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
>> 409 else { norm = G4ThreeVector(0,0, 1.0) ; }
>> 410 }
>> 411 }
>> 412 return norm;
312 } 413 }
313 414
314 ////////////////////////////////////////////// << 415 ///////////////////////////////////////////////////////////////////////////
315 // 416 //
316 // Calculate distance to box from an outside p 417 // Calculate distance to box from an outside point
317 // - return kInfinity if no intersection << 418 // - return kInfinity if no intersection.
>> 419 //
>> 420 // ALGORITHM:
>> 421 //
>> 422 // Check that if point lies outside x/y/z extent of box, travel is towards
>> 423 // the box (ie. there is a possibility of an intersection)
318 // 424 //
>> 425 // Calculate pairs of minimum and maximum distances for x/y/z travel for
>> 426 // intersection with the box's x/y/z extent.
>> 427 // If there is a valid intersection, it is given by the maximum min distance
>> 428 // (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance
>> 429 // (ie. distance after which 1+ of x/y/z intersections not satisfied)
>> 430 //
>> 431 // NOTE:
>> 432 //
>> 433 // `Inside' safe - meaningful answers given if point is inside the exact
>> 434 // shape.
319 435
320 G4double G4Box::DistanceToIn(const G4ThreeVect 436 G4double G4Box::DistanceToIn(const G4ThreeVector& p,
321 const G4ThreeVect 437 const G4ThreeVector& v) const
322 { 438 {
323 // Check if point is on the surface and trav << 439 G4double safx, safy, safz ;
324 // << 440 G4double smin=0.0, sminy, sminz ; // , sminx ;
325 if ((std::abs(p.x()) - fDx) >= -delta && p.x << 441 G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ; // they always > 0
326 if ((std::abs(p.y()) - fDy) >= -delta && p.y << 442 G4double stmp ;
327 if ((std::abs(p.z()) - fDz) >= -delta && p.z << 443 G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
328 444
329 // Find intersection << 445 safx = std::fabs(p.x()) - fDx ; // minimum distance to x surface of shape
330 // << 446 safy = std::fabs(p.y()) - fDy ;
331 G4double invx = (v.x() == 0) ? DBL_MAX : -1. << 447 safz = std::fabs(p.z()) - fDz ;
332 G4double dx = std::copysign(fDx,invx); << 448
333 G4double txmin = (p.x() - dx)*invx; << 449 // Will we intersect?
334 G4double txmax = (p.x() + dx)*invx; << 450 // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent.
335 << 451 // If both p.x/y/z and v.x/y/z repectively are both positive/negative,
336 G4double invy = (v.y() == 0) ? DBL_MAX : -1. << 452 // travel is in a direction away from the shape.
337 G4double dy = std::copysign(fDy,invy); << 453
338 G4double tymin = std::max(txmin,(p.y() - dy) << 454 if ( ((p.x()*v.x() >= 0.0) && (safx > -delta))
339 G4double tymax = std::min(txmax,(p.y() + dy) << 455 || ((p.y()*v.y() >= 0.0) && (safy > -delta))
340 << 456 || ((p.z()*v.z() >= 0.0) && (safz > -delta)) )
341 G4double invz = (v.z() == 0) ? DBL_MAX : -1. << 457 {
342 G4double dz = std::copysign(fDz,invz); << 458 return kInfinity ; // travel away or parallel within tolerance
343 G4double tmin = std::max(tymin,(p.z() - dz)* << 459 }
344 G4double tmax = std::min(tymax,(p.z() + dz)* << 460
>> 461 // Compute min / max distances for x/y/z travel:
>> 462 // X Planes
>> 463
>> 464 if ( v.x() ) // != 0
>> 465 {
>> 466 stmp = 1.0/std::fabs(v.x()) ;
>> 467
>> 468 if (safx >= 0.0)
>> 469 {
>> 470 smin = safx*stmp ;
>> 471 smax = (fDx+std::fabs(p.x()))*stmp ;
>> 472 }
>> 473 else
>> 474 {
>> 475 if (v.x() < 0) { sOut = (fDx + p.x())*stmp ; }
>> 476 else { sOut = (fDx - p.x())*stmp ; }
>> 477 }
>> 478 }
345 479
346 if (tmax <= tmin + delta) return kInfinity; << 480 // Y Planes
347 return (tmin < delta) ? 0. : tmin; << 481
>> 482 if ( v.y() ) // != 0
>> 483 {
>> 484 stmp = 1.0/std::fabs(v.y()) ;
>> 485
>> 486 if (safy >= 0.0)
>> 487 {
>> 488 sminy = safy*stmp ;
>> 489 smaxy = (fDy+std::fabs(p.y()))*stmp ;
>> 490
>> 491 if (sminy > smin) { smin=sminy ; }
>> 492 if (smaxy < smax) { smax=smaxy ; }
>> 493
>> 494 if (smin >= (smax-delta))
>> 495 {
>> 496 return kInfinity ; // touch XY corner
>> 497 }
>> 498 }
>> 499 else
>> 500 {
>> 501 if (v.y() < 0) { sOuty = (fDy + p.y())*stmp ; }
>> 502 else { sOuty = (fDy - p.y())*stmp ; }
>> 503 if( sOuty < sOut ) { sOut = sOuty ; }
>> 504 }
>> 505 }
>> 506
>> 507 // Z planes
>> 508
>> 509 if ( v.z() ) // != 0
>> 510 {
>> 511 stmp = 1.0/std::fabs(v.z()) ;
>> 512
>> 513 if ( safz >= 0.0 )
>> 514 {
>> 515 sminz = safz*stmp ;
>> 516 smaxz = (fDz+std::fabs(p.z()))*stmp ;
>> 517
>> 518 if (sminz > smin) { smin = sminz ; }
>> 519 if (smaxz < smax) { smax = smaxz ; }
>> 520
>> 521 if (smin >= (smax-delta))
>> 522 {
>> 523 return kInfinity ; // touch ZX or ZY corners
>> 524 }
>> 525 }
>> 526 else
>> 527 {
>> 528 if (v.z() < 0) { sOutz = (fDz + p.z())*stmp ; }
>> 529 else { sOutz = (fDz - p.z())*stmp ; }
>> 530 if( sOutz < sOut ) { sOut = sOutz ; }
>> 531 }
>> 532 }
>> 533
>> 534 if (sOut <= (smin + delta)) // travel over edge
>> 535 {
>> 536 return kInfinity ;
>> 537 }
>> 538 if (smin < delta) { smin = 0.0 ; }
>> 539
>> 540 return smin ;
348 } 541 }
349 542
350 ////////////////////////////////////////////// 543 //////////////////////////////////////////////////////////////////////////
351 // << 544 //
352 // Appoximate distance to box. 545 // Appoximate distance to box.
353 // Returns largest perpendicular distance to t 546 // Returns largest perpendicular distance to the closest x/y/z sides of
354 // the box, which is the most fast estimation 547 // the box, which is the most fast estimation of the shortest distance to box
355 // - If inside return 0 548 // - If inside return 0
356 549
357 G4double G4Box::DistanceToIn(const G4ThreeVect 550 G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
358 { 551 {
359 G4double dist = std::max(std::max( << 552 G4double safex, safey, safez, safe = 0.0 ;
360 std::abs(p.x())-fDx, << 553
361 std::abs(p.y())-fDy), << 554 safex = std::fabs(p.x()) - fDx ;
362 std::abs(p.z())-fDz); << 555 safey = std::fabs(p.y()) - fDy ;
363 return (dist > 0) ? dist : 0.; << 556 safez = std::fabs(p.z()) - fDz ;
>> 557
>> 558 if (safex > safe) { safe = safex ; }
>> 559 if (safey > safe) { safe = safey ; }
>> 560 if (safez > safe) { safe = safez ; }
>> 561
>> 562 return safe ;
364 } 563 }
365 564
366 ////////////////////////////////////////////// << 565 /////////////////////////////////////////////////////////////////////////
367 // 566 //
368 // Calculate distance to surface of the box fr << 567 // Calculate distance to surface of box from inside
369 // find normal at exit point, if required << 568 // by calculating distances to box's x/y/z planes.
370 // - when leaving the surface, return 0 << 569 // Smallest distance is exact distance to exiting.
>> 570 // - Eliminate one side of each pair by considering direction of v
>> 571 // - when leaving a surface & v.close, return 0
371 572
372 G4double G4Box::DistanceToOut( const G4ThreeVe << 573 G4double G4Box::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
373 const G4ThreeVe <<
374 const G4bool ca 574 const G4bool calcNorm,
375 G4bool* validNo << 575 G4bool *validNorm,G4ThreeVector *n) const
376 { 576 {
377 // Check if point is on the surface and trav << 577 ESide side = kUndefined ;
378 // << 578 G4double pdist,stmp,snxt=kInfinity;
379 if ((std::abs(p.x()) - fDx) >= -delta && p.x << 579
>> 580 if (calcNorm) { *validNorm = true ; } // All normals are valid
>> 581
>> 582 if (v.x() > 0) // X planes
380 { 583 {
381 if (calcNorm) << 584 pdist = fDx - p.x() ;
>> 585
>> 586 if (pdist > delta)
382 { 587 {
383 *validNorm = true; << 588 snxt = pdist/v.x() ;
384 n->set((p.x() < 0) ? -1. : 1., 0., 0.); << 589 side = kPX ;
>> 590 }
>> 591 else
>> 592 {
>> 593 if (calcNorm) { *n = G4ThreeVector(1,0,0) ; }
>> 594 return snxt = 0 ;
385 } 595 }
386 return 0.; <<
387 } 596 }
388 if ((std::abs(p.y()) - fDy) >= -delta && p.y << 597 else if (v.x() < 0)
389 { 598 {
390 if (calcNorm) << 599 pdist = fDx + p.x() ;
>> 600
>> 601 if (pdist > delta)
391 { 602 {
392 *validNorm = true; << 603 snxt = -pdist/v.x() ;
393 n->set(0., (p.y() < 0) ? -1. : 1., 0.); << 604 side = kMX ;
>> 605 }
>> 606 else
>> 607 {
>> 608 if (calcNorm) { *n = G4ThreeVector(-1,0,0) ; }
>> 609 return snxt = 0 ;
394 } 610 }
395 return 0.; <<
396 } 611 }
397 if ((std::abs(p.z()) - fDz) >= -delta && p.z << 612
>> 613 if (v.y() > 0) // Y planes
398 { 614 {
399 if (calcNorm) << 615 pdist = fDy-p.y();
>> 616
>> 617 if (pdist > delta)
>> 618 {
>> 619 stmp = pdist/v.y();
>> 620
>> 621 if (stmp < snxt)
>> 622 {
>> 623 snxt = stmp;
>> 624 side = kPY;
>> 625 }
>> 626 }
>> 627 else
400 { 628 {
401 *validNorm = true; << 629 if (calcNorm) { *n = G4ThreeVector(0,1,0) ; }
402 n->set(0., 0., (p.z() < 0) ? -1. : 1.); << 630 return snxt = 0 ;
403 } 631 }
404 return 0.; <<
405 } 632 }
>> 633 else if (v.y() < 0)
>> 634 {
>> 635 pdist = fDy + p.y() ;
406 636
407 // Find intersection << 637 if (pdist > delta)
408 // << 638 {
409 G4double vx = v.x(); << 639 stmp = -pdist/v.y();
410 G4double tx = (vx == 0) ? DBL_MAX : (std::co <<
411 640
412 G4double vy = v.y(); << 641 if ( stmp < snxt )
413 G4double ty = (vy == 0) ? tx : (std::copysig << 642 {
414 G4double txy = std::min(tx,ty); << 643 snxt = stmp;
415 << 644 side = kMY;
416 G4double vz = v.z(); << 645 }
417 G4double tz = (vz == 0) ? txy : (std::copysi << 646 }
418 G4double tmax = std::min(txy,tz); << 647 else
>> 648 {
>> 649 if (calcNorm) { *n = G4ThreeVector(0,-1,0) ; }
>> 650 return snxt = 0 ;
>> 651 }
>> 652 }
419 653
420 // Set normal, if required, and return dista << 654 if (v.z() > 0) // Z planes
421 // << 655 {
422 if (calcNorm) << 656 pdist = fDz-p.z();
>> 657
>> 658 if ( pdist > delta )
>> 659 {
>> 660 stmp = pdist/v.z();
>> 661
>> 662 if ( stmp < snxt )
>> 663 {
>> 664 snxt = stmp;
>> 665 side = kPZ;
>> 666 }
>> 667 }
>> 668 else
>> 669 {
>> 670 if (calcNorm) { *n = G4ThreeVector(0,0,1) ; }
>> 671 return snxt = 0 ;
>> 672 }
>> 673 }
>> 674 else if (v.z() < 0)
423 { 675 {
424 *validNorm = true; << 676 pdist = fDz + p.z();
425 if (tmax == tx) n->set((v.x() < 0) ? << 677
426 else if (tmax == ty) n->set(0., (v.y() < 0 << 678 if ( pdist > delta )
427 else n->set(0., 0., (v.z() << 679 {
>> 680 stmp = -pdist/v.z();
>> 681
>> 682 if ( stmp < snxt )
>> 683 {
>> 684 snxt = stmp;
>> 685 side = kMZ;
>> 686 }
>> 687 }
>> 688 else
>> 689 {
>> 690 if (calcNorm) { *n = G4ThreeVector(0,0,-1) ; }
>> 691 return snxt = 0 ;
>> 692 }
428 } 693 }
429 return tmax; << 694
>> 695 if (calcNorm)
>> 696 {
>> 697 switch (side)
>> 698 {
>> 699 case kPX:
>> 700 *n=G4ThreeVector(1,0,0);
>> 701 break;
>> 702 case kMX:
>> 703 *n=G4ThreeVector(-1,0,0);
>> 704 break;
>> 705 case kPY:
>> 706 *n=G4ThreeVector(0,1,0);
>> 707 break;
>> 708 case kMY:
>> 709 *n=G4ThreeVector(0,-1,0);
>> 710 break;
>> 711 case kPZ:
>> 712 *n=G4ThreeVector(0,0,1);
>> 713 break;
>> 714 case kMZ:
>> 715 *n=G4ThreeVector(0,0,-1);
>> 716 break;
>> 717 default:
>> 718 G4cout << G4endl;
>> 719 DumpInfo();
>> 720 std::ostringstream message;
>> 721 G4int oldprc = message.precision(16);
>> 722 message << "Undefined side for valid surface normal to solid."
>> 723 << G4endl
>> 724 << "Position:" << G4endl << G4endl
>> 725 << "p.x() = " << p.x()/mm << " mm" << G4endl
>> 726 << "p.y() = " << p.y()/mm << " mm" << G4endl
>> 727 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
>> 728 << "Direction:" << G4endl << G4endl
>> 729 << "v.x() = " << v.x() << G4endl
>> 730 << "v.y() = " << v.y() << G4endl
>> 731 << "v.z() = " << v.z() << G4endl << G4endl
>> 732 << "Proposed distance :" << G4endl << G4endl
>> 733 << "snxt = " << snxt/mm << " mm" << G4endl;
>> 734 message.precision(oldprc);
>> 735 G4Exception("G4Box::DistanceToOut(p,v,..)", "GeomSolids1002",
>> 736 JustWarning, message);
>> 737 break;
>> 738 }
>> 739 }
>> 740 return snxt;
430 } 741 }
431 742
432 ////////////////////////////////////////////// 743 ////////////////////////////////////////////////////////////////////////////
433 // 744 //
434 // Calculate exact shortest distance to any bo 745 // Calculate exact shortest distance to any boundary from inside
435 // - if outside return 0 << 746 // - If outside return 0
436 747
437 G4double G4Box::DistanceToOut(const G4ThreeVec 748 G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
438 { 749 {
>> 750 G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0;
>> 751
439 #ifdef G4CSGDEBUG 752 #ifdef G4CSGDEBUG
440 if( Inside(p) == kOutside ) 753 if( Inside(p) == kOutside )
441 { 754 {
442 std::ostringstream message; << 755 G4int oldprc = G4cout.precision(16) ;
443 G4int oldprc = message.precision(16); << 756 G4cout << G4endl ;
444 message << "Point p is outside (!?) of sol << 757 DumpInfo();
445 message << "Position:\n"; << 758 G4cout << "Position:" << G4endl << G4endl ;
446 message << " p.x() = " << p.x()/mm << " << 759 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
447 message << " p.y() = " << p.y()/mm << " << 760 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
448 message << " p.z() = " << p.z()/mm << " << 761 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
449 G4cout.precision(oldprc); << 762 G4cout.precision(oldprc) ;
450 G4Exception("G4Box::DistanceToOut(p)", "Ge << 763 G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
451 JustWarning, message ); << 764 JustWarning, "Point p is outside !?" );
452 DumpInfo(); <<
453 } 765 }
454 #endif 766 #endif
455 G4double dist = std::min(std::min( <<
456 fDx-std::abs(p.x()), <<
457 fDy-std::abs(p.y())), <<
458 fDz-std::abs(p.z())); <<
459 return (dist > 0) ? dist : 0.; <<
460 } <<
461 767
462 ////////////////////////////////////////////// << 768 safx1 = fDx - p.x() ;
463 // << 769 safx2 = fDx + p.x() ;
464 // GetEntityType << 770 safy1 = fDy - p.y() ;
>> 771 safy2 = fDy + p.y() ;
>> 772 safz1 = fDz - p.z() ;
>> 773 safz2 = fDz + p.z() ;
>> 774
>> 775 // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
>> 776
>> 777 if (safx2 < safx1) { safe = safx2; }
>> 778 else { safe = safx1; }
>> 779 if (safy1 < safe) { safe = safy1; }
>> 780 if (safy2 < safe) { safe = safy2; }
>> 781 if (safz1 < safe) { safe = safz1; }
>> 782 if (safz2 < safe) { safe = safz2; }
465 783
466 G4GeometryType G4Box::GetEntityType() const << 784 if (safe < 0) { safe = 0 ; }
467 { << 785 return safe ;
468 return {"G4Box"}; <<
469 } 786 }
470 787
471 ////////////////////////////////////////////// 788 //////////////////////////////////////////////////////////////////////////
472 // 789 //
473 // IsFaceted << 790 // GetEntityType
474 791
475 G4bool G4Box::IsFaceted() const << 792 G4GeometryType G4Box::GetEntityType() const
476 { 793 {
477 return true; << 794 return G4String("G4Box");
478 } 795 }
479 796
480 ////////////////////////////////////////////// 797 //////////////////////////////////////////////////////////////////////////
481 // 798 //
482 // Stream object contents to an output stream 799 // Stream object contents to an output stream
483 800
484 std::ostream& G4Box::StreamInfo(std::ostream& 801 std::ostream& G4Box::StreamInfo(std::ostream& os) const
485 { 802 {
486 G4long oldprc = os.precision(16); << 803 G4int oldprc = os.precision(16);
487 os << "------------------------------------- 804 os << "-----------------------------------------------------------\n"
488 << " *** Dump for solid - " << GetName 805 << " *** Dump for solid - " << GetName() << " ***\n"
489 << " ================================= 806 << " ===================================================\n"
490 << "Solid type: G4Box\n" << 807 << " Solid type: G4Box\n"
491 << "Parameters: \n" << 808 << " Parameters: \n"
492 << " half length X: " << fDx/mm << " mm << 809 << " half length X: " << fDx/mm << " mm \n"
493 << " half length Y: " << fDy/mm << " mm << 810 << " half length Y: " << fDy/mm << " mm \n"
494 << " half length Z: " << fDz/mm << " mm << 811 << " half length Z: " << fDz/mm << " mm \n"
495 << "------------------------------------- 812 << "-----------------------------------------------------------\n";
496 os.precision(oldprc); 813 os.precision(oldprc);
>> 814
497 return os; 815 return os;
498 } 816 }
499 817
500 ////////////////////////////////////////////// << 818 /////////////////////////////////////////////////////////////////////////////////////
501 // 819 //
502 // Return a point randomly and uniformly selec << 820 // GetPointOnSurface
>> 821 //
>> 822 // Return a point (G4ThreeVector) randomly and uniformly selected
>> 823 // on the solid surface
503 824
504 G4ThreeVector G4Box::GetPointOnSurface() const 825 G4ThreeVector G4Box::GetPointOnSurface() const
505 { 826 {
506 G4double sxy = fDx*fDy, sxz = fDx*fDz, syz = << 827 G4double px, py, pz, select, sumS;
507 G4double select = (sxy + sxz + syz)*G4QuickR << 828 G4double Sxy = fDx*fDy, Sxz = fDx*fDz, Syz = fDy*fDz;
508 G4double u = 2.*G4QuickRand() - 1.; << 829
509 G4double v = 2.*G4QuickRand() - 1.; << 830 sumS = Sxy + Sxz + Syz;
510 << 831 select = sumS*G4UniformRand();
511 if (select < sxy) << 832
512 return { u*fDx, v*fDy, ((select < 0.5*sxy) << 833 if( select < Sxy )
513 else if (select < sxy + sxz) << 834 {
514 return { u*fDx, ((select < sxy + 0.5*sxz) << 835 px = -fDx +2*fDx*G4UniformRand();
515 else << 836 py = -fDy +2*fDy*G4UniformRand();
516 return { ((select < sxy + sxz + 0.5*syz) ? << 837
>> 838 if(G4UniformRand() > 0.5) { pz = fDz; }
>> 839 else { pz = -fDz; }
>> 840 }
>> 841 else if ( ( select - Sxy ) < Sxz )
>> 842 {
>> 843 px = -fDx +2*fDx*G4UniformRand();
>> 844 pz = -fDz +2*fDz*G4UniformRand();
>> 845
>> 846 if(G4UniformRand() > 0.5) { py = fDy; }
>> 847 else { py = -fDy; }
>> 848 }
>> 849 else
>> 850 {
>> 851 py = -fDy +2*fDy*G4UniformRand();
>> 852 pz = -fDz +2*fDz*G4UniformRand();
>> 853
>> 854 if(G4UniformRand() > 0.5) { px = fDx; }
>> 855 else { px = -fDx; }
>> 856 }
>> 857 return G4ThreeVector(px,py,pz);
517 } 858 }
518 859
519 ////////////////////////////////////////////// 860 //////////////////////////////////////////////////////////////////////////
520 // 861 //
521 // Make a clone of the object 862 // Make a clone of the object
522 // 863 //
523 G4VSolid* G4Box::Clone() const 864 G4VSolid* G4Box::Clone() const
524 { 865 {
525 return new G4Box(*this); 866 return new G4Box(*this);
526 } 867 }
527 868
528 ////////////////////////////////////////////// 869 //////////////////////////////////////////////////////////////////////////
529 // 870 //
530 // Methods for visualisation 871 // Methods for visualisation
531 872
532 void G4Box::DescribeYourselfTo (G4VGraphicsSce << 873 void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const
533 { 874 {
534 scene.AddSolid (*this); 875 scene.AddSolid (*this);
535 } 876 }
536 877
537 G4VisExtent G4Box::GetExtent() const << 878 G4VisExtent G4Box::GetExtent() const
538 { 879 {
539 return { -fDx, fDx, -fDy, fDy, -fDz, fDz }; << 880 return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz);
540 } 881 }
541 882
542 G4Polyhedron* G4Box::CreatePolyhedron () const << 883 G4Polyhedron* G4Box::CreatePolyhedron () const
543 { 884 {
544 return new G4PolyhedronBox (fDx, fDy, fDz); 885 return new G4PolyhedronBox (fDx, fDy, fDz);
545 } 886 }
546 #endif 887 #endif
547 888