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