Geant4_10
UTubs.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 // UTubs
12 //
13 // 19.10.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "UUtils.hh"
18 #include <string>
19 #include <cmath>
20 #include <sstream>
21 #include "UTubs.hh"
22 
23 using namespace std;
24 
26 //
27 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
28 // - note if pdphi>2PI then reset to 2PI
29 
30 UTubs::UTubs(const std::string& pName,
31  double pRMin, double pRMax,
32  double pDz,
33  double pSPhi, double pDPhi)
34  : VUSolid(pName.c_str()), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
35 {
36  fCubicVolume=0.;
37  fSurfaceArea=0.;
40 
41  if (pDz <= 0) // Check z-len
42  {
43  std::ostringstream message;
44  message << "Negative Z half-length (" << pDz << ") in solid: " << GetName();
45  UUtils::Exception("UTubs::UTubs()", "GeomSolids0002", FatalErrorInArguments, 1, message.str().c_str());
46  }
47  if ((pRMin >= pRMax) || (pRMin < 0)) // Check radii
48  {
49  std::ostringstream message;
50  message << "Invalid values for radii in solid: " << GetName()
51  << std::endl
52  << " pRMin = " << pRMin << ", pRMax = " << pRMax;
53  UUtils::Exception("UTubs::UTubs()", "GeomSolids0002", FatalErrorInArguments, 1, message.str().c_str());
54  }
55 
56  // Check angles
57  //
58  CheckPhiAngles(pSPhi, pDPhi);
59 }
60 
62 //
63 // Fake default constructor - sets only member data and allocates memory
64 // for usage restricted to object persistency.
65 //
67  : VUSolid(""), fCubicVolume(0.), fSurfaceArea(0.),
68  kRadTolerance(0.), kAngTolerance(0.),
69  fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
70  fSinCPhi(0.), fCosCPhi(0.), fCosHDPhiOT(0.), fCosHDPhiIT(0.),
71  fSinSPhi(0.), fCosSPhi(0.), fSinEPhi(0.), fCosEPhi(0.),
72  fSinSPhiDPhi(0.), fCosSPhiDPhi(0.),
73  fPhiFullTube(false)
74 {
75 }
76 
78 //
79 // Destructor
80 
82 {
83 }
84 
86 //
87 // Copy constructor
88 
89 UTubs::UTubs(const UTubs& rhs)
90  : VUSolid(rhs),
91  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
92  kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
93  fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
94  fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
95  fSinCPhi(rhs.fSinCPhi), fCosCPhi(rhs.fSinCPhi),
96  fCosHDPhiOT(rhs.fCosHDPhiOT), fCosHDPhiIT(rhs.fCosHDPhiOT),
97  fSinSPhi(rhs.fSinSPhi), fCosSPhi(rhs.fCosSPhi),
98  fSinEPhi(rhs.fSinEPhi), fCosEPhi(rhs.fCosEPhi),
99  fSinSPhiDPhi(rhs.fSinSPhiDPhi), fCosSPhiDPhi(rhs.fCosSPhiDPhi),
100  fPhiFullTube(rhs.fPhiFullTube)
101 {
102 }
103 
105 //
106 // Assignment operator
107 
109 {
110  // Check assignment to self
111  //
112  if (this == &rhs)
113  {
114  return *this;
115  }
116 
117  // Copy base class data
118  //
119  VUSolid::operator=(rhs);
120 
121  // Copy data
122  //
127  fRMin = rhs.fRMin;
128  fRMax = rhs.fRMax;
129  fDz = rhs.fDz;
130  fSPhi = rhs.fSPhi;
131  fDPhi = rhs.fDPhi;
132  fSinCPhi = rhs.fSinCPhi;
133  fCosCPhi = rhs.fSinCPhi;
134  fCosHDPhiOT = rhs.fCosHDPhiOT;
135  fCosHDPhiIT = rhs.fCosHDPhiOT;
136  fSinSPhi = rhs.fSinSPhi;
137  fCosSPhi = rhs.fCosSPhi;
138  fSinEPhi = rhs.fSinEPhi;
139  fCosEPhi = rhs.fCosEPhi;
142 
144 
145  return *this;
146 }
147 
148 
150 //
151 // Return whether point inside/outside/on surface
152 
154 {
155  double r2, pPhi, tolRMin, tolRMax;
157  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
158  static const double halfRadTolerance = kRadTolerance * 0.5;
159  static const double halfAngTolerance = kAngTolerance * 0.5;
160 
161  if (std::fabs(p.z) <= fDz - halfCarTolerance)
162  {
163  r2 = p.x * p.x + p.y * p.y;
164 
165  if (fRMin)
166  {
167  tolRMin = fRMin + halfRadTolerance;
168  }
169  else
170  {
171  tolRMin = 0;
172  }
173 
174  tolRMax = fRMax - halfRadTolerance;
175 
176  if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
177  {
178  if (fPhiFullTube)
179  {
180  in = eInside;
181  }
182  else
183  {
184  // Try inner tolerant phi boundaries (=>inside)
185  // if not inside, try outer tolerant phi boundaries
186 
187  if ((tolRMin == 0) && (std::fabs(p.x) <= halfCarTolerance)
188  && (std::fabs(p.y) <= halfCarTolerance))
189  {
190  in = eSurface;
191  }
192  else
193  {
194  pPhi = std::atan2(p.y, p.x);
195  if (pPhi < -halfAngTolerance)
196  {
197  pPhi += 2 * UUtils::kPi; // 0<=pPhi<2UUtils::kPi
198  }
199 
200  if (fSPhi >= 0)
201  {
202  if ((std::fabs(pPhi) < halfAngTolerance)
203  && (std::fabs(fSPhi + fDPhi - 2 * UUtils::kPi) < halfAngTolerance))
204  {
205  pPhi += 2 * UUtils::kPi; // 0 <= pPhi < 2UUtils::kPi
206  }
207  if ((pPhi >= fSPhi + halfAngTolerance)
208  && (pPhi <= fSPhi + fDPhi - halfAngTolerance))
209  {
210  in = eInside;
211  }
212  else if ((pPhi >= fSPhi - halfAngTolerance)
213  && (pPhi <= fSPhi + fDPhi + halfAngTolerance))
214  {
215  in = eSurface;
216  }
217  }
218  else // fSPhi < 0
219  {
220  if ((pPhi <= fSPhi + 2 * UUtils::kPi - halfAngTolerance)
221  && (pPhi >= fSPhi + fDPhi + halfAngTolerance))
222  {
223  ; //eOutside
224  }
225  else if ((pPhi <= fSPhi + 2 * UUtils::kPi + halfAngTolerance)
226  && (pPhi >= fSPhi + fDPhi - halfAngTolerance))
227  {
228  in = eSurface;
229  }
230  else
231  {
232  in = eInside;
233  }
234  }
235  }
236  }
237  }
238  else // Try generous boundaries
239  {
240  tolRMin = fRMin - halfRadTolerance;
241  tolRMax = fRMax + halfRadTolerance;
242 
243  if (tolRMin < 0)
244  {
245  tolRMin = 0;
246  }
247 
248  if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
249  {
250  if (fPhiFullTube || (r2 <= halfRadTolerance * halfRadTolerance))
251  {
252  // Continuous in phi or on z-axis
253  in = eSurface;
254  }
255  else // Try outer tolerant phi boundaries only
256  {
257  pPhi = std::atan2(p.y, p.x);
258 
259  if (pPhi < -halfAngTolerance)
260  {
261  pPhi += 2 * UUtils::kPi; // 0<=pPhi<2UUtils::kPi
262  }
263  if (fSPhi >= 0)
264  {
265  if ((std::fabs(pPhi) < halfAngTolerance)
266  && (std::fabs(fSPhi + fDPhi - 2 * UUtils::kPi) < halfAngTolerance))
267  {
268  pPhi += 2 * UUtils::kPi; // 0 <= pPhi < 2UUtils::kPi
269  }
270  if ((pPhi >= fSPhi - halfAngTolerance)
271  && (pPhi <= fSPhi + fDPhi + halfAngTolerance))
272  {
273  in = eSurface;
274  }
275  }
276  else // fSPhi < 0
277  {
278  if ((pPhi <= fSPhi + 2 * UUtils::kPi - halfAngTolerance)
279  && (pPhi >= fSPhi + fDPhi + halfAngTolerance))
280  {
281  ; // eOutside
282  }
283  else
284  {
285  in = eSurface;
286  }
287  }
288  }
289  }
290  }
291  }
292  else if (std::fabs(p.z) <= fDz + halfCarTolerance)
293  {
294  // Check within tolerant r limits
295  r2 = p.x * p.x + p.y * p.y;
296  tolRMin = fRMin - halfRadTolerance;
297  tolRMax = fRMax + halfRadTolerance;
298 
299  if (tolRMin < 0)
300  {
301  tolRMin = 0;
302  }
303 
304  if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
305  {
306  if (fPhiFullTube || (r2 <= halfRadTolerance * halfRadTolerance))
307  {
308  // Continuous in phi or on z-axis
309  in = eSurface;
310  }
311  else // Try outer tolerant phi boundaries
312  {
313  pPhi = std::atan2(p.y, p.x);
314 
315  if (pPhi < -halfAngTolerance)
316  {
317  pPhi += 2 * UUtils::kPi; // 0<=pPhi<2UUtils::kPi
318  }
319  if (fSPhi >= 0)
320  {
321  if ((std::fabs(pPhi) < halfAngTolerance)
322  && (std::fabs(fSPhi + fDPhi - 2 * UUtils::kPi) < halfAngTolerance))
323  {
324  pPhi += 2 * UUtils::kPi; // 0 <= pPhi < 2UUtils::kPi
325  }
326  if ((pPhi >= fSPhi - halfAngTolerance)
327  && (pPhi <= fSPhi + fDPhi + halfAngTolerance))
328  {
329  in = eSurface;
330  }
331  }
332  else // fSPhi < 0
333  {
334  if ((pPhi <= fSPhi + 2 * UUtils::kPi - halfAngTolerance)
335  && (pPhi >= fSPhi + fDPhi + halfAngTolerance))
336  {
337  ;
338  }
339  else
340  {
341  in = eSurface;
342  }
343  }
344  }
345  }
346  }
347  return in;
348 }
349 
351 //
352 // Return Unit normal of surface closest to p
353 // - note if point on z axis, ignore phi divided sides
354 // - unsafe if point close to z axis a rmin=0 - no explicit checks
355 
356 bool UTubs::Normal(const UVector3& p, UVector3& n) const
357 {
358  int noSurfaces = 0;
359  double rho, pPhi;
360  double distZ, distRMin, distRMax;
361  double distSPhi = UUtils::kInfinity, distEPhi = UUtils::kInfinity;
362 
363  static const double halfCarTolerance = 0.5 * VUSolid::Tolerance();
364  static const double halfAngTolerance = 0.5 * kAngTolerance;
365 
366  UVector3 norm, sumnorm(0., 0., 0.);
367  UVector3 nZ = UVector3(0, 0, 1.0);
368  UVector3 nR, nPs, nPe;
369 
370  rho = std::sqrt(p.x * p.x + p.y * p.y);
371 
372  distRMin = std::fabs(rho - fRMin);
373  distRMax = std::fabs(rho - fRMax);
374  distZ = std::fabs(std::fabs(p.z) - fDz);
375 
376  if (!fPhiFullTube) // Protected against (0,0,z)
377  {
378  if (rho > halfCarTolerance)
379  {
380  pPhi = std::atan2(p.y, p.x);
381 
382  if (pPhi < fSPhi - halfCarTolerance)
383  {
384  pPhi += 2 * UUtils::kPi;
385  }
386  else if (pPhi > fSPhi + fDPhi + halfCarTolerance)
387  {
388  pPhi -= 2 * UUtils::kPi;
389  }
390 
391  distSPhi = std::fabs(pPhi - fSPhi);
392  distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
393  }
394  else if (!fRMin)
395  {
396  distSPhi = 0.;
397  distEPhi = 0.;
398  }
399  nPs = UVector3(fSinSPhi, -fCosSPhi, 0);
400  nPe = UVector3(-fSinSPhiDPhi /* std::sin(fSPhi+fDPhi)*/, fCosSPhiDPhi, 0);
401  }
402  if (rho > halfCarTolerance)
403  {
404  nR = UVector3(p.x / rho, p.y / rho, 0);
405  }
406 
407  if (distRMax <= halfCarTolerance)
408  {
409  noSurfaces ++;
410  sumnorm += nR;
411  }
412  if (fRMin && (distRMin <= halfCarTolerance))
413  {
414  noSurfaces ++;
415  sumnorm -= nR;
416  }
417  if (fDPhi < 2 * UUtils::kPi)
418  {
419  if (distSPhi <= halfAngTolerance)
420  {
421  noSurfaces ++;
422  sumnorm += nPs;
423  }
424  if (distEPhi <= halfAngTolerance)
425  {
426  noSurfaces ++;
427  sumnorm += nPe;
428  }
429  }
430  if (distZ <= halfCarTolerance)
431  {
432  noSurfaces ++;
433  if (p.z >= 0.)
434  {
435  sumnorm += nZ;
436  }
437  else
438  {
439  sumnorm -= nZ;
440  }
441  }
442  if (noSurfaces == 0)
443  {
444 #ifdef UDEBUG
445  UUtils::Exception("UTubs::SurfaceNormal(p)", "GeomSolids1002",
446  Warning, 1, "Point p is not on surface !?");
447  int oldprc = cout.precision(20);
448  cout << "UTubs::SN ( " << p.x << ", " << p.y << ", " << p.z << " ); "
449  << std::endl << std::endl;
450  cout.precision(oldprc);
451 #endif
452  norm = ApproxSurfaceNormal(p);
453  }
454  else if (noSurfaces == 1)
455  {
456  norm = sumnorm;
457  }
458  else
459  {
460  norm = sumnorm.Unit();
461  }
462 
463  n = norm;
464 
465  return noSurfaces; // TODO: return true or false on validity
466 }
467 
469 //
470 // Algorithm for SurfaceNormal() following the original specification
471 // for points not on the surface
472 
474 {
475  ENorm side;
476  UVector3 norm;
477  double rho, phi;
478  double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin;
479 
480  rho = std::sqrt(p.x * p.x + p.y * p.y);
481 
482  distRMin = std::fabs(rho - fRMin);
483  distRMax = std::fabs(rho - fRMax);
484  distZ = std::fabs(std::fabs(p.z) - fDz);
485 
486  if (distRMin < distRMax) // First minimum
487  {
488  if (distZ < distRMin)
489  {
490  distMin = distZ;
491  side = kNZ;
492  }
493  else
494  {
495  distMin = distRMin;
496  side = kNRMin ;
497  }
498  }
499  else
500  {
501  if (distZ < distRMax)
502  {
503  distMin = distZ;
504  side = kNZ ;
505  }
506  else
507  {
508  distMin = distRMax;
509  side = kNRMax ;
510  }
511  }
512  if (!fPhiFullTube && rho) // Protected against (0,0,z)
513  {
514  phi = std::atan2(p.y, p.x);
515 
516  if (phi < 0)
517  {
518  phi += 2 * UUtils::kPi;
519  }
520 
521  if (fSPhi < 0)
522  {
523  distSPhi = std::fabs(phi - (fSPhi + 2 * UUtils::kPi)) * rho;
524  }
525  else
526  {
527  distSPhi = std::fabs(phi - fSPhi) * rho;
528  }
529  distEPhi = std::fabs(phi - fSPhi - fDPhi) * rho;
530 
531  if (distSPhi < distEPhi) // Find new minimum
532  {
533  if (distSPhi < distMin)
534  {
535  side = kNSPhi;
536  }
537  }
538  else
539  {
540  if (distEPhi < distMin)
541  {
542  side = kNEPhi;
543  }
544  }
545  }
546  switch (side)
547  {
548  case kNRMin : // Inner radius
549  {
550  norm = UVector3(-p.x / rho, -p.y / rho, 0);
551  break;
552  }
553  case kNRMax : // Outer radius
554  {
555  norm = UVector3(p.x / rho, p.y / rho, 0);
556  break;
557  }
558  case kNZ : // + or - dz
559  {
560  if (p.z > 0)
561  {
562  norm = UVector3(0, 0, 1);
563  }
564  else
565  {
566  norm = UVector3(0, 0, -1);
567  }
568  break;
569  }
570  case kNSPhi:
571  {
572  norm = UVector3(fSinSPhi, -fCosSPhi, 0);
573  break;
574  }
575  case kNEPhi:
576  {
577  norm = UVector3(-fSinSPhiDPhi, fCosSPhiDPhi, 0);
578  break;
579  }
580  default: // Should never reach this case ...
581  {
582  // DumpInfo();
583  UUtils::Exception("UTubs::ApproxSurfaceNormal()",
584  "GeomSolids1002", Warning, 1,
585  "Undefined side for valid surface normal to solid.");
586  break;
587  }
588  }
589  return norm;
590 }
591 
593 //
594 //
595 // Calculate distance to shape from outside, along normalised vector
596 // - return UUtils::kInfinity if no intersection, or intersection distance <= tolerance
597 //
598 // - Compute the intersection with the z planes
599 // - if at valid r, phi, return
600 //
601 // -> If point is outer outer radius, compute intersection with rmax
602 // - if at valid phi,z return
603 //
604 // -> Compute intersection with inner radius, taking largest +ve root
605 // - if valid (in z,phi), save intersction
606 //
607 // -> If phi segmented, compute intersections with phi half planes
608 // - return smallest of valid phi intersections and
609 // inner radius intersection
610 //
611 // NOTE:
612 // - 'if valid' implies tolerant checking of intersection points
613 
614 double UTubs::DistanceToIn(const UVector3& p, const UVector3& v, double) const
615 {
616  double snxt = UUtils::kInfinity; // snxt = default return value
617  double tolORMin2, tolIRMax2; // 'generous' radii squared
618  double tolORMax2, tolIRMin2, tolODz, tolIDz;
619  const double dRmax = 100.*fRMax;
620 
621  static const double halfCarTolerance = 0.5 * VUSolid::Tolerance();
622  static const double halfRadTolerance = 0.5 * kRadTolerance;
623 
624  // Intersection point variables
625  //
626  double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp;
627  double t1, t2, t3, b, c, d; // Quadratic solver variables
628 
629  // Calculate tolerant rmin and rmax
630 
631  if (fRMin > kRadTolerance)
632  {
633  tolORMin2 = (fRMin - halfRadTolerance) * (fRMin - halfRadTolerance);
634  tolIRMin2 = (fRMin + halfRadTolerance) * (fRMin + halfRadTolerance);
635  }
636  else
637  {
638  tolORMin2 = 0.0;
639  tolIRMin2 = 0.0;
640  }
641  tolORMax2 = (fRMax + halfRadTolerance) * (fRMax + halfRadTolerance);
642  tolIRMax2 = (fRMax - halfRadTolerance) * (fRMax - halfRadTolerance);
643 
644  // Intersection with Z surfaces
645 
646  tolIDz = fDz - halfCarTolerance;
647  tolODz = fDz + halfCarTolerance;
648 
649  if (std::fabs(p.z) >= tolIDz)
650  {
651  if (p.z * v.z < 0) // at +Z going in -Z or visa versa
652  {
653  sd = (std::fabs(p.z) - fDz) / std::fabs(v.z); // Z intersect distance
654 
655  if (sd < 0.0)
656  {
657  sd = 0.0;
658  }
659 
660  xi = p.x + sd * v.x; // Intersection coords
661  yi = p.y + sd * v.y;
662  rho2 = xi * xi + yi * yi;
663 
664  // Check validity of intersection
665 
666  if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
667  {
668  if (!fPhiFullTube && rho2)
669  {
670  // Psi = angle made with central (average) phi of shape
671  //
672  inum = xi * fCosCPhi + yi * fSinCPhi;
673  iden = std::sqrt(rho2);
674  cosPsi = inum / iden;
675  if (cosPsi >= fCosHDPhiIT)
676  {
677  return sd;
678  }
679  }
680  else
681  {
682  return sd;
683  }
684  }
685  }
686  else
687  {
688  if (snxt < halfCarTolerance)
689  {
690  snxt = 0;
691  }
692  return snxt; // On/outside extent, and heading away
693  // -> cannot intersect
694  }
695  }
696 
697  // -> Can not intersect z surfaces
698  //
699  // Intersection with rmax (possible return) and rmin (must also check phi)
700  //
701  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
702  //
703  // Intersects with x^2+y^2=R^2
704  //
705  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
706  // t1 t2 t3
707 
708  t1 = 1.0 - v.z * v.z;
709  t2 = p.x * v.x + p.y * v.y;
710  t3 = p.x * p.x + p.y * p.y;
711 
712  if (t1 > 0) // Check not || to z axis
713  {
714  b = t2 / t1;
715  c = t3 - fRMax * fRMax;
716  if ((t3 >= tolORMax2) && (t2 < 0)) // This also handles the tangent case
717  {
718  // Try outer cylinder intersection
719  // c=(t3-fRMax*fRMax)/t1;
720 
721  c /= t1;
722  d = b * b - c;
723 
724  if (d >= 0) // If real root
725  {
726  sd = c / (-b + std::sqrt(d));
727  if (sd >= 0) // If 'forwards'
728  {
729  if (sd > dRmax) // Avoid rounding errors due to precision issues on
730  {
731  // 64 bits systems. Split long distances and recompute
732  double fTerm = sd - std::fmod(sd, dRmax);
733  sd = fTerm + DistanceToIn(p + fTerm * v, v);
734  }
735  // Check z intersection
736  //
737  zi = p.z + sd * v.z;
738  if (std::fabs(zi) <= tolODz)
739  {
740  // Z ok. Check phi intersection if reqd
741  //
742  if (fPhiFullTube)
743  {
744  return sd;
745  }
746  else
747  {
748  xi = p.x + sd * v.x;
749  yi = p.y + sd * v.y;
750  cosPsi = (xi * fCosCPhi + yi * fSinCPhi) / fRMax;
751  if (cosPsi >= fCosHDPhiIT)
752  {
753  return sd;
754  }
755  }
756  } // end if std::fabs(zi)
757  } // end if (sd>=0)
758  } // end if (d>=0)
759  } // end if (r>=fRMax)
760  else
761  {
762  // Inside outer radius :
763  // check not inside, and heading through tubs (-> 0 to in)
764 
765  if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z) <= tolIDz))
766  {
767  // Inside both radii, delta r -ve, inside z extent
768 
769  if (!fPhiFullTube)
770  {
771  inum = p.x * fCosCPhi + p.y * fSinCPhi;
772  iden = std::sqrt(t3);
773  cosPsi = inum / iden;
774  if (cosPsi >= fCosHDPhiIT)
775  {
776  // In the old version, the small negative tangent for the point
777  // on surface was not taken in account, and returning 0.0 ...
778  // New version: check the tangent for the point on surface and
779  // if no intersection, return UUtils::kInfinity, if intersection instead
780  // return sd.
781  //
782  c = t3 - fRMax * fRMax;
783  if (c <= 0.0)
784  {
785  return 0.0;
786  }
787  else
788  {
789  c = c / t1;
790  d = b * b - c;
791  if (d >= 0.0)
792  {
793  snxt = c / (-b + std::sqrt(d)); // using safe solution
794  // for quadratic equation
795  if (snxt < halfCarTolerance)
796  {
797  snxt = 0;
798  }
799  return snxt;
800  }
801  else
802  {
803  return UUtils::kInfinity;
804  }
805  }
806  }
807  }
808  else
809  {
810  // In the old version, the small negative tangent for the point
811  // on surface was not taken in account, and returning 0.0 ...
812  // New version: check the tangent for the point on surface and
813  // if no intersection, return UUtils::kInfinity, if intersection instead
814  // return sd.
815  //
816  c = t3 - fRMax * fRMax;
817  if (c <= 0.0)
818  {
819  return 0.0;
820  }
821  else
822  {
823  c = c / t1;
824  d = b * b - c;
825  if (d >= 0.0)
826  {
827  snxt = c / (-b + std::sqrt(d)); // using safe solution
828  // for quadratic equation
829  if (snxt < halfCarTolerance)
830  {
831  snxt = 0;
832  }
833  return snxt;
834  }
835  else
836  {
837  return UUtils::kInfinity;
838  }
839  }
840  } // end if (!fPhiFullTube)
841  } // end if (t3>tolIRMin2)
842  } // end if (Inside Outer Radius)
843  if (fRMin) // Try inner cylinder intersection
844  {
845  c = (t3 - fRMin * fRMin) / t1;
846  d = b * b - c;
847  if (d >= 0.0) // If real root
848  {
849  // Always want 2nd root - we are outside and know rmax Hit was bad
850  // - If on surface of rmin also need farthest root
851 
852  sd = (b > 0.) ? c / (-b - std::sqrt(d)) : (-b + std::sqrt(d));
853  if (sd >= -halfCarTolerance) // check forwards
854  {
855  // Check z intersection
856  //
857  if (sd < 0.0)
858  {
859  sd = 0.0;
860  }
861  if (sd > dRmax) // Avoid rounding errors due to precision issues seen
862  {
863  // 64 bits systems. Split long distances and recompute
864  double fTerm = sd - std::fmod(sd, dRmax);
865  sd = fTerm + DistanceToIn(p + fTerm * v, v);
866  }
867  zi = p.z + sd * v.z;
868  if (std::fabs(zi) <= tolODz)
869  {
870  // Z ok. Check phi
871  //
872  if (fPhiFullTube)
873  {
874  return sd;
875  }
876  else
877  {
878  xi = p.x + sd * v.x;
879  yi = p.y + sd * v.y;
880  cosPsi = (xi * fCosCPhi + yi * fSinCPhi) / fRMin;
881  if (cosPsi >= fCosHDPhiIT)
882  {
883  // Good inner radius isect
884  // - but earlier phi isect still possible
885 
886  snxt = sd;
887  }
888  }
889  } // end if std::fabs(zi)
890  } // end if (sd>=0)
891  } // end if (d>=0)
892  } // end if (fRMin)
893  }
894 
895  // Phi segment intersection
896  //
897  // o Tolerant of points inside phi planes by up to VUSolid::Tolerance()*0.5
898  //
899  // o NOTE: Large duplication of code between sphi & ephi checks
900  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
901  // intersection check <=0 -> >=0
902  // -> use some form of loop Construct ?
903  //
904  if (!fPhiFullTube)
905  {
906  // First phi surface (Starting phi)
907  //
908  Comp = v.x * fSinSPhi - v.y * fCosSPhi;
909 
910  if (Comp < 0) // Component in outwards normal dirn
911  {
912  Dist = (p.y * fCosSPhi - p.x * fSinSPhi);
913 
914  if (Dist < halfCarTolerance)
915  {
916  sd = Dist / Comp;
917 
918  if (sd < snxt)
919  {
920  if (sd < 0)
921  {
922  sd = 0.0;
923  }
924  zi = p.z + sd * v.z;
925  if (std::fabs(zi) <= tolODz)
926  {
927  xi = p.x + sd * v.x;
928  yi = p.y + sd * v.y;
929  rho2 = xi * xi + yi * yi;
930 
931  if (((rho2 >= tolIRMin2) && (rho2 <= tolIRMax2))
932  || ((rho2 > tolORMin2) && (rho2 < tolIRMin2)
933  && (v.y * fCosSPhi - v.x * fSinSPhi > 0)
934  && (v.x * fCosSPhi + v.y * fSinSPhi >= 0))
935  || ((rho2 > tolIRMax2) && (rho2 < tolORMax2)
936  && (v.y * fCosSPhi - v.x * fSinSPhi > 0)
937  && (v.x * fCosSPhi + v.y * fSinSPhi < 0)))
938  {
939  // z and r intersections good
940  // - check intersecting with correct half-plane
941  //
942  if ((yi * fCosCPhi - xi * fSinCPhi) <= halfCarTolerance)
943  {
944  snxt = sd;
945  }
946  }
947  }
948  }
949  }
950  }
951 
952  // Second phi surface (Ending phi)
953 
954  Comp = -(v.x * fSinEPhi - v.y * fCosEPhi);
955 
956  if (Comp < 0) // Component in outwards normal dirn
957  {
958  Dist = -(p.y * fCosEPhi - p.x * fSinEPhi);
959 
960  if (Dist < halfCarTolerance)
961  {
962  sd = Dist / Comp;
963 
964  if (sd < snxt)
965  {
966  if (sd < 0)
967  {
968  sd = 0;
969  }
970  zi = p.z + sd * v.z;
971  if (std::fabs(zi) <= tolODz)
972  {
973  xi = p.x + sd * v.x;
974  yi = p.y + sd * v.y;
975  rho2 = xi * xi + yi * yi;
976  if (((rho2 >= tolIRMin2) && (rho2 <= tolIRMax2))
977  || ((rho2 > tolORMin2) && (rho2 < tolIRMin2)
978  && (v.x * fSinEPhi - v.y * fCosEPhi > 0)
979  && (v.x * fCosEPhi + v.y * fSinEPhi >= 0))
980  || ((rho2 > tolIRMax2) && (rho2 < tolORMax2)
981  && (v.x * fSinEPhi - v.y * fCosEPhi > 0)
982  && (v.x * fCosEPhi + v.y * fSinEPhi < 0)))
983  {
984  // z and r intersections good
985  // - check intersecting with correct half-plane
986  //
987  if ((yi * fCosCPhi - xi * fSinCPhi) >= 0)
988  {
989  snxt = sd;
990  }
991  } //?? >=-halfCarTolerance
992  }
993  }
994  }
995  } // Comp < 0
996  } // !fPhiFullTube
997  if (snxt < halfCarTolerance)
998  {
999  snxt = 0;
1000  }
1001  return snxt;
1002 }
1003 
1005 //
1006 // Calculate distance to shape from outside, along normalised vector
1007 // - return UUtils::kInfinity if no intersection, or intersection distance <= tolerance
1008 //
1009 // - Compute the intersection with the z planes
1010 // - if at valid r, phi, return
1011 //
1012 // -> If point is outer outer radius, compute intersection with rmax
1013 // - if at valid phi,z return
1014 //
1015 // -> Compute intersection with inner radius, taking largest +ve root
1016 // - if valid (in z,phi), save intersction
1017 //
1018 // -> If phi segmented, compute intersections with phi half planes
1019 // - return smallest of valid phi intersections and
1020 // inner radius intersection
1021 //
1022 // NOTE:
1023 // - Precalculations for phi trigonometry are Done `just in time'
1024 // - `if valid' implies tolerant checking of intersection points
1025 // Calculate distance (<= actual) to closest surface of shape from outside
1026 // - Calculate distance to z, radial planes
1027 // - Only to phi planes if outside phi extent
1028 // - Return 0 if point inside
1029 
1030 double UTubs::SafetyFromOutside(const UVector3& p, bool) const
1031 {
1032  double safe = 0.0, rho, safe1, safe2, safe3;
1033  double safePhi, cosPsi;
1034 
1035  rho = std::sqrt(p.x * p.x + p.y * p.y);
1036  safe1 = fRMin - rho;
1037  safe2 = rho - fRMax;
1038  safe3 = std::fabs(p.z) - fDz;
1039 
1040  if (safe1 > safe2)
1041  {
1042  safe = safe1;
1043  }
1044  else
1045  {
1046  safe = safe2;
1047  }
1048  if (safe3 > safe)
1049  {
1050  safe = safe3;
1051  }
1052 
1053  if ((!fPhiFullTube) && (rho))
1054  {
1055  // Psi=angle from central phi to point
1056  //
1057  cosPsi = (p.x * fCosCPhi + p.y * fSinCPhi) / rho;
1058 
1059  if (cosPsi < std::cos(fDPhi * 0.5))
1060  {
1061  // Point lies outside phi range
1062 
1063  if ((p.y * fCosCPhi - p.x * fSinCPhi) <= 0)
1064  {
1065  safePhi = std::fabs(p.x * fSinSPhi - p.y * fCosSPhi);
1066  }
1067  else
1068  {
1069  safePhi = std::fabs(p.x * fSinEPhi - p.y * fCosEPhi);
1070  }
1071  if (safePhi > safe)
1072  {
1073  safe = safePhi;
1074  }
1075  }
1076  }
1077  if (safe < 0)
1078  {
1079  safe = 0;
1080  }
1081  return safe;
1082 }
1083 
1085 //
1086 // Calculate distance to surface of shape from `inside', allowing for tolerance
1087 // - Only Calc rmax intersection if no valid rmin intersection
1088 
1089 // double UTubs::DistanceToOut( const UVector3& p, const UVector3& v, const bool calcNorm, bool *validNorm, UVector3 *n ) const
1090 double UTubs::DistanceToOut(const UVector3& p, const UVector3& v, UVector3& n, bool& validNorm, double) const
1091 {
1092  ESide side = kNull , sider = kNull, sidephi = kNull;
1093  double snxt, srd = UUtils::kInfinity, sphi = UUtils::kInfinity, pdist;
1094  double deltaR, t1, t2, t3, b, c, d2, roMin2;
1095 
1096  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
1097  static const double halfAngTolerance = kAngTolerance * 0.5;
1098 
1099  // Vars for phi intersection:
1100 
1101  double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2;
1102 
1103  // Z plane intersection
1104 
1105  if (v.z > 0)
1106  {
1107  pdist = fDz - p.z;
1108  if (pdist > halfCarTolerance)
1109  {
1110  snxt = pdist / v.z;
1111  side = kPZ;
1112  }
1113  else
1114  {
1115  // if (calcNorm)
1116  {
1117  n = UVector3(0, 0, 1);
1118  validNorm = true;
1119  }
1120  return snxt = 0;
1121  }
1122  }
1123  else if (v.z < 0)
1124  {
1125  pdist = fDz + p.z;
1126 
1127  if (pdist > halfCarTolerance)
1128  {
1129  snxt = -pdist / v.z;
1130  side = kMZ;
1131  }
1132  else
1133  {
1134  // if (calcNorm)
1135  {
1136  n = UVector3(0, 0, -1);
1137  validNorm = true;
1138  }
1139  return snxt = 0.0;
1140  }
1141  }
1142  else
1143  {
1144  snxt = UUtils::kInfinity; // Travel perpendicular to z axis
1145  side = kNull;
1146  }
1147 
1148  // Radial Intersections
1149  //
1150  // Find intersection with cylinders at rmax/rmin
1151  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1152  //
1153  // Intersects with x^2+y^2=R^2
1154  //
1155  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1156  //
1157  // t1 t2 t3
1158 
1159  t1 = 1.0 - v.z * v.z; // since v normalised
1160  t2 = p.x * v.x + p.y * v.y;
1161  t3 = p.x * p.x + p.y * p.y;
1162 
1163  if (snxt > 10 * (fDz + fRMax))
1164  {
1165  roi2 = 2 * fRMax * fRMax;
1166  }
1167  else
1168  {
1169  roi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3; // radius^2 on +-fDz
1170  }
1171 
1172  if (t1 > 0) // Check not parallel
1173  {
1174  // Calculate srd, r exit distance
1175 
1176  if ((t2 >= 0.0) && (roi2 > fRMax * (fRMax + kRadTolerance)))
1177  {
1178  // Delta r not negative => leaving via rmax
1179 
1180  deltaR = t3 - fRMax * fRMax;
1181 
1182  // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1183  // - avoid sqrt for efficiency
1184 
1185  if (deltaR < -kRadTolerance * fRMax)
1186  {
1187  b = t2 / t1;
1188  c = deltaR / t1;
1189  d2 = b * b - c;
1190  if (d2 >= 0)
1191  {
1192  srd = c / (-b - std::sqrt(d2));
1193  }
1194  else
1195  {
1196  srd = 0.;
1197  }
1198  sider = kRMax;
1199  }
1200  else
1201  {
1202  // On tolerant boundary & heading outwards (or perpendicular to)
1203  // outer radial surface -> leaving immediately
1204 
1205  // if ( calcNorm )
1206  {
1207  n = UVector3(p.x / fRMax, p.y / fRMax, 0);
1208  validNorm = true;
1209  }
1210  return snxt = 0; // Leaving by rmax immediately
1211  }
1212  }
1213  else if (t2 < 0.) // i.e. t2 < 0; Possible rmin intersection
1214  {
1215  roMin2 = t3 - t2 * t2 / t1; // min ro2 of the plane of movement
1216 
1217  if (fRMin && (roMin2 < fRMin * (fRMin - kRadTolerance)))
1218  {
1219  deltaR = t3 - fRMin * fRMin;
1220  b = t2 / t1;
1221  c = deltaR / t1;
1222  d2 = b * b - c;
1223 
1224  if (d2 >= 0) // Leaving via rmin
1225  {
1226  // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1227  // - avoid sqrt for efficiency
1228 
1229  if (deltaR > kRadTolerance * fRMin)
1230  {
1231  srd = c / (-b + std::sqrt(d2));
1232  sider = kRMin;
1233  }
1234  else
1235  {
1236  //if ( calcNorm )
1237  {
1238  validNorm = false; // Concave side
1239  }
1240  return snxt = 0.0;
1241  }
1242  }
1243  else // No rmin intersect -> must be rmax intersect
1244  {
1245  deltaR = t3 - fRMax * fRMax;
1246  c = deltaR / t1;
1247  d2 = b * b - c;
1248  if (d2 >= 0.)
1249  {
1250  srd = -b + std::sqrt(d2);
1251  sider = kRMax;
1252  }
1253  else // Case: On the border+t2<kRadTolerance
1254  // (v is perpendicular to the surface)
1255  {
1256  // if (calcNorm)
1257  {
1258  n = UVector3(p.x / fRMax, p.y / fRMax, 0);
1259  validNorm = true;
1260  }
1261  return snxt = 0.0;
1262  }
1263  }
1264  }
1265  else if (roi2 > fRMax * (fRMax + kRadTolerance))
1266  // No rmin intersect -> must be rmax intersect
1267  {
1268  deltaR = t3 - fRMax * fRMax;
1269  b = t2 / t1;
1270  c = deltaR / t1;
1271  d2 = b * b - c;
1272  if (d2 >= 0)
1273  {
1274  srd = -b + std::sqrt(d2);
1275  sider = kRMax;
1276  }
1277  else // Case: On the border+t2<kRadTolerance
1278  // (v is perpendicular to the surface)
1279  {
1280  // if (calcNorm)
1281  {
1282  n = UVector3(p.x / fRMax, p.y / fRMax, 0);
1283  validNorm = true;
1284  }
1285  return snxt = 0.0;
1286  }
1287  }
1288  }
1289 
1290  // Phi Intersection
1291 
1292  if (!fPhiFullTube)
1293  {
1294  // add angle calculation with correction
1295  // of the difference in domain of atan2 and Sphi
1296  //
1297  vphi = std::atan2(v.y, v.x);
1298 
1299  if (vphi < fSPhi - halfAngTolerance)
1300  {
1301  vphi += 2 * UUtils::kPi;
1302  }
1303  else if (vphi > fSPhi + fDPhi + halfAngTolerance)
1304  {
1305  vphi -= 2 * UUtils::kPi;
1306  }
1307 
1308 
1309  if (p.x || p.y) // Check if on z axis (rho not needed later)
1310  {
1311  // pDist -ve when inside
1312 
1313  pDistS = p.x * fSinSPhi - p.y * fCosSPhi;
1314  pDistE = -p.x * fSinEPhi + p.y * fCosEPhi;
1315 
1316  // Comp -ve when in direction of outwards normal
1317 
1318  compS = -fSinSPhi * v.x + fCosSPhi * v.y;
1319  compE = fSinEPhi * v.x - fCosEPhi * v.y;
1320 
1321  sidephi = kNull;
1322 
1323  if (((fDPhi <= UUtils::kPi) && ((pDistS <= halfCarTolerance)
1324  && (pDistE <= halfCarTolerance)))
1325  || ((fDPhi > UUtils::kPi) && !((pDistS > halfCarTolerance)
1326  && (pDistE > halfCarTolerance))))
1327  {
1328  // Inside both phi *full* planes
1329 
1330  if (compS < 0)
1331  {
1332  sphi = pDistS / compS;
1333 
1334  if (sphi >= -halfCarTolerance)
1335  {
1336  xi = p.x + sphi * v.x;
1337  yi = p.y + sphi * v.y;
1338 
1339  // Check intersecting with correct half-plane
1340  // (if not -> no intersect)
1341  //
1342  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
1343  {
1344  sidephi = kSPhi;
1345  if (((fSPhi - halfAngTolerance) <= vphi)
1346  && ((fSPhi + fDPhi + halfAngTolerance) >= vphi))
1347  {
1348  sphi = UUtils::kInfinity;
1349  }
1350  }
1351  else if (yi * fCosCPhi - xi * fSinCPhi >= 0)
1352  {
1353  sphi = UUtils::kInfinity;
1354  }
1355  else
1356  {
1357  sidephi = kSPhi;
1358  if (pDistS > -halfCarTolerance)
1359  {
1360  sphi = 0.0; // Leave by sphi immediately
1361  }
1362  }
1363  }
1364  else
1365  {
1366  sphi = UUtils::kInfinity;
1367  }
1368  }
1369  else
1370  {
1371  sphi = UUtils::kInfinity;
1372  }
1373 
1374  if (compE < 0)
1375  {
1376  sphi2 = pDistE / compE;
1377 
1378  // Only check further if < starting phi intersection
1379  //
1380  if ((sphi2 > -halfCarTolerance) && (sphi2 < sphi))
1381  {
1382  xi = p.x + sphi2 * v.x;
1383  yi = p.y + sphi2 * v.y;
1384 
1385  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
1386  {
1387  // Leaving via ending phi
1388  //
1389  if (!((fSPhi - halfAngTolerance <= vphi)
1390  && (fSPhi + fDPhi + halfAngTolerance >= vphi)))
1391  {
1392  sidephi = kEPhi;
1393  if (pDistE <= -halfCarTolerance)
1394  {
1395  sphi = sphi2;
1396  }
1397  else
1398  {
1399  sphi = 0.0;
1400  }
1401  }
1402  }
1403  else // Check intersecting with correct half-plane
1404 
1405  if ((yi * fCosCPhi - xi * fSinCPhi) >= 0)
1406  {
1407  // Leaving via ending phi
1408  //
1409  sidephi = kEPhi;
1410  if (pDistE <= -halfCarTolerance)
1411  {
1412  sphi = sphi2;
1413  }
1414  else
1415  {
1416  sphi = 0.0;
1417  }
1418  }
1419  }
1420  }
1421  }
1422  else
1423  {
1424  sphi = UUtils::kInfinity;
1425  }
1426  }
1427  else
1428  {
1429  // On z axis + travel not || to z axis -> if phi of vector direction
1430  // within phi of shape, Step limited by rmax, else Step =0
1431 
1432  if ((fSPhi - halfAngTolerance <= vphi)
1433  && (vphi <= fSPhi + fDPhi + halfAngTolerance))
1434  {
1435  sphi = UUtils::kInfinity;
1436  }
1437  else
1438  {
1439  sidephi = kSPhi; // arbitrary
1440  sphi = 0.0;
1441  }
1442  }
1443  if (sphi < snxt) // Order intersecttions
1444  {
1445  snxt = sphi;
1446  side = sidephi;
1447  }
1448  }
1449  if (srd < snxt) // Order intersections
1450  {
1451  snxt = srd;
1452  side = sider;
1453  }
1454  }
1455  // if (calcNorm)
1456  {
1457  switch (side)
1458  {
1459  case kRMax:
1460  // Note: returned vector not normalised
1461  // (divide by fRMax for Unit vector)
1462  //
1463  xi = p.x + snxt * v.x;
1464  yi = p.y + snxt * v.y;
1465  n = UVector3(xi / fRMax, yi / fRMax, 0);
1466  validNorm = true;
1467  break;
1468 
1469  case kRMin:
1470  validNorm = false; // Rmin is inconvex
1471  break;
1472 
1473  case kSPhi:
1474  if (fDPhi <= UUtils::kPi)
1475  {
1476  n = UVector3(fSinSPhi, -fCosSPhi, 0);
1477  validNorm = true;
1478  }
1479  else
1480  {
1481  validNorm = false;
1482  }
1483  break;
1484 
1485  case kEPhi:
1486  if (fDPhi <= UUtils::kPi)
1487  {
1488  n = UVector3(-fSinEPhi, fCosEPhi, 0);
1489  validNorm = true;
1490  }
1491  else
1492  {
1493  validNorm = false;
1494  }
1495  break;
1496 
1497  case kPZ:
1498  n = UVector3(0, 0, 1);
1499  validNorm = true;
1500  break;
1501 
1502  case kMZ:
1503  n = UVector3(0, 0, -1);
1504  validNorm = true;
1505  break;
1506 
1507  default:
1508  cout << std::endl;
1509  // DumpInfo();
1510  std::ostringstream message;
1511  int oldprc = message.precision(16);
1512  message << "Undefined side for valid surface normal to solid."
1513  << std::endl
1514  << "Position:" << std::endl << std::endl
1515  << "p.x = " << p.x << " mm" << std::endl
1516  << "p.y = " << p.y << " mm" << std::endl
1517  << "p.z = " << p.z << " mm" << std::endl << std::endl
1518  << "Direction:" << std::endl << std::endl
1519  << "v.x = " << v.x << std::endl
1520  << "v.y = " << v.y << std::endl
1521  << "v.z = " << v.z << std::endl << std::endl
1522  << "Proposed distance :" << std::endl << std::endl
1523  << "snxt = " << snxt << " mm" << std::endl;
1524  message.precision(oldprc);
1525  UUtils::Exception("UTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1526  Warning, 1, message.str().c_str());
1527  break;
1528  }
1529  }
1530  if (snxt < halfCarTolerance)
1531  {
1532  snxt = 0;
1533  }
1534 
1535  return snxt;
1536 }
1537 
1538 
1540 //
1541 // Calculate distance (<=actual) to closest surface of shape from inside
1542 
1543 double UTubs::SafetyFromInside(const UVector3& p, bool) const
1544 {
1545  double safe = 0.0, rho, safeR1, safeR2, safeZ, safePhi;
1546  rho = std::sqrt(p.x * p.x + p.y * p.y);
1547 
1548 #ifdef UDEBUG
1549  if (Inside(p) == eOutside)
1550  {
1551  int oldprc = cout.precision(16);
1552  cout << std::endl;
1553  DumpInfo();
1554  cout << "Position:" << std::endl << std::endl;
1555  cout << "p.x = " << p.x << " mm" << std::endl;
1556  cout << "p.y = " << p.y << " mm" << std::endl;
1557  cout << "p.z = " << p.z << " mm" << std::endl << std::endl;
1558  cout.precision(oldprc);
1559  UUtils::Exception("UTubs::DistanceToOut(p)", "GeomSolids1002",
1560  Warning, 1, "Point p is outside !?");
1561  }
1562 #endif
1563 
1564  if (fRMin)
1565  {
1566  safeR1 = rho - fRMin;
1567  safeR2 = fRMax - rho;
1568 
1569  if (safeR1 < safeR2)
1570  {
1571  safe = safeR1;
1572  }
1573  else
1574  {
1575  safe = safeR2;
1576  }
1577  }
1578  else
1579  {
1580  safe = fRMax - rho;
1581  }
1582  safeZ = fDz - std::fabs(p.z);
1583 
1584  if (safeZ < safe)
1585  {
1586  safe = safeZ;
1587  }
1588 
1589  // Check if phi divided, Calc distances closest phi plane
1590  //
1591  if (!fPhiFullTube)
1592  {
1593  if (p.y * fCosCPhi - p.x * fSinCPhi <= 0)
1594  {
1595  safePhi = -(p.x * fSinCPhi - p.y * fCosSPhi);
1596  }
1597  else
1598  {
1599  safePhi = (p.x * fSinEPhi - p.y * fCosEPhi);
1600  }
1601  if (safePhi < safe)
1602  {
1603  safe = safePhi;
1604  }
1605  }
1606  if (safe < 0)
1607  {
1608  safe = 0;
1609  }
1610 
1611  return safe;
1612 }
1613 
1614 
1616 //
1617 // Stream object contents to an output stream
1618 
1620 {
1621  return std::string("Tubs");
1622 }
1623 
1625 //
1626 // Make a clone of the object
1627 //
1629 {
1630  return new UTubs(*this);
1631 }
1632 
1634 //
1635 // Stream object contents to an output stream
1636 
1637 std::ostream& UTubs::StreamInfo(std::ostream& os) const
1638 {
1639  int oldprc = os.precision(16);
1640  os << "-----------------------------------------------------------\n"
1641  << " *** Dump for solid - " << GetName() << " ***\n"
1642  << " ===================================================\n"
1643  << " Solid type: UTubs\n"
1644  << " Parameters: \n"
1645  << " inner radius : " << fRMin << " mm \n"
1646  << " outer radius : " << fRMax << " mm \n"
1647  << " half length Z: " << fDz << " mm \n"
1648  << " starting phi : " << fSPhi / (UUtils::kPi / 180.0) << " degrees \n"
1649  << " delta phi : " << fDPhi / (UUtils::kPi / 180.0) << " degrees \n"
1650  << "-----------------------------------------------------------\n";
1651  os.precision(oldprc);
1652 
1653  return os;
1654 }
1655 
1657 //
1658 // GetPointOnSurface
1659 
1661 {
1662  double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1663  aOne, aTwo, aThr, aFou;
1664  double rRand;
1665 
1666  aOne = 2.*fDz * fDPhi * fRMax;
1667  aTwo = 2.*fDz * fDPhi * fRMin;
1668  aThr = 0.5 * fDPhi * (fRMax * fRMax - fRMin * fRMin);
1669  aFou = 2.*fDz * (fRMax - fRMin);
1670 
1671  phi = UUtils::Random(fSPhi, fSPhi + fDPhi);
1672  cosphi = std::cos(phi);
1673  sinphi = std::sin(phi);
1674 
1675  rRand = UUtils::GetRadiusInRing(fRMin, fRMax);
1676 
1677  if ((fSPhi == 0) && (fDPhi == 2 * UUtils::kPi))
1678  {
1679  aFou = 0;
1680  }
1681 
1682  chose = UUtils::Random(0., aOne + aTwo + 2.*aThr + 2.*aFou);
1683 
1684  if ((chose >= 0) && (chose < aOne))
1685  {
1686  xRand = fRMax * cosphi;
1687  yRand = fRMax * sinphi;
1688  zRand = UUtils::Random(-1.*fDz, fDz);
1689  return UVector3(xRand, yRand, zRand);
1690  }
1691  else if ((chose >= aOne) && (chose < aOne + aTwo))
1692  {
1693  xRand = fRMin * cosphi;
1694  yRand = fRMin * sinphi;
1695  zRand = UUtils::Random(-1.*fDz, fDz);
1696  return UVector3(xRand, yRand, zRand);
1697  }
1698  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr))
1699  {
1700  xRand = rRand * cosphi;
1701  yRand = rRand * sinphi;
1702  zRand = fDz;
1703  return UVector3(xRand, yRand, zRand);
1704  }
1705  else if ((chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr))
1706  {
1707  xRand = rRand * cosphi;
1708  yRand = rRand * sinphi;
1709  zRand = -1.*fDz;
1710  return UVector3(xRand, yRand, zRand);
1711  }
1712  else if ((chose >= aOne + aTwo + 2.*aThr)
1713  && (chose < aOne + aTwo + 2.*aThr + aFou))
1714  {
1715  xRand = rRand * fCosSPhi;
1716  yRand = rRand * fSinSPhi;
1717  zRand = UUtils::Random(-1.*fDz, fDz);
1718  return UVector3(xRand, yRand, zRand);
1719  }
1720  else
1721  {
1722  xRand = rRand * fCosSPhiDPhi;
1723  yRand = rRand * fSinSPhiDPhi;
1724  zRand = UUtils::Random(-1.*fDz, fDz);
1725  return UVector3(xRand, yRand, zRand);
1726  }
1727 }
1728 
1729 void UTubs::Extent(UVector3& aMin, UVector3& aMax) const
1730 {
1731  aMin = UVector3(-fRMax, -fRMax, -fDz);
1732  aMax = UVector3(fRMax, fRMax, fDz);
1733 }
1734 
1735 void UTubs::GetParametersList(int, double* aArray) const
1736 {
1737  aArray[0] = GetInnerRadius();
1738  aArray[1] = GetOuterRadius();
1739  aArray[2] = GetZHalfLength();
1740  aArray[3] = GetStartPhiAngle();
1741  aArray[4] = GetDeltaPhiAngle();
1742 }
double fCosSPhiDPhi
Definition: UTubs.hh:174
Definition: UTubs.hh:47
double fCosEPhi
Definition: UTubs.hh:174
double fSinSPhiDPhi
Definition: UTubs.hh:174
double GetStartPhiAngle() const
double fCosHDPhiIT
Definition: UTubs.hh:174
TTree * t1
Definition: plottest35.C:26
static double frTolerance
Definition: VUSolid.hh:31
double GetZHalfLength() const
double SafetyFromOutside(const UVector3 &p, bool precise=false) const
Definition: UTubs.cc:1030
double fSPhi
Definition: UTubs.hh:170
double SafetyFromInside(const UVector3 &p, bool precise=false) const
Definition: UTubs.cc:1543
Float_t d
Definition: plot.C:237
virtual void GetParametersList(int, double *) const
Definition: UTubs.cc:1735
const std::string & GetName() const
Definition: VUSolid.hh:103
ifstream in
Definition: comparison.C:7
double GetInnerRadius() const
const char * p
Definition: xmltok.h:285
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: UTubs.cc:153
Float_t norm
double fSinCPhi
Definition: UTubs.hh:174
VUSolid * Clone() const
Definition: UTubs.cc:1628
static double Tolerance()
Definition: VUSolid.hh:127
double x
Definition: UVector3.hh:136
double fDPhi
Definition: UTubs.hh:170
void CheckPhiAngles(double sPhi, double dPhi)
double fDz
Definition: UTubs.hh:170
tuple b
Definition: test.py:12
Char_t n[5]
double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &n, bool &validNorm, double aPstep=UUtils::kInfinity) const
Definition: UTubs.cc:1090
double GetDeltaPhiAngle() const
UGeometryType GetEntityType() const
Definition: UTubs.cc:1619
double fCosSPhi
Definition: UTubs.hh:174
bool Normal(const UVector3 &p, UVector3 &normal) const
Definition: UTubs.cc:356
virtual UVector3 ApproxSurfaceNormal(const UVector3 &p) const
Definition: UTubs.cc:473
UTubs & operator=(const UTubs &rhs)
Definition: UTubs.cc:108
double fRMin
Definition: UTubs.hh:170
double fSinSPhi
Definition: UTubs.hh:174
static double faTolerance
Definition: VUSolid.hh:32
bool fPhiFullTube
Definition: UTubs.hh:179
double kAngTolerance
Definition: UTubs.hh:166
double fCubicVolume
Definition: UTubs.hh:157
ENorm
Definition: UTubs.hh:164
EnumInside
Definition: VUSolid.hh:23
UVector3 GetPointOnSurface() const
Definition: UTubs.cc:1660
tuple v
Definition: test.py:18
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UTubs.cc:1729
virtual ~UTubs()
Definition: UTubs.cc:81
double GetRadiusInRing(double rmin, double rmax)
Definition: UUtils.hh:160
std::ostream & StreamInfo(std::ostream &os) const
Definition: UTubs.cc:1637
double fRMax
Definition: UTubs.hh:170
UVector3 Unit() const
Definition: UVector3.cc:80
double GetOuterRadius() const
std::string UGeometryType
Definition: UTypes.hh:70
double kRadTolerance
Definition: UTubs.hh:166
TTree * t2
Definition: plottest35.C:36
UTubs()
Definition: UTubs.cc:66
ESide
Definition: UTubs.hh:160
double fCosCPhi
Definition: UTubs.hh:174
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
double fSinEPhi
Definition: UTubs.hh:174
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UTubs.cc:614
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double y
Definition: UVector3.hh:137
double fSurfaceArea
Definition: UTubs.hh:157
double fCosHDPhiOT
Definition: UTubs.hh:174