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