| 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 69788 2013-05-15 12:06:57Z 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
31 // ------------------------------------------- 40 // --------------------------------------------------------------------
32 41
33 #include "G4Box.hh" 42 #include "G4Box.hh"
34 43
35 #if !defined(G4GEOM_USE_UBOX) <<
36 <<
37 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
38 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 46 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh" << 47 #include "Randomize.hh"
41 #include "G4QuickRand.hh" <<
42 48
43 #include "G4VPVParameterisation.hh" 49 #include "G4VPVParameterisation.hh"
44 50
45 #include "G4VGraphicsScene.hh" 51 #include "G4VGraphicsScene.hh"
>> 52 #include "G4Polyhedron.hh"
>> 53 #include "G4NURBS.hh"
>> 54 #include "G4NURBSbox.hh"
46 #include "G4VisExtent.hh" 55 #include "G4VisExtent.hh"
47 56
48 ////////////////////////////////////////////// 57 ////////////////////////////////////////////////////////////////////////
49 // 58 //
50 // Constructor - check & set half widths 59 // Constructor - check & set half widths
51 60
52 G4Box::G4Box(const G4String& pName, 61 G4Box::G4Box(const G4String& pName,
53 G4double pX, 62 G4double pX,
54 G4double pY, 63 G4double pY,
55 G4double pZ) 64 G4double pZ)
56 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(p 65 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
57 { 66 {
58 delta = 0.5*kCarTolerance; << 67 if ( (pX < 2*kCarTolerance)
59 if (pX < 2*kCarTolerance || << 68 && (pY < 2*kCarTolerance)
60 pY < 2*kCarTolerance || << 69 && (pZ < 2*kCarTolerance) ) // limit to thickness of surfaces
61 pZ < 2*kCarTolerance) // limit to thick <<
62 { 70 {
63 std::ostringstream message; 71 std::ostringstream message;
64 message << "Dimensions too small for Solid 72 message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
65 << " hX, hY, hZ = " << pX << " 73 << " hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
66 G4Exception("G4Box::G4Box()", "GeomSolids0 74 G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
67 } 75 }
68 } 76 }
69 77
70 ////////////////////////////////////////////// 78 //////////////////////////////////////////////////////////////////////////
71 // 79 //
72 // Fake default constructor - sets only member 80 // Fake default constructor - sets only member data and allocates memory
73 // for usage restri << 81 // for usage restricted to object persistency.
74 82
75 G4Box::G4Box( __void__& a ) 83 G4Box::G4Box( __void__& a )
76 : G4CSGSolid(a) << 84 : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.)
77 { 85 {
78 } 86 }
79 87
80 ////////////////////////////////////////////// 88 //////////////////////////////////////////////////////////////////////////
81 // 89 //
82 // Destructor 90 // Destructor
83 91
84 G4Box::~G4Box() = default; << 92 G4Box::~G4Box()
>> 93 {
>> 94 }
85 95
86 ////////////////////////////////////////////// 96 //////////////////////////////////////////////////////////////////////////
87 // 97 //
88 // Copy constructor 98 // Copy constructor
89 99
90 G4Box::G4Box(const G4Box&) = default; << 100 G4Box::G4Box(const G4Box& rhs)
>> 101 : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz)
>> 102 {
>> 103 }
91 104
92 ////////////////////////////////////////////// 105 //////////////////////////////////////////////////////////////////////////
93 // 106 //
94 // Assignment operator 107 // Assignment operator
95 108
96 G4Box& G4Box::operator = (const G4Box& rhs) << 109 G4Box& G4Box::operator = (const G4Box& rhs)
97 { 110 {
98 // Check assignment to self 111 // Check assignment to self
99 // 112 //
100 if (this == &rhs) { return *this; } 113 if (this == &rhs) { return *this; }
101 114
102 // Copy base class data 115 // Copy base class data
103 // 116 //
104 G4CSGSolid::operator=(rhs); 117 G4CSGSolid::operator=(rhs);
105 118
106 // Copy data 119 // Copy data
107 // 120 //
108 fDx = rhs.fDx; 121 fDx = rhs.fDx;
109 fDy = rhs.fDy; 122 fDy = rhs.fDy;
110 fDz = rhs.fDz; 123 fDz = rhs.fDz;
111 delta = rhs.delta; <<
112 124
113 return *this; 125 return *this;
114 } 126 }
115 127
116 ////////////////////////////////////////////// << 128 //////////////////////////////////////////////////////////////////////////////
117 // <<
118 // Set X dimension <<
119 129
120 void G4Box::SetXHalfLength(G4double dx) 130 void G4Box::SetXHalfLength(G4double dx)
121 { 131 {
122 if(dx > 2*kCarTolerance) // limit to thickn 132 if(dx > 2*kCarTolerance) // limit to thickness of surfaces
123 { 133 {
124 fDx = dx; 134 fDx = dx;
125 } 135 }
126 else 136 else
127 { 137 {
128 std::ostringstream message; 138 std::ostringstream message;
129 message << "Dimension X too small for soli 139 message << "Dimension X too small for solid: " << GetName() << "!"
130 << G4endl 140 << G4endl
131 << " hX = " << dx; 141 << " hX = " << dx;
132 G4Exception("G4Box::SetXHalfLength()", "Ge 142 G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
133 FatalException, message); 143 FatalException, message);
134 } 144 }
135 fCubicVolume = 0.; << 145 fCubicVolume= 0.;
136 fSurfaceArea = 0.; << 146 fSurfaceArea= 0.;
137 fRebuildPolyhedron = true; << 147 fpPolyhedron = 0;
138 } << 148 }
139 <<
140 ////////////////////////////////////////////// <<
141 // <<
142 // Set Y dimension <<
143 149
144 void G4Box::SetYHalfLength(G4double dy) << 150 void G4Box::SetYHalfLength(G4double dy)
145 { 151 {
146 if(dy > 2*kCarTolerance) // limit to thickn 152 if(dy > 2*kCarTolerance) // limit to thickness of surfaces
147 { 153 {
148 fDy = dy; 154 fDy = dy;
149 } 155 }
150 else 156 else
151 { 157 {
152 std::ostringstream message; 158 std::ostringstream message;
153 message << "Dimension Y too small for soli << 159 message << "Dimension Y too small for solid: " << GetName() << "!"
>> 160 << G4endl
154 << " hY = " << dy; 161 << " hY = " << dy;
155 G4Exception("G4Box::SetYHalfLength()", "Ge 162 G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
156 FatalException, message); 163 FatalException, message);
157 } 164 }
158 fCubicVolume = 0.; << 165 fCubicVolume= 0.;
159 fSurfaceArea = 0.; << 166 fSurfaceArea= 0.;
160 fRebuildPolyhedron = true; << 167 fpPolyhedron = 0;
161 } << 168 }
162 169
163 ////////////////////////////////////////////// << 170 void G4Box::SetZHalfLength(G4double dz)
164 // <<
165 // Set Z dimension <<
166 <<
167 void G4Box::SetZHalfLength(G4double dz) <<
168 { 171 {
169 if(dz > 2*kCarTolerance) // limit to thickn 172 if(dz > 2*kCarTolerance) // limit to thickness of surfaces
170 { 173 {
171 fDz = dz; 174 fDz = dz;
172 } 175 }
173 else 176 else
174 { 177 {
175 std::ostringstream message; 178 std::ostringstream message;
176 message << "Dimension Z too small for soli << 179 message << "Dimension Z too small for solid: " << GetName() << "!"
>> 180 << G4endl
177 << " hZ = " << dz; 181 << " hZ = " << dz;
178 G4Exception("G4Box::SetZHalfLength()", "Ge 182 G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
179 FatalException, message); 183 FatalException, message);
180 } 184 }
181 fCubicVolume = 0.; << 185 fCubicVolume= 0.;
182 fSurfaceArea = 0.; << 186 fSurfaceArea= 0.;
183 fRebuildPolyhedron = true; << 187 fpPolyhedron = 0;
184 } << 188 }
185 189
186 ////////////////////////////////////////////// << 190 ////////////////////////////////////////////////////////////////////////
187 // 191 //
188 // Dispatch to parameterisation for replicatio 192 // Dispatch to parameterisation for replication mechanism dimension
189 // computation & modification. 193 // computation & modification.
190 194
191 void G4Box::ComputeDimensions(G4VPVParameteris 195 void G4Box::ComputeDimensions(G4VPVParameterisation* p,
192 const G4int n, 196 const G4int n,
193 const G4VPhysica 197 const G4VPhysicalVolume* pRep)
194 { 198 {
195 p->ComputeDimensions(*this,n,pRep); 199 p->ComputeDimensions(*this,n,pRep);
196 } 200 }
197 201
198 ////////////////////////////////////////////// 202 //////////////////////////////////////////////////////////////////////////
199 // 203 //
200 // Get bounding box <<
201 <<
202 void G4Box::BoundingLimits(G4ThreeVector& pMin <<
203 { <<
204 pMin.set(-fDx,-fDy,-fDz); <<
205 pMax.set( fDx, fDy, fDz); <<
206 <<
207 // Check correctness of the bounding box <<
208 // <<
209 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
210 { <<
211 std::ostringstream message; <<
212 message << "Bad bounding box (min >= max) <<
213 << GetName() << " !" <<
214 << "\npMin = " << pMin <<
215 << "\npMax = " << pMax; <<
216 G4Exception("G4Box::BoundingLimits()", "Ge <<
217 DumpInfo(); <<
218 } <<
219 } <<
220 <<
221 ////////////////////////////////////////////// <<
222 // <<
223 // Calculate extent under transform and specif 204 // Calculate extent under transform and specified limit
224 205
225 G4bool G4Box::CalculateExtent(const EAxis pAxi 206 G4bool G4Box::CalculateExtent(const EAxis pAxis,
226 const G4VoxelLim 207 const G4VoxelLimits& pVoxelLimit,
227 const G4AffineTr 208 const G4AffineTransform& pTransform,
228 G4double& 209 G4double& pMin, G4double& pMax) const
229 { 210 {
230 G4ThreeVector bmin, bmax; << 211 if (!pTransform.IsRotated())
>> 212 {
>> 213 // Special case handling for unrotated boxes
>> 214 // Compute x/y/z mins and maxs respecting limits, with early returns
>> 215 // if outside limits. Then switch() on pAxis
>> 216
>> 217 G4double xoffset,xMin,xMax;
>> 218 G4double yoffset,yMin,yMax;
>> 219 G4double zoffset,zMin,zMax;
>> 220
>> 221 xoffset = pTransform.NetTranslation().x() ;
>> 222 xMin = xoffset - fDx ;
>> 223 xMax = xoffset + fDx ;
231 224
232 // Get bounding box << 225 if (pVoxelLimit.IsXLimited())
233 BoundingLimits(bmin,bmax); << 226 {
>> 227 if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) ||
>> 228 (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; }
>> 229 else
>> 230 {
>> 231 xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
>> 232 xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
>> 233 }
>> 234 }
>> 235 yoffset = pTransform.NetTranslation().y() ;
>> 236 yMin = yoffset - fDy ;
>> 237 yMax = yoffset + fDy ;
234 238
235 // Find extent << 239 if (pVoxelLimit.IsYLimited())
236 G4BoundingEnvelope bbox(bmin,bmax); << 240 {
237 return bbox.CalculateExtent(pAxis,pVoxelLimi << 241 if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) ||
238 } << 242 (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; }
>> 243 else
>> 244 {
>> 245 yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
>> 246 yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
>> 247 }
>> 248 }
>> 249 zoffset = pTransform.NetTranslation().z() ;
>> 250 zMin = zoffset - fDz ;
>> 251 zMax = zoffset + fDz ;
239 252
240 ////////////////////////////////////////////// << 253 if (pVoxelLimit.IsZLimited())
>> 254 {
>> 255 if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) ||
>> 256 (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; }
>> 257 else
>> 258 {
>> 259 zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
>> 260 zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
>> 261 }
>> 262 }
>> 263 switch (pAxis)
>> 264 {
>> 265 case kXAxis:
>> 266 pMin = xMin ;
>> 267 pMax = xMax ;
>> 268 break ;
>> 269 case kYAxis:
>> 270 pMin=yMin;
>> 271 pMax=yMax;
>> 272 break;
>> 273 case kZAxis:
>> 274 pMin=zMin;
>> 275 pMax=zMax;
>> 276 break;
>> 277 default:
>> 278 break;
>> 279 }
>> 280 pMin -= kCarTolerance ;
>> 281 pMax += kCarTolerance ;
>> 282
>> 283 return true;
>> 284 }
>> 285 else // General rotated case - create and clip mesh to boundaries
>> 286 {
>> 287 G4bool existsAfterClip = false ;
>> 288 G4ThreeVectorList* vertices ;
>> 289
>> 290 pMin = +kInfinity ;
>> 291 pMax = -kInfinity ;
>> 292
>> 293 // Calculate rotated vertex coordinates
>> 294
>> 295 vertices = CreateRotatedVertices(pTransform) ;
>> 296 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
>> 297 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
>> 298 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
>> 299
>> 300 if (pVoxelLimit.IsLimited(pAxis) == false)
>> 301 {
>> 302 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
>> 303 {
>> 304 existsAfterClip = true ;
>> 305
>> 306 // Add 2*tolerance to avoid precision troubles
>> 307
>> 308 pMin -= kCarTolerance;
>> 309 pMax += kCarTolerance;
>> 310 }
>> 311 }
>> 312 else
>> 313 {
>> 314 G4ThreeVector clipCentre(
>> 315 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 316 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 317 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
>> 318
>> 319 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
>> 320 {
>> 321 existsAfterClip = true ;
>> 322
>> 323
>> 324 // Check to see if endpoints are in the solid
>> 325
>> 326 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
>> 327
>> 328 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
>> 329 {
>> 330 pMin = pVoxelLimit.GetMinExtent(pAxis);
>> 331 }
>> 332 else
>> 333 {
>> 334 pMin -= kCarTolerance;
>> 335 }
>> 336 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
>> 337
>> 338 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
>> 339 {
>> 340 pMax = pVoxelLimit.GetMaxExtent(pAxis);
>> 341 }
>> 342 else
>> 343 {
>> 344 pMax += kCarTolerance;
>> 345 }
>> 346 }
>> 347
>> 348 // Check for case where completely enveloping clipping volume
>> 349 // If point inside then we are confident that the solid completely
>> 350 // envelopes the clipping volume. Hence set min/max extents according
>> 351 // to clipping volume extents along the specified axis.
>> 352
>> 353 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
>> 354 != kOutside)
>> 355 {
>> 356 existsAfterClip = true ;
>> 357 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
>> 358 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
>> 359 }
>> 360 }
>> 361 delete vertices;
>> 362 return existsAfterClip;
>> 363 }
>> 364 }
>> 365
>> 366 /////////////////////////////////////////////////////////////////////////
241 // 367 //
242 // Return whether point inside/outside/on surf 368 // Return whether point inside/outside/on surface, using tolerance
243 369
244 EInside G4Box::Inside(const G4ThreeVector& p) 370 EInside G4Box::Inside(const G4ThreeVector& p) const
245 { 371 {
246 G4double dist = std::max(std::max( << 372 static const G4double delta=0.5*kCarTolerance;
247 std::abs(p.x())-fDx, << 373 EInside in = kOutside ;
248 std::abs(p.y())-fDy), << 374 G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
249 std::abs(p.z())-fDz); << 375
250 return (dist > delta) ? kOutside : << 376 if ( q.x() <= (fDx - delta) )
251 ((dist > -delta) ? kSurface : kInside); << 377 {
>> 378 if (q.y() <= (fDy - delta) )
>> 379 {
>> 380 if ( q.z() <= (fDz - delta) ) { in = kInside ; }
>> 381 else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
>> 382 }
>> 383 else if ( q.y() <= (fDy + delta) )
>> 384 {
>> 385 if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
>> 386 }
>> 387 }
>> 388 else if ( q.x() <= (fDx + delta) )
>> 389 {
>> 390 if ( q.y() <= (fDy + delta) )
>> 391 {
>> 392 if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
>> 393 }
>> 394 }
>> 395 return in ;
252 } 396 }
253 397
254 ////////////////////////////////////////////// << 398 ///////////////////////////////////////////////////////////////////////
255 // 399 //
256 // Detect the side(s) and return corresponding << 400 // Calculate side nearest to p, and return normal
>> 401 // If two sides are equidistant, normal of first side (x/y/z)
>> 402 // encountered returned
257 403
258 G4ThreeVector G4Box::SurfaceNormal( const G4Th 404 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
259 { 405 {
260 G4ThreeVector norm(0,0,0); << 406 G4double distx, disty, distz ;
261 G4double px = p.x(); << 407 G4ThreeVector norm(0.,0.,0.);
262 if (std::abs(std::abs(px) - fDx) <= delta) n << 408
263 G4double py = p.y(); << 409 // Calculate distances as if in 1st octant
264 if (std::abs(std::abs(py) - fDy) <= delta) n << 410
265 G4double pz = p.z(); << 411 distx = std::fabs(std::fabs(p.x()) - fDx) ;
266 if (std::abs(std::abs(pz) - fDz) <= delta) n << 412 disty = std::fabs(std::fabs(p.y()) - fDy) ;
267 << 413 distz = std::fabs(std::fabs(p.z()) - fDz) ;
268 G4double nside = norm.mag2(); // number of s << 414
269 if (nside == 1) << 415 // New code for particle on surface including edges and corners with specific
270 return norm; << 416 // normals
271 else if (nside > 1) << 417
272 return norm.unit(); // edge or corner << 418 static const G4double delta = 0.5*kCarTolerance;
>> 419 const G4ThreeVector nX = G4ThreeVector( 1.0, 0,0 );
>> 420 const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0 );
>> 421 const G4ThreeVector nY = G4ThreeVector( 0, 1.0,0 );
>> 422 const G4ThreeVector nmY = G4ThreeVector( 0,-1.0,0 );
>> 423 const G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
>> 424 const G4ThreeVector nmZ = G4ThreeVector( 0, 0,- 1.0);
>> 425
>> 426 G4ThreeVector normX(0.,0.,0.), normY(0.,0.,0.), normZ(0.,0.,0.);
>> 427 G4ThreeVector sumnorm(0., 0., 0.);
>> 428 G4int noSurfaces=0;
>> 429
>> 430 if (distx <= delta) // on X/mX surface and around
>> 431 {
>> 432 noSurfaces ++;
>> 433 if ( p.x() >= 0. ) { normX= nX ; } // on +X surface : (1,0,0)
>> 434 else { normX= nmX; } // (-1,0,0)
>> 435 sumnorm= normX;
>> 436 }
>> 437
>> 438 if (disty <= delta) // on one of the +Y or -Y surfaces
>> 439 {
>> 440 noSurfaces ++;
>> 441 if ( p.y() >= 0. ) { normY= nY; } // on +Y surface
>> 442 else { normY= nmY; }
>> 443 sumnorm += normY;
>> 444 }
>> 445
>> 446 if (distz <= delta) // on one of the +Z or -Z surfaces
>> 447 {
>> 448 noSurfaces ++;
>> 449 if ( p.z() >= 0. ) { normZ= nZ; } // on +Z surface
>> 450 else { normZ= nmZ; }
>> 451 sumnorm += normZ;
>> 452 }
>> 453
>> 454 static const G4double invSqrt2 = 1.0 / std::sqrt(2.0);
>> 455 static const G4double invSqrt3 = 1.0 / std::sqrt(3.0);
>> 456
>> 457 if( noSurfaces > 0 )
>> 458 {
>> 459 if( noSurfaces == 1 )
>> 460 {
>> 461 norm= sumnorm;
>> 462 }
>> 463 else
>> 464 {
>> 465 // norm = sumnorm . unit();
>> 466 if( noSurfaces == 2 )
>> 467 {
>> 468 // 2 surfaces -> on edge
>> 469 norm = invSqrt2 * sumnorm;
>> 470 }
>> 471 else
>> 472 {
>> 473 // 3 surfaces (on corner)
>> 474 norm = invSqrt3 * sumnorm;
>> 475 }
>> 476 }
>> 477 }
273 else 478 else
274 { 479 {
275 // Point is not on the surface <<
276 // <<
277 #ifdef G4CSGDEBUG 480 #ifdef G4CSGDEBUG
278 std::ostringstream message; << 481 G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning,
279 G4int oldprc = message.precision(16); << 482 "Point p is not on surface !?" );
280 message << "Point p is not on surface (!?) << 483 #endif
281 << GetName() << G4endl; << 484 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 } 485 }
>> 486
>> 487 return norm;
293 } 488 }
294 489
295 ////////////////////////////////////////////// 490 //////////////////////////////////////////////////////////////////////////
296 // 491 //
297 // Algorithm for SurfaceNormal() following the 492 // Algorithm for SurfaceNormal() following the original specification
298 // for points not on the surface 493 // for points not on the surface
299 494
300 G4ThreeVector G4Box::ApproxSurfaceNormal(const << 495 G4ThreeVector G4Box::ApproxSurfaceNormal( const G4ThreeVector& p ) const
301 { 496 {
302 G4double distx = std::abs(p.x()) - fDx; << 497 G4double distx, disty, distz ;
303 G4double disty = std::abs(p.y()) - fDy; << 498 G4ThreeVector norm(0.,0.,0.);
304 G4double distz = std::abs(p.z()) - fDz; << 499
305 << 500 // Calculate distances as if in 1st octant
306 if (distx >= disty && distx >= distz) << 501
307 return {std::copysign(1.,p.x()), 0., 0.}; << 502 distx = std::fabs(std::fabs(p.x()) - fDx) ;
308 if (disty >= distx && disty >= distz) << 503 disty = std::fabs(std::fabs(p.y()) - fDy) ;
309 return {0., std::copysign(1.,p.y()), 0.}; << 504 distz = std::fabs(std::fabs(p.z()) - fDz) ;
>> 505
>> 506 if ( distx <= disty )
>> 507 {
>> 508 if ( distx <= distz ) // Closest to X
>> 509 {
>> 510 if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
>> 511 else { norm = G4ThreeVector( 1.0,0,0) ; }
>> 512 }
>> 513 else // Closest to Z
>> 514 {
>> 515 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
>> 516 else { norm = G4ThreeVector(0,0, 1.0) ; }
>> 517 }
>> 518 }
310 else 519 else
311 return {0., 0., std::copysign(1.,p.z())}; << 520 {
>> 521 if ( disty <= distz ) // Closest to Y
>> 522 {
>> 523 if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
>> 524 else { norm = G4ThreeVector(0, 1.0,0) ; }
>> 525 }
>> 526 else // Closest to Z
>> 527 {
>> 528 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
>> 529 else { norm = G4ThreeVector(0,0, 1.0) ; }
>> 530 }
>> 531 }
>> 532 return norm;
312 } 533 }
313 534
314 ////////////////////////////////////////////// << 535 ///////////////////////////////////////////////////////////////////////////
315 // 536 //
316 // Calculate distance to box from an outside p 537 // Calculate distance to box from an outside point
317 // - return kInfinity if no intersection << 538 // - return kInfinity if no intersection.
>> 539 //
>> 540 // ALGORITHM:
>> 541 //
>> 542 // Check that if point lies outside x/y/z extent of box, travel is towards
>> 543 // the box (ie. there is a possibility of an intersection)
>> 544 //
>> 545 // Calculate pairs of minimum and maximum distances for x/y/z travel for
>> 546 // intersection with the box's x/y/z extent.
>> 547 // If there is a valid intersection, it is given by the maximum min distance
>> 548 // (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance
>> 549 // (ie. distance after which 1+ of x/y/z intersections not satisfied)
318 // 550 //
>> 551 // NOTE:
>> 552 //
>> 553 // `Inside' safe - meaningful answers given if point is inside the exact
>> 554 // shape.
319 555
320 G4double G4Box::DistanceToIn(const G4ThreeVect 556 G4double G4Box::DistanceToIn(const G4ThreeVector& p,
321 const G4ThreeVect 557 const G4ThreeVector& v) const
322 { 558 {
323 // Check if point is on the surface and trav << 559 G4double safx, safy, safz ;
324 // << 560 G4double smin=0.0, sminy, sminz ; // , sminx ;
325 if ((std::abs(p.x()) - fDx) >= -delta && p.x << 561 G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ; // they always > 0
326 if ((std::abs(p.y()) - fDy) >= -delta && p.y << 562 G4double stmp ;
327 if ((std::abs(p.z()) - fDz) >= -delta && p.z << 563 G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
328 << 564
329 // Find intersection << 565 static const G4double delta = 0.5*kCarTolerance;
330 // << 566
331 G4double invx = (v.x() == 0) ? DBL_MAX : -1. << 567 safx = std::fabs(p.x()) - fDx ; // minimum distance to x surface of shape
332 G4double dx = std::copysign(fDx,invx); << 568 safy = std::fabs(p.y()) - fDy ;
333 G4double txmin = (p.x() - dx)*invx; << 569 safz = std::fabs(p.z()) - fDz ;
334 G4double txmax = (p.x() + dx)*invx; << 570
335 << 571 // Will we intersect?
336 G4double invy = (v.y() == 0) ? DBL_MAX : -1. << 572 // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent.
337 G4double dy = std::copysign(fDy,invy); << 573 // If both p.x/y/z and v.x/y/z repectively are both positive/negative,
338 G4double tymin = std::max(txmin,(p.y() - dy) << 574 // travel is in a direction away from the shape.
339 G4double tymax = std::min(txmax,(p.y() + dy) << 575
340 << 576 if ( ((p.x()*v.x() >= 0.0) && (safx > -delta))
341 G4double invz = (v.z() == 0) ? DBL_MAX : -1. << 577 || ((p.y()*v.y() >= 0.0) && (safy > -delta))
342 G4double dz = std::copysign(fDz,invz); << 578 || ((p.z()*v.z() >= 0.0) && (safz > -delta)) )
343 G4double tmin = std::max(tymin,(p.z() - dz)* << 579 {
344 G4double tmax = std::min(tymax,(p.z() + dz)* << 580 return kInfinity ; // travel away or parallel within tolerance
>> 581 }
>> 582
>> 583 // Compute min / max distances for x/y/z travel:
>> 584 // X Planes
>> 585
>> 586 if ( v.x() ) // != 0
>> 587 {
>> 588 stmp = 1.0/std::fabs(v.x()) ;
345 589
346 if (tmax <= tmin + delta) return kInfinity; << 590 if (safx >= 0.0)
347 return (tmin < delta) ? 0. : tmin; << 591 {
>> 592 smin = safx*stmp ;
>> 593 smax = (fDx+std::fabs(p.x()))*stmp ;
>> 594 }
>> 595 else
>> 596 {
>> 597 if (v.x() < 0) { sOut = (fDx + p.x())*stmp ; }
>> 598 else { sOut = (fDx - p.x())*stmp ; }
>> 599 }
>> 600 }
>> 601
>> 602 // Y Planes
>> 603
>> 604 if ( v.y() ) // != 0
>> 605 {
>> 606 stmp = 1.0/std::fabs(v.y()) ;
>> 607
>> 608 if (safy >= 0.0)
>> 609 {
>> 610 sminy = safy*stmp ;
>> 611 smaxy = (fDy+std::fabs(p.y()))*stmp ;
>> 612
>> 613 if (sminy > smin) { smin=sminy ; }
>> 614 if (smaxy < smax) { smax=smaxy ; }
>> 615
>> 616 if (smin >= (smax-delta))
>> 617 {
>> 618 return kInfinity ; // touch XY corner
>> 619 }
>> 620 }
>> 621 else
>> 622 {
>> 623 if (v.y() < 0) { sOuty = (fDy + p.y())*stmp ; }
>> 624 else { sOuty = (fDy - p.y())*stmp ; }
>> 625 if( sOuty < sOut ) { sOut = sOuty ; }
>> 626 }
>> 627 }
>> 628
>> 629 // Z planes
>> 630
>> 631 if ( v.z() ) // != 0
>> 632 {
>> 633 stmp = 1.0/std::fabs(v.z()) ;
>> 634
>> 635 if ( safz >= 0.0 )
>> 636 {
>> 637 sminz = safz*stmp ;
>> 638 smaxz = (fDz+std::fabs(p.z()))*stmp ;
>> 639
>> 640 if (sminz > smin) { smin = sminz ; }
>> 641 if (smaxz < smax) { smax = smaxz ; }
>> 642
>> 643 if (smin >= (smax-delta))
>> 644 {
>> 645 return kInfinity ; // touch ZX or ZY corners
>> 646 }
>> 647 }
>> 648 else
>> 649 {
>> 650 if (v.z() < 0) { sOutz = (fDz + p.z())*stmp ; }
>> 651 else { sOutz = (fDz - p.z())*stmp ; }
>> 652 if( sOutz < sOut ) { sOut = sOutz ; }
>> 653 }
>> 654 }
>> 655
>> 656 if (sOut <= (smin + delta)) // travel over edge
>> 657 {
>> 658 return kInfinity ;
>> 659 }
>> 660 if (smin < delta) { smin = 0.0 ; }
>> 661
>> 662 return smin ;
348 } 663 }
349 664
350 ////////////////////////////////////////////// 665 //////////////////////////////////////////////////////////////////////////
351 // << 666 //
352 // Appoximate distance to box. 667 // Appoximate distance to box.
353 // Returns largest perpendicular distance to t 668 // Returns largest perpendicular distance to the closest x/y/z sides of
354 // the box, which is the most fast estimation 669 // the box, which is the most fast estimation of the shortest distance to box
355 // - If inside return 0 670 // - If inside return 0
356 671
357 G4double G4Box::DistanceToIn(const G4ThreeVect 672 G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
358 { 673 {
359 G4double dist = std::max(std::max( << 674 G4double safex, safey, safez, safe = 0.0 ;
360 std::abs(p.x())-fDx, << 675
361 std::abs(p.y())-fDy), << 676 safex = std::fabs(p.x()) - fDx ;
362 std::abs(p.z())-fDz); << 677 safey = std::fabs(p.y()) - fDy ;
363 return (dist > 0) ? dist : 0.; << 678 safez = std::fabs(p.z()) - fDz ;
>> 679
>> 680 if (safex > safe) { safe = safex ; }
>> 681 if (safey > safe) { safe = safey ; }
>> 682 if (safez > safe) { safe = safez ; }
>> 683
>> 684 return safe ;
364 } 685 }
365 686
366 ////////////////////////////////////////////// << 687 /////////////////////////////////////////////////////////////////////////
367 // 688 //
368 // Calculate distance to surface of the box fr << 689 // Calculate distance to surface of box from inside
369 // find normal at exit point, if required << 690 // by calculating distances to box's x/y/z planes.
370 // - when leaving the surface, return 0 << 691 // Smallest distance is exact distance to exiting.
>> 692 // - Eliminate one side of each pair by considering direction of v
>> 693 // - when leaving a surface & v.close, return 0
371 694
372 G4double G4Box::DistanceToOut( const G4ThreeVe << 695 G4double G4Box::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
373 const G4ThreeVe <<
374 const G4bool ca 696 const G4bool calcNorm,
375 G4bool* validNo << 697 G4bool *validNorm,G4ThreeVector *n) const
376 { 698 {
377 // Check if point is on the surface and trav << 699 ESide side = kUndefined ;
378 // << 700 G4double pdist,stmp,snxt=kInfinity;
379 if ((std::abs(p.x()) - fDx) >= -delta && p.x << 701
>> 702 static const G4double delta = 0.5*kCarTolerance;
>> 703
>> 704 if (calcNorm) { *validNorm = true ; } // All normals are valid
>> 705
>> 706 if (v.x() > 0) // X planes
380 { 707 {
381 if (calcNorm) << 708 pdist = fDx - p.x() ;
>> 709
>> 710 if (pdist > delta)
382 { 711 {
383 *validNorm = true; << 712 snxt = pdist/v.x() ;
384 n->set((p.x() < 0) ? -1. : 1., 0., 0.); << 713 side = kPX ;
>> 714 }
>> 715 else
>> 716 {
>> 717 if (calcNorm) { *n = G4ThreeVector(1,0,0) ; }
>> 718 return snxt = 0 ;
385 } 719 }
386 return 0.; <<
387 } 720 }
388 if ((std::abs(p.y()) - fDy) >= -delta && p.y << 721 else if (v.x() < 0)
389 { 722 {
390 if (calcNorm) << 723 pdist = fDx + p.x() ;
>> 724
>> 725 if (pdist > delta)
391 { 726 {
392 *validNorm = true; << 727 snxt = -pdist/v.x() ;
393 n->set(0., (p.y() < 0) ? -1. : 1., 0.); << 728 side = kMX ;
>> 729 }
>> 730 else
>> 731 {
>> 732 if (calcNorm) { *n = G4ThreeVector(-1,0,0) ; }
>> 733 return snxt = 0 ;
394 } 734 }
395 return 0.; <<
396 } 735 }
397 if ((std::abs(p.z()) - fDz) >= -delta && p.z << 736
>> 737 if (v.y() > 0) // Y planes
398 { 738 {
399 if (calcNorm) << 739 pdist = fDy-p.y();
>> 740
>> 741 if (pdist > delta)
>> 742 {
>> 743 stmp = pdist/v.y();
>> 744
>> 745 if (stmp < snxt)
>> 746 {
>> 747 snxt = stmp;
>> 748 side = kPY;
>> 749 }
>> 750 }
>> 751 else
400 { 752 {
401 *validNorm = true; << 753 if (calcNorm) { *n = G4ThreeVector(0,1,0) ; }
402 n->set(0., 0., (p.z() < 0) ? -1. : 1.); << 754 return snxt = 0 ;
403 } 755 }
404 return 0.; <<
405 } 756 }
>> 757 else if (v.y() < 0)
>> 758 {
>> 759 pdist = fDy + p.y() ;
406 760
407 // Find intersection << 761 if (pdist > delta)
408 // << 762 {
409 G4double vx = v.x(); << 763 stmp = -pdist/v.y();
410 G4double tx = (vx == 0) ? DBL_MAX : (std::co <<
411 764
412 G4double vy = v.y(); << 765 if ( stmp < snxt )
413 G4double ty = (vy == 0) ? tx : (std::copysig << 766 {
414 G4double txy = std::min(tx,ty); << 767 snxt = stmp;
>> 768 side = kMY;
>> 769 }
>> 770 }
>> 771 else
>> 772 {
>> 773 if (calcNorm) { *n = G4ThreeVector(0,-1,0) ; }
>> 774 return snxt = 0 ;
>> 775 }
>> 776 }
415 777
416 G4double vz = v.z(); << 778 if (v.z() > 0) // Z planes
417 G4double tz = (vz == 0) ? txy : (std::copysi << 779 {
418 G4double tmax = std::min(txy,tz); << 780 pdist = fDz-p.z();
419 781
420 // Set normal, if required, and return dista << 782 if ( pdist > delta )
421 // << 783 {
422 if (calcNorm) << 784 stmp = pdist/v.z();
>> 785
>> 786 if ( stmp < snxt )
>> 787 {
>> 788 snxt = stmp;
>> 789 side = kPZ;
>> 790 }
>> 791 }
>> 792 else
>> 793 {
>> 794 if (calcNorm) { *n = G4ThreeVector(0,0,1) ; }
>> 795 return snxt = 0 ;
>> 796 }
>> 797 }
>> 798 else if (v.z() < 0)
423 { 799 {
424 *validNorm = true; << 800 pdist = fDz + p.z();
425 if (tmax == tx) n->set((v.x() < 0) ? << 801
426 else if (tmax == ty) n->set(0., (v.y() < 0 << 802 if ( pdist > delta )
427 else n->set(0., 0., (v.z() << 803 {
>> 804 stmp = -pdist/v.z();
>> 805
>> 806 if ( stmp < snxt )
>> 807 {
>> 808 snxt = stmp;
>> 809 side = kMZ;
>> 810 }
>> 811 }
>> 812 else
>> 813 {
>> 814 if (calcNorm) { *n = G4ThreeVector(0,0,-1) ; }
>> 815 return snxt = 0 ;
>> 816 }
428 } 817 }
429 return tmax; << 818
>> 819 if (calcNorm)
>> 820 {
>> 821 switch (side)
>> 822 {
>> 823 case kPX:
>> 824 *n=G4ThreeVector(1,0,0);
>> 825 break;
>> 826 case kMX:
>> 827 *n=G4ThreeVector(-1,0,0);
>> 828 break;
>> 829 case kPY:
>> 830 *n=G4ThreeVector(0,1,0);
>> 831 break;
>> 832 case kMY:
>> 833 *n=G4ThreeVector(0,-1,0);
>> 834 break;
>> 835 case kPZ:
>> 836 *n=G4ThreeVector(0,0,1);
>> 837 break;
>> 838 case kMZ:
>> 839 *n=G4ThreeVector(0,0,-1);
>> 840 break;
>> 841 default:
>> 842 G4cout << G4endl;
>> 843 DumpInfo();
>> 844 std::ostringstream message;
>> 845 G4int oldprc = message.precision(16);
>> 846 message << "Undefined side for valid surface normal to solid."
>> 847 << G4endl
>> 848 << "Position:" << G4endl << G4endl
>> 849 << "p.x() = " << p.x()/mm << " mm" << G4endl
>> 850 << "p.y() = " << p.y()/mm << " mm" << G4endl
>> 851 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
>> 852 << "Direction:" << G4endl << G4endl
>> 853 << "v.x() = " << v.x() << G4endl
>> 854 << "v.y() = " << v.y() << G4endl
>> 855 << "v.z() = " << v.z() << G4endl << G4endl
>> 856 << "Proposed distance :" << G4endl << G4endl
>> 857 << "snxt = " << snxt/mm << " mm" << G4endl;
>> 858 message.precision(oldprc);
>> 859 G4Exception("G4Box::DistanceToOut(p,v,..)", "GeomSolids1002",
>> 860 JustWarning, message);
>> 861 break;
>> 862 }
>> 863 }
>> 864 return snxt;
430 } 865 }
431 866
432 ////////////////////////////////////////////// 867 ////////////////////////////////////////////////////////////////////////////
433 // 868 //
434 // Calculate exact shortest distance to any bo 869 // Calculate exact shortest distance to any boundary from inside
435 // - if outside return 0 << 870 // - If outside return 0
436 871
437 G4double G4Box::DistanceToOut(const G4ThreeVec 872 G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
438 { 873 {
>> 874 G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0;
>> 875
439 #ifdef G4CSGDEBUG 876 #ifdef G4CSGDEBUG
440 if( Inside(p) == kOutside ) 877 if( Inside(p) == kOutside )
441 { 878 {
442 std::ostringstream message; << 879 G4int oldprc = G4cout.precision(16) ;
443 G4int oldprc = message.precision(16); << 880 G4cout << G4endl ;
444 message << "Point p is outside (!?) of sol << 881 DumpInfo();
445 message << "Position:\n"; << 882 G4cout << "Position:" << G4endl << G4endl ;
446 message << " p.x() = " << p.x()/mm << " << 883 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
447 message << " p.y() = " << p.y()/mm << " << 884 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
448 message << " p.z() = " << p.z()/mm << " << 885 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
449 G4cout.precision(oldprc); << 886 G4cout.precision(oldprc) ;
450 G4Exception("G4Box::DistanceToOut(p)", "Ge << 887 G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
451 JustWarning, message ); << 888 JustWarning, "Point p is outside !?" );
452 DumpInfo(); <<
453 } 889 }
454 #endif 890 #endif
455 G4double dist = std::min(std::min( << 891
456 fDx-std::abs(p.x()), << 892 safx1 = fDx - p.x() ;
457 fDy-std::abs(p.y())), << 893 safx2 = fDx + p.x() ;
458 fDz-std::abs(p.z())); << 894 safy1 = fDy - p.y() ;
459 return (dist > 0) ? dist : 0.; << 895 safy2 = fDy + p.y() ;
>> 896 safz1 = fDz - p.z() ;
>> 897 safz2 = fDz + p.z() ;
>> 898
>> 899 // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
>> 900
>> 901 if (safx2 < safx1) { safe = safx2; }
>> 902 else { safe = safx1; }
>> 903 if (safy1 < safe) { safe = safy1; }
>> 904 if (safy2 < safe) { safe = safy2; }
>> 905 if (safz1 < safe) { safe = safz1; }
>> 906 if (safz2 < safe) { safe = safz2; }
>> 907
>> 908 if (safe < 0) { safe = 0 ; }
>> 909 return safe ;
460 } 910 }
461 911
462 ////////////////////////////////////////////// << 912 ////////////////////////////////////////////////////////////////////////
463 // 913 //
464 // GetEntityType << 914 // Create a List containing the transformed vertices
>> 915 // Ordering [0-3] -fDz cross section
>> 916 // [4-7] +fDz cross section such that [0] is below [4],
>> 917 // [1] below [5] etc.
>> 918 // Note:
>> 919 // Caller has deletion resposibility
465 920
466 G4GeometryType G4Box::GetEntityType() const << 921 G4ThreeVectorList*
>> 922 G4Box::CreateRotatedVertices(const G4AffineTransform& pTransform) const
467 { 923 {
468 return {"G4Box"}; << 924 G4ThreeVectorList* vertices = new G4ThreeVectorList();
>> 925
>> 926 if (vertices)
>> 927 {
>> 928 vertices->reserve(8);
>> 929 G4ThreeVector vertex0(-fDx,-fDy,-fDz) ;
>> 930 G4ThreeVector vertex1(fDx,-fDy,-fDz) ;
>> 931 G4ThreeVector vertex2(fDx,fDy,-fDz) ;
>> 932 G4ThreeVector vertex3(-fDx,fDy,-fDz) ;
>> 933 G4ThreeVector vertex4(-fDx,-fDy,fDz) ;
>> 934 G4ThreeVector vertex5(fDx,-fDy,fDz) ;
>> 935 G4ThreeVector vertex6(fDx,fDy,fDz) ;
>> 936 G4ThreeVector vertex7(-fDx,fDy,fDz) ;
>> 937
>> 938 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 939 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 940 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 941 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 942 vertices->push_back(pTransform.TransformPoint(vertex4));
>> 943 vertices->push_back(pTransform.TransformPoint(vertex5));
>> 944 vertices->push_back(pTransform.TransformPoint(vertex6));
>> 945 vertices->push_back(pTransform.TransformPoint(vertex7));
>> 946 }
>> 947 else
>> 948 {
>> 949 DumpInfo();
>> 950 G4Exception("G4Box::CreateRotatedVertices()",
>> 951 "GeomSolids0003", FatalException,
>> 952 "Error in allocation of vertices. Out of memory !");
>> 953 }
>> 954 return vertices;
469 } 955 }
470 956
471 ////////////////////////////////////////////// 957 //////////////////////////////////////////////////////////////////////////
472 // 958 //
473 // IsFaceted << 959 // GetEntityType
474 960
475 G4bool G4Box::IsFaceted() const << 961 G4GeometryType G4Box::GetEntityType() const
476 { 962 {
477 return true; << 963 return G4String("G4Box");
478 } 964 }
479 965
480 ////////////////////////////////////////////// 966 //////////////////////////////////////////////////////////////////////////
481 // 967 //
482 // Stream object contents to an output stream 968 // Stream object contents to an output stream
483 969
484 std::ostream& G4Box::StreamInfo(std::ostream& 970 std::ostream& G4Box::StreamInfo(std::ostream& os) const
485 { 971 {
486 G4long oldprc = os.precision(16); << 972 G4int oldprc = os.precision(16);
487 os << "------------------------------------- 973 os << "-----------------------------------------------------------\n"
488 << " *** Dump for solid - " << GetName 974 << " *** Dump for solid - " << GetName() << " ***\n"
489 << " ================================= 975 << " ===================================================\n"
490 << "Solid type: G4Box\n" << 976 << " Solid type: G4Box\n"
491 << "Parameters: \n" << 977 << " Parameters: \n"
492 << " half length X: " << fDx/mm << " mm << 978 << " half length X: " << fDx/mm << " mm \n"
493 << " half length Y: " << fDy/mm << " mm << 979 << " half length Y: " << fDy/mm << " mm \n"
494 << " half length Z: " << fDz/mm << " mm << 980 << " half length Z: " << fDz/mm << " mm \n"
495 << "------------------------------------- 981 << "-----------------------------------------------------------\n";
496 os.precision(oldprc); 982 os.precision(oldprc);
>> 983
497 return os; 984 return os;
498 } 985 }
499 986
500 ////////////////////////////////////////////// << 987 /////////////////////////////////////////////////////////////////////////////////////
>> 988 //
>> 989 // GetPointOnSurface
501 // 990 //
502 // Return a point randomly and uniformly selec << 991 // Return a point (G4ThreeVector) randomly and uniformly selected
>> 992 // on the solid surface
503 993
504 G4ThreeVector G4Box::GetPointOnSurface() const 994 G4ThreeVector G4Box::GetPointOnSurface() const
505 { 995 {
506 G4double sxy = fDx*fDy, sxz = fDx*fDz, syz = << 996 G4double px, py, pz, select, sumS;
507 G4double select = (sxy + sxz + syz)*G4QuickR << 997 G4double Sxy = fDx*fDy, Sxz = fDx*fDz, Syz = fDy*fDz;
508 G4double u = 2.*G4QuickRand() - 1.; << 998
509 G4double v = 2.*G4QuickRand() - 1.; << 999 sumS = Sxy + Sxz + Syz;
510 << 1000 select = sumS*G4UniformRand();
511 if (select < sxy) << 1001
512 return { u*fDx, v*fDy, ((select < 0.5*sxy) << 1002 if( select < Sxy )
513 else if (select < sxy + sxz) << 1003 {
514 return { u*fDx, ((select < sxy + 0.5*sxz) << 1004 px = -fDx +2*fDx*G4UniformRand();
515 else << 1005 py = -fDy +2*fDy*G4UniformRand();
516 return { ((select < sxy + sxz + 0.5*syz) ? << 1006
>> 1007 if(G4UniformRand() > 0.5) { pz = fDz; }
>> 1008 else { pz = -fDz; }
>> 1009 }
>> 1010 else if ( ( select - Sxy ) < Sxz )
>> 1011 {
>> 1012 px = -fDx +2*fDx*G4UniformRand();
>> 1013 pz = -fDz +2*fDz*G4UniformRand();
>> 1014
>> 1015 if(G4UniformRand() > 0.5) { py = fDy; }
>> 1016 else { py = -fDy; }
>> 1017 }
>> 1018 else
>> 1019 {
>> 1020 py = -fDy +2*fDy*G4UniformRand();
>> 1021 pz = -fDz +2*fDz*G4UniformRand();
>> 1022
>> 1023 if(G4UniformRand() > 0.5) { px = fDx; }
>> 1024 else { px = -fDx; }
>> 1025 }
>> 1026 return G4ThreeVector(px,py,pz);
517 } 1027 }
518 1028
519 ////////////////////////////////////////////// 1029 //////////////////////////////////////////////////////////////////////////
520 // 1030 //
521 // Make a clone of the object 1031 // Make a clone of the object
522 // 1032 //
523 G4VSolid* G4Box::Clone() const 1033 G4VSolid* G4Box::Clone() const
524 { 1034 {
525 return new G4Box(*this); 1035 return new G4Box(*this);
526 } 1036 }
527 1037
528 ////////////////////////////////////////////// 1038 //////////////////////////////////////////////////////////////////////////
529 // 1039 //
530 // Methods for visualisation 1040 // Methods for visualisation
531 1041
532 void G4Box::DescribeYourselfTo (G4VGraphicsSce << 1042 void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const
533 { 1043 {
534 scene.AddSolid (*this); 1044 scene.AddSolid (*this);
535 } 1045 }
536 1046
537 G4VisExtent G4Box::GetExtent() const << 1047 G4VisExtent G4Box::GetExtent() const
538 { 1048 {
539 return { -fDx, fDx, -fDy, fDy, -fDz, fDz }; << 1049 return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz);
540 } 1050 }
541 1051
542 G4Polyhedron* G4Box::CreatePolyhedron () const << 1052 G4Polyhedron* G4Box::CreatePolyhedron () const
543 { 1053 {
544 return new G4PolyhedronBox (fDx, fDy, fDz); 1054 return new G4PolyhedronBox (fDx, fDy, fDz);
545 } 1055 }
546 #endif << 1056
>> 1057 G4NURBS* G4Box::CreateNURBS () const
>> 1058 {
>> 1059 return new G4NURBSbox (fDx, fDy, fDz);
>> 1060 }
547 1061