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