Geant4_10
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;
90  fCubicVolume = rhs.fCubicVolume;
91  fSurfaceArea = rhs.fSurfaceArea;
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  // aNormal.x = UUtils::Sign(1.0, aPoint.x);
146  return (aPoint.x * aDirection.x > 0) ? UUtils::kInfinity : 0.0;
147  }
148  if (safy > -delta)
149  {
150  // aNormal.y = UUtils::Sign(1.0, aPoint.y);
151  return (aPoint.y * aDirection.y > 0) ? UUtils::kInfinity : 0.0;
152  }
153  if (safz > -delta)
154  {
155  // aNormal.z = UUtils::Sign(1.0, aPoint.z);
156  return (aPoint.z * aDirection.z > 0) ? UUtils::kInfinity : 0.0;
157  }
158  // Point actually "deep" inside, return zero distance, normal un-defined
159  return 0.0;
160  }
161  // The point is really outside. Only axis with positive safety to be
162  // considered. Early exit on each axis if point and direction components
163  // have the same .sign.
164  double dist = 0.0;
165  double coordinate = 0.0;
166  if (safx > 0)
167  {
168  if (aPoint.x * aDirection.x >= 0) return UUtils::kInfinity;
169  dist = safx / std::abs(aDirection.x);
170  coordinate = aPoint.y + dist * aDirection.y;
171  if (std::abs(coordinate) < fDy)
172  {
173  coordinate = aPoint.z + dist * aDirection.z;
174  if (std::abs(coordinate) < fDz)
175  {
176  // aNormal.x = UUtils::Sign(1.0, aPoint.x);
177  if (dist < 0.5 * delta) dist = 0.;
178  return dist;
179  }
180  }
181  }
182  if (safy > 0)
183  {
184  if (aPoint.y * aDirection.y >= 0) return UUtils::kInfinity;
185  dist = safy / std::abs(aDirection.y);
186  coordinate = aPoint.x + dist * aDirection.x;
187  if (std::abs(coordinate) < fDx)
188  {
189  coordinate = aPoint.z + dist * aDirection.z;
190  if (std::abs(coordinate) < fDz)
191  {
192  // aNormal.y = UUtils::Sign(1.0, aPoint.y);
193  if (dist < 0.5 * delta) dist = 0.;
194  return dist;
195  }
196  }
197  }
198  if (safz > 0)
199  {
200  if (aPoint.z * aDirection.z >= 0) return UUtils::kInfinity;
201  dist = safz / std::abs(aDirection.z);
202  coordinate = aPoint.x + dist * aDirection.x;
203  if (std::abs(coordinate) < fDx)
204  {
205  coordinate = aPoint.y + dist * aDirection.y;
206  if (std::abs(coordinate) < fDy)
207  {
208  // aNormal.z = UUtils::Sign(1.0, aPoint.z);
209  if (dist < 0.5 * delta) dist = 0.;
210  return dist;
211  }
212  }
213  }
214  return UUtils::kInfinity;
215 }
216 
217 //______________________________________________________________________________
218 double UBox::DistanceToOut(const UVector3& aPoint, const UVector3& aDirection,
219  UVector3& aNormal,
220  bool& convex,
221  double /*aPstep*/) const
222 {
223  // Computes distance from a point presumably intside the solid to the solid
224  // surface. Ignores first surface along each axis systematically (for points
225  // inside or outside. Early returns zero in case the second surface is behind
226  // the starting point.
227  // o The proposed step is ignored.
228  // o The normal vector to the crossed surface is always filled.
229  double smin = UUtils::kInfinity;
230  double snxt, signDir;
231  convex = true; // Box is convex (even if the starting point is outside)
232  // Check always the "away" surface along direction on axis. This responds
233  // corectly even for points outside the solid (no need for tolerance check)
234  if (aDirection.x != 0.0)
235  {
236  signDir = UUtils::Sign(1.0, aDirection.x);
237  aNormal.Set(signDir, 0., 0.);
238  snxt = (-aPoint.x + signDir * fDx) / aDirection.x;
239  if (snxt <= 0) return 0.0; // point outside moving outwards
240  smin = snxt;
241  }
242 
243  if (aDirection.y != 0.0)
244  {
245  signDir = UUtils::Sign(1.0, aDirection.y);
246  snxt = (-aPoint.y + signDir * fDy) / aDirection.y;
247  if (snxt <= 0)
248  {
249  aNormal.Set(0., signDir, 0.);
250  return 0.0; // point outside moving outwards
251  }
252  if (snxt < smin)
253  {
254  smin = snxt;
255  aNormal.Set(0., signDir, 0.);
256  }
257  }
258 
259  if (aDirection.z != 0.0)
260  {
261  signDir = UUtils::Sign(1.0, aDirection.z);
262  snxt = (-aPoint.z + signDir * fDz) / aDirection.z;
263  if (snxt <= 0)
264  {
265  aNormal.Set(0., 0., signDir);
266  return 0.0; // point outside moving outwards
267  }
268  if (snxt < smin)
269  {
270  smin = snxt;
271  aNormal.Set(0., 0., signDir);
272  }
273  }
274  if (smin < 0.5 * VUSolid::fgTolerance) smin = 0.;
275  return smin;
276 }
277 
278 //______________________________________________________________________________
279 double UBox::SafetyFromInside(const UVector3& aPoint,
280  bool /*aAccurate*/) const
281 {
282  // Estimates the isotropic safety from a point inside the current solid to any
283  // of its surfaces. The algorithm may be accurate or should provide a fast
284  // underestimate.
285  double safe, safy, safz;
286  safe = fDx - std::abs(aPoint.x);
287  safy = fDy - std::abs(aPoint.y);
288  if (safy < safe) safe = safy;
289  safz = fDz - std::abs(aPoint.z);
290  if (safz < safe) safe = safz;
291  return std::max(0.0, safe);
292 }
293 
294 //______________________________________________________________________________
295 double UBox::SafetyFromOutside(const UVector3& aPoint,
296  bool aAccurate) const
297 {
298  // Estimates the isotropic safety from a point outside the current solid to any
299  // of its surfaces. The algorithm may be accurate or should provide a fast
300  // underestimate.
301  double safe, safx, safy, safz;
302  safe = safx = -fDx + std::abs(aPoint.x);
303  safy = -fDy + std::abs(aPoint.y);
304  if (safy > safe) safe = safy;
305  safz = -fDz + std::abs(aPoint.z);
306  if (safz > safe) safe = safz;
307  if (safe < 0.0) return 0.0; // point is inside
308  if (!aAccurate) return safe;
309  double safsq = 0.0;
310  int count = 0;
311  if (safx > 0)
312  {
313  safsq += safx * safx;
314  count++;
315  }
316  if (safy > 0)
317  {
318  safsq += safy * safy;
319  count++;
320  }
321  if (safz > 0)
322  {
323  safsq += safz * safz;
324  count++;
325  }
326  if (count == 1) return safe;
327  return std::sqrt(safsq);
328 }
329 
330 //______________________________________________________________________________
331 bool UBox::Normal(const UVector3& aPoint, UVector3& aNormal) const
332 {
333  // Computes the normal on a surface and returns it as a unit vector
334  // In case a point is further than tolerance_normal from a surface, set validNormal=false
335  // Must return a valid vector. (even if the point is not on the surface.)
336  //
337  // On an edge or corner, provide an average normal of all facets within tolerance
338  // NOTE: the tolerance value used in here is not yet the global surface
339  // tolerance - we will have to revise this value - TODO
340  static const double delta = 100.*VUSolid::fgTolerance;
341  static const double kInvSqrt2 = 1. / std::sqrt(2.);
342  static const double kInvSqrt3 = 1. / std::sqrt(3.);
343  aNormal.Set(0.);
344  UVector3 crt_normal, min_normal;
345  int nsurf = 0;
346  double safx = std::abs(std::abs(aPoint.x) - fDx);
347  double safmin = safx;
348  crt_normal.Set(UUtils::Sign(1., aPoint.x), 0., 0.);
349  min_normal = crt_normal;
350  if (safx < delta)
351  {
352  nsurf++;
353  aNormal += crt_normal;
354  }
355  double safy = std::abs(std::abs(aPoint.y) - fDy);
356  crt_normal.Set(0., UUtils::Sign(1., aPoint.y), 0.);
357  if (safy < delta)
358  {
359  nsurf++;
360  aNormal += crt_normal;
361  }
362  if (safy < safmin)
363  {
364  min_normal = crt_normal;
365  safmin = safy;
366  }
367  double safz = std::abs(std::abs(aPoint.z) - fDz);
368  crt_normal.Set(0., 0., UUtils::Sign(1., aPoint.z));
369  if (safz < delta)
370  {
371  nsurf++;
372  aNormal += crt_normal;
373  }
374  if (safz < safmin)
375  {
376  min_normal = crt_normal;
377  safmin = safz;
378  }
379 
380  bool valid = true;
381  switch (nsurf)
382  {
383  case 0:
384  aNormal = min_normal;
385  valid = false;
386  break;
387  case 1:
388  break;
389  case 2:
390  aNormal *= kInvSqrt2;
391  break;
392  case 3:
393  aNormal *= kInvSqrt3;
394  };
395  return valid;
396 }
397 
398 
399 //______________________________________________________________________________
400 void UBox::Extent(UVector3& aMin, UVector3& aMax) const
401 {
402  // Returns the full 3D cartesian extent of the solid.
403  aMin.x = -fDx;
404  aMax.x = fDx;
405  aMin.y = -fDy;
406  aMax.y = fDy;
407  aMin.z = -fDz;
408  aMax.z = fDz;
409 }
410 
412 //
413 // GetPointOnSurface
414 //
415 // Return a point (UVector3) randomly and uniformly selected
416 // on the solid surface
417 
419 {
420  double px, py, pz, select, sumS;
421  double Sxy = fDx * fDy, Sxz = fDx * fDz, Syz = fDy * fDz;
422 
423  sumS = Sxy + Sxz + Syz;
424  select = sumS * UUtils::Random();
425 
426  if (select < Sxy)
427  {
428  px = -fDx + 2 * fDx * UUtils::Random();
429  py = -fDy + 2 * fDy * UUtils::Random();
430 
431  if (UUtils::Random() > 0.5)
432  {
433  pz = fDz;
434  }
435  else
436  {
437  pz = -fDz;
438  }
439  }
440  else if ((select - Sxy) < Sxz)
441  {
442  px = -fDx + 2 * fDx * UUtils::Random();
443  pz = -fDz + 2 * fDz * UUtils::Random();
444 
445  if (UUtils::Random() > 0.5)
446  {
447  py = fDy;
448  }
449  else
450  {
451  py = -fDy;
452  }
453  }
454  else
455  {
456  py = -fDy + 2 * fDy * UUtils::Random();
457  pz = -fDz + 2 * fDz * UUtils::Random();
458 
459  if (UUtils::Random() > 0.5)
460  {
461  px = fDx;
462  }
463  else
464  {
465  px = -fDx;
466  }
467  }
468  return UVector3(px, py, pz);
469 }
470 
471 std::ostream& UBox::StreamInfo(std::ostream& os) const
472 {
473  int oldprc = os.precision(16);
474  os << "-----------------------------------------------------------\n"
475  << " *** Dump for solid - " << GetName() << " ***\n"
476  << " ===================================================\n"
477  << " Solid type: UBox\n"
478  << " Parameters: \n"
479  << " half length X: " << fDx << " mm \n"
480  << " half length Y: " << fDy << " mm \n"
481  << " half length Z: " << fDz << " mm \n"
482  << "-----------------------------------------------------------\n";
483  os.precision(oldprc);
484 
485  return os;
486 }
487 
488 void UBox::SetXHalfLength(double dx)
489 {
490  if (dx > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
491  {
492  fDx = dx;
493  }
494  else
495  {
496  std::ostringstream message;
497  message << "Dimension X too small for solid: " << GetName() << "!"
498  << std::endl
499  << " hX = " << dx;
500  UUtils::Exception("UBox::SetXHalfLength()", "GeomSolids0002",
501  FatalErrorInArguments, 1, message.str().c_str());
502  }
503  fCubicVolume = 0.;
504  fSurfaceArea = 0.;
505 
506 }
507 
508 void UBox::SetYHalfLength(double dy)
509 {
510  if (dy > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
511  {
512  fDy = dy;
513  }
514  else
515  {
516  std::ostringstream message;
517  message << "Dimension Y too small for solid: " << GetName() << "!"
518  << std::endl
519  << " hY = " << dy;
520  UUtils::Exception("UBox::SetYHalfLength()", "GeomSolids0002",
521  FatalErrorInArguments, 1, message.str().c_str());
522  }
523  fCubicVolume = 0.;
524  fSurfaceArea = 0.;
525 
526 }
527 
528 void UBox::SetZHalfLength(double dz)
529 {
530  if (dz > 2 * VUSolid::fgTolerance) // limit to thickness of surfaces
531  {
532  fDz = dz;
533  }
534  else
535  {
536  std::ostringstream message;
537  message << "Dimension Z too small for solid: " << GetName() << "!"
538  << std::endl
539  << " hZ = " << dz;
540  UUtils::Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
541  FatalErrorInArguments, 1, message.str().c_str());
542  }
543  fCubicVolume = 0.;
544  fSurfaceArea = 0.;
545 
546 }
548 {
549  return "Box";
550 }
551 
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Definition: UBox.cc:331
UVector3 GetPointOnSurface() const
Definition: UBox.cc:418
std::ostream & StreamInfo(std::ostream &os) const
Definition: UBox.cc:471
const std::string & GetName() const
Definition: VUSolid.hh:103
static double fgTolerance
Definition: VUSolid.hh:30
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UBox.cc:218
UGeometryType GetEntityType() const
Definition: UBox.cc:547
const XML_Char * name
Definition: expat.h:151
double x
Definition: UVector3.hh:136
Definition: UBox.hh:33
virtual ~UBox()
Definition: UBox.cc:57
UBox()
Definition: UBox.hh:37
if(nlines<=0)
EnumInside
Definition: VUSolid.hh:23
void SetXHalfLength(double dx)
Definition: UBox.cc:488
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:279
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
void Set(double dx, double dy, double dz)
Definition: UBox.cc:43
std::string UGeometryType
Definition: UTypes.hh:70
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UBox.cc:295
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:177
double z
Definition: UVector3.hh:138
void SetYHalfLength(double dy)
Definition: UBox.cc:508
short Sign(short a, short b)
Definition: UUtils.hh:184
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UBox.cc:400
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:528