| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation for G4Box class
27 //
28 // 30.06.95 - P.Kent: First version
29 // 20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
30 // 18.04.17 - E.Tcherniaev: complete revision, speed-up
31 // --------------------------------------------------------------------
32
33 #include "G4Box.hh"
34
35 #if !defined(G4GEOM_USE_UBOX)
36
37 #include "G4SystemOfUnits.hh"
38 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh"
41 #include "G4QuickRand.hh"
42
43 #include "G4VPVParameterisation.hh"
44
45 #include "G4VGraphicsScene.hh"
46 #include "G4VisExtent.hh"
47
48 ////////////////////////////////////////////////////////////////////////
49 //
50 // Constructor - check & set half widths
51
52 G4Box::G4Box(const G4String& pName,
53 G4double pX,
54 G4double pY,
55 G4double pZ)
56 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
57 {
58 delta = 0.5*kCarTolerance;
59 if (pX < 2*kCarTolerance ||
60 pY < 2*kCarTolerance ||
61 pZ < 2*kCarTolerance) // limit to thickness of surfaces
62 {
63 std::ostringstream message;
64 message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
65 << " hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
66 G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
67 }
68 }
69
70 //////////////////////////////////////////////////////////////////////////
71 //
72 // Fake default constructor - sets only member data and allocates memory
73 // for usage restricted to object persistency
74
75 G4Box::G4Box( __void__& a )
76 : G4CSGSolid(a)
77 {
78 }
79
80 //////////////////////////////////////////////////////////////////////////
81 //
82 // Destructor
83
84 G4Box::~G4Box() = default;
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 {
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 }
115
116 //////////////////////////////////////////////////////////////////////////
117 //
118 // Set X dimension
119
120 void G4Box::SetXHalfLength(G4double dx)
121 {
122 if(dx > 2*kCarTolerance) // limit to thickness of surfaces
123 {
124 fDx = dx;
125 }
126 else
127 {
128 std::ostringstream message;
129 message << "Dimension X too small for solid: " << GetName() << "!"
130 << G4endl
131 << " hX = " << dx;
132 G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
133 FatalException, message);
134 }
135 fCubicVolume = 0.;
136 fSurfaceArea = 0.;
137 fRebuildPolyhedron = true;
138 }
139
140 //////////////////////////////////////////////////////////////////////////
141 //
142 // Set Y dimension
143
144 void G4Box::SetYHalfLength(G4double dy)
145 {
146 if(dy > 2*kCarTolerance) // limit to thickness of surfaces
147 {
148 fDy = dy;
149 }
150 else
151 {
152 std::ostringstream message;
153 message << "Dimension Y too small for solid: " << GetName() << "!\n"
154 << " hY = " << dy;
155 G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
156 FatalException, message);
157 }
158 fCubicVolume = 0.;
159 fSurfaceArea = 0.;
160 fRebuildPolyhedron = true;
161 }
162
163 //////////////////////////////////////////////////////////////////////////
164 //
165 // Set Z dimension
166
167 void G4Box::SetZHalfLength(G4double dz)
168 {
169 if(dz > 2*kCarTolerance) // limit to thickness of surfaces
170 {
171 fDz = dz;
172 }
173 else
174 {
175 std::ostringstream message;
176 message << "Dimension Z too small for solid: " << GetName() << "!\n"
177 << " hZ = " << dz;
178 G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
179 FatalException, message);
180 }
181 fCubicVolume = 0.;
182 fSurfaceArea = 0.;
183 fRebuildPolyhedron = true;
184 }
185
186 //////////////////////////////////////////////////////////////////////////
187 //
188 // Dispatch to parameterisation for replication mechanism dimension
189 // computation & modification.
190
191 void G4Box::ComputeDimensions(G4VPVParameterisation* p,
192 const G4int n,
193 const G4VPhysicalVolume* pRep)
194 {
195 p->ComputeDimensions(*this,n,pRep);
196 }
197
198 //////////////////////////////////////////////////////////////////////////
199 //
200 // Get bounding box
201
202 void G4Box::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
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.y() || pMin.z() >= pMax.z())
210 {
211 std::ostringstream message;
212 message << "Bad bounding box (min >= max) for solid: "
213 << GetName() << " !"
214 << "\npMin = " << pMin
215 << "\npMax = " << pMax;
216 G4Exception("G4Box::BoundingLimits()", "GeomMgt0001", JustWarning, message);
217 DumpInfo();
218 }
219 }
220
221 //////////////////////////////////////////////////////////////////////////
222 //
223 // Calculate extent under transform and specified limit
224
225 G4bool G4Box::CalculateExtent(const EAxis pAxis,
226 const G4VoxelLimits& pVoxelLimit,
227 const G4AffineTransform& pTransform,
228 G4double& pMin, G4double& pMax) const
229 {
230 G4ThreeVector bmin, bmax;
231
232 // Get bounding box
233 BoundingLimits(bmin,bmax);
234
235 // Find extent
236 G4BoundingEnvelope bbox(bmin,bmax);
237 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
238 }
239
240 //////////////////////////////////////////////////////////////////////////
241 //
242 // Return whether point inside/outside/on surface, using tolerance
243
244 EInside G4Box::Inside(const G4ThreeVector& p) const
245 {
246 G4double dist = std::max(std::max(
247 std::abs(p.x())-fDx,
248 std::abs(p.y())-fDy),
249 std::abs(p.z())-fDz);
250 return (dist > delta) ? kOutside :
251 ((dist > -delta) ? kSurface : kInside);
252 }
253
254 //////////////////////////////////////////////////////////////////////////
255 //
256 // Detect the side(s) and return corresponding normal
257
258 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
259 {
260 G4ThreeVector norm(0,0,0);
261 G4double px = p.x();
262 if (std::abs(std::abs(px) - fDx) <= delta) norm.setX(px < 0 ? -1. : 1.);
263 G4double py = p.y();
264 if (std::abs(std::abs(py) - fDy) <= delta) norm.setY(py < 0 ? -1. : 1.);
265 G4double pz = p.z();
266 if (std::abs(std::abs(pz) - fDz) <= delta) norm.setZ(pz < 0 ? -1. : 1.);
267
268 G4double nside = norm.mag2(); // number of sides = magnitude squared
269 if (nside == 1)
270 return norm;
271 else if (nside > 1)
272 return norm.unit(); // edge or corner
273 else
274 {
275 // Point is not on the surface
276 //
277 #ifdef G4CSGDEBUG
278 std::ostringstream message;
279 G4int oldprc = message.precision(16);
280 message << "Point p is not on surface (!?) of solid: "
281 << GetName() << G4endl;
282 message << "Position:\n";
283 message << " p.x() = " << p.x()/mm << " mm\n";
284 message << " p.y() = " << p.y()/mm << " mm\n";
285 message << " p.z() = " << p.z()/mm << " mm";
286 G4cout.precision(oldprc);
287 G4Exception("G4Box::SurfaceNormal(p)", "GeomSolids1002",
288 JustWarning, message );
289 DumpInfo();
290 #endif
291 return ApproxSurfaceNormal(p);
292 }
293 }
294
295 //////////////////////////////////////////////////////////////////////////
296 //
297 // Algorithm for SurfaceNormal() following the original specification
298 // for points not on the surface
299
300 G4ThreeVector G4Box::ApproxSurfaceNormal(const G4ThreeVector& p) const
301 {
302 G4double distx = std::abs(p.x()) - fDx;
303 G4double disty = std::abs(p.y()) - fDy;
304 G4double distz = std::abs(p.z()) - fDz;
305
306 if (distx >= disty && distx >= distz)
307 return {std::copysign(1.,p.x()), 0., 0.};
308 if (disty >= distx && disty >= distz)
309 return {0., std::copysign(1.,p.y()), 0.};
310 else
311 return {0., 0., std::copysign(1.,p.z())};
312 }
313
314 //////////////////////////////////////////////////////////////////////////
315 //
316 // Calculate distance to box from an outside point
317 // - return kInfinity if no intersection
318 //
319
320 G4double G4Box::DistanceToIn(const G4ThreeVector& p,
321 const G4ThreeVector& v) const
322 {
323 // Check if point is on the surface and traveling away
324 //
325 if ((std::abs(p.x()) - fDx) >= -delta && p.x()*v.x() >= 0) return kInfinity;
326 if ((std::abs(p.y()) - fDy) >= -delta && p.y()*v.y() >= 0) return kInfinity;
327 if ((std::abs(p.z()) - fDz) >= -delta && p.z()*v.z() >= 0) return kInfinity;
328
329 // Find intersection
330 //
331 G4double invx = (v.x() == 0) ? DBL_MAX : -1./v.x();
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./v.y();
337 G4double dy = std::copysign(fDy,invy);
338 G4double tymin = std::max(txmin,(p.y() - dy)*invy);
339 G4double tymax = std::min(txmax,(p.y() + dy)*invy);
340
341 G4double invz = (v.z() == 0) ? DBL_MAX : -1./v.z();
342 G4double dz = std::copysign(fDz,invz);
343 G4double tmin = std::max(tymin,(p.z() - dz)*invz);
344 G4double tmax = std::min(tymax,(p.z() + dz)*invz);
345
346 if (tmax <= tmin + delta) return kInfinity; // touch or no hit
347 return (tmin < delta) ? 0. : tmin;
348 }
349
350 //////////////////////////////////////////////////////////////////////////
351 //
352 // Appoximate distance to box.
353 // Returns largest perpendicular distance to the closest x/y/z sides of
354 // the box, which is the most fast estimation of the shortest distance to box
355 // - If inside return 0
356
357 G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
358 {
359 G4double dist = std::max(std::max(
360 std::abs(p.x())-fDx,
361 std::abs(p.y())-fDy),
362 std::abs(p.z())-fDz);
363 return (dist > 0) ? dist : 0.;
364 }
365
366 //////////////////////////////////////////////////////////////////////////
367 //
368 // Calculate distance to surface of the box from inside and
369 // find normal at exit point, if required
370 // - when leaving the surface, return 0
371
372 G4double G4Box::DistanceToOut( const G4ThreeVector& p,
373 const G4ThreeVector& v,
374 const G4bool calcNorm,
375 G4bool* validNorm, G4ThreeVector* n) const
376 {
377 // Check if point is on the surface and traveling away
378 //
379 if ((std::abs(p.x()) - fDx) >= -delta && p.x()*v.x() > 0)
380 {
381 if (calcNorm)
382 {
383 *validNorm = true;
384 n->set((p.x() < 0) ? -1. : 1., 0., 0.);
385 }
386 return 0.;
387 }
388 if ((std::abs(p.y()) - fDy) >= -delta && p.y()*v.y() > 0)
389 {
390 if (calcNorm)
391 {
392 *validNorm = true;
393 n->set(0., (p.y() < 0) ? -1. : 1., 0.);
394 }
395 return 0.;
396 }
397 if ((std::abs(p.z()) - fDz) >= -delta && p.z()*v.z() > 0)
398 {
399 if (calcNorm)
400 {
401 *validNorm = true;
402 n->set(0., 0., (p.z() < 0) ? -1. : 1.);
403 }
404 return 0.;
405 }
406
407 // Find intersection
408 //
409 G4double vx = v.x();
410 G4double tx = (vx == 0) ? DBL_MAX : (std::copysign(fDx,vx) - p.x())/vx;
411
412 G4double vy = v.y();
413 G4double ty = (vy == 0) ? tx : (std::copysign(fDy,vy) - p.y())/vy;
414 G4double txy = std::min(tx,ty);
415
416 G4double vz = v.z();
417 G4double tz = (vz == 0) ? txy : (std::copysign(fDz,vz) - p.z())/vz;
418 G4double tmax = std::min(txy,tz);
419
420 // Set normal, if required, and return distance
421 //
422 if (calcNorm)
423 {
424 *validNorm = true;
425 if (tmax == tx) n->set((v.x() < 0) ? -1. : 1., 0., 0.);
426 else if (tmax == ty) n->set(0., (v.y() < 0) ? -1. : 1., 0.);
427 else n->set(0., 0., (v.z() < 0) ? -1. : 1.);
428 }
429 return tmax;
430 }
431
432 ////////////////////////////////////////////////////////////////////////////
433 //
434 // Calculate exact shortest distance to any boundary from inside
435 // - if outside return 0
436
437 G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
438 {
439 #ifdef G4CSGDEBUG
440 if( Inside(p) == kOutside )
441 {
442 std::ostringstream message;
443 G4int oldprc = message.precision(16);
444 message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
445 message << "Position:\n";
446 message << " p.x() = " << p.x()/mm << " mm\n";
447 message << " p.y() = " << p.y()/mm << " mm\n";
448 message << " p.z() = " << p.z()/mm << " mm";
449 G4cout.precision(oldprc);
450 G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
451 JustWarning, message );
452 DumpInfo();
453 }
454 #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
462 //////////////////////////////////////////////////////////////////////////
463 //
464 // GetEntityType
465
466 G4GeometryType G4Box::GetEntityType() const
467 {
468 return {"G4Box"};
469 }
470
471 //////////////////////////////////////////////////////////////////////////
472 //
473 // IsFaceted
474
475 G4bool G4Box::IsFaceted() const
476 {
477 return true;
478 }
479
480 //////////////////////////////////////////////////////////////////////////
481 //
482 // Stream object contents to an output stream
483
484 std::ostream& G4Box::StreamInfo(std::ostream& os) const
485 {
486 G4long oldprc = os.precision(16);
487 os << "-----------------------------------------------------------\n"
488 << " *** Dump for solid - " << GetName() << " ***\n"
489 << " ===================================================\n"
490 << "Solid type: G4Box\n"
491 << "Parameters: \n"
492 << " half length X: " << fDx/mm << " mm \n"
493 << " half length Y: " << fDy/mm << " mm \n"
494 << " half length Z: " << fDz/mm << " mm \n"
495 << "-----------------------------------------------------------\n";
496 os.precision(oldprc);
497 return os;
498 }
499
500 //////////////////////////////////////////////////////////////////////////
501 //
502 // Return a point randomly and uniformly selected on the surface
503
504 G4ThreeVector G4Box::GetPointOnSurface() const
505 {
506 G4double sxy = fDx*fDy, sxz = fDx*fDz, syz = fDy*fDz;
507 G4double select = (sxy + sxz + syz)*G4QuickRand();
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) ? -fDz : fDz) };
513 else if (select < sxy + sxz)
514 return { u*fDx, ((select < sxy + 0.5*sxz) ? -fDy : fDy), v*fDz };
515 else
516 return { ((select < sxy + sxz + 0.5*syz) ? -fDx : fDx), u*fDy, v*fDz };
517 }
518
519 //////////////////////////////////////////////////////////////////////////
520 //
521 // Make a clone of the object
522 //
523 G4VSolid* G4Box::Clone() const
524 {
525 return new G4Box(*this);
526 }
527
528 //////////////////////////////////////////////////////////////////////////
529 //
530 // Methods for visualisation
531
532 void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const
533 {
534 scene.AddSolid (*this);
535 }
536
537 G4VisExtent G4Box::GetExtent() const
538 {
539 return { -fDx, fDx, -fDy, fDy, -fDz, fDz };
540 }
541
542 G4Polyhedron* G4Box::CreatePolyhedron () const
543 {
544 return new G4PolyhedronBox (fDx, fDy, fDz);
545 }
546 #endif
547