Geant4  10.01.p02
UOrb.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 // UOrb
12 //
13 // 19.10.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include <cmath>
18 #include <iostream>
19 
20 #include "UOrb.hh"
21 #include "UUtils.hh"
22 
23 using namespace std;
24 
25 //______________________________________________________________________________
26 UOrb::UOrb(const std::string& name, double r)
27  : VUSolid(name), fR(r), fCubicVolume(0), fSurfaceArea(0)
28 {
29  const double epsilon = 2.e-11; // relative tolerance of fR
30 
31  // Check radius
32  //
33  if (r < 10 * VUSolid::fgTolerance) // cartesian tolerance
34  {
35  UUtils::Exception("UOrb::UOrb()", "GeomSolids0002",
36  UFatalErrorInArguments, 1, "Invalid radius > 10*kCarTolerance.");
37  }
38  // VUSolid::fRTolerance is radial tolerance (note: half of G4 tolerance)
39  fRTolerance = max(VUSolid::frTolerance, epsilon * r);
40 }
41 
42 //______________________________________________________________________________
53 // ok
55 {
56  double rad2 = p.x() * p.x() + p.y() * p.y() + p.z() * p.z();
57  // if (false) double rad = sqrt(rad2);
58 
59  double tolRMax = fR - fRTolerance * 0.5;
60 
61  // Check radial surface
62  double tolRMax2 = tolRMax * tolRMax;
63  if (rad2 <= tolRMax2)
64  return eInside;
65  else
66  {
67  tolRMax = fR + fRTolerance * 0.5;
68  tolRMax2 = tolRMax * tolRMax;
69  if (rad2 <= tolRMax2)
70  return eSurface;
71  else
72  return eOutside;
73  }
74 }
75 
76 
77 /*
78 * Computes distance from a point presumably outside the solid to the solid
79 * surface. Ignores first surface if the point is actually inside. Early return
80 * infinity in case the safety to any surface is found greater than the proposed
81 * step aPstep.
82 * The normal vector to the crossed surface is filled only in case the Orb is
83 * crossed, otherwise aNormal.IsNull() is true.
84 */
85 double UOrb::DistanceToIn(const UVector3& p,
86  const UVector3& v,
87  // UVector3 &aNormal,
88  double /*aPstep*/) const
89 {
90  double snxt = UUtils::kInfinity; // snxt = default return value
91 
92  double rad, pDotV3d; // , tolORMax2, tolIRMax2;
93  double c, d2, s = UUtils::kInfinity;
94 
95  const double dRmax = 100.*fR;
96 
97  // General Precalcs
98 
99  rad = sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z());
100  pDotV3d = p.x() * v.x() + p.y() * v.y() + p.z() * v.z();
101 
102  // Radial Precalcs
103 
104  // tolORMax2 = (fR+fRTolerance*0.5)*(fR+fRTolerance*0.5);
105  // tolIRMax2 = (fR-fRTolerance*0.5)*(fR-fRTolerance*0.5);
106 
107  // Outer spherical shell intersection
108  // - Only if outside tolerant fR
109  // - Check for if inside and outer UOrb heading through solid (-> 0)
110  // - No intersect -> no intersection with UOrb
111  //
112  // Shell eqn: x^2+y^2+z^2 = RSPH^2
113  //
114  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
115  //
116  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
117  // => rad2 +2s(pDotV3d) +s^2 =R^2
118  //
119  // => s=-pDotV3d+-sqrt(pDotV3d^2-(rad2-R^2))
120 
121  c = (rad - fR) * (rad + fR); // c = (rad2-R^2))
122 
123  if (rad > fR - fRTolerance * 0.5) // not inside in terms of Inside(p)
124  {
125  if (c > fRTolerance * fR)
126  {
127  // If outside tolerant boundary of outer UOrb in terms of c
128  // [ should be sqrt(rad2) - fR > fRTolerance*0.5 ]
129 
130  d2 = pDotV3d * pDotV3d - c;
131 
132  if (d2 >= 0)
133  {
134  s = -pDotV3d - sqrt(d2); // ok! = [ ( -2 p dot v) +- sqrt [(-2p dot v)2 - 4*(rad - fR)*(rad + fR)] ] / 2
135  // pDotV3d must be positive always, if not use alternative http://en.wikipedia.org/wiki/Quadratic_equation#Alternative_quadratic_formula
136  if (s >= 0)
137  {
138  if (s > dRmax) // Avoid rounding errors due to precision issues seen on
139  {
140  // 64 bits systems. Split long distances and recompute
141  double fTerm = s - fmod(s, dRmax);
142  s = fTerm + DistanceToIn(p + fTerm * v, v);
143  }
144  return snxt = s;
145  }
146  }
147  else // No intersection with UOrb
148  {
149  return snxt = UUtils::kInfinity;
150  }
151  }
152  else // not outside in terms of c
153  {
154  if (c > -fRTolerance * fR) // on surface
155  {
156  d2 = pDotV3d * pDotV3d - c;
157  if ((d2 < fRTolerance * fR) || (pDotV3d >= 0)) // pDotV3d = cos si >= 0
158  {
159  return snxt = UUtils::kInfinity;
160  }
161  else
162  {
163  return snxt = 0.;
164  }
165  }
166  }
167  }
168 #ifdef UDEBUG
169  else // inside ???
170  {
171  UUtils::Exception("UOrb::DistanceToIn(p,v)", "GeomSolids1002", UWarning, 1, "Point p is inside !?");
172  }
173 #endif
174 
175  return snxt;
176 }
177 
179 {
180  double distanceIn = DistanceToIn(p, v);
181  UVector3 shift = distanceIn * v;
182  UVector3 surfacePoint = p + shift;
184  ((UOrb&)*this).Normal(surfacePoint, normal);
185  double dot = normal.Dot(v);
186  if (dot > 0) return 0;
187  else
188  {
189  bool convex;
190  double distanceOut = DistanceToOut(surfacePoint, v, n, convex);
191  return distanceIn + distanceOut;
192  }
193 }
194 
195 /*
196 * Computes distance from a point presumably intside the solid to the solid
197 * surface. Ignores first surface along each axis systematically (for points
198 * inside or outside. Early returns zero in case the second surface is behind
199 * the starting point.
200 * o The proposed step is ignored.
201 * o The normal vector to the crossed surface is always filled.
202 * ______________________________________________________________________________
203 */
204 double UOrb::DistanceToOut(const UVector3& p, const UVector3& v,
205  UVector3& n, bool& convex, double /*aPstep*/) const
206 {
207  double snxt = 0; // snxt: distance to next surface, is default return value
208  bool notOutside = false;
209  convex = true; // orb is always convex, if we leave surface of Orb, we will neber bump on the orb again ...
210 
211  double rad2, pDotV3d;
212  double xi, yi, zi; // Intersection point
213  double c, d2;
214 
215  rad2 = p.x() * p.x() + p.y() * p.y() + p.z() * p.z();
216  pDotV3d = p.x() * v.x() + p.y() * v.y() + p.z() * v.z();
217 
218  // Radial Intersection from UOrb::DistanceToIn
219  //
220  // Outer spherical shell intersection
221  // - Only if outside tolerant fR
222  // - Check for if inside and outer UOrb heading through solid (-> 0)
223  // - No intersect -> no intersection with UOrb
224  //
225  // Shell eqn: x^2+y^2+z^2=RSPH^2
226  //
227  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
228  //
229  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
230  // => rad2 +2s(pDotV3d) +s^2 =R^2
231  //
232  // => s=-pDotV3d+-sqrt(pDotV3d^2-(rad2-R^2))
233 
234  const double rPlus = fR + fRTolerance;
235  double rad = sqrt(rad2);
236 
237  if (rad <= rPlus)
238  {
239  c = (rad - fR) * (rad + fR); // rad2 - fR2
240 
241  if (c < fRTolerance * fR)
242  {
243  // Within tolerant Outer radius
244  //
245  // The test is
246  // rad - fR < 0.5*fRTolerance
247  // => rad < fR + 0.5*kRadTol
248  // => rad2 < (fR + 0.5*kRadTol)^2
249  // => rad2 < fR^2 + 2.*0.5*fR*kRadTol + 0.25*kRadTol*kRadTol
250  // => rad2 - fR^2 <~ fR*kRadTol
251 
252  d2 = pDotV3d * pDotV3d - c;
253 
254  if ((c > -2 * fRTolerance * fR) && // => point is on tolerant surface (i.e. within +- tolerance)
255  ((pDotV3d >= 0) || (d2 < 0))) // if (pDotV3d >= 0 ) => leaving outside from Rmax; i.e. from surface
256  // not re-entering
257  // if (d2 < 0) => it means the point is already outside
258  {
259  // if(calcNorm) // NOTE: we do not have this variable, calcNorm is true always
260  {
261  // *validNorm = true; // NOTE: we do not have this variable, probably always true
262  n = UVector3(p.x() / fR, p.y() / fR, p.z() / fR);
263  }
264  return snxt = 0;
265  }
266  else
267  {
268  // we are inside, with + version of quadratic eq. solution we calculate solution for distance
269  snxt = -pDotV3d + sqrt(d2); // second root since inside Rmax
270  // the solution is safe because pDotV3d is negative
271  // c alternative formula, see http://en.wikipedia.org/wiki/Quadratic_equation#Alternative_quadratic_formula
272  // is not neccessary in this case
273  notOutside = true;
274  }
275  }
276  }
277  else // p is outside ???
278  {
279  // Rule 2: DistanceToOut
280  // Surface points = 0
281  // Outside points: If pointing outwards (dot product with normal positive) return 0, otherwise ignore first surface, another surface should always be on the direction line, use instead distance to this one.
282 
283  // double res = DistanceToOutForOutsidePoints(p, v, n);
284 
285  cout.precision(16);
286  cout << endl;
287  // DumpInfo();
288  cout << "Position:" << endl << endl;
289  cout << "p.x() = " << p.x() << endl;
290  cout << "p.y() = " << p.y() << endl;
291  cout << "p.z() = " << p.z() << endl << endl;
292  cout << "Rp = " << sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()) << endl << endl;
293  cout << "Direction:" << endl << endl;
294  cout << "v.x() = " << v.x() << endl;
295  cout << "v.y() = " << v.y() << endl;
296  cout << "v.z() = " << v.z() << endl << endl;
297  cout << "Proposed distance :" << endl << endl;
298  cout << "snxt = " << snxt << endl << endl;
299 
300 
301 
302  cout.precision(6);
303  UUtils::Exception("UOrb::DistanceToOut(p,v,..)", "GeomSolids1002",
304  UWarning, 1, "Logic error: snxt = kInfinity ???");
305 
306  }
307 
308  // if (calcNorm) // Output switch operator
309  {
310  if (notOutside)
311  {
312  xi = p.x() + snxt * v.x(); // we move to the point on surface, then return normal at that point which for orb, see method bool UOrb::Normal( const UVector3& p, UVector3 &n)
313  yi = p.y() + snxt * v.y();
314  zi = p.z() + snxt * v.z();
315  n = UVector3(xi / fR, yi / fR, zi / fR); // we return normalized vector
316  }
317  else
318  {
319 
320  cout.precision(16);
321  cout << endl;
322  // DumpInfo();
323  cout << "Position:" << endl << endl;
324  cout << "p.x() = " << p.x() << " mm" << endl;
325  cout << "p.y() = " << p.y() << " mm" << endl;
326  cout << "p.z() = " << p.z() << " mm" << endl << endl;
327  cout << "Direction:" << endl << endl;
328  cout << "v.x() = " << v.x() << endl;
329  cout << "v.y() = " << v.y() << endl;
330  cout << "v.z() = " << v.z() << endl << endl;
331  cout << "Proposed distance :" << endl << endl;
332  cout << "snxt = " << snxt << " mm" << endl << endl;
333  cout.precision(6);
334 
335 
336  UUtils::Exception("UOrb::DistanceToOut(p,v,..)", "GeomSolids1002",
337  UWarning, 1, "Undefined side for valid surface normal to solid.");
338  }
339  }
340  return snxt;
341 }
342 
343 /*
344 * Estimates the isotropic safety from a point inside the current solid to any
345 * of its surfaces. The algorithm may be accurate or should provide a fast
346 * underestimate.
347 * ______________________________________________________________________________
348 * Note: In geant4, these methods are DistanceToOut, without given direction
349 * Note: ??? Should not Return 0 anymore if point outside, just the value
350 * OK
351 */
352 double UOrb::SafetyFromInside(const UVector3& p, bool /*aAccurate*/) const
353 {
354 
356  //
357  // Calculate distance (<=actual) to closest surface of shape from inside
358 
359  double safe = 0.0, rad = sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z());
360 
361 #ifdef UDEBUG
362  if (Inside(p) == kOutside)
363  {
364  // int oldprc = cout.precision(16);
365  cout << endl;
366  // DumpInfo();
367  cout << "Position:" << endl << endl;
368  cout << "p.x = " << p.x() << endl;
369  cout << "p.y = " << p.y() << endl;
370  cout << "p.z = " << p.z() << endl << endl;
371  // cout.precision(oldprc);
372  UUtils::Exception("UOrb::DistanceToOut(p)", "GeomSolids1002", UWarning, 1,
373  "Point p is outside !?");
374  }
375 #endif
376 
377  safe = fR - rad;
378  if (safe < 0.) safe = 0.;
379  return safe;
380 }
381 
382 
383 /*
384 * Estimates the isotropic safety from a point outside the current solid to any
385 * of its surfaces. The algorithm may be accurate or should provide a fast
386 * underestimate.
387 * Note: In geant4, this method is equivalent to DistanceToIn, without given direction
388 * ______________________________________________________________________________
389 *
390 * Calculate distance (<= actual) to closest surface of shape from outside
391 * - Calculate distance to radial plane
392 * - Return 0 if point inside
393 * OK
394 */
395 double UOrb::SafetyFromOutside(const UVector3& p, bool /*aAccurate*/) const
396 {
397  double safe = 0.0;
398  double rad = sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z());
399  safe = rad - fR;
400  if (safe < 0)
401  {
402  safe = 0.;
403  }
404  return safe;
405 }
406 
407 
423 bool UOrb::Normal(const UVector3& p, UVector3& n) const
424 {
425  double rad2 = p.x() * p.x() + p.y() * p.y() + p.z() * p.z();
426  double rad = sqrt(rad2);
427 
428  n = UVector3(p.x() / rad, p.y() / rad, p.z() / rad);
429 
430  double tolRMaxP = fR + fRTolerance;
431  double tolRMaxM = fR - fRTolerance;
432 
433  // Check radial surface
434  bool result = ((rad2 <= tolRMaxP * tolRMaxP) && (rad2 >= tolRMaxM * tolRMaxM)); // means we are on surface
435  return result;
436 }
437 
438 
445 void UOrb::Extent(UVector3& aMin, UVector3& aMax) const
446 {
447  aMin.Set(-fR);
448  aMax.Set(fR);
449 }
450 
451 
453 //
454 // Stream object contents to an output stream
455 
456 std::ostream& UOrb::StreamInfo(std::ostream& os) const
457 {
458  int oldprc = os.precision(16);
459  os << "-----------------------------------------------------------\n"
460  << " *** Dump for solid - " << GetName() << " ***\n"
461  << " ===================================================\n"
462  << " Solid type: UOrb\n"
463  << " Parameters: \n"
464 
465  << " outer radius: " << fR << " mm \n"
466  << "-----------------------------------------------------------\n";
467  os.precision(oldprc);
468 
469  return os;
470 }
471 
473 //
474 // GetPointOnSurface
475 
477 {
478  // generate a random number from zero to 2UUtils::kPi...
479  //
480  double phi = UUtils::Random(0., 2.*UUtils::kPi);
481  double cosphi = std::cos(phi);
482  double sinphi = std::sin(phi);
483 
484  // generate a random point uniform in area
485  double costheta = UUtils::Random(-1., 1.);
486  double sintheta = std::sqrt(1. - UUtils::sqr(costheta));
487 
488  return UVector3(fR * sintheta * cosphi, fR * sintheta * sinphi, fR * costheta);
489 }
490 
492 {
493  return new UOrb(*this);
494 }
495 
496 // Copy constructor
497 
498 UOrb::UOrb(const UOrb& rhs)
499  : VUSolid(rhs), fR(rhs.fR), fRTolerance(rhs.fRTolerance), fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea)
500 {
501 }
502 
504 //
505 // Assignment operator
506 
508 {
509  // Check assignment to self
510  //
511  if (this == &rhs)
512  {
513  return *this;
514  }
515 
516  // Copy base class data
517  //
518  VUSolid::operator=(rhs);
519 
520  // Copy data
521  //
522  fR = rhs.fR;
523  fRTolerance = rhs.fRTolerance;
526  return *this;
527 }
529 //
530 // Get Parameters List for visualisation
531 
532 void UOrb::GetParametersList(int, double* aArray)const
533 {
534  aArray[0] = GetRadius();
535 }
537 //
538 // Get Entity Type
539 
541 {
542  return "Orb";
543 }
UOrb & operator=(const UOrb &rhs)
Definition: UOrb.cc:507
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Return unit normal of surface closest to p.
Definition: UOrb.cc:423
double & z()
Definition: UVector3.hh:352
double & y()
Definition: UVector3.hh:348
static double frTolerance
Definition: VUSolid.hh:31
UOrb()
Definition: UOrb.hh:32
const std::string & GetName() const
Definition: VUSolid.hh:103
VUSolid * Clone() const
Definition: UOrb.cc:491
static double fgTolerance
Definition: VUSolid.hh:30
G4String name
Definition: TRTMaterials.hh:40
std::ostream & StreamInfo(std::ostream &os) const
Definition: UOrb.cc:456
UGeometryType GetEntityType() const
Definition: UOrb.cc:540
double GetRadius() const
Definition: UOrb.hh:93
double DistanceToOutForOutsidePoints(const UVector3 &p, const UVector3 &v, UVector3 &n) const
Definition: UOrb.cc:178
UVector3 GetPointOnSurface() const
Definition: UOrb.cc:476
double fCubicVolume
Definition: UOrb.hh:86
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
static const double s
Definition: G4SIunits.hh:150
static const double kInfinity
Definition: UUtils.hh:54
Definition: UOrb.hh:28
double & x()
Definition: UVector3.hh:344
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UOrb.cc:352
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UOrb.cc:85
void GetParametersList(int, double *) const
Definition: UOrb.cc:532
EnumInside
Definition: VUSolid.hh:23
const G4int n
double Dot(const UVector3 &) const
Definition: UVector3.hh:276
void Extent(UVector3 &aMin, UVector3 &aMax) const
Returns extent of the solid along a given cartesian axis OK.
Definition: UOrb.cc:445
static const double rad
Definition: G4SIunits.hh:130
double fSurfaceArea
Definition: UOrb.hh:87
T sqr(const T &x)
Definition: UUtils.hh:138
double fR
Definition: UOrb.hh:84
double fRTolerance
Definition: UOrb.hh:85
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
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
static const double kPi
Definition: UUtils.hh:49
std::string UGeometryType
Definition: UTypes.hh:39
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UOrb.cc:395
static const G4double d2
EnumInside Inside(const UVector3 &aPo6int) const
Return whether point inside/outside/on surface Split into radius checks.
Definition: UOrb.cc:54
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UOrb.cc:204