Geant4_10
USphere.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 // USphere
12 //
13 // 19.10.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "USphere.hh"
18 
19 using namespace std;
20 
21 // Private enum: Not for external use - used by distanceToOut
22 
23 enum ESide {kNull, kRMin, kRMax, kSPhi, kEPhi, kSTheta, kETheta};
24 
25 // used by normal
26 
27 enum ENorm {kNRMin, kNRMax, kNSPhi, kNEPhi, kNSTheta, kNETheta};
28 
30 //
31 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
32 // - note if pDPhi>2PI then reset to 2PI
33 
34 USphere::USphere(const std::string& pName,
35  double pRmin, double pRmax,
36  double pSPhi, double pDPhi,
37  double pSTheta, double pDTheta)
38  : VUSolid(pName), fCubicVolume(0.),
39  fSurfaceArea(0.),fEpsilon(2.e-11),
40  fFullPhiSphere(true), fFullThetaSphere(true)
41 {
42  kAngTolerance = faTolerance;
43 
44  // Check radii and Set radial tolerances
45 
46  kRadTolerance = frTolerance;
47  if ((pRmin >= pRmax) || (pRmax < 1.1 * kRadTolerance) || (pRmin < 0))
48  {
49  std::ostringstream message;
50  message << "Invalid radii for Solid: " << GetName() << std::endl
51  << " pRmin = " << pRmin << ", pRmax = " << pRmax;
52  UUtils::Exception("USphere::USphere()", "GeomSolids0002",
53  FatalErrorInArguments, 1, message.str().c_str());
54  }
55  fRmin = pRmin;
56  fRmax = pRmax;
57  fRminTolerance = (fRmin) ? std::max(kRadTolerance, fEpsilon * fRmin) : 0;
58  kTolerance = std::max(kRadTolerance, fEpsilon * fRmax);
59 
60  // Check angles
61 
62  CheckPhiAngles(pSPhi, pDPhi);
63  CheckThetaAngles(pSTheta, pDTheta);
64 }
65 
67 //
68 // Destructor
69 
70 USphere::~USphere()
71 {
72 }
73 
75 //
76 // Copy constructor
77 
78 USphere::USphere(const USphere& rhs)
79  : VUSolid(rhs), fCubicVolume(0.),fSurfaceArea(0.),
80  fRminTolerance(rhs.fRminTolerance),
81  kTolerance(rhs.kTolerance), kAngTolerance(rhs.kAngTolerance),
82  kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
83  fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
84  fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
85  sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
86  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
87  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
88  sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
89  hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
90  sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
91  sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
92  tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
93  tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
94  fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
95  fFullSphere(rhs.fFullSphere)
96 {
97 }
98 
100 //
101 // Assignment operator
102 
103 USphere& USphere::operator = (const USphere& rhs)
104 {
105  // Check assignment to self
106  //
107  if (this == &rhs)
108  {
109  return *this;
110  }
111 
112  // Copy base class data
113  //
114  VUSolid::operator=(rhs);
115 
116  // Copy data
117  //
118  fCubicVolume = rhs.fCubicVolume;
119  fSurfaceArea = rhs.fSurfaceArea;
120  fRminTolerance = rhs.fRminTolerance;
121  kTolerance = rhs.kTolerance;
122  kAngTolerance = rhs.kAngTolerance;
123  kRadTolerance = rhs.kRadTolerance;
124  fEpsilon = rhs.fEpsilon;
125  fRmin = rhs.fRmin;
126  fRmax = rhs.fRmax;
127  fSPhi = rhs.fSPhi;
128  fDPhi = rhs.fDPhi;
129  fSTheta = rhs.fSTheta;
130  fDTheta = rhs.fDTheta;
131  sinCPhi = rhs.sinCPhi;
132  cosCPhi = rhs.cosCPhi;
133  cosHDPhiOT = rhs.cosHDPhiOT;
134  cosHDPhiIT = rhs.cosHDPhiIT;
135  sinSPhi = rhs.sinSPhi;
136  cosSPhi = rhs.cosSPhi;
137  sinEPhi = rhs.sinEPhi;
138  cosEPhi = rhs.cosEPhi;
139  hDPhi = rhs.hDPhi;
140  cPhi = rhs.cPhi;
141  ePhi = rhs.ePhi;
142  sinSTheta = rhs.sinSTheta;
143  cosSTheta = rhs.cosSTheta;
144  sinETheta = rhs.sinETheta;
145  cosETheta = rhs.cosETheta;
146  tanSTheta = rhs.tanSTheta;
147  tanSTheta2 = rhs.tanSTheta2;
148  tanETheta = rhs.tanETheta;
149  tanETheta2 = rhs.tanETheta2;
150  eTheta = rhs.eTheta;
151  fFullPhiSphere = rhs.fFullPhiSphere;
152  fFullThetaSphere = rhs.fFullThetaSphere;
153  fFullSphere = rhs.fFullSphere;
154 
155  return *this;
156 }
157 
159 //
160 // Return whether point inside/outside/on surface
161 // Split into radius, phi, theta checks
162 // Each check modifies 'in', or returns as approprate
163 
164 VUSolid::EnumInside USphere::Inside(const UVector3& p) const
165 {
166  double rho, rho2, rad2, tolRMin, tolRMax;
167  double pPhi, pTheta;
168  VUSolid::EnumInside in = eOutside;
169  static const double halfAngTolerance = kAngTolerance * 0.5;
170  const double halfTolerance = kTolerance * 0.5;
171  const double halfRminTolerance = fRminTolerance * 0.5;
172  const double rMaxMinus = fRmax - halfTolerance;
173  const double rMinPlus = (fRmin > 0) ? fRmin + halfRminTolerance : 0;
174 
175  rho2 = p.x * p.x + p.y * p.y;
176  rad2 = rho2 + p.z * p.z;
177 
178  // Check radial surfaces. Sets 'in'
179 
180  tolRMin = rMinPlus;
181  tolRMax = rMaxMinus;
182 
183  if ((rad2 <= rMaxMinus * rMaxMinus) && (rad2 >= rMinPlus * rMinPlus))
184  {
185  in = eInside;
186  }
187  else
188  {
189  tolRMax = fRmax + halfTolerance; // outside case
190  tolRMin = std::max(fRmin - halfRminTolerance, 0.); // outside case
191  if ((rad2 <= tolRMax * tolRMax) && (rad2 >= tolRMin * tolRMin))
192  {
193  in = eSurface;
194  }
195  else
196  {
197  return in = eOutside;
198  }
199  }
200 
201  // Phi boundaries : Do not check if it has no phi boundary!
202 
203  if (!fFullPhiSphere && rho2) // [fDPhi < 2*UUtils::kPi] and [p.x or p.y]
204  {
205  pPhi = std::atan2(p.y, p.x);
206 
207  if (pPhi < fSPhi - halfAngTolerance)
208  {
209  pPhi += 2 * UUtils::kPi;
210  }
211  else if (pPhi > ePhi + halfAngTolerance)
212  {
213  pPhi -= 2 * UUtils::kPi;
214  }
215 
216  if ((pPhi < fSPhi - halfAngTolerance)
217  || (pPhi > ePhi + halfAngTolerance))
218  {
219  return in = eOutside;
220  }
221 
222  else if (in == eInside) // else it's eSurface anyway already
223  {
224  if ((pPhi < fSPhi + halfAngTolerance)
225  || (pPhi > ePhi - halfAngTolerance))
226  {
227  in = eSurface;
228  }
229  }
230  }
231 
232  // Theta bondaries
233 
234  if ((rho2 || p.z) && (!fFullThetaSphere))
235  {
236  rho = std::sqrt(rho2);
237  pTheta = std::atan2(rho, p.z);
238 
239  if (in == eInside)
240  {
241  if ((pTheta < fSTheta + halfAngTolerance)
242  || (pTheta > eTheta - halfAngTolerance))
243  {
244  if ((pTheta >= fSTheta - halfAngTolerance)
245  && (pTheta <= eTheta + halfAngTolerance))
246  {
247  in = eSurface;
248  }
249  else
250  {
251  in = eOutside;
252  }
253  }
254  }
255  else
256  {
257  if ((pTheta < fSTheta - halfAngTolerance)
258  || (pTheta > eTheta + halfAngTolerance))
259  {
260  in = eOutside;
261  }
262  }
263  }
264  return in;
265 }
266 
268 //
269 // Return Unit normal of surface closest to p
270 // - note if point on z axis, ignore phi divided sides
271 // - unsafe if point close to z axis a rmin=0 - no explicit checks
272 
273 bool USphere::Normal(const UVector3& p, UVector3& n) const
274 {
275  int noSurfaces = 0;
276  double rho, rho2, radius, pTheta, pPhi = 0.;
277  double distRMin = UUtils::Infinity();
278  double distSPhi = UUtils::Infinity(), distEPhi = UUtils::Infinity();
279  double distSTheta = UUtils::Infinity(), distETheta = UUtils::Infinity();
280  UVector3 nR, nPs, nPe, nTs, nTe, nZ(0., 0., 1.);
281  UVector3 norm, sumnorm(0., 0., 0.);
282 
283  static const double halfCarTolerance = 0.5 * VUSolid::Tolerance();
284  static const double halfAngTolerance = 0.5 * kAngTolerance;
285 
286  rho2 = p.x * p.x + p.y * p.y;
287  radius = std::sqrt(rho2 + p.z * p.z);
288  rho = std::sqrt(rho2);
289 
290  double distRMax = std::fabs(radius - fRmax);
291  if (fRmin) distRMin = std::fabs(radius - fRmin);
292 
293  if (rho && !fFullSphere)
294  {
295  pPhi = std::atan2(p.y, p.x);
296 
297  if (pPhi < fSPhi - halfAngTolerance)
298  {
299  pPhi += 2 * UUtils::kPi;
300  }
301  else if (pPhi > ePhi + halfAngTolerance)
302  {
303  pPhi -= 2 * UUtils::kPi;
304  }
305  }
306  if (!fFullPhiSphere)
307  {
308  if (rho)
309  {
310  distSPhi = std::fabs(pPhi - fSPhi);
311  distEPhi = std::fabs(pPhi - ePhi);
312  }
313  else if (!fRmin)
314  {
315  distSPhi = 0.;
316  distEPhi = 0.;
317  }
318  nPs = UVector3(sinSPhi, -cosSPhi, 0);
319  nPe = UVector3(-sinEPhi, cosEPhi, 0);
320  }
321  if (!fFullThetaSphere)
322  {
323  if (rho)
324  {
325  pTheta = std::atan2(rho, p.z);
326  distSTheta = std::fabs(pTheta - fSTheta);
327  distETheta = std::fabs(pTheta - eTheta);
328 
329  nTs = UVector3(-cosSTheta * p.x / rho,
330  -cosSTheta * p.y / rho,
331  sinSTheta);
332 
333  nTe = UVector3(cosETheta * p.x / rho,
334  cosETheta * p.y / rho,
335  -sinETheta);
336  }
337  else if (!fRmin)
338  {
339  if (fSTheta)
340  {
341  distSTheta = 0.;
342  nTs = UVector3(0., 0., -1.);
343  }
344  if (eTheta < UUtils::kPi)
345  {
346  distETheta = 0.;
347  nTe = UVector3(0., 0., 1.);
348  }
349  }
350  }
351  if (radius)
352  {
353  nR = UVector3(p.x / radius, p.y / radius, p.z / radius);
354  }
355 
356  if (distRMax <= halfCarTolerance)
357  {
358  noSurfaces ++;
359  sumnorm += nR;
360  }
361  if (fRmin && (distRMin <= halfCarTolerance))
362  {
363  noSurfaces ++;
364  sumnorm -= nR;
365  }
366  if (!fFullPhiSphere)
367  {
368  if (distSPhi <= halfAngTolerance)
369  {
370  noSurfaces ++;
371  sumnorm += nPs;
372  }
373  if (distEPhi <= halfAngTolerance)
374  {
375  noSurfaces ++;
376  sumnorm += nPe;
377  }
378  }
379  if (!fFullThetaSphere)
380  {
381  if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
382  {
383  noSurfaces ++;
384  if ((radius <= halfCarTolerance) && fFullPhiSphere)
385  {
386  sumnorm += nZ;
387  }
388  else
389  {
390  sumnorm += nTs;
391  }
392  }
393  if ((distETheta <= halfAngTolerance) && (eTheta < UUtils::kPi))
394  {
395  noSurfaces ++;
396  if ((radius <= halfCarTolerance) && fFullPhiSphere)
397  {
398  sumnorm -= nZ;
399  }
400  else
401  {
402  sumnorm += nTe;
403  }
404  if (sumnorm.z == 0.)
405  {
406  sumnorm += nZ;
407  }
408  }
409  }
410  if (noSurfaces == 0)
411  {
412 #ifdef UDEBUG
413  UUtils::Exception("USphere::SurfaceNormal(p)", "GeomSolids1002",
414  Warning, 1, "Point p is not on surface !?");
415 #endif
416  norm = ApproxSurfaceNormal(p);
417  }
418  else if (noSurfaces == 1)
419  {
420  norm = sumnorm;
421  }
422  else
423  {
424  norm = sumnorm.Unit();
425  }
426  n = norm;
427  return (noSurfaces > 0);
428 }
429 
430 
432 //
433 // Algorithm for SurfaceNormal() following the original specification
434 // for points not on the surface
435 
436 UVector3 USphere::ApproxSurfaceNormal(const UVector3& p) const
437 {
438  ENorm side;
439  UVector3 norm;
440  double rho, rho2, radius, pPhi, pTheta;
441  double distRMin, distRMax, distSPhi, distEPhi,
442  distSTheta, distETheta, distMin;
443 
444  rho2 = p.x * p.x + p.y * p.y;
445  radius = std::sqrt(rho2 + p.z * p.z);
446  rho = std::sqrt(rho2);
447 
448  //
449  // Distance to r shells
450  //
451 
452  distRMax = std::fabs(radius - fRmax);
453  if (fRmin)
454  {
455  distRMin = std::fabs(radius - fRmin);
456 
457  if (distRMin < distRMax)
458  {
459  distMin = distRMin;
460  side = kNRMin;
461  }
462  else
463  {
464  distMin = distRMax;
465  side = kNRMax;
466  }
467  }
468  else
469  {
470  distMin = distRMax;
471  side = kNRMax;
472  }
473  //
474  // Distance to phi planes
475  //
476  // Protected against (0,0,z)
477 
478  pPhi = std::atan2(p.y, p.x);
479  if (pPhi < 0)
480  {
481  pPhi += 2 * UUtils::kPi;
482  }
483 
484  if (!fFullPhiSphere && rho)
485  {
486  if (fSPhi < 0)
487  {
488  distSPhi = std::fabs(pPhi - (fSPhi + 2 * UUtils::kPi)) * rho;
489  }
490  else
491  {
492  distSPhi = std::fabs(pPhi - fSPhi) * rho;
493  }
494 
495  distEPhi = std::fabs(pPhi - fSPhi - fDPhi) * rho;
496 
497  // Find new minimum
498  //
499  if (distSPhi < distEPhi)
500  {
501  if (distSPhi < distMin)
502  {
503  distMin = distSPhi;
504  side = kNSPhi;
505  }
506  }
507  else
508  {
509  if (distEPhi < distMin)
510  {
511  distMin = distEPhi;
512  side = kNEPhi;
513  }
514  }
515  }
516 
517  //
518  // Distance to theta planes
519  //
520 
521  if (!fFullThetaSphere && radius)
522  {
523  pTheta = std::atan2(rho, p.z);
524  distSTheta = std::fabs(pTheta - fSTheta) * radius;
525  distETheta = std::fabs(pTheta - fSTheta - fDTheta) * radius;
526 
527  // Find new minimum
528  //
529  if (distSTheta < distETheta)
530  {
531  if (distSTheta < distMin)
532  {
533  distMin = distSTheta;
534  side = kNSTheta;
535  }
536  }
537  else
538  {
539  if (distETheta < distMin)
540  {
541  distMin = distETheta;
542  side = kNETheta;
543  }
544  }
545  }
546 
547  switch (side)
548  {
549  case kNRMin: // Inner radius
550  norm = UVector3(-p.x / radius, -p.y / radius, -p.z / radius);
551  break;
552  case kNRMax: // Outer radius
553  norm = UVector3(p.x / radius, p.y / radius, p.z / radius);
554  break;
555  case kNSPhi:
556  norm = UVector3(sinSPhi, -cosSPhi, 0);
557  break;
558  case kNEPhi:
559  norm = UVector3(-sinEPhi, cosEPhi, 0);
560  break;
561  case kNSTheta:
562  norm = UVector3(-cosSTheta * std::cos(pPhi),
563  -cosSTheta * std::sin(pPhi),
564  sinSTheta);
565  break;
566  case kNETheta:
567  norm = UVector3(cosETheta * std::cos(pPhi),
568  cosETheta * std::sin(pPhi),
569  sinETheta);
570  break;
571  default: // Should never reach this case ...
572 
573  UUtils::Exception("USphere::ApproxSurfaceNormal()",
574  "GeomSolids1002", Warning, 1,
575  "Undefined side for valid surface normal to solid.");
576  break;
577  }
578 
579  return norm;
580 }
581 
583 //
584 // Calculate distance to shape from outside, along normalised vector
585 // - return UUtils::Infinity() if no intersection, or intersection distance <= tolerance
586 //
587 // -> If point is outside outer radius, compute intersection with rmax
588 // - if no intersection return
589 // - if valid phi,theta return intersection Dist
590 //
591 // -> If shell, compute intersection with inner radius, taking largest +ve root
592 // - if valid phi,theta, save intersection
593 //
594 // -> If phi segmented, compute intersection with phi half planes
595 // - if valid intersection(r,theta), return smallest intersection of
596 // inner shell & phi intersection
597 //
598 // -> If theta segmented, compute intersection with theta cones
599 // - if valid intersection(r,phi), return smallest intersection of
600 // inner shell & theta intersection
601 //
602 //
603 // NOTE:
604 // - `if valid' (above) implies tolerant checking of intersection points
605 //
606 // OPT:
607 // Move tolIO/ORmin/RMax2 precalcs to where they are needed -
608 // not required for most cases.
609 // Avoid atan2 for non theta cut USphere.
610 
611 double USphere::DistanceToIn(const UVector3& p, const UVector3& v, double /*aPstep*/) const
612 {
613  double snxt = UUtils::Infinity(); // snxt = default return value
614  double rho2, rad2, pDotV2d, pDotV3d, pTheta;
615  double tolSTheta = 0., tolETheta = 0.;
616  const double dRmax = 100.*fRmax;
617 
618  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
619  static const double halfAngTolerance = kAngTolerance * 0.5;
620  const double halfTolerance = kTolerance * 0.5;
621  const double halfRminTolerance = fRminTolerance * 0.5;
622  const double tolORMin2 = (fRmin > halfRminTolerance)
623  ? (fRmin - halfRminTolerance) * (fRmin - halfRminTolerance) : 0;
624  const double tolIRMin2 =
625  (fRmin + halfRminTolerance) * (fRmin + halfRminTolerance);
626  const double tolORMax2 =
627  (fRmax + halfTolerance) * (fRmax + halfTolerance);
628  const double tolIRMax2 =
629  (fRmax - halfTolerance) * (fRmax - halfTolerance);
630 
631  // Intersection point
632  //
633  double xi, yi, zi, rhoi, rhoi2, radi2, iTheta;
634 
635  // Phi intersection
636  //
637  double Comp;
638 
639  // Phi precalcs
640  //
641  double Dist, cosPsi;
642 
643  // Theta precalcs
644  //
645  double dist2STheta, dist2ETheta;
646  double t1, t2, b, c, d2, d, sd = UUtils::Infinity();
647 
648  // General Precalcs
649  //
650  rho2 = p.x * p.x + p.y * p.y;
651  rad2 = rho2 + p.z * p.z;
652  pTheta = std::atan2(std::sqrt(rho2), p.z);
653 
654  pDotV2d = p.x * v.x + p.y * v.y;
655  pDotV3d = pDotV2d + p.z * v.z;
656 
657  // Theta precalcs
658  //
659  if (!fFullThetaSphere)
660  {
661  tolSTheta = fSTheta - halfAngTolerance;
662  tolETheta = eTheta + halfAngTolerance;
663  }
664 
665  // Outer spherical shell intersection
666  // - Only if outside tolerant fRmax
667  // - Check for if inside and outer USphere heading through solid (-> 0)
668  // - No intersect -> no intersection with USphere
669  //
670  // Shell eqn: x^2+y^2+z^2=RSPH^2
671  //
672  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
673  //
674  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
675  // => rad2 +2sd(pDotV3d) +sd^2 =R^2
676  //
677  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
678 
679  c = rad2 - fRmax * fRmax;
680 
681  if (c > kTolerance * fRmax)
682  {
683  // If outside tolerant boundary of outer USphere
684  // [should be std::sqrt(rad2)-fRmax > halfTolerance]
685 
686  d2 = pDotV3d * pDotV3d - c;
687 
688  if (d2 >= 0)
689  {
690  sd = -pDotV3d - std::sqrt(d2);
691 
692  if (sd >= 0)
693  {
694  if (sd > dRmax) // Avoid rounding errors due to precision issues seen on
695  {
696  // 64 bits systems. Split long distances and recompute
697  double fTerm = sd - std::fmod(sd, dRmax);
698  sd = fTerm + DistanceToIn(p + fTerm * v, v);
699  }
700  xi = p.x + sd * v.x;
701  yi = p.y + sd * v.y;
702  rhoi = std::sqrt(xi * xi + yi * yi);
703 
704  if (!fFullPhiSphere && rhoi) // Check phi intersection
705  {
706  cosPsi = (xi * cosCPhi + yi * sinCPhi) / rhoi;
707 
708  if (cosPsi >= cosHDPhiOT)
709  {
710  if (!fFullThetaSphere) // Check theta intersection
711  {
712  zi = p.z + sd * v.z;
713 
714  // rhoi & zi can never both be 0
715  // (=>intersect at origin =>fRmax=0)
716  //
717  iTheta = std::atan2(rhoi, zi);
718  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
719  {
720  return snxt = sd;
721  }
722  }
723  else
724  {
725  return snxt = sd;
726  }
727  }
728  }
729  else
730  {
731  if (!fFullThetaSphere) // Check theta intersection
732  {
733  zi = p.z + sd * v.z;
734 
735  // rhoi & zi can never both be 0
736  // (=>intersect at origin => fRmax=0 !)
737  //
738  iTheta = std::atan2(rhoi, zi);
739  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
740  {
741  return snxt = sd;
742  }
743  }
744  else
745  {
746  return snxt = sd;
747  }
748  }
749  }
750  }
751  else // No intersection with USphere
752  {
753  return snxt = UUtils::Infinity();
754  }
755  }
756  else
757  {
758  // Inside outer radius
759  // check not inside, and heading through USphere (-> 0 to in)
760 
761  d2 = pDotV3d * pDotV3d - c;
762 
763  if ((rad2 > tolIRMax2)
764  && ((d2 >= kTolerance * fRmax) && (pDotV3d < 0)))
765  {
766  if (!fFullPhiSphere)
767  {
768  // Use inner phi tolerant boundary -> if on tolerant
769  // phi boundaries, phi intersect code handles leaving/entering checks
770 
771  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
772 
773  if (cosPsi >= cosHDPhiIT)
774  {
775  // inside radii, delta r -ve, inside phi
776 
777  if (!fFullThetaSphere)
778  {
779  if ((pTheta >= tolSTheta + kAngTolerance)
780  && (pTheta <= tolETheta - kAngTolerance))
781  {
782  return snxt = 0;
783  }
784  }
785  else // strictly inside Theta in both cases
786  {
787  return snxt = 0;
788  }
789  }
790  }
791  else
792  {
793  if (!fFullThetaSphere)
794  {
795  if ((pTheta >= tolSTheta + kAngTolerance)
796  && (pTheta <= tolETheta - kAngTolerance))
797  {
798  return snxt = 0;
799  }
800  }
801  else // strictly inside Theta in both cases
802  {
803  return snxt = 0;
804  }
805  }
806  }
807  }
808 
809  // Inner spherical shell intersection
810  // - Always farthest root, because would have passed through outer
811  // surface first.
812  // - Tolerant check if travelling through solid
813 
814  if (fRmin)
815  {
816  c = rad2 - fRmin * fRmin;
817  d2 = pDotV3d * pDotV3d - c;
818 
819  // Within tolerance inner radius of inner USphere
820  // Check for immediate entry/already inside and travelling outwards
821 
822  if ((c > -halfRminTolerance) && (rad2 < tolIRMin2)
823  && ((d2 < fRmin * VUSolid::Tolerance()) || (pDotV3d >= 0)))
824  {
825  if (!fFullPhiSphere)
826  {
827  // Use inner phi tolerant boundary -> if on tolerant
828  // phi boundaries, phi intersect code handles leaving/entering checks
829 
830  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
831  if (cosPsi >= cosHDPhiIT)
832  {
833  // inside radii, delta r -ve, inside phi
834  //
835  if (!fFullThetaSphere)
836  {
837  if ((pTheta >= tolSTheta + kAngTolerance)
838  && (pTheta <= tolETheta - kAngTolerance))
839  {
840  return snxt = 0;
841  }
842  }
843  else
844  {
845  return snxt = 0;
846  }
847  }
848  }
849  else
850  {
851  if (!fFullThetaSphere)
852  {
853  if ((pTheta >= tolSTheta + kAngTolerance)
854  && (pTheta <= tolETheta - kAngTolerance))
855  {
856  return snxt = 0;
857  }
858  }
859  else
860  {
861  return snxt = 0;
862  }
863  }
864  }
865  else // Not special tolerant case
866  {
867  if (d2 >= 0)
868  {
869  sd = -pDotV3d + std::sqrt(d2);
870  if (sd >= halfRminTolerance) // It was >= 0 ??
871  {
872  xi = p.x + sd * v.x;
873  yi = p.y + sd * v.y;
874  rhoi = std::sqrt(xi * xi + yi * yi);
875 
876  if (!fFullPhiSphere && rhoi) // Check phi intersection
877  {
878  cosPsi = (xi * cosCPhi + yi * sinCPhi) / rhoi;
879 
880  if (cosPsi >= cosHDPhiOT)
881  {
882  if (!fFullThetaSphere) // Check theta intersection
883  {
884  zi = p.z + sd * v.z;
885 
886  // rhoi & zi can never both be 0
887  // (=>intersect at origin =>fRmax=0)
888  //
889  iTheta = std::atan2(rhoi, zi);
890  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
891  {
892  snxt = sd;
893  }
894  }
895  else
896  {
897  snxt = sd;
898  }
899  }
900  }
901  else
902  {
903  if (!fFullThetaSphere) // Check theta intersection
904  {
905  zi = p.z + sd * v.z;
906 
907  // rhoi & zi can never both be 0
908  // (=>intersect at origin => fRmax=0 !)
909  //
910  iTheta = std::atan2(rhoi, zi);
911  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
912  {
913  snxt = sd;
914  }
915  }
916  else
917  {
918  snxt = sd;
919  }
920  }
921  }
922  }
923  }
924  }
925 
926  // Phi segment intersection
927  //
928  // o Tolerant of points inside phi planes by up to VUSolid::Tolerance()*0.5
929  //
930  // o NOTE: Large duplication of code between sphi & ephi checks
931  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
932  // intersection check <=0 -> >=0
933  // -> Should use some form of loop Construct
934  //
935  if (!fFullPhiSphere)
936  {
937  // First phi surface ('S'tarting phi)
938  // Comp = Component in outwards normal dirn
939  //
940  Comp = v.x * sinSPhi - v.y * cosSPhi;
941 
942  if (Comp < 0)
943  {
944  Dist = p.y * cosSPhi - p.x * sinSPhi;
945 
946  if (Dist < halfCarTolerance)
947  {
948  sd = Dist / Comp;
949 
950  if (sd < snxt)
951  {
952  if (sd > 0)
953  {
954  xi = p.x + sd * v.x;
955  yi = p.y + sd * v.y;
956  zi = p.z + sd * v.z;
957  rhoi2 = xi * xi + yi * yi ;
958  radi2 = rhoi2 + zi * zi ;
959  }
960  else
961  {
962  sd = 0;
963  xi = p.x;
964  yi = p.y;
965  zi = p.z;
966  rhoi2 = rho2;
967  radi2 = rad2;
968  }
969  if ((radi2 <= tolORMax2)
970  && (radi2 >= tolORMin2)
971  && ((yi * cosCPhi - xi * sinCPhi) <= 0))
972  {
973  // Check theta intersection
974  // rhoi & zi can never both be 0
975  // (=>intersect at origin =>fRmax=0)
976  //
977  if (!fFullThetaSphere)
978  {
979  iTheta = std::atan2(std::sqrt(rhoi2), zi);
980  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
981  {
982  // r and theta intersections good
983  // - check intersecting with correct half-plane
984 
985  if ((yi * cosCPhi - xi * sinCPhi) <= 0)
986  {
987  snxt = sd;
988  }
989  }
990  }
991  else
992  {
993  snxt = sd;
994  }
995  }
996  }
997  }
998  }
999 
1000  // Second phi surface ('E'nding phi)
1001  // Component in outwards normal dirn
1002 
1003  Comp = -(v.x * sinEPhi - v.y * cosEPhi);
1004 
1005  if (Comp < 0)
1006  {
1007  Dist = -(p.y * cosEPhi - p.x * sinEPhi);
1008  if (Dist < halfCarTolerance)
1009  {
1010  sd = Dist / Comp;
1011 
1012  if (sd < snxt)
1013  {
1014  if (sd > 0)
1015  {
1016  xi = p.x + sd * v.x;
1017  yi = p.y + sd * v.y;
1018  zi = p.z + sd * v.z;
1019  rhoi2 = xi * xi + yi * yi ;
1020  radi2 = rhoi2 + zi * zi ;
1021  }
1022  else
1023  {
1024  sd = 0 ;
1025  xi = p.x;
1026  yi = p.y;
1027  zi = p.z;
1028  rhoi2 = rho2 ;
1029  radi2 = rad2 ;
1030  }
1031  if ((radi2 <= tolORMax2)
1032  && (radi2 >= tolORMin2)
1033  && ((yi * cosCPhi - xi * sinCPhi) >= 0))
1034  {
1035  // Check theta intersection
1036  // rhoi & zi can never both be 0
1037  // (=>intersect at origin =>fRmax=0)
1038  //
1039  if (!fFullThetaSphere)
1040  {
1041  iTheta = std::atan2(std::sqrt(rhoi2), zi);
1042  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
1043  {
1044  // r and theta intersections good
1045  // - check intersecting with correct half-plane
1046 
1047  if ((yi * cosCPhi - xi * sinCPhi) >= 0)
1048  {
1049  snxt = sd;
1050  }
1051  }
1052  }
1053  else
1054  {
1055  snxt = sd;
1056  }
1057  }
1058  }
1059  }
1060  }
1061  }
1062 
1063  // Theta segment intersection
1064 
1065  if (!fFullThetaSphere)
1066  {
1067 
1068  // Intersection with theta surfaces
1069  // Known failure cases:
1070  // o Inside tolerance of stheta surface, skim
1071  // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1072  //
1073  // To solve: Check 2nd root of etheta surface in addition to stheta
1074  //
1075  // o start/end theta is exactly UUtils::kPi/2
1076  // Intersections with cones
1077  //
1078  // Cone equation: x^2+y^2=z^2tan^2(t)
1079  //
1080  // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1081  //
1082  // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1083  // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1084  //
1085  // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1086 
1087  if (fSTheta)
1088  {
1089  dist2STheta = rho2 - p.z * p.z * tanSTheta2;
1090  }
1091  else
1092  {
1093  dist2STheta = UUtils::Infinity();
1094  }
1095  if (eTheta < UUtils::kPi)
1096  {
1097  dist2ETheta = rho2 - p.z * p.z * tanETheta2;
1098  }
1099  else
1100  {
1101  dist2ETheta = UUtils::Infinity();
1102  }
1103  if (pTheta < tolSTheta)
1104  {
1105  // Inside (theta<stheta-tol) stheta cone
1106  // First root of stheta cone, second if first root -ve
1107 
1108  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1109  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1110  if (t1)
1111  {
1112  b = t2 / t1;
1113  c = dist2STheta / t1;
1114  d2 = b * b - c;
1115 
1116  if (d2 >= 0)
1117  {
1118  d = std::sqrt(d2);
1119  sd = -b - d; // First root
1120  zi = p.z + sd * v.z;
1121 
1122  if ((sd < 0) || (zi * (fSTheta - UUtils::kPi / 2) > 0))
1123  {
1124  sd = -b + d; // Second root
1125  }
1126  if ((sd >= 0) && (sd < snxt))
1127  {
1128  xi = p.x + sd * v.x;
1129  yi = p.y + sd * v.y;
1130  zi = p.z + sd * v.z;
1131  rhoi2 = xi * xi + yi * yi;
1132  radi2 = rhoi2 + zi * zi;
1133  if ((radi2 <= tolORMax2)
1134  && (radi2 >= tolORMin2)
1135  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1136  {
1137  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1138  {
1139  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1140  if (cosPsi >= cosHDPhiOT)
1141  {
1142  snxt = sd;
1143  }
1144  }
1145  else
1146  {
1147  snxt = sd;
1148  }
1149  }
1150  }
1151  }
1152  }
1153 
1154  // Possible intersection with ETheta cone.
1155  // Second >= 0 root should be considered
1156 
1157  if (eTheta < UUtils::kPi)
1158  {
1159  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1160  t2 = pDotV2d - p.z * v.z * tanETheta2;
1161  if (t1)
1162  {
1163  b = t2 / t1;
1164  c = dist2ETheta / t1;
1165  d2 = b * b - c;
1166 
1167  if (d2 >= 0)
1168  {
1169  d = std::sqrt(d2);
1170  sd = -b + d; // Second root
1171 
1172  if ((sd >= 0) && (sd < snxt))
1173  {
1174  xi = p.x + sd * v.x;
1175  yi = p.y + sd * v.y;
1176  zi = p.z + sd * v.z;
1177  rhoi2 = xi * xi + yi * yi;
1178  radi2 = rhoi2 + zi * zi;
1179 
1180  if ((radi2 <= tolORMax2)
1181  && (radi2 >= tolORMin2)
1182  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1183  {
1184  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1185  {
1186  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1187  if (cosPsi >= cosHDPhiOT)
1188  {
1189  snxt = sd;
1190  }
1191  }
1192  else
1193  {
1194  snxt = sd;
1195  }
1196  }
1197  }
1198  }
1199  }
1200  }
1201  }
1202  else if (pTheta > tolETheta)
1203  {
1204  // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
1205  // Inside (theta > etheta+tol) e-theta cone
1206  // First root of etheta cone, second if first root 'imaginary'
1207 
1208  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1209  t2 = pDotV2d - p.z * v.z * tanETheta2;
1210  if (t1)
1211  {
1212  b = t2 / t1;
1213  c = dist2ETheta / t1;
1214  d2 = b * b - c;
1215 
1216  if (d2 >= 0)
1217  {
1218  d = std::sqrt(d2);
1219  sd = -b - d; // First root
1220  zi = p.z + sd * v.z;
1221 
1222  if ((sd < 0) || (zi * (eTheta - UUtils::kPi / 2) > 0))
1223  {
1224  sd = -b + d; // second root
1225  }
1226  if ((sd >= 0) && (sd < snxt))
1227  {
1228  xi = p.x + sd * v.x;
1229  yi = p.y + sd * v.y;
1230  zi = p.z + sd * v.z;
1231  rhoi2 = xi * xi + yi * yi;
1232  radi2 = rhoi2 + zi * zi;
1233 
1234  if ((radi2 <= tolORMax2)
1235  && (radi2 >= tolORMin2)
1236  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1237  {
1238  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1239  {
1240  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1241  if (cosPsi >= cosHDPhiOT)
1242  {
1243  snxt = sd;
1244  }
1245  }
1246  else
1247  {
1248  snxt = sd;
1249  }
1250  }
1251  }
1252  }
1253  }
1254 
1255  // Possible intersection with STheta cone.
1256  // Second >= 0 root should be considered
1257 
1258  if (fSTheta)
1259  {
1260  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1261  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1262  if (t1)
1263  {
1264  b = t2 / t1;
1265  c = dist2STheta / t1;
1266  d2 = b * b - c;
1267 
1268  if (d2 >= 0)
1269  {
1270  d = std::sqrt(d2);
1271  sd = -b + d; // Second root
1272 
1273  if ((sd >= 0) && (sd < snxt))
1274  {
1275  xi = p.x + sd * v.x;
1276  yi = p.y + sd * v.y;
1277  zi = p.z + sd * v.z;
1278  rhoi2 = xi * xi + yi * yi;
1279  radi2 = rhoi2 + zi * zi;
1280 
1281  if ((radi2 <= tolORMax2)
1282  && (radi2 >= tolORMin2)
1283  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1284  {
1285  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1286  {
1287  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1288  if (cosPsi >= cosHDPhiOT)
1289  {
1290  snxt = sd;
1291  }
1292  }
1293  else
1294  {
1295  snxt = sd;
1296  }
1297  }
1298  }
1299  }
1300  }
1301  }
1302  }
1303  else if ((pTheta < tolSTheta + kAngTolerance)
1304  && (fSTheta > halfAngTolerance))
1305  {
1306  // In tolerance of stheta
1307  // If entering through solid [r,phi] => 0 to in
1308  // else try 2nd root
1309 
1310  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1311  if ((t2 >= 0 && tolIRMin2 < rad2 && rad2 < tolIRMax2 && fSTheta < UUtils::kPi / 2)
1312  || (t2 < 0 && tolIRMin2 < rad2 && rad2 < tolIRMax2 && fSTheta > UUtils::kPi / 2)
1313  || (v.z < 0 && tolIRMin2 < rad2 && rad2 < tolIRMax2 && fSTheta == UUtils::kPi / 2))
1314  {
1315  if (!fFullPhiSphere && rho2) // Check phi intersection
1316  {
1317  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
1318  if (cosPsi >= cosHDPhiIT)
1319  {
1320  return 0;
1321  }
1322  }
1323  else
1324  {
1325  return 0;
1326  }
1327  }
1328 
1329  // Not entering immediately/travelling through
1330 
1331  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1332  if (t1)
1333  {
1334  b = t2 / t1;
1335  c = dist2STheta / t1;
1336  d2 = b * b - c;
1337 
1338  if (d2 >= 0)
1339  {
1340  d = std::sqrt(d2);
1341  sd = -b + d;
1342  if ((sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < UUtils::kPi / 2))
1343  {
1344  // ^^^^^^^^^^^^^^^^^^^^^ shouldn't it be >=0 instead ?
1345  xi = p.x + sd * v.x;
1346  yi = p.y + sd * v.y;
1347  zi = p.z + sd * v.z;
1348  rhoi2 = xi * xi + yi * yi;
1349  radi2 = rhoi2 + zi * zi;
1350 
1351  if ((radi2 <= tolORMax2)
1352  && (radi2 >= tolORMin2)
1353  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1354  {
1355  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1356  {
1357  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1358  if (cosPsi >= cosHDPhiOT)
1359  {
1360  snxt = sd;
1361  }
1362  }
1363  else
1364  {
1365  snxt = sd;
1366  }
1367  }
1368  }
1369  }
1370  }
1371  }
1372  else if ((pTheta > tolETheta - kAngTolerance) && (eTheta < UUtils::kPi - kAngTolerance))
1373  {
1374 
1375  // In tolerance of etheta
1376  // If entering through solid [r,phi] => 0 to in
1377  // else try 2nd root
1378 
1379  t2 = pDotV2d - p.z * v.z * tanETheta2;
1380 
1381  if (((t2 < 0) && (eTheta < UUtils::kPi / 2)
1382  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1383  || ((t2 >= 0) && (eTheta > UUtils::kPi / 2)
1384  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1385  || ((v.z > 0) && (eTheta == UUtils::kPi / 2)
1386  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)))
1387  {
1388  if (!fFullPhiSphere && rho2) // Check phi intersection
1389  {
1390  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
1391  if (cosPsi >= cosHDPhiIT)
1392  {
1393  return 0;
1394  }
1395  }
1396  else
1397  {
1398  return 0;
1399  }
1400  }
1401 
1402  // Not entering immediately/travelling through
1403 
1404  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1405  if (t1)
1406  {
1407  b = t2 / t1;
1408  c = dist2ETheta / t1;
1409  d2 = b * b - c;
1410 
1411  if (d2 >= 0)
1412  {
1413  d = std::sqrt(d2);
1414  sd = -b + d;
1415 
1416  if ((sd >= halfCarTolerance)
1417  && (sd < snxt) && (eTheta > UUtils::kPi / 2))
1418  {
1419  xi = p.x + sd * v.x;
1420  yi = p.y + sd * v.y;
1421  zi = p.z + sd * v.z;
1422  rhoi2 = xi * xi + yi * yi;
1423  radi2 = rhoi2 + zi * zi;
1424 
1425  if ((radi2 <= tolORMax2)
1426  && (radi2 >= tolORMin2)
1427  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1428  {
1429  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1430  {
1431  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1432  if (cosPsi >= cosHDPhiOT)
1433  {
1434  snxt = sd;
1435  }
1436  }
1437  else
1438  {
1439  snxt = sd;
1440  }
1441  }
1442  }
1443  }
1444  }
1445  }
1446  else
1447  {
1448  // stheta+tol<theta<etheta-tol
1449  // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1450 
1451  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1452  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1453  if (t1)
1454  {
1455  b = t2 / t1;
1456  c = dist2STheta / t1;
1457  d2 = b * b - c;
1458 
1459  if (d2 >= 0)
1460  {
1461  d = std::sqrt(d2);
1462  sd = -b + d; // second root
1463 
1464  if ((sd >= 0) && (sd < snxt))
1465  {
1466  xi = p.x + sd * v.x;
1467  yi = p.y + sd * v.y;
1468  zi = p.z + sd * v.z;
1469  rhoi2 = xi * xi + yi * yi;
1470  radi2 = rhoi2 + zi * zi;
1471 
1472  if ((radi2 <= tolORMax2)
1473  && (radi2 >= tolORMin2)
1474  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1475  {
1476  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1477  {
1478  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1479  if (cosPsi >= cosHDPhiOT)
1480  {
1481  snxt = sd;
1482  }
1483  }
1484  else
1485  {
1486  snxt = sd;
1487  }
1488  }
1489  }
1490  }
1491  }
1492  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1493  t2 = pDotV2d - p.z * v.z * tanETheta2;
1494  if (t1)
1495  {
1496  b = t2 / t1;
1497  c = dist2ETheta / t1;
1498  d2 = b * b - c;
1499 
1500  if (d2 >= 0)
1501  {
1502  d = std::sqrt(d2);
1503  sd = -b + d; // second root
1504 
1505  if ((sd >= 0) && (sd < snxt))
1506  {
1507  xi = p.x + sd * v.x;
1508  yi = p.y + sd * v.y;
1509  zi = p.z + sd * v.z;
1510  rhoi2 = xi * xi + yi * yi;
1511  radi2 = rhoi2 + zi * zi;
1512 
1513  if ((radi2 <= tolORMax2)
1514  && (radi2 >= tolORMin2)
1515  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1516  {
1517  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1518  {
1519  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1520  if (cosPsi >= cosHDPhiOT)
1521  {
1522  snxt = sd;
1523  }
1524  }
1525  else
1526  {
1527  snxt = sd;
1528  }
1529  }
1530  }
1531  }
1532  }
1533  }
1534  }
1535  return snxt;
1536 }
1537 
1539 //
1540 // Calculate distance (<= actual) to closest surface of shape from outside
1541 // - Calculate distance to radial planes
1542 // - Only to phi planes if outside phi extent
1543 // - Only to theta planes if outside theta extent
1544 // - Return 0 if point inside
1545 
1546 double USphere::SafetyFromOutside(const UVector3& p, bool /*aAccurate*/) const
1547 {
1548  double safe = 0.0, safeRMin, safeRMax, safePhi, safeTheta;
1549  double rho2, rds, rho;
1550  double cosPsi;
1551  double pTheta, dTheta1, dTheta2;
1552  rho2 = p.x * p.x + p.y * p.y;
1553  rds = std::sqrt(rho2 + p.z * p.z);
1554  rho = std::sqrt(rho2);
1555 
1556  //
1557  // Distance to r shells
1558  //
1559  if (fRmin)
1560  {
1561  safeRMin = fRmin - rds;
1562  safeRMax = rds - fRmax;
1563  if (safeRMin > safeRMax)
1564  {
1565  safe = safeRMin;
1566  }
1567  else
1568  {
1569  safe = safeRMax;
1570  }
1571  }
1572  else
1573  {
1574  safe = rds - fRmax;
1575  }
1576 
1577  //
1578  // Distance to phi extent
1579  //
1580  if (!fFullPhiSphere && rho)
1581  {
1582  // Psi=angle from central phi to point
1583  //
1584  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / rho;
1585  if (cosPsi < std::cos(hDPhi))
1586  {
1587  // Point lies outside phi range
1588  //
1589  if ((p.y * cosCPhi - p.x * sinCPhi) <= 0)
1590  {
1591  safePhi = std::fabs(p.x * sinSPhi - p.y * cosSPhi);
1592  }
1593  else
1594  {
1595  safePhi = std::fabs(p.x * sinEPhi - p.y * cosEPhi);
1596  }
1597  if (safePhi > safe)
1598  {
1599  safe = safePhi;
1600  }
1601  }
1602  }
1603  //
1604  // Distance to Theta extent
1605  //
1606  if ((rds != 0.0) && (!fFullThetaSphere))
1607  {
1608  pTheta = std::acos(p.z / rds);
1609  if (pTheta < 0)
1610  {
1611  pTheta += UUtils::kPi;
1612  }
1613  dTheta1 = fSTheta - pTheta;
1614  dTheta2 = pTheta - eTheta;
1615  if (dTheta1 > dTheta2)
1616  {
1617  if (dTheta1 >= 0) // WHY ???????????
1618  {
1619  safeTheta = rds * std::sin(dTheta1);
1620  if (safe <= safeTheta)
1621  {
1622  safe = safeTheta;
1623  }
1624  }
1625  }
1626  else
1627  {
1628  if (dTheta2 >= 0)
1629  {
1630  safeTheta = rds * std::sin(dTheta2);
1631  if (safe <= safeTheta)
1632  {
1633  safe = safeTheta;
1634  }
1635  }
1636  }
1637  }
1638 
1639  if (safe < 0)
1640  {
1641  safe = 0;
1642  }
1643  return safe;
1644 }
1645 
1647 //
1648 // Calculate distance to surface of shape from 'inside', allowing for tolerance
1649 // - Only Calc rmax intersection if no valid rmin intersection
1650 
1651 double USphere::DistanceToOut(const UVector3& p, const UVector3& v, UVector3& n, bool& validNorm, double /*aPstep*/) const
1652 {
1653  double snxt = UUtils::Infinity(); // snxt is default return value
1654  double sphi = UUtils::Infinity(), stheta = UUtils::Infinity();
1655  ESide side = kNull, sidephi = kNull, sidetheta = kNull;
1656 
1657  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
1658  static const double halfAngTolerance = kAngTolerance * 0.5;
1659  const double halfTolerance = kTolerance * 0.5;
1660  const double halfRminTolerance = fRminTolerance * 0.5;
1661  const double rMaxPlus = fRmax + halfTolerance;
1662  const double Rmin_minus = (fRmin) ? fRmin - halfRminTolerance : 0;
1663  double t1, t2;
1664  double b, c, d;
1665 
1666  // Variables for phi intersection:
1667 
1668  double pDistS, compS, pDistE, compE, sphi2, vphi;
1669 
1670  double rho2, rad2, pDotV2d, pDotV3d;
1671 
1672  double xi, yi, zi; // Intersection point
1673 
1674  // Theta precals
1675  //
1676  double rhoSecTheta;
1677  double dist2STheta, dist2ETheta, distTheta;
1678  double d2, sd;
1679 
1680  // General Precalcs
1681  //
1682  rho2 = p.x * p.x + p.y * p.y;
1683  rad2 = rho2 + p.z * p.z;
1684 
1685  pDotV2d = p.x * v.x + p.y * v.y;
1686  pDotV3d = pDotV2d + p.z * v.z;
1687 
1688  // Radial Intersections from USphere::DistanceToIn
1689  //
1690  // Outer spherical shell intersection
1691  // - Only if outside tolerant fRmax
1692  // - Check for if inside and outer USphere heading through solid (-> 0)
1693  // - No intersect -> no intersection with USphere
1694  //
1695  // Shell eqn: x^2+y^2+z^2=RSPH^2
1696  //
1697  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
1698  //
1699  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
1700  // => rad2 +2sd(pDotV3d) +sd^2 =R^2
1701  //
1702  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
1703 
1704  if ((rad2 <= rMaxPlus * rMaxPlus) && (rad2 >= Rmin_minus * Rmin_minus))
1705  {
1706  c = rad2 - fRmax * fRmax;
1707 
1708  if (c < kTolerance * fRmax)
1709  {
1710  // Within tolerant Outer radius
1711  //
1712  // The test is
1713  // rad - fRmax < 0.5*kRadTolerance
1714  // => rad < fRmax + 0.5*kRadTol
1715  // => rad2 < (fRmax + 0.5*kRadTol)^2
1716  // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
1717  // => rad2 - fRmax^2 <~ fRmax*kRadTol
1718 
1719  d2 = pDotV3d * pDotV3d - c;
1720 
1721  if ((c > - kTolerance * fRmax) // on tolerant surface
1722  && ((pDotV3d >= 0) || (d2 < 0))) // leaving outside from Rmax
1723  // not re-entering
1724  {
1725  validNorm = true;
1726  n = UVector3(p.x / fRmax, p.y / fRmax, p.z / fRmax);
1727  return snxt = 0;
1728  }
1729  else
1730  {
1731  snxt = -pDotV3d + std::sqrt(d2); // second root since inside Rmax
1732  side = kRMax;
1733  }
1734  }
1735 
1736  // Inner spherical shell intersection:
1737  // Always first >=0 root, because would have passed
1738  // from outside of Rmin surface .
1739 
1740  if (fRmin)
1741  {
1742  c = rad2 - fRmin * fRmin;
1743  d2 = pDotV3d * pDotV3d - c;
1744 
1745  if (c > - fRminTolerance * fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
1746  {
1747  if ((c < fRminTolerance * fRmin) // leaving from Rmin
1748  && (d2 >= fRminTolerance * fRmin) && (pDotV3d < 0))
1749  {
1750  validNorm = false; // Rmin surface is concave
1751  return snxt = 0;
1752  }
1753  else
1754  {
1755  if (d2 >= 0.)
1756  {
1757  sd = -pDotV3d - std::sqrt(d2);
1758 
1759  if (sd >= 0.) // Always intersect Rmin first
1760  {
1761  snxt = sd;
1762  side = kRMin;
1763  }
1764  }
1765  }
1766  }
1767  }
1768  }
1769 
1770  // Theta segment intersection
1771 
1772  if (!fFullThetaSphere)
1773  {
1774  // Intersection with theta surfaces
1775  //
1776  // Known failure cases:
1777  // o Inside tolerance of stheta surface, skim
1778  // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1779  //
1780  // To solve: Check 2nd root of etheta surface in addition to stheta
1781  //
1782  // o start/end theta is exactly UUtils::kPi/2
1783  //
1784  // Intersections with cones
1785  //
1786  // Cone equation: x^2+y^2=z^2tan^2(t)
1787  //
1788  // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1789  //
1790  // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1791  // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1792  //
1793  // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1794  //
1795 
1796  if (fSTheta) // intersection with first cons
1797  {
1798  if (std::fabs(tanSTheta) > 5. / kAngTolerance) // kons is plane z=0
1799  {
1800  if (v.z > 0.)
1801  {
1802  if (std::fabs(p.z) <= halfTolerance)
1803  {
1804  validNorm = true;
1805  n = UVector3(0., 0., 1.);
1806  return snxt = 0;
1807  }
1808  stheta = -p.z / v.z;
1809  sidetheta = kSTheta;
1810  }
1811  }
1812  else // kons is not plane
1813  {
1814  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1815  t2 = pDotV2d - p.z * v.z * tanSTheta2; // ~vDotN if p on cons
1816  dist2STheta = rho2 - p.z * p.z * tanSTheta2; // t3
1817 
1818  distTheta = std::sqrt(rho2) - p.z * tanSTheta;
1819 
1820  if (std::fabs(t1) < halfAngTolerance) // 1st order equation,
1821  {
1822  // v parallel to kons
1823  if (v.z > 0.)
1824  {
1825  if (std::fabs(distTheta) < halfTolerance) // p on surface
1826  {
1827  if ((fSTheta < UUtils::kPi / 2) && (p.z > 0.))
1828  {
1829  validNorm = false;
1830  return snxt = 0.;
1831  }
1832  else if ((fSTheta > UUtils::kPi / 2) && (p.z <= 0))
1833  {
1834  validNorm = true;
1835  if (rho2)
1836  {
1837  rhoSecTheta = std::sqrt(rho2 * (1 + tanSTheta2));
1838 
1839  n = UVector3(p.x / rhoSecTheta,
1840  p.y / rhoSecTheta,
1841  std::sin(fSTheta));
1842  }
1843  else n = UVector3(0., 0., 1.);
1844  return snxt = 0.;
1845  }
1846  }
1847  stheta = -0.5 * dist2STheta / t2;
1848  sidetheta = kSTheta;
1849  }
1850  } // 2nd order equation, 1st root of fSTheta cone,
1851  else // 2nd if 1st root -ve
1852  {
1853  if (std::fabs(distTheta) < halfTolerance)
1854  {
1855  if ((fSTheta > UUtils::kPi / 2) && (t2 >= 0.)) // leave
1856  {
1857  validNorm = true;
1858  if (rho2)
1859  {
1860  rhoSecTheta = std::sqrt(rho2 * (1 + tanSTheta2));
1861 
1862  n = UVector3(p.x / rhoSecTheta,
1863  p.y / rhoSecTheta,
1864  std::sin(fSTheta));
1865  }
1866  else
1867  {
1868  n = UVector3(0., 0., 1.);
1869  }
1870  return snxt = 0.;
1871  }
1872  else if ((fSTheta < UUtils::kPi / 2) && (t2 < 0.) && (p.z >= 0.)) // leave
1873  {
1874  validNorm = false;
1875  return snxt = 0.;
1876  }
1877  }
1878  b = t2 / t1;
1879  c = dist2STheta / t1;
1880  d2 = b * b - c;
1881 
1882  if (d2 >= 0.)
1883  {
1884  d = std::sqrt(d2);
1885 
1886  if (fSTheta > UUtils::kPi / 2)
1887  {
1888  sd = -b - d; // First root
1889 
1890  if (((std::fabs(sd) < halfTolerance) && (t2 < 0.))
1891  || (sd < 0.) || ((sd > 0.) && (p.z + sd * v.z > 0.)))
1892  {
1893  sd = -b + d; // 2nd root
1894  }
1895  if ((sd > halfTolerance) && (p.z + sd * v.z <= 0.))
1896  {
1897  stheta = sd;
1898  sidetheta = kSTheta;
1899  }
1900  }
1901  else // sTheta < UUtils::kPi/2, concave surface, no normal
1902  {
1903  sd = -b - d; // First root
1904 
1905  if (((std::fabs(sd) < halfTolerance) && (t2 >= 0.))
1906  || (sd < 0.) || ((sd > 0.) && (p.z + sd * v.z < 0.)))
1907  {
1908  sd = -b + d; // 2nd root
1909  }
1910  if ((sd > halfTolerance) && (p.z + sd * v.z >= 0.))
1911  {
1912  stheta = sd;
1913  sidetheta = kSTheta;
1914  }
1915  }
1916  }
1917  }
1918  }
1919  }
1920  if (eTheta < UUtils::kPi) // intersection with second cons
1921  {
1922  if (std::fabs(tanETheta) > 5. / kAngTolerance) // kons is plane z=0
1923  {
1924  if (v.z < 0.)
1925  {
1926  if (std::fabs(p.z) <= halfTolerance)
1927  {
1928  validNorm = true;
1929  n = UVector3(0., 0., -1.);
1930  return snxt = 0;
1931  }
1932  sd = -p.z / v.z;
1933 
1934  if (sd < stheta)
1935  {
1936  stheta = sd;
1937  sidetheta = kETheta;
1938  }
1939  }
1940  }
1941  else // kons is not plane
1942  {
1943  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1944  t2 = pDotV2d - p.z * v.z * tanETheta2; // ~vDotN if p on cons
1945  dist2ETheta = rho2 - p.z * p.z * tanETheta2; // t3
1946 
1947  distTheta = std::sqrt(rho2) - p.z * tanETheta;
1948 
1949  if (std::fabs(t1) < halfAngTolerance) // 1st order equation,
1950  {
1951  // v parallel to kons
1952  if (v.z < 0.)
1953  {
1954  if (std::fabs(distTheta) < halfTolerance) // p on surface
1955  {
1956  if ((eTheta > UUtils::kPi / 2) && (p.z < 0.))
1957  {
1958  validNorm = false;
1959  return snxt = 0.;
1960  }
1961  else if ((eTheta < UUtils::kPi / 2) && (p.z >= 0))
1962  {
1963  validNorm = true;
1964  if (rho2)
1965  {
1966  rhoSecTheta = std::sqrt(rho2 * (1 + tanETheta2));
1967  n = UVector3(p.x / rhoSecTheta,
1968  p.y / rhoSecTheta,
1969  -sinETheta);
1970  }
1971  else
1972  {
1973  n = UVector3(0., 0., -1.);
1974  }
1975  return snxt = 0.;
1976  }
1977  }
1978  sd = -0.5 * dist2ETheta / t2;
1979 
1980  if (sd < stheta)
1981  {
1982  stheta = sd;
1983  sidetheta = kETheta;
1984  }
1985  }
1986  } // 2nd order equation, 1st root of fSTheta cone
1987  else // 2nd if 1st root -ve
1988  {
1989  if (std::fabs(distTheta) < halfTolerance)
1990  {
1991  if ((eTheta < UUtils::kPi / 2) && (t2 >= 0.)) // leave
1992  {
1993  validNorm = true;
1994  if (rho2)
1995  {
1996  rhoSecTheta = std::sqrt(rho2 * (1 + tanETheta2));
1997  n = UVector3(p.x / rhoSecTheta,
1998  p.y / rhoSecTheta,
1999  -sinETheta);
2000  }
2001  else n = UVector3(0., 0., -1.);
2002  return snxt = 0.;
2003  }
2004  else if ((eTheta > UUtils::kPi / 2)
2005  && (t2 < 0.) && (p.z <= 0.)) // leave
2006  {
2007  validNorm = false;
2008  return snxt = 0.;
2009  }
2010  }
2011  b = t2 / t1;
2012  c = dist2ETheta / t1;
2013  d2 = b * b - c;
2014 
2015  if (d2 >= 0.)
2016  {
2017  d = std::sqrt(d2);
2018 
2019  if (eTheta < UUtils::kPi / 2)
2020  {
2021  sd = -b - d; // First root
2022 
2023  if (((std::fabs(sd) < halfTolerance) && (t2 < 0.))
2024  || (sd < 0.))
2025  {
2026  sd = -b + d; // 2nd root
2027  }
2028  if (sd > halfTolerance)
2029  {
2030  if (sd < stheta)
2031  {
2032  stheta = sd;
2033  sidetheta = kETheta;
2034  }
2035  }
2036  }
2037  else // sTheta+fDTheta > UUtils::kPi/2, concave surface, no normal
2038  {
2039  sd = -b - d; // First root
2040 
2041  if (((std::fabs(sd) < halfTolerance) && (t2 >= 0.))
2042  || (sd < 0.) || ((sd > 0.) && (p.z + sd * v.z > 0.)))
2043  {
2044  sd = -b + d; // 2nd root
2045  }
2046  if ((sd > halfTolerance) && (p.z + sd * v.z <= 0.))
2047  {
2048  if (sd < stheta)
2049  {
2050  stheta = sd;
2051  sidetheta = kETheta;
2052  }
2053  }
2054  }
2055  }
2056  }
2057  }
2058  }
2059 
2060  } // end theta intersections
2061 
2062  // Phi Intersection
2063 
2064  if (!fFullPhiSphere)
2065  {
2066  if (p.x || p.y) // Check if on z axis (rho not needed later)
2067  {
2068  // pDist -ve when inside
2069 
2070  pDistS = p.x * sinSPhi - p.y * cosSPhi;
2071  pDistE = -p.x * sinEPhi + p.y * cosEPhi;
2072 
2073  // Comp -ve when in direction of outwards normal
2074 
2075  compS = -sinSPhi * v.x + cosSPhi * v.y;
2076  compE = sinEPhi * v.x - cosEPhi * v.y;
2077  sidephi = kNull;
2078 
2079  if ((pDistS <= 0) && (pDistE <= 0))
2080  {
2081  // Inside both phi *full* planes
2082 
2083  if (compS < 0)
2084  {
2085  sphi = pDistS / compS;
2086  xi = p.x + sphi * v.x;
2087  yi = p.y + sphi * v.y;
2088 
2089  // Check intersection with correct half-plane (if not -> no intersect)
2090  //
2091  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2092  {
2093  vphi = std::atan2(v.y, v.x);
2094  sidephi = kSPhi;
2095  if (((fSPhi - halfAngTolerance) <= vphi)
2096  && ((ePhi + halfAngTolerance) >= vphi))
2097  {
2098  sphi = UUtils::Infinity();
2099  }
2100  }
2101  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
2102  {
2103  sphi = UUtils::Infinity();
2104  }
2105  else
2106  {
2107  sidephi = kSPhi;
2108  if (pDistS > -halfCarTolerance)
2109  {
2110  sphi = 0; // Leave by sphi
2111  }
2112  }
2113  }
2114  else
2115  {
2116  sphi = UUtils::Infinity();
2117  }
2118 
2119  if (compE < 0)
2120  {
2121  sphi2 = pDistE / compE;
2122  if (sphi2 < sphi) // Only check further if < starting phi intersection
2123  {
2124  xi = p.x + sphi2 * v.x;
2125  yi = p.y + sphi2 * v.y;
2126 
2127  // Check intersection with correct half-plane
2128  //
2129  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2130  {
2131  // Leaving via ending phi
2132  //
2133  vphi = std::atan2(v.y, v.x);
2134 
2135  if (!((fSPhi - halfAngTolerance <= vphi)
2136  && (fSPhi + fDPhi + halfAngTolerance >= vphi)))
2137  {
2138  sidephi = kEPhi;
2139  if (pDistE <= -halfCarTolerance)
2140  {
2141  sphi = sphi2;
2142  }
2143  else
2144  {
2145  sphi = 0.0;
2146  }
2147  }
2148  }
2149  else if ((yi * cosCPhi - xi * sinCPhi) >= 0) // Leaving via ending phi
2150  {
2151  sidephi = kEPhi;
2152  if (pDistE <= -halfCarTolerance)
2153  {
2154  sphi = sphi2;
2155  }
2156  else
2157  {
2158  sphi = 0;
2159  }
2160  }
2161  }
2162  }
2163  }
2164  else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
2165  {
2166  if (pDistS <= pDistE)
2167  {
2168  sidephi = kSPhi;
2169  }
2170  else
2171  {
2172  sidephi = kEPhi;
2173  }
2174  if (fDPhi > UUtils::kPi)
2175  {
2176  if ((compS < 0) && (compE < 0))
2177  {
2178  sphi = 0;
2179  }
2180  else
2181  {
2182  sphi = UUtils::Infinity();
2183  }
2184  }
2185  else
2186  {
2187  // if towards both >=0 then once inside (after error)
2188  // will remain inside
2189 
2190  if ((compS >= 0) && (compE >= 0))
2191  {
2192  sphi = UUtils::Infinity();
2193  }
2194  else
2195  {
2196  sphi = 0;
2197  }
2198  }
2199  }
2200  else if ((pDistS > 0) && (pDistE < 0))
2201  {
2202  // Outside full starting plane, inside full ending plane
2203 
2204  if (fDPhi > UUtils::kPi)
2205  {
2206  if (compE < 0)
2207  {
2208  sphi = pDistE / compE;
2209  xi = p.x + sphi * v.x;
2210  yi = p.y + sphi * v.y;
2211 
2212  // Check intersection in correct half-plane
2213  // (if not -> not leaving phi extent)
2214  //
2215  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2216  {
2217  vphi = std::atan2(v.y, v.x);
2218  sidephi = kSPhi;
2219  if (((fSPhi - halfAngTolerance) <= vphi)
2220  && ((ePhi + halfAngTolerance) >= vphi))
2221  {
2222  sphi = UUtils::Infinity();
2223  }
2224  }
2225  else if ((yi * cosCPhi - xi * sinCPhi) <= 0)
2226  {
2227  sphi = UUtils::Infinity();
2228  }
2229  else // Leaving via Ending phi
2230  {
2231  sidephi = kEPhi;
2232  if (pDistE > -halfCarTolerance)
2233  {
2234  sphi = 0.;
2235  }
2236  }
2237  }
2238  else
2239  {
2240  sphi = UUtils::Infinity();
2241  }
2242  }
2243  else
2244  {
2245  if (compS >= 0)
2246  {
2247  if (compE < 0)
2248  {
2249  sphi = pDistE / compE;
2250  xi = p.x + sphi * v.x;
2251  yi = p.y + sphi * v.y;
2252 
2253  // Check intersection in correct half-plane
2254  // (if not -> remain in extent)
2255  //
2256  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2257  {
2258  vphi = std::atan2(v.y, v.x);
2259  sidephi = kSPhi;
2260  if (((fSPhi - halfAngTolerance) <= vphi)
2261  && ((ePhi + halfAngTolerance) >= vphi))
2262  {
2263  sphi = UUtils::Infinity();
2264  }
2265  }
2266  else if ((yi * cosCPhi - xi * sinCPhi) <= 0)
2267  {
2268  sphi = UUtils::Infinity();
2269  }
2270  else // otherwise leaving via Ending phi
2271  {
2272  sidephi = kEPhi;
2273  }
2274  }
2275  else sphi = UUtils::Infinity();
2276  }
2277  else // leaving immediately by starting phi
2278  {
2279  sidephi = kSPhi;
2280  sphi = 0;
2281  }
2282  }
2283  }
2284  else
2285  {
2286  // Must be pDistS < 0 && pDistE > 0
2287  // Inside full starting plane, outside full ending plane
2288 
2289  if (fDPhi > UUtils::kPi)
2290  {
2291  if (compS < 0)
2292  {
2293  sphi = pDistS / compS;
2294  xi = p.x + sphi * v.x;
2295  yi = p.y + sphi * v.y;
2296 
2297  // Check intersection in correct half-plane
2298  // (if not -> not leaving phi extent)
2299  //
2300  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2301  {
2302  vphi = std::atan2(v.y, v.x);
2303  sidephi = kSPhi;
2304  if (((fSPhi - halfAngTolerance) <= vphi)
2305  && ((ePhi + halfAngTolerance) >= vphi))
2306  {
2307  sphi = UUtils::Infinity();
2308  }
2309  }
2310  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
2311  {
2312  sphi = UUtils::Infinity();
2313  }
2314  else // Leaving via Starting phi
2315  {
2316  sidephi = kSPhi;
2317  if (pDistS > -halfCarTolerance)
2318  {
2319  sphi = 0;
2320  }
2321  }
2322  }
2323  else
2324  {
2325  sphi = UUtils::Infinity();
2326  }
2327  }
2328  else
2329  {
2330  if (compE >= 0)
2331  {
2332  if (compS < 0)
2333  {
2334  sphi = pDistS / compS;
2335  xi = p.x + sphi * v.x;
2336  yi = p.y + sphi * v.y;
2337 
2338  // Check intersection in correct half-plane
2339  // (if not -> remain in extent)
2340  //
2341  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2342  {
2343  vphi = std::atan2(v.y, v.x);
2344  sidephi = kSPhi;
2345  if (((fSPhi - halfAngTolerance) <= vphi)
2346  && ((ePhi + halfAngTolerance) >= vphi))
2347  {
2348  sphi = UUtils::Infinity();
2349  }
2350  }
2351  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
2352  {
2353  sphi = UUtils::Infinity();
2354  }
2355  else // otherwise leaving via Starting phi
2356  {
2357  sidephi = kSPhi;
2358  }
2359  }
2360  else
2361  {
2362  sphi = UUtils::Infinity();
2363  }
2364  }
2365  else // leaving immediately by ending
2366  {
2367  sidephi = kEPhi;
2368  sphi = 0 ;
2369  }
2370  }
2371  }
2372  }
2373  else
2374  {
2375  // On z axis + travel not || to z axis -> if phi of vector direction
2376  // within phi of shape, Step limited by rmax, else Step =0
2377 
2378  if (v.x || v.y)
2379  {
2380  vphi = std::atan2(v.y, v.x);
2381  if ((fSPhi - halfAngTolerance < vphi) && (vphi < ePhi + halfAngTolerance))
2382  {
2383  sphi = UUtils::Infinity();
2384  }
2385  else
2386  {
2387  sidephi = kSPhi; // arbitrary
2388  sphi = 0;
2389  }
2390  }
2391  else // travel along z - no phi intersection
2392  {
2393  sphi = UUtils::Infinity();
2394  }
2395  }
2396  if (sphi < snxt) // Order intersecttions
2397  {
2398  snxt = sphi;
2399  side = sidephi;
2400  }
2401  }
2402  if (stheta < snxt) // Order intersections
2403  {
2404  snxt = stheta;
2405  side = sidetheta;
2406  }
2407 
2408  switch (side)
2409  {
2410  case kRMax:
2411  xi = p.x + snxt * v.x;
2412  yi = p.y + snxt * v.y;
2413  zi = p.z + snxt * v.z;
2414  n = UVector3(xi / fRmax, yi / fRmax, zi / fRmax);
2415  validNorm = true;
2416  break;
2417 
2418  case kRMin:
2419  validNorm = false; // Rmin is concave
2420  break;
2421 
2422  case kSPhi:
2423  if (fDPhi <= UUtils::kPi) // Normal to Phi-
2424  {
2425  n = UVector3(sinSPhi, -cosSPhi, 0);
2426  validNorm = true;
2427  }
2428  else
2429  {
2430  validNorm = false;
2431  }
2432  break;
2433 
2434  case kEPhi:
2435  if (fDPhi <= UUtils::kPi) // Normal to Phi+
2436  {
2437  n = UVector3(-sinEPhi, cosEPhi, 0);
2438  validNorm = true;
2439  }
2440  else
2441  {
2442  validNorm = false;
2443  }
2444  break;
2445 
2446  case kSTheta:
2447  if (fSTheta == UUtils::kPi / 2)
2448  {
2449  n = UVector3(0., 0., 1.);
2450  validNorm = true;
2451  }
2452  else if (fSTheta > UUtils::kPi / 2)
2453  {
2454  xi = p.x + snxt * v.x;
2455  yi = p.y + snxt * v.y;
2456  rho2 = xi * xi + yi * yi;
2457  if (rho2)
2458  {
2459  rhoSecTheta = std::sqrt(rho2 * (1 + tanSTheta2));
2460  n = UVector3(xi / rhoSecTheta, yi / rhoSecTheta,
2461  -tanSTheta / std::sqrt(1 + tanSTheta2));
2462  }
2463  else
2464  {
2465  n = UVector3(0., 0., 1.);
2466  }
2467  validNorm = true;
2468  }
2469  else
2470  {
2471  validNorm = false; // Concave STheta cone
2472  }
2473  break;
2474 
2475  case kETheta:
2476  if (eTheta == UUtils::kPi / 2)
2477  {
2478  n = UVector3(0., 0., -1.);
2479  validNorm = true;
2480  }
2481  else if (eTheta < UUtils::kPi / 2)
2482  {
2483  xi = p.x + snxt * v.x;
2484  yi = p.y + snxt * v.y;
2485  rho2 = xi * xi + yi * yi;
2486  if (rho2)
2487  {
2488  rhoSecTheta = std::sqrt(rho2 * (1 + tanETheta2));
2489  n = UVector3(xi / rhoSecTheta, yi / rhoSecTheta,
2490  -tanETheta / std::sqrt(1 + tanETheta2));
2491  }
2492  else
2493  {
2494  n = UVector3(0., 0., -1.);
2495  }
2496  validNorm = true;
2497  }
2498  else
2499  {
2500  validNorm = false; // Concave ETheta cone
2501  }
2502  break;
2503 
2504  default:
2505  cout << std::endl;
2506 
2507  std::ostringstream message;
2508  int oldprc = message.precision(16);
2509  message << "Undefined side for valid surface normal to solid."
2510  << std::endl
2511  << "Position:" << std::endl << std::endl
2512  << "p.x = " << p.x << " mm" << std::endl
2513  << "p.y = " << p.y << " mm" << std::endl
2514  << "p.z = " << p.z << " mm" << std::endl << std::endl
2515  << "Direction:" << std::endl << std::endl
2516  << "v.x = " << v.x << std::endl
2517  << "v.y = " << v.y << std::endl
2518  << "v.z = " << v.z << std::endl << std::endl
2519  << "Proposed distance :" << std::endl << std::endl
2520  << "snxt = " << snxt << " mm" << std::endl;
2521  message.precision(oldprc);
2522  UUtils::Exception("USphere::DistanceToOut(p,v,..)",
2523  "GeomSolids1002", Warning, 1, message.str().c_str());
2524  break;
2525  }
2526  if (snxt == UUtils::Infinity())
2527  {
2528  cout << std::endl;
2529 
2530  std::ostringstream message;
2531  int oldprc = message.precision(16);
2532  message << "Logic error: snxt = UUtils::Infinity() ???" << std::endl
2533  << "Position:" << std::endl << std::endl
2534  << "p.x = " << p.x << " mm" << std::endl
2535  << "p.y = " << p.y << " mm" << std::endl
2536  << "p.z = " << p.z << " mm" << std::endl << std::endl
2537  << "Rp = " << std::sqrt(p.x * p.x + p.y * p.y + p.z * p.z)
2538  << " mm" << std::endl << std::endl
2539  << "Direction:" << std::endl << std::endl
2540  << "v.x = " << v.x << std::endl
2541  << "v.y = " << v.y << std::endl
2542  << "v.z = " << v.z << std::endl << std::endl
2543  << "Proposed distance :" << std::endl << std::endl
2544  << "snxt = " << snxt << " mm" << std::endl;
2545  message.precision(oldprc);
2546  UUtils::Exception("USphere::DistanceToOut(p,v,..)",
2547  "GeomSolids1002", Warning, 1, message.str().c_str());
2548  }
2549 
2550  return snxt;
2551 }
2552 
2554 //
2555 // Calculate distance (<=actual) to closest surface of shape from inside
2556 
2557 double USphere::SafetyFromInside(const UVector3& p, bool /*aAccurate*/) const
2558 {
2559  double safe = 0.0, safeRMin, safeRMax, safePhi, safeTheta;
2560  double rho2, rds, rho;
2561  double pTheta, dTheta1, dTheta2;
2562  rho2 = p.x * p.x + p.y * p.y;
2563  rds = std::sqrt(rho2 + p.z * p.z);
2564  rho = std::sqrt(rho2);
2565 
2566 #ifdef UDEBUG
2567  if (Inside(p) == eOutside)
2568  {
2569  int old_prc = cout.precision(16);
2570  cout << std::endl;
2571 
2572  cout << "Position:" << std::endl << std::endl;
2573  cout << "p.x = " << p.x << " mm" << std::endl;
2574  cout << "p.y = " << p.y << " mm" << std::endl;
2575  cout << "p.z = " << p.z << " mm" << std::endl << std::endl;
2576  cout.precision(old_prc);
2577  UUtils::Exception("USphere::DistanceToOut(p)",
2578  "GeomSolids1002", Warning, 1, "Point p is outside !?");
2579  }
2580 #endif
2581 
2582  //
2583  // Distance to r shells
2584  //
2585  if (fRmin)
2586  {
2587  safeRMin = rds - fRmin;
2588  safeRMax = fRmax - rds;
2589  if (safeRMin < safeRMax)
2590  {
2591  safe = safeRMin;
2592  }
2593  else
2594  {
2595  safe = safeRMax;
2596  }
2597  }
2598  else
2599  {
2600  safe = fRmax - rds;
2601  }
2602 
2603  //
2604  // Distance to phi extent
2605  //
2606  if (!fFullPhiSphere && rho)
2607  {
2608  if ((p.y * cosCPhi - p.x * sinCPhi) <= 0)
2609  {
2610  safePhi = -(p.x * sinSPhi - p.y * cosSPhi);
2611  }
2612  else
2613  {
2614  safePhi = (p.x * sinEPhi - p.y * cosEPhi);
2615  }
2616  if (safePhi < safe)
2617  {
2618  safe = safePhi;
2619  }
2620  }
2621 
2622  //
2623  // Distance to Theta extent
2624  //
2625  if (rds)
2626  {
2627  pTheta = std::acos(p.z / rds);
2628  if (pTheta < 0)
2629  {
2630  pTheta += UUtils::kPi;
2631  }
2632  dTheta1 = pTheta - fSTheta;
2633  dTheta2 = eTheta - pTheta;
2634  if (dTheta1 < dTheta2)
2635  {
2636  safeTheta = rds * std::sin(dTheta1);
2637  }
2638  else
2639  {
2640  safeTheta = rds * std::sin(dTheta2);
2641  }
2642  if (safe > safeTheta)
2643  {
2644  safe = safeTheta;
2645  }
2646  }
2647 
2648  if (safe < 0)
2649  {
2650  safe = 0;
2651  }
2652  return safe;
2653 }
2654 
2656 //
2657 // Create a List containing the transformed vertices
2658 // Ordering [0-3] -fDz Cross section
2659 // [4-7] +fDz Cross section such that [0] is below [4],
2660 // [1] below [5] etc.
2661 // Note:
2662 // Caller has deletion resposibility
2663 // Potential improvement: For last slice, use actual ending angle
2664 // to avoid rounding error problems.
2665 
2666 /*
2667 UVector3List*
2668 USphere::CreateRotatedVertices(const UAffineTransform& pTransform,
2669  int& noPolygonVertices) const
2670 {
2671  UVector3List *vertices;
2672  UVector3 vertex;
2673  double meshAnglePhi,meshRMax,crossAnglePhi,
2674  coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
2675  double meshTheta,crossTheta,startTheta;
2676  double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
2677  int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;
2678 
2679  // Phi Cross sections
2680 
2681  noPhiCrossSections = int(fDPhi/UUtils::kMeshAngleDefault)+1;
2682 
2683  if (noPhiCrossSections<UUtils::kMinMeshSections)
2684  {
2685  noPhiCrossSections=UUtils::kMinMeshSections;
2686  }
2687  else if (noPhiCrossSections>UUtils::kMaxMeshSections)
2688  {
2689  noPhiCrossSections=UUtils::kMaxMeshSections;
2690  }
2691  meshAnglePhi=fDPhi/(noPhiCrossSections-1);
2692 
2693  // If complete in phi, Set start angle such that mesh will be at fRMax
2694  // on the x axis. Will give better extent calculations when not rotated.
2695 
2696  if (fFullPhiSphere)
2697  {
2698  sAnglePhi = -meshAnglePhi*0.5;
2699  }
2700  else
2701  {
2702  sAnglePhi=fSPhi;
2703  }
2704 
2705  // Theta Cross sections
2706 
2707  noThetaSections = int(fDTheta/UUtils::kMeshAngleDefault)+1;
2708 
2709  if (noThetaSections<UUtils::kMinMeshSections)
2710  {
2711  noThetaSections=UUtils::kMinMeshSections;
2712  }
2713  else if (noThetaSections>UUtils::kMaxMeshSections)
2714  {
2715  noThetaSections=UUtils::kMaxMeshSections;
2716  }
2717  meshTheta=fDTheta/(noThetaSections-1);
2718 
2719  // If complete in Theta, Set start angle such that mesh will be at fRMax
2720  // on the z axis. Will give better extent calculations when not rotated.
2721 
2722  if (fFullThetaSphere)
2723  {
2724  startTheta = -meshTheta*0.5;
2725  }
2726  else
2727  {
2728  startTheta=fSTheta;
2729  }
2730 
2731  meshRMax = (meshAnglePhi >= meshTheta) ?
2732  fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5);
2733  double* cosCrossTheta = new double[noThetaSections];
2734  double* sinCrossTheta = new double[noThetaSections];
2735  vertices=new UVector3List();
2736  if (vertices && cosCrossTheta && sinCrossTheta)
2737  {
2738  vertices->reserve(noPhiCrossSections*(noThetaSections*2));
2739  for (crossSectionPhi=0;
2740  crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
2741  {
2742  crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
2743  coscrossAnglePhi=std::cos(crossAnglePhi);
2744  sincrossAnglePhi=std::sin(crossAnglePhi);
2745  for (crossSectionTheta=0;
2746  crossSectionTheta<noThetaSections;crossSectionTheta++)
2747  {
2748  // Compute coordinates of Cross section at section crossSectionPhi
2749  //
2750  crossTheta=startTheta+crossSectionTheta*meshTheta;
2751  cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
2752  sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
2753 
2754  rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
2755  rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
2756  rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
2757 
2758  vertex=UVector3(rMinX,rMinY,rMinZ);
2759  vertices->push_back(pTransform.TransformPoint(vertex));
2760 
2761  } // Theta forward
2762 
2763  for (crossSectionTheta=noThetaSections-1;
2764  crossSectionTheta>=0; crossSectionTheta--)
2765  {
2766  rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
2767  rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
2768  rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
2769 
2770  vertex=UVector3(rMaxX,rMaxY,rMaxZ);
2771  vertices->push_back(pTransform.TransformPoint(vertex));
2772 
2773  } // Theta back
2774  } // Phi
2775  noPolygonVertices = noThetaSections*2;
2776  }
2777  else
2778  {
2779 
2780  UUtils::Exception("USphere::CreateRotatedVertices()",
2781  "GeomSolids0003", FatalError,1,
2782  "Error in allocation of vertices. Out of memory !");
2783  }
2784 
2785  delete [] cosCrossTheta;
2786  delete [] sinCrossTheta;
2787 
2788  return vertices;
2789 }
2790 */
2791 
2793 //
2794 // UEntityType
2795 
2796 UGeometryType USphere::GetEntityType() const
2797 {
2798  return std::string("Sphere");
2799 }
2800 
2802 //
2803 // Make a clone of the object
2804 //
2805 VUSolid* USphere::Clone() const
2806 {
2807  return new USphere(*this);
2808 }
2809 
2811 //
2812 // Stream object contents to an output stream
2813 
2814 std::ostream& USphere::StreamInfo(std::ostream& os) const
2815 {
2816  int oldprc = os.precision(16);
2817  os << "-----------------------------------------------------------\n"
2818  << " *** Dump for solid - " << GetName() << " ***\n"
2819  << " ===================================================\n"
2820  << " Solid type: USphere\n"
2821  << " Parameters: \n"
2822  << " inner radius: " << fRmin << " mm \n"
2823  << " outer radius: " << fRmax << " mm \n"
2824  << " starting phi of segment : " << fSPhi / (UUtils::kPi / 180.0) << " degrees \n"
2825  << " delta phi of segment : " << fDPhi / (UUtils::kPi / 180.0) << " degrees \n"
2826  << " starting theta of segment: " << fSTheta / (UUtils::kPi / 180.0) << " degrees \n"
2827  << " delta theta of segment : " << fDTheta / (UUtils::kPi / 180.0) << " degrees \n"
2828  << "-----------------------------------------------------------\n";
2829  os.precision(oldprc);
2830 
2831  return os;
2832 }
2833 
2835 //
2836 // GetPointOnSurface
2837 
2838 UVector3 USphere::GetPointOnSurface() const
2839 {
2840  double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
2841  double height1, height2, slant1, slant2, costheta, sintheta, rRand;
2842 
2843  height1 = (fRmax - fRmin) * cosSTheta;
2844  height2 = (fRmax - fRmin) * cosETheta;
2845  slant1 = std::sqrt(UUtils::sqr((fRmax - fRmin) * sinSTheta) + height1 * height1);
2846  slant2 = std::sqrt(UUtils::sqr((fRmax - fRmin) * sinETheta) + height2 * height2);
2847  rRand = UUtils::GetRadiusInRing(fRmin, fRmax);
2848 
2849  aOne = fRmax * fRmax * fDPhi * (cosSTheta - cosETheta);
2850  aTwo = fRmin * fRmin * fDPhi * (cosSTheta - cosETheta);
2851  aThr = fDPhi * ((fRmax + fRmin) * sinSTheta) * slant1;
2852  aFou = fDPhi * ((fRmax + fRmin) * sinETheta) * slant2;
2853  aFiv = 0.5 * fDTheta * (fRmax * fRmax - fRmin * fRmin);
2854 
2855  phi = UUtils::Random(fSPhi, ePhi);
2856  cosphi = std::cos(phi);
2857  sinphi = std::sin(phi);
2858  costheta = UUtils::Random(cosETheta, cosSTheta);
2859  sintheta = std::sqrt(1. - UUtils::sqr(costheta));
2860 
2861  if (fFullPhiSphere)
2862  {
2863  aFiv = 0;
2864  }
2865  if (fSTheta == 0)
2866  {
2867  aThr = 0;
2868  }
2869  if (eTheta == UUtils::kPi)
2870  {
2871  aFou = 0;
2872  }
2873  if (fSTheta == UUtils::kPi / 2)
2874  {
2875  aThr = UUtils::kPi * (fRmax * fRmax - fRmin * fRmin);
2876  }
2877  if (eTheta == UUtils::kPi / 2)
2878  {
2879  aFou = UUtils::kPi * (fRmax * fRmax - fRmin * fRmin);
2880  }
2881 
2882  chose = UUtils::Random(0., aOne + aTwo + aThr + aFou + 2.*aFiv);
2883  if ((chose >= 0.) && (chose < aOne))
2884  {
2885  return UVector3(fRmax * sintheta * cosphi,
2886  fRmax * sintheta * sinphi, fRmax * costheta);
2887  }
2888  else if ((chose >= aOne) && (chose < aOne + aTwo))
2889  {
2890  return UVector3(fRmin * sintheta * cosphi,
2891  fRmin * sintheta * sinphi, fRmin * costheta);
2892  }
2893  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr))
2894  {
2895  if (fSTheta != UUtils::kPi / 2)
2896  {
2897  zRand = UUtils::Random(fRmin * cosSTheta, fRmax * cosSTheta);
2898  return UVector3(tanSTheta * zRand * cosphi,
2899  tanSTheta * zRand * sinphi, zRand);
2900  }
2901  else
2902  {
2903  return UVector3(rRand * cosphi, rRand * sinphi, 0.);
2904  }
2905  }
2906  else if ((chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + aThr + aFou))
2907  {
2908  if (eTheta != UUtils::kPi / 2)
2909  {
2910  zRand = UUtils::Random(fRmin * cosETheta, fRmax * cosETheta);
2911  return UVector3(tanETheta * zRand * cosphi,
2912  tanETheta * zRand * sinphi, zRand);
2913  }
2914  else
2915  {
2916  return UVector3(rRand * cosphi, rRand * sinphi, 0.);
2917  }
2918  }
2919  else if ((chose >= aOne + aTwo + aThr + aFou) && (chose < aOne + aTwo + aThr + aFou + aFiv))
2920  {
2921  return UVector3(rRand * sintheta * cosSPhi,
2922  rRand * sintheta * sinSPhi, rRand * costheta);
2923  }
2924  else
2925  {
2926  return UVector3(rRand * sintheta * cosEPhi,
2927  rRand * sintheta * sinEPhi, rRand * costheta);
2928  }
2929 }
2930 
2932 //
2933 // GetSurfaceArea
2934 
2935 double USphere::SurfaceArea()
2936 {
2937  if (fSurfaceArea != 0.)
2938  {
2939  ;
2940  }
2941  else
2942  {
2943  double Rsq = fRmax * fRmax;
2944  double rsq = fRmin * fRmin;
2945 
2946  fSurfaceArea = fDPhi * (rsq + Rsq) * (cosSTheta - cosETheta);
2947  if (!fFullPhiSphere)
2948  {
2949  fSurfaceArea = fSurfaceArea + fDTheta * (Rsq - rsq);
2950  }
2951  if (fSTheta > 0)
2952  {
2953  double acos1 = std::acos(std::pow(sinSTheta, 2) * std::cos(fDPhi)
2954  + std::pow(cosSTheta, 2));
2955  if (fDPhi > UUtils::kPi)
2956  {
2957  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * (2 * UUtils::kPi - acos1);
2958  }
2959  else
2960  {
2961  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * acos1;
2962  }
2963  }
2964  if (eTheta < UUtils::kPi)
2965  {
2966  double acos2 = std::acos(std::pow(sinETheta, 2) * std::cos(fDPhi)
2967  + std::pow(cosETheta, 2));
2968  if (fDPhi > UUtils::kPi)
2969  {
2970  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * (2 * UUtils::kPi - acos2);
2971  }
2972  else
2973  {
2974  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * acos2;
2975  }
2976  }
2977  }
2978  return fSurfaceArea;
2979 }
2980 
2981 
2982 
2983 void USphere::Extent(UVector3& aMin, UVector3& aMax) const
2984 {
2985  aMin.Set(-fRmax);
2986  aMax.Set(fRmax);
2987 }
2988 
2989 void USphere::GetParametersList(int, double* aArray) const
2990 {
2991  aArray[0] = GetInnerRadius();
2992  aArray[1] = GetOuterRadius();
2993  aArray[2] = GetStartPhiAngle();
2994  aArray[3] = GetDeltaPhiAngle();
2995  aArray[4] = GetStartThetaAngle();
2996  aArray[5] = GetDeltaThetaAngle();
2997 }
kETheta
Definition: G4Sphere.cc:79
t1
TTree * t1
Definition: plottest35.C:26
VUSolid::frTolerance
static double frTolerance
Definition: VUSolid.hh:31
USphere::~USphere
~USphere()
Definition: USphere.cc:70
d
Float_t d
Definition: plot.C:237
USphere::GetInnerRadius
double GetInnerRadius() const
Definition: USphere.hh:212
kRMin
Definition: G4Cons.cc:68
FatalErrorInArguments
Definition: UUtils.hh:36
VUSolid::GetName
const std::string & GetName() const
Definition: VUSolid.hh:103
in
ifstream in
Definition: comparison.C:7
kNSTheta
Definition: G4Sphere.cc:83
p
const char * p
Definition: xmltok.h:285
USphere::Inside
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: USphere.cc:164
USphere::Extent
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: USphere.cc:2983
USphere::DistanceToOut
double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &n, bool &validNorm, double aPstep=UUtils::kInfinity) const
Definition: USphere.cc:1651
norm
Float_t norm
Definition: comparison_ascii.C:71
Warning
Definition: UUtils.hh:36
USphere::GetDeltaThetaAngle
double GetDeltaThetaAngle() const
Definition: USphere.hh:241
USphere::operator=
USphere & operator=(const USphere &rhs)
Definition: USphere.cc:103
USphere::GetOuterRadius
double GetOuterRadius() const
Definition: USphere.hh:218
UUtils::Infinity
double Infinity()
Definition: UUtils.hh:177
USphere::DistanceToIn
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: USphere.cc:611
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
USphere::GetPointOnSurface
UVector3 GetPointOnSurface() const
Definition: USphere.cc:2838
kNull
Definition: G4Cons.cc:68
ENorm
ENorm
Definition: G4Cons.cc:72
USphere::SafetyFromOutside
double SafetyFromOutside(const UVector3 &p, bool aAccurate=false) const
Definition: USphere.cc:1546
test.b
tuple b
Definition: test.py:12
n
Char_t n[5]
Definition: comparison_ascii.C:55
USphere::GetEntityType
UGeometryType GetEntityType() const
Definition: USphere.cc:2796
kSPhi
Definition: G4Cons.cc:68
VUSolid::faTolerance
static double faTolerance
Definition: VUSolid.hh:32
VUSolid::EnumInside
EnumInside
Definition: VUSolid.hh:23
kSTheta
Definition: G4Sphere.cc:79
test.v
tuple v
Definition: test.py:18
UUtils::sqr
T sqr(const T &x)
Definition: UUtils.hh:142
UUtils::GetRadiusInRing
double GetRadiusInRing(double rmin, double rmax)
Definition: UUtils.hh:160
USphere::USphere
USphere(const std::string &pName, double pRmin, double pRmax, double pSPhi, double pDPhi, double pSTheta, double pDTheta)
Definition: USphere.cc:34
UVector3::Set
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:245
USphere::GetParametersList
void GetParametersList(int, double *) const
Definition: USphere.cc:2989
kNETheta
Definition: G4Sphere.cc:83
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
USphere::GetStartThetaAngle
double GetStartThetaAngle() const
Definition: USphere.hh:236
USphere::Clone
VUSolid * Clone() const
Definition: USphere.cc:2805
UVector3::Unit
UVector3 Unit() const
Definition: UVector3.cc:80
USphere::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: USphere.cc:2814
USphere.hh
UGeometryType
std::string UGeometryType
Definition: UTypes.hh:70
t2
TTree * t2
Definition: plottest35.C:36
USphere::Normal
bool Normal(const UVector3 &p, UVector3 &n) const
Definition: USphere.cc:273
kNEPhi
Definition: G4Cons.cc:72
USphere::GetStartPhiAngle
double GetStartPhiAngle() const
Definition: USphere.hh:224
ESide
ESide
Definition: G4Cons.cc:68
UUtils::Exception
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:177
UVector3::z
double z
Definition: UVector3.hh:138
kRMax
Definition: G4Cons.cc:68
test.c
tuple c
Definition: test.py:13
UUtils::Random
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
kNRMax
Definition: G4Cons.cc:72
USphere::GetDeltaPhiAngle
double GetDeltaPhiAngle() const
Definition: USphere.hh:230
UVector3::y
double y
Definition: UVector3.hh:137
kEPhi
Definition: G4Cons.cc:68
USphere::SafetyFromInside
double SafetyFromInside(const UVector3 &p, bool aAccurate=false) const
Definition: USphere.cc:2557
USphere::SurfaceArea
double SurfaceArea()
Definition: USphere.cc:2935
kNRMin
Definition: G4Cons.cc:72
kNSPhi
Definition: G4Cons.cc:72
Generated on Sat Dec 14 2013 14:34:14 for Geant4_10 by   doxygen 1.8.5