Geant4  10.01
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 aAccurate) const
1031 {
1032  double safe = 0.0, rho, safe1, safe2, safe3;
1033  double safePhi;
1034  bool outside;
1035 
1036  rho = std::sqrt(p.x * p.x + p.y * p.y);
1037  safe1 = fRMin - rho;
1038  safe2 = rho - fRMax;
1039  safe3 = std::fabs(p.z) - fDz;
1040 
1041  if (safe1 > safe2)
1042  {
1043  safe = safe1;
1044  }
1045  else
1046  {
1047  safe = safe2;
1048  }
1049  if (safe3 > safe)
1050  {
1051  safe = safe3;
1052  }
1053 
1054  if ((!fPhiFullTube) && (rho))
1055  {
1056  safePhi = SafetyToPhi(p,rho,outside);
1057  if ((outside) && (safePhi > safe))
1058  {
1059  safe = safePhi;
1060  }
1061  }
1062 
1063  if (safe < 0)
1064  {
1065  safe = 0; return safe; // point is Inside;
1066  }
1067  if (!aAccurate) return safe;
1068  double safsq = 0.0;
1069  int count = 0;
1070  if (safe1 > 0)
1071  {
1072  safsq += safe1 * safe1;
1073  count++;
1074  }
1075  if (safe2 > 0)
1076  {
1077  safsq += safe2 * safe2;
1078  count++;
1079  }
1080  if (safe3 > 0)
1081  {
1082  safsq += safe3 * safe3;
1083  count++;
1084  }
1085  if (count == 1) return safe;
1086  return std::sqrt(safsq);
1087 }
1088 
1090 //
1091 // Calculate distance to surface of shape from `inside', allowing for tolerance
1092 // - Only Calc rmax intersection if no valid rmin intersection
1093 
1094 // double UTubs::DistanceToOut( const UVector3& p, const UVector3& v, const bool calcNorm, bool *validNorm, UVector3 *n ) const
1095 double UTubs::DistanceToOut(const UVector3& p, const UVector3& v, UVector3& n, bool& validNorm, double) const
1096 {
1097  ESide side = kNull , sider = kNull, sidephi = kNull;
1098  double snxt, srd = UUtils::kInfinity, sphi = UUtils::kInfinity, pdist;
1099  double deltaR, t1, t2, t3, b, c, d2, roMin2;
1100 
1101  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
1102  static const double halfAngTolerance = kAngTolerance * 0.5;
1103 
1104  // Vars for phi intersection:
1105 
1106  double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2;
1107 
1108  // Z plane intersection
1109 
1110  if (v.z > 0)
1111  {
1112  pdist = fDz - p.z;
1113  if (pdist > halfCarTolerance)
1114  {
1115  snxt = pdist / v.z;
1116  side = kPZ;
1117  }
1118  else
1119  {
1120  n = UVector3(0, 0, 1);
1121  validNorm = true;
1122  return snxt = 0;
1123  }
1124  }
1125  else if (v.z < 0)
1126  {
1127  pdist = fDz + p.z;
1128 
1129  if (pdist > halfCarTolerance)
1130  {
1131  snxt = -pdist / v.z;
1132  side = kMZ;
1133  }
1134  else
1135  {
1136  n = UVector3(0, 0, -1);
1137  validNorm = true;
1138  return snxt = 0.0;
1139  }
1140  }
1141  else
1142  {
1143  snxt = UUtils::kInfinity; // Travel perpendicular to z axis
1144  side = kNull;
1145  }
1146 
1147  // Radial Intersections
1148  //
1149  // Find intersection with cylinders at rmax/rmin
1150  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1151  //
1152  // Intersects with x^2+y^2=R^2
1153  //
1154  // 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
1155  //
1156  // t1 t2 t3
1157 
1158  t1 = 1.0 - v.z * v.z; // since v normalised
1159  t2 = p.x * v.x + p.y * v.y;
1160  t3 = p.x * p.x + p.y * p.y;
1161 
1162  if (snxt > 10 * (fDz + fRMax))
1163  {
1164  roi2 = 2 * fRMax * fRMax;
1165  }
1166  else
1167  {
1168  roi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3; // radius^2 on +-fDz
1169  }
1170 
1171  if (t1 > 0) // Check not parallel
1172  {
1173  // Calculate srd, r exit distance
1174 
1175  if ((t2 >= 0.0) && (roi2 > fRMax * (fRMax + kRadTolerance)))
1176  {
1177  // Delta r not negative => leaving via rmax
1178 
1179  deltaR = t3 - fRMax * fRMax;
1180 
1181  // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1182  // - avoid sqrt for efficiency
1183 
1184  if (deltaR < -kRadTolerance * fRMax)
1185  {
1186  b = t2 / t1;
1187  c = deltaR / t1;
1188  d2 = b * b - c;
1189  if (d2 >= 0)
1190  {
1191  srd = c / (-b - std::sqrt(d2));
1192  }
1193  else
1194  {
1195  srd = 0.;
1196  }
1197  sider = kRMax;
1198  }
1199  else
1200  {
1201  // On tolerant boundary & heading outwards (or perpendicular to)
1202  // outer radial surface -> leaving immediately
1203  n = UVector3(p.x / fRMax, p.y / fRMax, 0);
1204  validNorm = true;
1205  return snxt = 0; // Leaving by rmax immediately
1206  }
1207  }
1208  else if (t2 < 0.) // i.e. t2 < 0; Possible rmin intersection
1209  {
1210  roMin2 = t3 - t2 * t2 / t1; // min ro2 of the plane of movement
1211 
1212  if (fRMin && (roMin2 < fRMin * (fRMin - kRadTolerance)))
1213  {
1214  deltaR = t3 - fRMin * fRMin;
1215  b = t2 / t1;
1216  c = deltaR / t1;
1217  d2 = b * b - c;
1218 
1219  if (d2 >= 0) // Leaving via rmin
1220  {
1221  // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1222  // - avoid sqrt for efficiency
1223 
1224  if (deltaR > kRadTolerance * fRMin)
1225  {
1226  srd = c / (-b + std::sqrt(d2));
1227  sider = kRMin;
1228  }
1229  else
1230  {
1231  validNorm = false; // Concave side
1232  n = UVector3(-p.x / fRMin, -p.y / fRMin, 0);
1233  return snxt = 0.0;
1234  }
1235  }
1236  else // No rmin intersect -> must be rmax intersect
1237  {
1238  deltaR = t3 - fRMax * fRMax;
1239  c = deltaR / t1;
1240  d2 = b * b - c;
1241  if (d2 >= 0.)
1242  {
1243  srd = -b + std::sqrt(d2);
1244  sider = kRMax;
1245  }
1246  else // Case: On the border+t2<kRadTolerance
1247  // (v is perpendicular to the surface)
1248  {
1249  n = UVector3(p.x / fRMax, p.y / fRMax, 0);
1250  validNorm = true;
1251  return snxt = 0.0;
1252  }
1253  }
1254  }
1255  else if (roi2 > fRMax * (fRMax + kRadTolerance))
1256  // No rmin intersect -> must be rmax intersect
1257  {
1258  deltaR = t3 - fRMax * fRMax;
1259  b = t2 / t1;
1260  c = deltaR / t1;
1261  d2 = b * b - c;
1262  if (d2 >= 0)
1263  {
1264  srd = -b + std::sqrt(d2);
1265  sider = kRMax;
1266  }
1267  else // Case: On the border+t2<kRadTolerance
1268  // (v is perpendicular to the surface)
1269  {
1270  n = UVector3(p.x / fRMax, p.y / fRMax, 0);
1271  validNorm = true;
1272  return snxt = 0.0;
1273  }
1274  }
1275  }
1276 
1277  // Phi Intersection
1278 
1279  if (!fPhiFullTube)
1280  {
1281  // add angle calculation with correction
1282  // of the difference in domain of atan2 and Sphi
1283  //
1284  vphi = std::atan2(v.y, v.x);
1285 
1286  if (vphi < fSPhi - halfAngTolerance)
1287  {
1288  vphi += 2 * UUtils::kPi;
1289  }
1290  else if (vphi > fSPhi + fDPhi + halfAngTolerance)
1291  {
1292  vphi -= 2 * UUtils::kPi;
1293  }
1294 
1295 
1296  if (p.x || p.y) // Check if on z axis (rho not needed later)
1297  {
1298  // pDist -ve when inside
1299 
1300  pDistS = p.x * fSinSPhi - p.y * fCosSPhi;
1301  pDistE = -p.x * fSinEPhi + p.y * fCosEPhi;
1302 
1303  // Comp -ve when in direction of outwards normal
1304 
1305  compS = -fSinSPhi * v.x + fCosSPhi * v.y;
1306  compE = fSinEPhi * v.x - fCosEPhi * v.y;
1307 
1308  sidephi = kNull;
1309 
1310  if (((fDPhi <= UUtils::kPi) && ((pDistS <= halfCarTolerance)
1311  && (pDistE <= halfCarTolerance)))
1312  || ((fDPhi > UUtils::kPi) && !((pDistS > halfCarTolerance)
1313  && (pDistE > halfCarTolerance))))
1314  {
1315  // Inside both phi *full* planes
1316 
1317  if (compS < 0)
1318  {
1319  sphi = pDistS / compS;
1320 
1321  if (sphi >= -halfCarTolerance)
1322  {
1323  xi = p.x + sphi * v.x;
1324  yi = p.y + sphi * v.y;
1325 
1326  // Check intersecting with correct half-plane
1327  // (if not -> no intersect)
1328  //
1329  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
1330  {
1331  sidephi = kSPhi;
1332  if (((fSPhi - halfAngTolerance) <= vphi)
1333  && ((fSPhi + fDPhi + halfAngTolerance) >= vphi))
1334  {
1335  sphi = UUtils::kInfinity;
1336  }
1337  }
1338  else if (yi * fCosCPhi - xi * fSinCPhi >= 0)
1339  {
1340  sphi = UUtils::kInfinity;
1341  }
1342  else
1343  {
1344  sidephi = kSPhi;
1345  if (pDistS > -halfCarTolerance)
1346  {
1347  sphi = 0.0; // Leave by sphi immediately
1348  }
1349  }
1350  }
1351  else
1352  {
1353  sphi = UUtils::kInfinity;
1354  }
1355  }
1356  else
1357  {
1358  sphi = UUtils::kInfinity;
1359  }
1360 
1361  if (compE < 0)
1362  {
1363  sphi2 = pDistE / compE;
1364 
1365  // Only check further if < starting phi intersection
1366  //
1367  if ((sphi2 > -halfCarTolerance) && (sphi2 < sphi))
1368  {
1369  xi = p.x + sphi2 * v.x;
1370  yi = p.y + sphi2 * v.y;
1371 
1372  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
1373  {
1374  // Leaving via ending phi
1375  //
1376  if (!((fSPhi - halfAngTolerance <= vphi)
1377  && (fSPhi + fDPhi + halfAngTolerance >= vphi)))
1378  {
1379  sidephi = kEPhi;
1380  if (pDistE <= -halfCarTolerance)
1381  {
1382  sphi = sphi2;
1383  }
1384  else
1385  {
1386  sphi = 0.0;
1387  }
1388  }
1389  }
1390  else // Check intersecting with correct half-plane
1391 
1392  if ((yi * fCosCPhi - xi * fSinCPhi) >= 0)
1393  {
1394  // Leaving via ending phi
1395  //
1396  sidephi = kEPhi;
1397  if (pDistE <= -halfCarTolerance)
1398  {
1399  sphi = sphi2;
1400  }
1401  else
1402  {
1403  sphi = 0.0;
1404  }
1405  }
1406  }
1407  }
1408  }
1409  else
1410  {
1411  sphi = UUtils::kInfinity;
1412  }
1413  }
1414  else
1415  {
1416  // On z axis + travel not || to z axis -> if phi of vector direction
1417  // within phi of shape, Step limited by rmax, else Step =0
1418 
1419  if ((fSPhi - halfAngTolerance <= vphi)
1420  && (vphi <= fSPhi + fDPhi + halfAngTolerance))
1421  {
1422  sphi = UUtils::kInfinity;
1423  }
1424  else
1425  {
1426  sidephi = kSPhi; // arbitrary
1427  sphi = 0.0;
1428  }
1429  }
1430  if (sphi < snxt) // Order intersecttions
1431  {
1432  snxt = sphi;
1433  side = sidephi;
1434  }
1435  }
1436  if (srd < snxt) // Order intersections
1437  {
1438  snxt = srd;
1439  side = sider;
1440  }
1441  }
1442  // if (calcNorm)
1443  {
1444  switch (side)
1445  {
1446  case kRMax:
1447  // Note: returned vector not normalised
1448  // (divide by fRMax for Unit vector)
1449  //
1450  xi = p.x + snxt * v.x;
1451  yi = p.y + snxt * v.y;
1452  n = UVector3(xi / fRMax, yi / fRMax, 0);
1453  validNorm = true;
1454  break;
1455 
1456  case kRMin:
1457  xi = p.x + snxt * v.x;
1458  yi = p.y + snxt * v.y;
1459  n = UVector3(-xi / fRMin, -yi / fRMin, 0);
1460  validNorm = false; // Rmin is inconvex
1461  break;
1462 
1463  case kSPhi:
1464  if (fDPhi <= UUtils::kPi)
1465  {
1466  n = UVector3(fSinSPhi, -fCosSPhi, 0);
1467  validNorm = true;
1468  }
1469  else
1470  {
1471  n = UVector3(fSinSPhi, -fCosSPhi, 0);
1472  validNorm = false;
1473  }
1474  break;
1475 
1476  case kEPhi:
1477  if (fDPhi <= UUtils::kPi)
1478  {
1479  n = UVector3(-fSinEPhi, fCosEPhi, 0);
1480  validNorm = true;
1481  }
1482  else
1483  {
1484  n = UVector3(-fSinEPhi, fCosEPhi, 0);
1485  validNorm = false;
1486  }
1487  break;
1488 
1489  case kPZ:
1490  n = UVector3(0, 0, 1);
1491  validNorm = true;
1492  break;
1493 
1494  case kMZ:
1495  n = UVector3(0, 0, -1);
1496  validNorm = true;
1497  break;
1498 
1499  default:
1500  cout << std::endl;
1501  // DumpInfo();
1502  std::ostringstream message;
1503  int oldprc = message.precision(16);
1504  message << "Undefined side for valid surface normal to solid."
1505  << std::endl
1506  << "Position:" << std::endl << std::endl
1507  << "p.x = " << p.x << " mm" << std::endl
1508  << "p.y = " << p.y << " mm" << std::endl
1509  << "p.z = " << p.z << " mm" << std::endl << std::endl
1510  << "Direction:" << std::endl << std::endl
1511  << "v.x = " << v.x << std::endl
1512  << "v.y = " << v.y << std::endl
1513  << "v.z = " << v.z << std::endl << std::endl
1514  << "Proposed distance :" << std::endl << std::endl
1515  << "snxt = " << snxt << " mm" << std::endl;
1516  message.precision(oldprc);
1517  UUtils::Exception("UTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1518  Warning, 1, message.str().c_str());
1519  break;
1520  }
1521  }
1522  if (snxt < halfCarTolerance)
1523  {
1524  snxt = 0;
1525  }
1526 
1527  return snxt;
1528 }
1529 
1531 //
1532 // Calculate distance (<=actual) to closest surface of shape from inside
1533 
1534 double UTubs::SafetyFromInside(const UVector3& p, bool) const
1535 {
1536  double safe = 0.0, rho, safeZ;
1537  rho = std::sqrt(p.x * p.x + p.y * p.y);
1538 
1539 #ifdef UDEBUG
1540  if (Inside(p) == eOutside)
1541  {
1542  int oldprc = cout.precision(16);
1543  cout << std::endl;
1544  DumpInfo();
1545  cout << "Position:" << std::endl << std::endl;
1546  cout << "p.x = " << p.x << " mm" << std::endl;
1547  cout << "p.y = " << p.y << " mm" << std::endl;
1548  cout << "p.z = " << p.z << " mm" << std::endl << std::endl;
1549  cout.precision(oldprc);
1550  UUtils::Exception("UTubs::DistanceToOut(p)", "GeomSolids1002",
1551  Warning, 1, "Point p is outside !?");
1552  }
1553 #endif
1554  safe = SafetyFromInsideR(p, rho);
1555  safeZ = fDz - std::fabs(p.z);
1556 
1557  if (safeZ < safe)
1558  {
1559  safe = safeZ;
1560  }
1561 
1562  if (safe < 0)
1563  {
1564  safe = 0;
1565  }
1566 
1567  return safe;
1568 }
1569 
1571 //
1572 // Stream object contents to an output stream
1573 
1575 {
1576  return std::string("Tubs");
1577 }
1578 
1580 //
1581 // Make a clone of the object
1582 //
1584 {
1585  return new UTubs(*this);
1586 }
1587 
1589 //
1590 // Stream object contents to an output stream
1591 
1592 std::ostream& UTubs::StreamInfo(std::ostream& os) const
1593 {
1594  int oldprc = os.precision(16);
1595  os << "-----------------------------------------------------------\n"
1596  << " *** Dump for solid - " << GetName() << " ***\n"
1597  << " ===================================================\n"
1598  << " Solid type: UTubs\n"
1599  << " Parameters: \n"
1600  << " inner radius : " << fRMin << " mm \n"
1601  << " outer radius : " << fRMax << " mm \n"
1602  << " half length Z: " << fDz << " mm \n"
1603  << " starting phi : " << fSPhi / (UUtils::kPi / 180.0) << " degrees \n"
1604  << " delta phi : " << fDPhi / (UUtils::kPi / 180.0) << " degrees \n"
1605  << "-----------------------------------------------------------\n";
1606  os.precision(oldprc);
1607 
1608  return os;
1609 }
1610 
1612 //
1613 // GetPointOnSurface
1614 
1616 {
1617  double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1618  aOne, aTwo, aThr, aFou;
1619  double rRand;
1620 
1621  aOne = 2.*fDz * fDPhi * fRMax;
1622  aTwo = 2.*fDz * fDPhi * fRMin;
1623  aThr = 0.5 * fDPhi * (fRMax * fRMax - fRMin * fRMin);
1624  aFou = 2.*fDz * (fRMax - fRMin);
1625 
1626  phi = UUtils::Random(fSPhi, fSPhi + fDPhi);
1627  cosphi = std::cos(phi);
1628  sinphi = std::sin(phi);
1629 
1630  rRand = UUtils::GetRadiusInRing(fRMin, fRMax);
1631 
1632  if ((fSPhi == 0) && (fDPhi == 2 * UUtils::kPi))
1633  {
1634  aFou = 0;
1635  }
1636 
1637  chose = UUtils::Random(0., aOne + aTwo + 2.*aThr + 2.*aFou);
1638 
1639  if ((chose >= 0) && (chose < aOne))
1640  {
1641  xRand = fRMax * cosphi;
1642  yRand = fRMax * sinphi;
1643  zRand = UUtils::Random(-1.*fDz, fDz);
1644  return UVector3(xRand, yRand, zRand);
1645  }
1646  else if ((chose >= aOne) && (chose < aOne + aTwo))
1647  {
1648  xRand = fRMin * cosphi;
1649  yRand = fRMin * sinphi;
1650  zRand = UUtils::Random(-1.*fDz, fDz);
1651  return UVector3(xRand, yRand, zRand);
1652  }
1653  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr))
1654  {
1655  xRand = rRand * cosphi;
1656  yRand = rRand * sinphi;
1657  zRand = fDz;
1658  return UVector3(xRand, yRand, zRand);
1659  }
1660  else if ((chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr))
1661  {
1662  xRand = rRand * cosphi;
1663  yRand = rRand * sinphi;
1664  zRand = -1.*fDz;
1665  return UVector3(xRand, yRand, zRand);
1666  }
1667  else if ((chose >= aOne + aTwo + 2.*aThr)
1668  && (chose < aOne + aTwo + 2.*aThr + aFou))
1669  {
1670  xRand = rRand * fCosSPhi;
1671  yRand = rRand * fSinSPhi;
1672  zRand = UUtils::Random(-1.*fDz, fDz);
1673  return UVector3(xRand, yRand, zRand);
1674  }
1675  else
1676  {
1677  xRand = rRand * fCosSPhiDPhi;
1678  yRand = rRand * fSinSPhiDPhi;
1679  zRand = UUtils::Random(-1.*fDz, fDz);
1680  return UVector3(xRand, yRand, zRand);
1681  }
1682 }
1683 
1684 void UTubs::Extent(UVector3& aMin, UVector3& aMax) const
1685 {
1686  aMin = UVector3(-fRMax, -fRMax, -fDz);
1687  aMax = UVector3(fRMax, fRMax, fDz);
1688 }
1689 
1690 void UTubs::GetParametersList(int, double* aArray) const
1691 {
1692  aArray[0] = GetInnerRadius();
1693  aArray[1] = GetOuterRadius();
1694  aArray[2] = GetZHalfLength();
1695  aArray[3] = GetStartPhiAngle();
1696  aArray[4] = GetDeltaPhiAngle();
1697 }
double fCosSPhiDPhi
Definition: UTubs.hh:180
Definition: UTubs.hh:47
double fCosEPhi
Definition: UTubs.hh:180
double fSinSPhiDPhi
Definition: UTubs.hh:180
double GetStartPhiAngle() const
double fCosHDPhiIT
Definition: UTubs.hh:180
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:176
double SafetyFromInside(const UVector3 &p, bool precise=false) const
Definition: UTubs.cc:1534
virtual void GetParametersList(int, double *) const
Definition: UTubs.cc:1690
const std::string & GetName() const
Definition: VUSolid.hh:103
double GetInnerRadius() const
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: UTubs.cc:153
double fSinCPhi
Definition: UTubs.hh:180
VUSolid * Clone() const
Definition: UTubs.cc:1583
static double Tolerance()
Definition: VUSolid.hh:127
double x
Definition: UVector3.hh:136
double fDPhi
Definition: UTubs.hh:176
void CheckPhiAngles(double sPhi, double dPhi)
double fDz
Definition: UTubs.hh:176
static const double kInfinity
Definition: UUtils.hh:53
double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &n, bool &validNorm, double aPstep=UUtils::kInfinity) const
Definition: UTubs.cc:1095
double GetDeltaPhiAngle() const
UGeometryType GetEntityType() const
Definition: UTubs.cc:1574
double fCosSPhi
Definition: UTubs.hh:180
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:176
double fSinSPhi
Definition: UTubs.hh:180
static double faTolerance
Definition: VUSolid.hh:32
bool fPhiFullTube
Definition: UTubs.hh:185
double kAngTolerance
Definition: UTubs.hh:172
double fCubicVolume
Definition: UTubs.hh:163
ENorm
Definition: UTubs.hh:170
EnumInside
Definition: VUSolid.hh:23
const G4int n
UVector3 GetPointOnSurface() const
Definition: UTubs.cc:1615
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UTubs.cc:1684
virtual ~UTubs()
Definition: UTubs.cc:81
double GetRadiusInRing(double rmin, double rmax)
Definition: UUtils.hh:155
std::ostream & StreamInfo(std::ostream &os) const
Definition: UTubs.cc:1592
double fRMax
Definition: UTubs.hh:176
static const double kPi
Definition: UUtils.hh:48
UVector3 Unit() const
Definition: UVector3.cc:80
double GetOuterRadius() const
std::string UGeometryType
Definition: UTypes.hh:39
double kRadTolerance
Definition: UTubs.hh:172
UTubs()
Definition: UTubs.cc:66
ESide
Definition: UTubs.hh:166
double fCosCPhi
Definition: UTubs.hh:180
double SafetyToPhi(const UVector3 &p, const double rho, bool &outside) const
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
double z
Definition: UVector3.hh:138
static const G4double d2
double fSinEPhi
Definition: UTubs.hh:180
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 SafetyFromInsideR(const UVector3 &p, const double rho, bool precise=false) const
double y
Definition: UVector3.hh:137
double fSurfaceArea
Definition: UTubs.hh:163
double fCosHDPhiOT
Definition: UTubs.hh:180