Geant4  10.00.p03
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()", "UGeomSolids", FatalErrorInArguments, 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 
468 std::ostream& UBox::StreamInfo(std::ostream& os) const
469 {
470  int oldprc = os.precision(16);
471  os << "-----------------------------------------------------------\n"
472  << " *** Dump for solid - " << GetName() << " ***\n"
473  << " ===================================================\n"
474  << " Solid type: UBox\n"
475  << " Parameters: \n"
476  << " half length X: " << fDx << " mm \n"
477  << " half length Y: " << fDy << " mm \n"
478  << " half length Z: " << fDz << " mm \n"
479  << "-----------------------------------------------------------\n";
480  os.precision(oldprc);
481 
482  return os;
483 }
484 
485 void UBox::SetXHalfLength(double dx)
486 {
487  if (dx > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
488  {
489  fDx = dx;
490  }
491  else
492  {
493  std::ostringstream message;
494  message << "Dimension X too small for solid: " << GetName() << "!"
495  << std::endl
496  << " hX = " << dx;
497  UUtils::Exception("UBox::SetXHalfLength()", "GeomSolids0002",
498  FatalErrorInArguments, 1, message.str().c_str());
499  }
500  fCubicVolume = 0.;
501  fSurfaceArea = 0.;
502 
503 }
504 
505 void UBox::SetYHalfLength(double dy)
506 {
507  if (dy > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
508  {
509  fDy = dy;
510  }
511  else
512  {
513  std::ostringstream message;
514  message << "Dimension Y too small for solid: " << GetName() << "!"
515  << std::endl
516  << " hY = " << dy;
517  UUtils::Exception("UBox::SetYHalfLength()", "GeomSolids0002",
518  FatalErrorInArguments, 1, message.str().c_str());
519  }
520  fCubicVolume = 0.;
521  fSurfaceArea = 0.;
522 
523 }
524 
525 void UBox::SetZHalfLength(double dz)
526 {
527  if (dz > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
528  {
529  fDz = dz;
530  }
531  else
532  {
533  std::ostringstream message;
534  message << "Dimension Z too small for solid: " << GetName() << "!"
535  << std::endl
536  << " hZ = " << dz;
537  UUtils::Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
538  FatalErrorInArguments, 1, message.str().c_str());
539  }
540  fCubicVolume = 0.;
541  fSurfaceArea = 0.;
542 
543 }
544 
546 {
547  return "Box";
548 }
549 
551 {
552  return new UBox(*this);
553 }
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:468
double fDz
Definition: UBox.hh:106
const std::string & GetName() const
Definition: VUSolid.hh:103
static double fgTolerance
Definition: VUSolid.hh:30
double fDx
Definition: UBox.hh:104
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:545
double x
Definition: UVector3.hh:136
Definition: UBox.hh:33
virtual ~UBox()
Definition: UBox.cc:57
UBox()
Definition: UBox.hh:37
static const double kInfinity
Definition: UUtils.hh:53
EnumInside
Definition: VUSolid.hh:23
void SetXHalfLength(double dx)
Definition: UBox.cc:485
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
VUSolid * Clone() const
Definition: UBox.cc:550
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:245
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:105
double fCubicVolume
Definition: UBox.hh:107
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 Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
double z
Definition: UVector3.hh:138
void SetYHalfLength(double dy)
Definition: UBox.cc:505
short Sign(short a, short b)
Definition: UUtils.hh:179
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
double y
Definition: UVector3.hh:137
void SetZHalfLength(double dz)
Definition: UBox.cc:525
double fSurfaceArea
Definition: UBox.hh:108