Geant4  10.01.p02
UBox.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UBox
12 //
13 // 10.06.11 J.Apostolakis, G.Cosmo, A.Gheata
14 // Created from original implementation in Geant4 and ROOT
15 // --------------------------------------------------------------------
16 
17 #include <iostream>
18 #include <sstream>
19 #include "UUtils.hh"
20 #include "UBox.hh"
21 
22 //______________________________________________________________________________
23 UBox::UBox(const std::string& name, double dx, double dy, double dz)
24  : VUSolid(name),
25  fDx(dx),
26  fDy(dy),
27  fDz(dz), fCubicVolume(0.), fSurfaceArea(0.)
28 {
29  // Named constructor
30 
31  if ((dx < 2 * VUSolid::fgTolerance)
32  || (dy < 2 * VUSolid::fgTolerance)
33  || (dz < 2 * VUSolid::fgTolerance)) // limit to thickness of surfaces
34  {
35  //std::ostringstream message;
36  std::ostringstream message;
37  message << "Dimensions too small for Solid: " << GetName() << "!" << std::endl
38  << " dx, dy, dz = " << dx << ", " << dy << ", " << dz;
39  UUtils::Exception("UBox::UBox()", "GeomSolids0002", UFatalErrorInArguments, 1, message.str().c_str());
40  }
41 }
42 
43 void UBox::Set(double dx, double dy, double dz)
44 {
45  fDx = dx;
46  fDy = dy;
47  fDz = dz;
48 }
49 
50 void UBox::Set(const UVector3& vec)
51 {
52  fDx = vec.x();
53  fDy = vec.y();
54  fDz = vec.z();
55 }
56 //Destructor
58 {
59 
60 }
61 // Copy constructor
62 
63 UBox::UBox(const UBox& rhs)
64  : VUSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
65  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea)
66 {
67 }
68 
69 //
70 // Assignment operator
71 
73 {
74  // Check assignment to self
75  //
76  if (this == &rhs)
77  {
78  return *this;
79  }
80 
81  // Copy base class data
82  //
83  VUSolid::operator=(rhs);
84 
85  // Copy data
86  //
87  fDx = rhs.fDx;
88  fDy = rhs.fDy;
89  fDz = rhs.fDz;
92 
93  return *this;
94 }
95 
96 //______________________________________________________________________________
98 {
99  // Classify point location with respect to solid:
100  // o eInside - inside the solid
101  // o eSurface - close to surface within tolerance
102  // o eOutside - outside the solid
103  static const double delta = VUSolid::fgTolerance;
104  // Early returns on outside condition on any axis. Check Z first for faster
105  // exclusion in phi symmetric geometries.
106  double ddz = std::abs(aPoint.z()) - fDz;
107  if (ddz > delta) return eOutside;
108  double ddx = std::abs(aPoint.x()) - fDx;
109  if (ddx > delta) return eOutside;
110  double ddy = std::abs(aPoint.y()) - fDy;
111  if (ddy > delta) return eOutside;
112  if (ddx > - delta || ddy > -delta || ddz > -delta) return eSurface;
113  return eInside;
114 }
115 
116 //______________________________________________________________________________
117 double UBox::DistanceToIn(const UVector3& aPoint,
118  const UVector3& aDirection,
119  // UVector3 &aNormal,
120  double aPstep) const
121 {
122  // Computes distance from a point presumably outside the solid to the solid
123  // surface. Ignores first surface if the point is actually inside. Early return
124  // infinity in case the safety to any surface is found greater than the proposed
125  // step aPstep.
126  // The normal vector to the crossed surface is filled only in case the box is
127  // crossed, otherwise aNormal.IsNull() is true.
128 
129  // Compute safety to the closest surface on each axis.
130  // Early exits if safety bigger than proposed step.
131  static const double delta = VUSolid::fgTolerance;
132  // aNormal.SetNull();
133  double safx = std::abs(aPoint.x()) - fDx;
134  double safy = std::abs(aPoint.y()) - fDy;
135  double safz = std::abs(aPoint.z()) - fDz;
136  if ((safx > aPstep) || (safy > aPstep) || (safz > aPstep))
137  return UUtils::kInfinity;
138  // Check numerical outside.
139  bool outside = (safx > 0) || (safy > 0) || (safz > 0);
140  if (!outside)
141  {
142  // If point close to this surface, check against the normal
143  if (safx > -delta)
144  {
145  if (aPoint.x() * aDirection.x() > 0) return UUtils::kInfinity ;
146  }
147  if (safy > -delta)
148  {
149  if (aPoint.y() * aDirection.y() > 0) return UUtils::kInfinity ;
150  }
151  if (safz > -delta)
152  {
153  if (aPoint.z() * aDirection.z() > 0) return UUtils::kInfinity ;
154  }
155  // Point actually "deep" inside, return zero distance, normal un-defined
156  return 0.0;
157  }
158  // The point is really outside. Only axis with positive safety to be
159  // considered. Early exit on each axis if point and direction components
160  // have the same .sign.
161  double dist = 0.0;
162  double coordinate = 0.0;
163  if (safx > 0)
164  {
165  if (aPoint.x() * aDirection.x() >= 0) return UUtils::kInfinity;
166  dist = safx / std::abs(aDirection.x());
167  coordinate = aPoint.y() + dist * aDirection.y();
168  if (std::abs(coordinate) < fDy)
169  {
170  coordinate = aPoint.z() + dist * aDirection.z();
171  if (std::abs(coordinate) < fDz)
172  {
173  // aNormal.x() = UUtils::Sign(1.0, aPoint.x());
174  if (dist < 0.5 * delta) dist = 0.;
175  return dist;
176  }
177  }
178  }
179  if (safy > 0)
180  {
181  if (aPoint.y() * aDirection.y() >= 0) return UUtils::kInfinity;
182  dist = safy / std::abs(aDirection.y());
183  coordinate = aPoint.x() + dist * aDirection.x();
184  if (std::abs(coordinate) < fDx)
185  {
186  coordinate = aPoint.z() + dist * aDirection.z();
187  if (std::abs(coordinate) < fDz)
188  {
189  // aNormal.y() = UUtils::Sign(1.0, aPoint.y());
190  if (dist < 0.5 * delta) dist = 0.;
191  return dist;
192  }
193  }
194  }
195  if (safz > 0)
196  {
197  if (aPoint.z() * aDirection.z() >= 0) return UUtils::kInfinity;
198  dist = safz / std::abs(aDirection.z());
199  coordinate = aPoint.x() + dist * aDirection.x();
200  if (std::abs(coordinate) < fDx)
201  {
202  coordinate = aPoint.y() + dist * aDirection.y();
203  if (std::abs(coordinate) < fDy)
204  {
205  // aNormal.z() = UUtils::Sign(1.0, aPoint.z());
206  if (dist < 0.5 * delta) dist = 0.;
207  return dist;
208  }
209  }
210  }
211  return UUtils::kInfinity;
212 }
213 
214 //______________________________________________________________________________
215 double UBox::DistanceToOut(const UVector3& aPoint, const UVector3& aDirection,
216  UVector3& aNormal,
217  bool& convex,
218  double /*aPstep*/) const
219 {
220  // Computes distance from a point presumably intside the solid to the solid
221  // surface. Ignores first surface along each axis systematically (for points
222  // inside or outside. Early returns zero in case the second surface is behind
223  // the starting point.
224  // o The proposed step is ignored.
225  // o The normal vector to the crossed surface is always filled.
226  double smin = UUtils::kInfinity;
227  double snxt, signDir;
228  convex = true; // Box is convex (even if the starting point is outside)
229  // Check always the "away" surface along direction on axis. This responds
230  // corectly even for points outside the solid (no need for tolerance check)
231  if (aDirection.x() != 0.0)
232  {
233  signDir = UUtils::Sign(1.0, aDirection.x());
234  aNormal.Set(signDir, 0., 0.);
235  snxt = (-aPoint.x() + signDir * fDx) / aDirection.x();
236  if (snxt <= 0) return 0.0; // point outside moving outwards
237  smin = snxt;
238  }
239 
240  if (aDirection.y() != 0.0)
241  {
242  signDir = UUtils::Sign(1.0, aDirection.y());
243  snxt = (-aPoint.y() + signDir * fDy) / aDirection.y();
244  if (snxt <= 0)
245  {
246  aNormal.Set(0., signDir, 0.);
247  return 0.0; // point outside moving outwards
248  }
249  if (snxt < smin)
250  {
251  smin = snxt;
252  aNormal.Set(0., signDir, 0.);
253  }
254  }
255 
256  if (aDirection.z() != 0.0)
257  {
258  signDir = UUtils::Sign(1.0, aDirection.z());
259  snxt = (-aPoint.z() + signDir * fDz) / aDirection.z();
260  if (snxt <= 0)
261  {
262  aNormal.Set(0., 0., signDir);
263  return 0.0; // point outside moving outwards
264  }
265  if (snxt < smin)
266  {
267  smin = snxt;
268  aNormal.Set(0., 0., signDir);
269  }
270  }
271  if (smin < 0.5 * VUSolid::fgTolerance) smin = 0.;
272  return smin;
273 }
274 
275 //______________________________________________________________________________
276 double UBox::SafetyFromInside(const UVector3& aPoint,
277  bool /*aAccurate*/) const
278 {
279  // Estimates the isotropic safety from a point inside the current solid to any
280  // of its surfaces. The algorithm may be accurate or should provide a fast
281  // underestimate.
282  double safe, safy, safz;
283  safe = fDx - std::abs(aPoint.x());
284  safy = fDy - std::abs(aPoint.y());
285  if (safy < safe) safe = safy;
286  safz = fDz - std::abs(aPoint.z());
287  if (safz < safe) safe = safz;
288  return std::max(0.0, safe);
289 }
290 
291 //______________________________________________________________________________
292 double UBox::SafetyFromOutside(const UVector3& aPoint,
293  bool aAccurate) const
294 {
295  // Estimates the isotropic safety from a point outside the current solid to any
296  // of its surfaces. The algorithm may be accurate or should provide a fast
297  // underestimate.
298  double safe, safx, safy, safz;
299  safe = safx = -fDx + std::abs(aPoint.x());
300  safy = -fDy + std::abs(aPoint.y());
301  if (safy > safe) safe = safy;
302  safz = -fDz + std::abs(aPoint.z());
303  if (safz > safe) safe = safz;
304  if (safe < 0.0) return 0.0; // point is inside
305  if (!aAccurate) return safe;
306  double safsq = 0.0;
307  int count = 0;
308  if (safx > 0)
309  {
310  safsq += safx * safx;
311  count++;
312  }
313  if (safy > 0)
314  {
315  safsq += safy * safy;
316  count++;
317  }
318  if (safz > 0)
319  {
320  safsq += safz * safz;
321  count++;
322  }
323  if (count == 1) return safe;
324  return std::sqrt(safsq);
325 }
326 
327 //______________________________________________________________________________
328 bool UBox::Normal(const UVector3& aPoint, UVector3& aNormal) const
329 {
330  // Computes the normal on a surface and returns it as a unit vector
331  // In case a point is further than tolerance_normal from a surface, set validNormal=false
332  // Must return a valid vector. (even if the point is not on the surface.)
333  //
334  // On an edge or corner, provide an average normal of all facets within tolerance
335  // NOTE: the tolerance value used in here is not yet the global surface
336  // tolerance - we will have to revise this value - TODO
337  static const double delta = 100.*VUSolid::fgTolerance;
338  static const double kInvSqrt2 = 1. / std::sqrt(2.);
339  static const double kInvSqrt3 = 1. / std::sqrt(3.);
340  aNormal.Set(0.);
341  UVector3 crt_normal, min_normal;
342  int nsurf = 0;
343  double safx = std::abs(std::abs(aPoint.x()) - fDx);
344  double safmin = safx;
345  crt_normal.Set(UUtils::Sign(1., aPoint.x()), 0., 0.);
346  min_normal = crt_normal;
347  if (safx < delta)
348  {
349  nsurf++;
350  aNormal += crt_normal;
351  }
352  double safy = std::abs(std::abs(aPoint.y()) - fDy);
353  crt_normal.Set(0., UUtils::Sign(1., aPoint.y()), 0.);
354  if (safy < delta)
355  {
356  nsurf++;
357  aNormal += crt_normal;
358  }
359  if (safy < safmin)
360  {
361  min_normal = crt_normal;
362  safmin = safy;
363  }
364  double safz = std::abs(std::abs(aPoint.z()) - fDz);
365  crt_normal.Set(0., 0., UUtils::Sign(1., aPoint.z()));
366  if (safz < delta)
367  {
368  nsurf++;
369  aNormal += crt_normal;
370  }
371  if (safz < safmin)
372  {
373  min_normal = crt_normal;
374  safmin = safz;
375  }
376 
377  bool valid = true;
378  switch (nsurf)
379  {
380  case 0:
381  aNormal = min_normal;
382  valid = false;
383  break;
384  case 1:
385  break;
386  case 2:
387  aNormal *= kInvSqrt2;
388  break;
389  case 3:
390  aNormal *= kInvSqrt3;
391  };
392  return valid;
393 }
394 
395 
396 //______________________________________________________________________________
397 void UBox::Extent(UVector3& aMin, UVector3& aMax) const
398 {
399  // Returns the full 3D cartesian extent of the solid.
400  aMin.x() = -fDx;
401  aMax.x() = fDx;
402  aMin.y() = -fDy;
403  aMax.y() = fDy;
404  aMin.z() = -fDz;
405  aMax.z() = fDz;
406 }
407 
409 //
410 // GetPointOnSurface
411 //
412 // Return a point (UVector3) randomly and uniformly selected
413 // on the solid surface
414 
416 {
417  double px, py, pz, select, sumS;
418  double Sxy = fDx * fDy, Sxz = fDx * fDz, Syz = fDy * fDz;
419 
420  sumS = Sxy + Sxz + Syz;
421  select = sumS * UUtils::Random();
422 
423  if (select < Sxy)
424  {
425  px = -fDx + 2 * fDx * UUtils::Random();
426  py = -fDy + 2 * fDy * UUtils::Random();
427 
428  if (UUtils::Random() > 0.5)
429  {
430  pz = fDz;
431  }
432  else
433  {
434  pz = -fDz;
435  }
436  }
437  else if ((select - Sxy) < Sxz)
438  {
439  px = -fDx + 2 * fDx * UUtils::Random();
440  pz = -fDz + 2 * fDz * UUtils::Random();
441 
442  if (UUtils::Random() > 0.5)
443  {
444  py = fDy;
445  }
446  else
447  {
448  py = -fDy;
449  }
450  }
451  else
452  {
453  py = -fDy + 2 * fDy * UUtils::Random();
454  pz = -fDz + 2 * fDz * UUtils::Random();
455 
456  if (UUtils::Random() > 0.5)
457  {
458  px = fDx;
459  }
460  else
461  {
462  px = -fDx;
463  }
464  }
465  return UVector3(px, py, pz);
466 }
467 
469 //
470 // GetPointOnEdge
471 //
472 // Return a point (UVector3) randomly and uniformly selected
473 // on the solid edge
474 
476 {
477  double select, sumL;
478  double Lx = 2 * fDx, Ly = 2 * fDy, Lz= 2 * fDz;
479 
480  sumL = Lx + Ly + Lz;
481  select = sumL * UUtils::Random();
482 
483  if (select < Lx)
484  {
485  select = UUtils::Random();
486  if (select < 0.25) return UVector3( -fDx + 2 * fDx* UUtils::Random(),-fDy,-fDz);
487  if (select < 0.5 ) return UVector3( -fDx + 2 * fDx* UUtils::Random(), fDy,-fDz);
488  if (select < 0.75) return UVector3( -fDx + 2 * fDx* UUtils::Random(),-fDy, fDz);
489  return UVector3( -fDx + 2 * fDx* UUtils::Random(),fDy,fDz);
490  }
491  else if ((select - Lx) < Ly)
492  {
493  select = UUtils::Random();
494  if (select < 0.25) return UVector3( -fDx,-fDy + 2 * fDy* UUtils::Random(),-fDz);
495  if (select < 0.5 ) return UVector3( fDx,-fDy + 2 * fDy* UUtils::Random(),-fDz);
496  if (select < 0.75) return UVector3( -fDx,-fDy + 2 * fDy* UUtils::Random(), fDz);
497  return UVector3( fDx,-fDy + 2 * fDy* UUtils::Random(),fDz);
498  }
499  else
500  {
501  select = UUtils::Random();
502  if (select < 0.25) return UVector3( -fDx,-fDy,-fDz + 2 * fDz* UUtils::Random());
503  if (select < 0.5 ) return UVector3( fDx,-fDy,-fDz + 2 * fDz* UUtils::Random());
504  if (select < 0.75) return UVector3( -fDx, fDy,-fDz + 2 * fDz* UUtils::Random());
505  return UVector3( fDx,fDy,-fDz + 2 * fDz* UUtils::Random());
506 
507  }
508  return 0;
509 }
510 
511 std::ostream& UBox::StreamInfo(std::ostream& os) const
512 {
513  int oldprc = os.precision(16);
514  os << "-----------------------------------------------------------\n"
515  << " *** Dump for solid - " << GetName() << " ***\n"
516  << " ===================================================\n"
517  << " Solid type: UBox\n"
518  << " Parameters: \n"
519  << " half length X: " << fDx << " mm \n"
520  << " half length Y: " << fDy << " mm \n"
521  << " half length Z: " << fDz << " mm \n"
522  << "-----------------------------------------------------------\n";
523  os.precision(oldprc);
524 
525  return os;
526 }
527 
528 void UBox::SetXHalfLength(double dx)
529 {
530  if (dx > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
531  {
532  fDx = dx;
533  }
534  else
535  {
536  std::ostringstream message;
537  message << "Dimension X too small for solid: " << GetName() << "!"
538  << std::endl
539  << " hX = " << dx;
540  UUtils::Exception("UBox::SetXHalfLength()", "GeomSolids0002",
541  UFatalErrorInArguments, 1, message.str().c_str());
542  }
543  fCubicVolume = 0.;
544  fSurfaceArea = 0.;
545 
546 }
547 
548 void UBox::SetYHalfLength(double dy)
549 {
550  if (dy > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
551  {
552  fDy = dy;
553  }
554  else
555  {
556  std::ostringstream message;
557  message << "Dimension Y too small for solid: " << GetName() << "!"
558  << std::endl
559  << " hY = " << dy;
560  UUtils::Exception("UBox::SetYHalfLength()", "GeomSolids0002",
561  UFatalErrorInArguments, 1, message.str().c_str());
562  }
563  fCubicVolume = 0.;
564  fSurfaceArea = 0.;
565 
566 }
567 
568 void UBox::SetZHalfLength(double dz)
569 {
570  if (dz > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
571  {
572  fDz = dz;
573  }
574  else
575  {
576  std::ostringstream message;
577  message << "Dimension Z too small for solid: " << GetName() << "!"
578  << std::endl
579  << " hZ = " << dz;
580  UUtils::Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
581  UFatalErrorInArguments, 1, message.str().c_str());
582  }
583  fCubicVolume = 0.;
584  fSurfaceArea = 0.;
585 
586 }
587 
589 {
590  return "Box";
591 }
592 
594 {
595  return new UBox(*this);
596 }
double & z()
Definition: UVector3.hh:352
double & y()
Definition: UVector3.hh:348
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Definition: UBox.cc:328
UVector3 GetPointOnSurface() const
Definition: UBox.cc:415
std::ostream & StreamInfo(std::ostream &os) const
Definition: UBox.cc:511
double fDz
Definition: UBox.hh:107
const std::string & GetName() const
Definition: VUSolid.hh:103
static double fgTolerance
Definition: VUSolid.hh:30
double fDx
Definition: UBox.hh:105
G4String name
Definition: TRTMaterials.hh:40
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UBox.cc:215
UGeometryType GetEntityType() const
Definition: UBox.cc:588
Definition: UBox.hh:33
virtual ~UBox()
Definition: UBox.cc:57
UBox()
Definition: UBox.hh:37
static const double kInfinity
Definition: UUtils.hh:54
double & x()
Definition: UVector3.hh:344
UVector3 GetPointOnEdge() const
Definition: UBox.cc:475
EnumInside
Definition: VUSolid.hh:23
void SetXHalfLength(double dx)
Definition: UBox.cc:528
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UBox.cc:117
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UBox.cc:276
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
VUSolid * Clone() const
Definition: UBox.cc:593
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:264
T max(const T t1, const T t2)
brief Return the largest of the two arguments
UBox & operator=(const UBox &rhs)
Definition: UBox.cc:72
double fDy
Definition: UBox.hh:106
double fCubicVolume
Definition: UBox.hh:108
void Set(double dx, double dy, double dz)
Definition: UBox.cc:43
std::string UGeometryType
Definition: UTypes.hh:39
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UBox.cc:292
EnumInside Inside(const UVector3 &aPoint) const
Definition: UBox.cc:97
void SetYHalfLength(double dy)
Definition: UBox.cc:548
short Sign(short a, short b)
Definition: UUtils.hh:180
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UBox.cc:397
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
void SetZHalfLength(double dz)
Definition: UBox.cc:568
double fSurfaceArea
Definition: UBox.hh:109