Geant4  10.00.p03
UCons.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 // UCons
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 "UCons.hh"
22 
23 using namespace std;
24 
26 //
27 // Private enum: Not for external use - used by distanceToOut
28 
30 
31 // used by normal
32 
34 
36 //
37 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
38 // - note if pDPhi>2PI then reset to 2PI
39 
40 UCons::UCons(const std::string& pName,
41  double pRmin1, double pRmax1,
42  double pRmin2, double pRmax2,
43  double pDz,
44  double pSPhi, double pDPhi)
45  : VUSolid(pName.c_str()), fRmin1(pRmin1), fRmin2(pRmin2),
46  fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
47 {
50  // Check z-len
51  //
52  if (pDz < 0)
53  {
54  std::ostringstream message;
55  message << "Invalid Z half-length for Solid: " << GetName() << std::endl
56  << " hZ = " << pDz;
57  UUtils::Exception("UCons::UCons()", "UGeomSolids", FatalErrorInArguments, 1, message.str().c_str());
58  }
59 
60  // Check radii
61  //
62  if (((pRmin1 >= pRmax1) || (pRmin2 >= pRmax2) || (pRmin1 < 0)) && (pRmin2 < 0))
63  {
64  std::ostringstream message;
65  message << "Invalid values of radii for Solid: " << GetName() << std::endl
66  << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
67  << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
68  UUtils::Exception("UCons::UCons()", "UGeomSolids", FatalErrorInArguments, 1, message.str().c_str());
69 
70  }
71  if ((pRmin1 == 0.0) && (pRmin2 > 0.0))
72  {
73  fRmin1 = 1e3 * kRadTolerance;
74  }
75  if ((pRmin2 == 0.0) && (pRmin1 > 0.0))
76  {
77  fRmin2 = 1e3 * kRadTolerance;
78  }
79 
80  // Check angles
81  //
82  CheckPhiAngles(pSPhi, pDPhi);
83 
84  Initialize();
85 }
86 
88 //
89 // Fake default constructor - sets only member data and allocates memory
90 // for usage restricted to object persistency.
91 //
92 UCons::UCons(/* __void__& a */)
93  : VUSolid(""), kRadTolerance(0.), kAngTolerance(0.),
94  fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
95  fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
96  cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
97  fPhiFullCone(false)
98 {
99  Initialize();
100 }
101 
103 //
104 // Destructor
105 
107 {
108 }
109 
111 //
112 // Copy constructor
113 
114 UCons::UCons(const UCons& rhs)
115  : VUSolid(rhs), kRadTolerance(rhs.kRadTolerance),
116  kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
117  fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
118  fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
119  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
120  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
121  cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone)
122 {
123  Initialize();
124 }
125 
127 //
128 // Assignment operator
129 
131 {
132  // Check assignment to self
133  //
134  if (this == &rhs)
135  {
136  return *this;
137  }
138 
139  // Copy base class data
140  //
141  VUSolid::operator=(rhs);
142 
143  // Copy data
144  //
147  fRmin1 = rhs.fRmin1;
148  fRmin2 = rhs.fRmin2;
149  fRmax1 = rhs.fRmax1;
150  fRmax2 = rhs.fRmax2;
151  fDz = rhs.fDz;
152  fSPhi = rhs.fSPhi;
153  fDPhi = rhs.fDPhi;
154  sinCPhi = rhs.sinCPhi;
155  cosCPhi = rhs.cosCPhi;
156  cosHDPhiOT = rhs.cosHDPhiOT;
157  cosHDPhiIT = rhs.cosHDPhiIT;
158  sinSPhi = rhs.sinSPhi;
159  cosSPhi = rhs.cosSPhi;
160  sinEPhi = rhs.sinEPhi;
161  cosEPhi = rhs.cosEPhi;
163 
164  Initialize();
165  return *this;
166 }
167 
169 //
170 // Return Unit normal of surface closest to p
171 // - note if point on z axis, ignore phi divided sides
172 // - unsafe if point close to z axis a rmin=0 - no explicit checks
173 
174 bool UCons::Normal(const UVector3& p, UVector3& n) const
175 {
176  int noSurfaces = 0;
177  double rho, pPhi;
178  double distZ, distRMin, distRMax;
179  double distSPhi = UUtils::kInfinity, distEPhi = UUtils::kInfinity;
180  double pRMin, widRMin;
181  double pRMax, widRMax;
182 
183  static const double delta = 0.5 * VUSolid::Tolerance();
184  static const double dAngle = 0.5 * kAngTolerance;
185 
186  UVector3 norm, sumnorm(0., 0., 0.), nZ = UVector3(0., 0., 1.);
187  UVector3 nR, nr(0., 0., 0.), nPs, nPe;
188 
189  distZ = std::fabs(std::fabs(p.z) - fDz);
190  rho = std::sqrt(p.x * p.x + p.y * p.y);
191 
192  pRMin = rho - p.z * tanRMin;
193  widRMin = fRmin2 - fDz * tanRMin;
194  distRMin = std::fabs(pRMin - widRMin) / secRMin;
195 
196  pRMax = rho - p.z * tanRMax;
197  widRMax = fRmax2 - fDz * tanRMax;
198  distRMax = std::fabs(pRMax - widRMax) / secRMax;
199 
200  if (!fPhiFullCone) // Protected against (0,0,z)
201  {
202  if (rho)
203  {
204  pPhi = std::atan2(p.y, p.x);
205 
206  if (pPhi < fSPhi - delta)
207  {
208  pPhi += 2 * UUtils::kPi;
209  }
210  else if (pPhi > fSPhi + fDPhi + delta)
211  {
212  pPhi -= 2 * UUtils::kPi;
213  }
214 
215  distSPhi = std::fabs(pPhi - fSPhi);
216  distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
217  }
218  else if (!(fRmin1) || !(fRmin2))
219  {
220  distSPhi = 0.;
221  distEPhi = 0.;
222  }
223  nPs = UVector3(std::sin(fSPhi), -std::cos(fSPhi), 0);
224  nPe = UVector3(-std::sin(fSPhi + fDPhi), std::cos(fSPhi + fDPhi), 0);
225  }
226  if (rho > delta)
227  {
228  nR = UVector3(p.x / rho / secRMax, p.y / rho / secRMax, -tanRMax / secRMax);
229  if (fRmin1 || fRmin2)
230  {
231  nr = UVector3(-p.x / rho / secRMin, -p.y / rho / secRMin, tanRMin / secRMin);
232  }
233  }
234 
235  if (distRMax <= delta)
236  {
237  noSurfaces ++;
238  sumnorm += nR;
239  }
240  if ((fRmin1 || fRmin2) && (distRMin <= delta))
241  {
242  noSurfaces ++;
243  sumnorm += nr;
244  }
245  if (!fPhiFullCone)
246  {
247  if (distSPhi <= dAngle)
248  {
249  noSurfaces ++;
250  sumnorm += nPs;
251  }
252  if (distEPhi <= dAngle)
253  {
254  noSurfaces ++;
255  sumnorm += nPe;
256  }
257  }
258  if (distZ <= delta)
259  {
260  noSurfaces ++;
261  if (p.z >= 0.)
262  {
263  sumnorm += nZ;
264  }
265  else
266  {
267  sumnorm -= nZ;
268  }
269  }
270  if (noSurfaces == 0)
271  {
272 #ifdef UDEBUG
273 
274  UUtils::Exception("UCons::SurfaceNormal(p)", "GeomSolids1002",
275  Warning, 1, "Point p is not on surface !?");
276 #endif
277  norm = ApproxSurfaceNormal(p);
278  }
279  else if (noSurfaces == 1)
280  {
281  norm = sumnorm;
282  }
283  else
284  {
285  norm = sumnorm.Unit();
286  }
287 
288  n = norm;
289 
290  return (bool) noSurfaces;
291 }
292 
294 //
295 // Algorithm for SurfaceNormal() following the original specification
296 // for points not on the surface
297 
299 {
300  ENorm side;
301  UVector3 norm;
302  double rho, phi;
303  double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin;
304  double pRMin, widRMin;
305  double pRMax, widRMax;
306 
307  distZ = std::fabs(std::fabs(p.z) - fDz);
308  rho = std::sqrt(p.x * p.x + p.y * p.y);
309 
310  pRMin = rho - p.z * tanRMin;
311  widRMin = fRmin2 - fDz * tanRMin;
312  distRMin = std::fabs(pRMin - widRMin) / secRMin;
313 
314  pRMax = rho - p.z * tanRMax;
315  widRMax = fRmax2 - fDz * tanRMax;
316  distRMax = std::fabs(pRMax - widRMax) / secRMax;
317 
318  if (distRMin < distRMax) // First minimum
319  {
320  if (distZ < distRMin)
321  {
322  distMin = distZ;
323  side = kNZ;
324  }
325  else
326  {
327  distMin = distRMin;
328  side = kNRMin;
329  }
330  }
331  else
332  {
333  if (distZ < distRMax)
334  {
335  distMin = distZ;
336  side = kNZ;
337  }
338  else
339  {
340  distMin = distRMax;
341  side = kNRMax;
342  }
343  }
344  if (!fPhiFullCone && rho) // Protected against (0,0,z)
345  {
346  phi = std::atan2(p.y, p.x);
347 
348  if (phi < 0)
349  {
350  phi += 2 * UUtils::kPi;
351  }
352 
353  if (fSPhi < 0)
354  {
355  distSPhi = std::fabs(phi - (fSPhi + 2 * UUtils::kPi)) * rho;
356  }
357  else
358  {
359  distSPhi = std::fabs(phi - fSPhi) * rho;
360  }
361 
362  distEPhi = std::fabs(phi - fSPhi - fDPhi) * rho;
363 
364  // Find new minimum
365 
366  if (distSPhi < distEPhi)
367  {
368  if (distSPhi < distMin)
369  {
370  side = kNSPhi;
371  }
372  }
373  else
374  {
375  if (distEPhi < distMin)
376  {
377  side = kNEPhi;
378  }
379  }
380  }
381  switch (side)
382  {
383  case kNRMin: // Inner radius
384  rho *= secRMin;
385  norm = UVector3(-p.x / rho, -p.y / rho, tanRMin / secRMin);
386  break;
387  case kNRMax: // Outer radius
388  rho *= secRMax;
389  norm = UVector3(p.x / rho, p.y / rho, -tanRMax / secRMax);
390  break;
391  case kNZ: // +/- dz
392  if (p.z > 0)
393  {
394  norm = UVector3(0, 0, 1);
395  }
396  else
397  {
398  norm = UVector3(0, 0, -1);
399  }
400  break;
401  case kNSPhi:
402  norm = UVector3(std::sin(fSPhi), -std::cos(fSPhi), 0);
403  break;
404  case kNEPhi:
405  norm = UVector3(-std::sin(fSPhi + fDPhi), std::cos(fSPhi + fDPhi), 0);
406  break;
407  default: // Should never reach this case...
408  UUtils::Exception("UCons::ApproxSurfaceNormal()",
409  "GeomSolids1002", Warning, 1,
410  "Undefined side for valid surface normal to solid.");
411  break;
412  }
413  return norm;
414 }
415 
417 //
418 // Calculate distance to shape from outside, along normalised vector
419 // - return UUtils::kInfinity if no intersection, or intersection distance <= tolerance
420 //
421 // - Compute the intersection with the z planes
422 // - if at valid r, phi, return
423 //
424 // -> If point is outside cone, compute intersection with rmax1*0.5
425 // - if at valid phi,z return
426 // - if inside outer cone, handle case when on tolerant outer cone
427 // boundary and heading inwards(->0 to in)
428 //
429 // -> Compute intersection with inner cone, taking largest +ve root
430 // - if valid (in z,phi), save intersction
431 //
432 // -> If phi segmented, compute intersections with phi half planes
433 // - return smallest of valid phi intersections and
434 // inner radius intersection
435 //
436 // NOTE:
437 // - `if valid' implies tolerant checking of intersection points
438 // - z, phi intersection from Tubs
439 
441  const UVector3& v, double /* aPstep */) const
442 {
443  double snxt = UUtils::kInfinity;
444  const double dRmax = 100 * std::max(fRmax1, fRmax2);
445  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
446  static const double halfRadTolerance = kRadTolerance * 0.5;
447 
448  double rMaxAv, rMaxOAv; // Data for cones
449  double rMinAv, rMinOAv;
450  double rout, rin;
451 
452  double tolORMin, tolORMin2, tolIRMin, tolIRMin2; // `generous' radii squared
453  double tolORMax2, tolIRMax, tolIRMax2;
454  double tolODz, tolIDz;
455 
456  double Dist, sd, xi, yi, zi, ri = 0., risec, rhoi2, cosPsi; // Intersection point vars
457 
458  double t1, t2, t3, b, c, d; // Quadratic solver variables
459  double nt1, nt2, nt3;
460  double Comp;
461 
462  UVector3 norm;
463 
464  // Cone Precalcs
465  rMinAv = (fRmin1 + fRmin2) * 0.5;
466 
467  if (rMinAv > halfRadTolerance)
468  {
469  rMinOAv = rMinAv - halfRadTolerance;
470  }
471  else
472  {
473  rMinOAv = 0.0;
474  }
475  rMaxAv = (fRmax1 + fRmax2) * 0.5;
476  rMaxOAv = rMaxAv + halfRadTolerance;
477 
478  // Intersection with z-surfaces
479 
480  tolIDz = fDz - halfCarTolerance;
481  tolODz = fDz + halfCarTolerance;
482 
483  if (std::fabs(p.z) >= tolIDz)
484  {
485  if (p.z * v.z < 0) // at +Z going in -Z or visa versa
486  {
487  sd = (std::fabs(p.z) - fDz) / std::fabs(v.z); // Z intersect distance
488 
489  if (sd < 0.0)
490  {
491  sd = 0.0; // negative dist -> zero
492  }
493 
494  xi = p.x + sd * v.x; // Intersection coords
495  yi = p.y + sd * v.y;
496  rhoi2 = xi * xi + yi * yi ;
497 
498  // Check validity of intersection
499  // Calculate (outer) tolerant radi^2 at intersecion
500 
501  if (v.z > 0)
502  {
503  tolORMin = fRmin1 - halfRadTolerance * secRMin;
504  tolIRMin = fRmin1 + halfRadTolerance * secRMin;
505  tolIRMax = fRmax1 - halfRadTolerance * secRMin;
506  tolORMax2 = (fRmax1 + halfRadTolerance * secRMax) *
507  (fRmax1 + halfRadTolerance * secRMax);
508  }
509  else
510  {
511  tolORMin = fRmin2 - halfRadTolerance * secRMin;
512  tolIRMin = fRmin2 + halfRadTolerance * secRMin;
513  tolIRMax = fRmax2 - halfRadTolerance * secRMin;
514  tolORMax2 = (fRmax2 + halfRadTolerance * secRMax) *
515  (fRmax2 + halfRadTolerance * secRMax);
516  }
517  if (tolORMin > 0)
518  {
519  tolORMin2 = tolORMin * tolORMin;
520  tolIRMin2 = tolIRMin * tolIRMin;
521  }
522  else
523  {
524  tolORMin2 = 0.0;
525  tolIRMin2 = 0.0;
526  }
527  if (tolIRMax > 0)
528  {
529  tolIRMax2 = tolIRMax * tolIRMax;
530  }
531  else
532  {
533  tolIRMax2 = 0.0;
534  }
535 
536  if ((tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2))
537  {
538  if (!fPhiFullCone && rhoi2)
539  {
540  // Psi = angle made with central (average) phi of shape
541 
542  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
543 
544  if (cosPsi >= cosHDPhiIT)
545  {
546  return sd;
547  }
548  }
549  else
550  {
551  return sd;
552  }
553  }
554  }
555  else // On/outside extent, and heading away -> cannot intersect
556  {
557  return snxt;
558  }
559  }
560 
561 // ----> Can not intersect z surfaces
562 
563 
564 // Intersection with outer cone (possible return) and
565 // inner cone (must also check phi)
566 //
567 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
568 //
569 // Intersects with x^2+y^2=(a*z+b)^2
570 //
571 // where a=tanRMax or tanRMin
572 // b=rMaxAv or rMinAv
573 //
574 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0;
575 // t1 t2 t3
576 //
577 // \--------u-------/ \-----------v----------/ \---------w--------/
578 //
579 
580  t1 = 1.0 - v.z * v.z;
581  t2 = p.x * v.x + p.y * v.y;
582  t3 = p.x * p.x + p.y * p.y;
583  rin = tanRMin * p.z + rMinAv;
584  rout = tanRMax * p.z + rMaxAv;
585 
586  // Outer Cone Intersection
587  // Must be outside/on outer cone for valid intersection
588 
589  nt1 = t1 - (tanRMax * v.z) * (tanRMax * v.z);
590  nt2 = t2 - tanRMax * v.z * rout;
591  nt3 = t3 - rout * rout;
592 
593  if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
594  {
595  b = nt2 / nt1;
596  c = nt3 / nt1;
597  d = b * b - c ;
598  if ((nt3 > rout * rout * kRadTolerance * kRadTolerance * secRMax * secRMax)
599  || (rout < 0))
600  {
601  // If outside real cone (should be rho-rout>kRadTolerance*0.5
602  // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
603 
604  if (d >= 0)
605  {
606 
607  if ((rout < 0) && (nt3 <= 0))
608  {
609  // Inside `shadow cone' with -ve radius
610  // -> 2nd root could be on real cone
611 
612  if (b > 0)
613  {
614  sd = c / (-b - std::sqrt(d));
615  }
616  else
617  {
618  sd = -b + std::sqrt(d);
619  }
620  }
621  else
622  {
623  if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
624  {
625  sd = c / (-b + std::sqrt(d));
626  }
627  else
628  {
629  if (c <= 0) // second >=0
630  {
631  sd = -b + std::sqrt(d);
632  }
633  else // both negative, travel away
634  {
635  return UUtils::kInfinity;
636  }
637  }
638  }
639  if (sd >= 0) // If 'forwards'. Check z intersection
640  {
641  if (sd > dRmax) // Avoid rounding errors due to precision issues on
642  {
643  // 64 bits systems. Split long distances and recompute
644  double fTerm = sd - std::fmod(sd, dRmax);
645  sd = fTerm + DistanceToIn(p + fTerm * v, v);
646  }
647  zi = p.z + sd * v.z;
648 
649  if (std::fabs(zi) <= tolODz)
650  {
651  // Z ok. Check phi intersection if reqd
652 
653  if (fPhiFullCone)
654  {
655  return sd;
656  }
657  else
658  {
659  xi = p.x + sd * v.x;
660  yi = p.y + sd * v.y;
661  ri = rMaxAv + zi * tanRMax;
662  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
663 
664  if (cosPsi >= cosHDPhiIT)
665  {
666  return sd;
667  }
668  }
669  }
670  } // end if (sd>0)
671  }
672  }
673  else
674  {
675  // Inside outer cone
676  // check not inside, and heading through UCons (-> 0 to in)
677 
678  if ((t3 > (rin + halfRadTolerance * secRMin)*
679  (rin + halfRadTolerance * secRMin))
680  && (nt2 < 0) && (d >= 0) && (std::fabs(p.z) <= tolIDz))
681  {
682  // Inside cones, delta r -ve, inside z extent
683  // Point is on the Surface => check Direction using Normal.Dot(v)
684 
685  xi = p.x;
686  yi = p.y ;
687  risec = std::sqrt(xi * xi + yi * yi) * secRMax;
688  norm = UVector3(xi / risec, yi / risec, -tanRMax / secRMax);
689  if (!fPhiFullCone)
690  {
691  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(t3);
692  if (cosPsi >= cosHDPhiIT)
693  {
694  if (norm.Dot(v) <= 0)
695  {
696  return 0.0;
697  }
698  }
699  }
700  else
701  {
702  if (norm.Dot(v) <= 0)
703  {
704  return 0.0;
705  }
706  }
707  }
708  }
709  }
710  else // Single root case
711  {
712  if (std::fabs(nt2) > kRadTolerance)
713  {
714  sd = -0.5 * nt3 / nt2;
715 
716  if (sd < 0)
717  {
718  return UUtils::kInfinity; // travel away
719  }
720  else // sd >= 0, If 'forwards'. Check z intersection
721  {
722  zi = p.z + sd * v.z;
723 
724  if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
725  {
726  // Z ok. Check phi intersection if reqd
727 
728  if (fPhiFullCone)
729  {
730  return sd;
731  }
732  else
733  {
734  xi = p.x + sd * v.x;
735  yi = p.y + sd * v.y;
736  ri = rMaxAv + zi * tanRMax;
737  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
738 
739  if (cosPsi >= cosHDPhiIT)
740  {
741  return sd;
742  }
743  }
744  }
745  }
746  }
747  else // travel || cone surface from its origin
748  {
749  sd = UUtils::kInfinity;
750  }
751  }
752 
753  // Inner Cone Intersection
754  // o Space is divided into 3 areas:
755  // 1) Radius greater than real inner cone & imaginary cone & outside
756  // tolerance
757  // 2) Radius less than inner or imaginary cone & outside tolarance
758  // 3) Within tolerance of real or imaginary cones
759  // - Extra checks needed for 3's intersections
760  // => lots of duplicated code
761 
762  if (rMinAv)
763  {
764  nt1 = t1 - (tanRMin * v.z) * (tanRMin * v.z);
765  nt2 = t2 - tanRMin * v.z * rin;
766  nt3 = t3 - rin * rin;
767 
768  if (nt1)
769  {
770  if (nt3 > rin * kRadTolerance * secRMin)
771  {
772  // At radius greater than real & imaginary cones
773  // -> 2nd root, with zi check
774 
775  b = nt2 / nt1;
776  c = nt3 / nt1;
777  d = b * b - c;
778  if (d >= 0) // > 0
779  {
780  if (b > 0)
781  {
782  sd = c / (-b - std::sqrt(d));
783  }
784  else
785  {
786  sd = -b + std::sqrt(d);
787  }
788 
789  if (sd >= 0) // > 0
790  {
791  if (sd > dRmax) // Avoid rounding errors due to precision issues on
792  {
793  // 64 bits systems. Split long distance and recompute
794  double fTerm = sd - std::fmod(sd, dRmax);
795  sd = fTerm + DistanceToIn(p + fTerm * v, v);
796  }
797  zi = p.z + sd * v.z;
798 
799  if (std::fabs(zi) <= tolODz)
800  {
801  if (!fPhiFullCone)
802  {
803  xi = p.x + sd * v.x;
804  yi = p.y + sd * v.y;
805  ri = rMinAv + zi * tanRMin;
806  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
807 
808  if (cosPsi >= cosHDPhiIT)
809  {
810  if (sd > halfRadTolerance)
811  {
812  snxt = sd;
813  }
814  else
815  {
816  // Calculate a normal vector in order to check Direction
817 
818  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
819  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
820  if (norm.Dot(v) <= 0)
821  {
822  snxt = sd;
823  }
824  }
825  }
826  }
827  else
828  {
829  if (sd > halfRadTolerance)
830  {
831  return sd;
832  }
833  else
834  {
835  // Calculate a normal vector in order to check Direction
836 
837  xi = p.x + sd * v.x;
838  yi = p.y + sd * v.y;
839  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
840  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
841  if (norm.Dot(v) <= 0)
842  {
843  return sd;
844  }
845  }
846  }
847  }
848  }
849  }
850  }
851  else if (nt3 < -rin * kRadTolerance * secRMin)
852  {
853  // Within radius of inner cone (real or imaginary)
854  // -> Try 2nd root, with checking intersection is with real cone
855  // -> If check fails, try 1st root, also checking intersection is
856  // on real cone
857 
858  b = nt2 / nt1;
859  c = nt3 / nt1;
860  d = b * b - c;
861 
862  if (d >= 0) // > 0
863  {
864  if (b > 0)
865  {
866  sd = c / (-b - std::sqrt(d));
867  }
868  else
869  {
870  sd = -b + std::sqrt(d);
871  }
872  zi = p.z + sd * v.z;
873  ri = rMinAv + zi * tanRMin;
874 
875  if (ri > 0)
876  {
877  if ((sd >= 0) && (std::fabs(zi) <= tolODz)) // sd > 0
878  {
879  if (sd > dRmax) // Avoid rounding errors due to precision issues
880  {
881  // seen on 64 bits systems. Split and recompute
882  double fTerm = sd - std::fmod(sd, dRmax);
883  sd = fTerm + DistanceToIn(p + fTerm * v, v);
884  }
885  if (!fPhiFullCone)
886  {
887  xi = p.x + sd * v.x;
888  yi = p.y + sd * v.y;
889  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
890 
891  if (cosPsi >= cosHDPhiOT)
892  {
893  if (sd > halfRadTolerance)
894  {
895  snxt = sd;
896  }
897  else
898  {
899  // Calculate a normal vector in order to check Direction
900 
901  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
902  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
903  if (norm.Dot(v) <= 0)
904  {
905  snxt = sd;
906  }
907  }
908  }
909  }
910  else
911  {
912  if (sd > halfRadTolerance)
913  {
914  return sd;
915  }
916  else
917  {
918  // Calculate a normal vector in order to check Direction
919 
920  xi = p.x + sd * v.x;
921  yi = p.y + sd * v.y;
922  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
923  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
924  if (norm.Dot(v) <= 0)
925  {
926  return sd;
927  }
928  }
929  }
930  }
931  }
932  else
933  {
934  if (b > 0)
935  {
936  sd = -b - std::sqrt(d);
937  }
938  else
939  {
940  sd = c / (-b + std::sqrt(d));
941  }
942  zi = p.z + sd * v.z;
943  ri = rMinAv + zi * tanRMin;
944 
945  if ((sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz)) // sd>0
946  {
947  if (sd > dRmax) // Avoid rounding errors due to precision issues
948  {
949  // seen on 64 bits systems. Split and recompute
950  double fTerm = sd - std::fmod(sd, dRmax);
951  sd = fTerm + DistanceToIn(p + fTerm * v, v);
952  }
953  if (!fPhiFullCone)
954  {
955  xi = p.x + sd * v.x;
956  yi = p.y + sd * v.y;
957  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
958 
959  if (cosPsi >= cosHDPhiIT)
960  {
961  if (sd > halfRadTolerance)
962  {
963  snxt = sd;
964  }
965  else
966  {
967  // Calculate a normal vector in order to check Direction
968 
969  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
970  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
971  if (norm.Dot(v) <= 0)
972  {
973  snxt = sd;
974  }
975  }
976  }
977  }
978  else
979  {
980  if (sd > halfRadTolerance)
981  {
982  return sd;
983  }
984  else
985  {
986  // Calculate a normal vector in order to check Direction
987 
988  xi = p.x + sd * v.x;
989  yi = p.y + sd * v.y;
990  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
991  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
992  if (norm.Dot(v) <= 0)
993  {
994  return sd;
995  }
996  }
997  }
998  }
999  }
1000  }
1001  }
1002  else
1003  {
1004  // Within kRadTol*0.5 of inner cone (real OR imaginary)
1005  // ----> Check not travelling through (=>0 to in)
1006  // ----> if not:
1007  // -2nd root with validity check
1008 
1009  if (std::fabs(p.z) <= tolODz)
1010  {
1011  if (nt2 > 0)
1012  {
1013  // Inside inner real cone, heading outwards, inside z range
1014 
1015  if (!fPhiFullCone)
1016  {
1017  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(t3);
1018 
1019  if (cosPsi >= cosHDPhiIT)
1020  {
1021  return 0.0;
1022  }
1023  }
1024  else
1025  {
1026  return 0.0;
1027  }
1028  }
1029  else
1030  {
1031  // Within z extent, but not travelling through
1032  // -> 2nd root or UUtils::kInfinity if 1st root on imaginary cone
1033 
1034  b = nt2 / nt1;
1035  c = nt3 / nt1;
1036  d = b * b - c;
1037 
1038  if (d >= 0) // > 0
1039  {
1040  if (b > 0)
1041  {
1042  sd = -b - std::sqrt(d);
1043  }
1044  else
1045  {
1046  sd = c / (-b + std::sqrt(d));
1047  }
1048  zi = p.z + sd * v.z;
1049  ri = rMinAv + zi * tanRMin;
1050 
1051  if (ri > 0) // 2nd root
1052  {
1053  if (b > 0)
1054  {
1055  sd = c / (-b - std::sqrt(d));
1056  }
1057  else
1058  {
1059  sd = -b + std::sqrt(d);
1060  }
1061 
1062  zi = p.z + sd * v.z;
1063 
1064  if ((sd >= 0) && (std::fabs(zi) <= tolODz)) // sd>0
1065  {
1066  if (sd > dRmax) // Avoid rounding errors due to precision issue
1067  {
1068  // seen on 64 bits systems. Split and recompute
1069  double fTerm = sd - std::fmod(sd, dRmax);
1070  sd = fTerm + DistanceToIn(p + fTerm * v, v);
1071  }
1072  if (!fPhiFullCone)
1073  {
1074  xi = p.x + sd * v.x;
1075  yi = p.y + sd * v.y;
1076  ri = rMinAv + zi * tanRMin;
1077  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
1078 
1079  if (cosPsi >= cosHDPhiIT)
1080  {
1081  snxt = sd;
1082  }
1083  }
1084  else
1085  {
1086  return sd;
1087  }
1088  }
1089  }
1090  else
1091  {
1092  return UUtils::kInfinity;
1093  }
1094  }
1095  }
1096  }
1097  else // 2nd root
1098  {
1099  b = nt2 / nt1;
1100  c = nt3 / nt1;
1101  d = b * b - c;
1102 
1103  if (d > 0)
1104  {
1105  if (b > 0)
1106  {
1107  sd = c / (-b - std::sqrt(d));
1108  }
1109  else
1110  {
1111  sd = -b + std::sqrt(d);
1112  }
1113  zi = p.z + sd * v.z;
1114 
1115  if ((sd >= 0) && (std::fabs(zi) <= tolODz)) // sd>0
1116  {
1117  if (sd > dRmax) // Avoid rounding errors due to precision issues
1118  {
1119  // seen on 64 bits systems. Split and recompute
1120  double fTerm = sd - std::fmod(sd, dRmax);
1121  sd = fTerm + DistanceToIn(p + fTerm * v, v);
1122  }
1123  if (!fPhiFullCone)
1124  {
1125  xi = p.x + sd * v.x;
1126  yi = p.y + sd * v.y;
1127  ri = rMinAv + zi * tanRMin;
1128  cosPsi = (xi * cosCPhi + yi * sinCPhi) / ri;
1129 
1130  if (cosPsi >= cosHDPhiIT)
1131  {
1132  snxt = sd;
1133  }
1134  }
1135  else
1136  {
1137  return sd;
1138  }
1139  }
1140  }
1141  }
1142  }
1143  }
1144  }
1145 
1146  // Phi segment intersection
1147  //
1148  // o Tolerant of points inside phi planes by up to VUSolid::Tolerance()*0.5
1149  //
1150  // o NOTE: Large duplication of code between sphi & ephi checks
1151  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1152  // intersection check <=0 -> >=0
1153  // -> Should use some form of loop Construct
1154 
1155  if (!fPhiFullCone)
1156  {
1157  // First phi surface (starting phi)
1158 
1159  Comp = v.x * sinSPhi - v.y * cosSPhi;
1160 
1161  if (Comp < 0) // Component in outwards normal dirn
1162  {
1163  Dist = (p.y * cosSPhi - p.x * sinSPhi);
1164 
1165  if (Dist < halfCarTolerance)
1166  {
1167  sd = Dist / Comp;
1168 
1169  if (sd < snxt)
1170  {
1171  if (sd < 0)
1172  {
1173  sd = 0.0;
1174  }
1175 
1176  zi = p.z + sd * v.z;
1177 
1178  if (std::fabs(zi) <= tolODz)
1179  {
1180  xi = p.x + sd * v.x;
1181  yi = p.y + sd * v.y;
1182  rhoi2 = xi * xi + yi * yi;
1183  tolORMin2 = (rMinOAv + zi * tanRMin) * (rMinOAv + zi * tanRMin);
1184  tolORMax2 = (rMaxOAv + zi * tanRMax) * (rMaxOAv + zi * tanRMax);
1185 
1186  if ((rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2))
1187  {
1188  // z and r intersections good - check intersecting with
1189  // correct half-plane
1190 
1191  if ((yi * cosCPhi - xi * sinCPhi) <= 0)
1192  {
1193  snxt = sd;
1194  }
1195  }
1196  }
1197  }
1198  }
1199  }
1200 
1201  // Second phi surface (Ending phi)
1202 
1203  Comp = -(v.x * sinEPhi - v.y * cosEPhi);
1204 
1205  if (Comp < 0) // Component in outwards normal dirn
1206  {
1207  Dist = -(p.y * cosEPhi - p.x * sinEPhi);
1208  if (Dist < halfCarTolerance)
1209  {
1210  sd = Dist / Comp;
1211 
1212  if (sd < snxt)
1213  {
1214  if (sd < 0)
1215  {
1216  sd = 0.0;
1217  }
1218 
1219  zi = p.z + sd * v.z;
1220 
1221  if (std::fabs(zi) <= tolODz)
1222  {
1223  xi = p.x + sd * v.x;
1224  yi = p.y + sd * v.y;
1225  rhoi2 = xi * xi + yi * yi;
1226  tolORMin2 = (rMinOAv + zi * tanRMin) * (rMinOAv + zi * tanRMin);
1227  tolORMax2 = (rMaxOAv + zi * tanRMax) * (rMaxOAv + zi * tanRMax);
1228 
1229  if ((rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2))
1230  {
1231  // z and r intersections good - check intersecting with
1232  // correct half-plane
1233 
1234  if ((yi * cosCPhi - xi * sinCPhi) >= 0.0)
1235  {
1236  snxt = sd;
1237  }
1238  }
1239  }
1240  }
1241  }
1242  }
1243  }
1244  if (snxt < halfCarTolerance)
1245  {
1246  snxt = 0.;
1247  }
1248 
1249  return snxt;
1250 }
1251 
1253 //
1254 // Calculate distance (<= actual) to closest surface of shape from outside
1255 // - Calculate distance to z, radial planes
1256 // - Only to phi planes if outside phi extent
1257 // - Return 0 if point inside
1258 
1259 double UCons::SafetyFromOutside(const UVector3& p, bool aAccurate) const
1260 {
1261  double safe = 0.0, rho, safeR1, safeR2, safeZ, safePhi;
1262  double pRMin, pRMax;
1263  bool outside;
1264 
1265  rho = std::sqrt(p.x * p.x + p.y * p.y);
1266  safeZ = std::fabs(p.z) - fDz;
1267  safeR1 = 0; safeR2 = 0;
1268 
1269  if (fRmin1 || fRmin2)
1270  {
1271  pRMin = tanRMin * p.z + (fRmin1 + fRmin2) * 0.5;
1272  safeR1 = (pRMin - rho) / secRMin;
1273 
1274  pRMax = tanRMax * p.z + (fRmax1 + fRmax2) * 0.5;
1275  safeR2 = (rho - pRMax) / secRMax;
1276 
1277  if (safeR1 > safeR2)
1278  {
1279  safe = safeR1;
1280  }
1281  else
1282  {
1283  safe = safeR2;
1284  }
1285  }
1286  else
1287  {
1288  pRMax = tanRMax * p.z + (fRmax1 + fRmax2) * 0.5;
1289  safe = (rho - pRMax) / secRMax;
1290  }
1291  if (safeZ > safe)
1292  {
1293  safe = safeZ;
1294  }
1295 
1296  if (!fPhiFullCone && rho)
1297  {
1298  // Psi=angle from central phi to point
1299  //
1300  safePhi = SafetyToPhi(p,rho,outside);
1301  if ((outside) && (safePhi > safe))
1302  {
1303  safe = safePhi;
1304  }
1305  }
1306  if (safe < 0.0)
1307  {
1308  safe = 0.0; return safe; //point is Inside
1309  }
1310  if (!aAccurate) return safe;
1311 
1312  double safsq = 0.0;
1313  int count = 0;
1314  if (safeR1 > 0)
1315  {
1316  safsq += safeR1 * safeR1;
1317  count++;
1318  }
1319  if (safeR2 > 0)
1320  {
1321  safsq += safeR2 * safeR2;
1322  count++;
1323  }
1324  if (safeZ > 0)
1325  {
1326  safsq += safeZ * safeZ;
1327  count++;
1328  }
1329  if (count == 1) return safe;
1330  return std::sqrt(safsq);
1331 }
1332 
1334 //
1335 // Calculate distance to surface of shape from 'inside', allowing for tolerance
1336 // - Only Calc rmax intersection if no valid rmin intersection
1337 
1338 
1339 
1340 
1341 // double DistanceToOut( const UVector3& p, const UVector3& v, bool calcNorm, bool aConvex, UVector3 *n ) const;
1342 
1344  const UVector3& v,
1345  UVector3& aNormalVector,
1346  bool& aConvex,
1347  double /* aPstep*/) const
1348 {
1349  ESide side = kNull, sider = kNull, sidephi = kNull;
1350 
1351  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
1352  static const double halfRadTolerance = kRadTolerance * 0.5;
1353  static const double halfAngTolerance = kAngTolerance * 0.5;
1354 
1355  double snxt, srd, sphi, pdist;
1356 
1357  double rMaxAv; // Data for outer cone
1358  double rMinAv; // Data for inner cone
1359 
1360  double t1, t2, t3, rout, rin, nt1, nt2, nt3;
1361  double b, c, d, sr2, sr3;
1362 
1363  // Vars for intersection within tolerance
1364 
1365  ESide sidetol = kNull;
1366  double slentol = UUtils::kInfinity;
1367 
1368  // Vars for phi intersection:
1369 
1370  double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi;
1371  double zi, ri, deltaRoi2;
1372 
1373  // Z plane intersection
1374 
1375  if (v.z > 0.0)
1376  {
1377  pdist = fDz - p.z;
1378 
1379  if (pdist > halfCarTolerance)
1380  {
1381  snxt = pdist / v.z;
1382  side = kPZ;
1383  }
1384  else
1385  {
1386  aNormalVector = UVector3(0, 0, 1);
1387  aConvex = true;
1388  return snxt = 0.0;
1389  }
1390  }
1391  else if (v.z < 0.0)
1392  {
1393  pdist = fDz + p.z;
1394 
1395  if (pdist > halfCarTolerance)
1396  {
1397  snxt = -pdist / v.z;
1398  side = kMZ;
1399  }
1400  else
1401  {
1402  aNormalVector = UVector3(0, 0, -1);
1403  aConvex = true;
1404  return snxt = 0.0;
1405  }
1406  }
1407  else // Travel perpendicular to z axis
1408  {
1409  snxt = UUtils::kInfinity;
1410  side = kNull;
1411  }
1412 
1413  // Radial Intersections
1414  //
1415  // Intersection with outer cone (possible return) and
1416  // inner cone (must also check phi)
1417  //
1418  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1419  //
1420  // Intersects with x^2+y^2=(a*z+b)^2
1421  //
1422  // where a=tanRMax or tanRMin
1423  // b=rMaxAv or rMinAv
1424  //
1425  // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0;
1426  // t1 t2 t3
1427  //
1428  // \--------u-------/ \-----------v----------/ \---------w--------/
1429 
1430  rMaxAv = (fRmax1 + fRmax2) * 0.5;
1431 
1432  t1 = 1.0 - v.z * v.z; // since v normalised
1433  t2 = p.x * v.x + p.y * v.y;
1434  t3 = p.x * p.x + p.y * p.y;
1435  rout = tanRMax * p.z + rMaxAv;
1436 
1437  nt1 = t1 - (tanRMax * v.z) * (tanRMax * v.z);
1438  nt2 = t2 - tanRMax * v.z * rout;
1439  nt3 = t3 - rout * rout;
1440 
1441  if (v.z > 0.0)
1442  {
1443  deltaRoi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3
1444  - fRmax2 * (fRmax2 + kRadTolerance * secRMax);
1445  }
1446  else if (v.z < 0.0)
1447  {
1448  deltaRoi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3
1449  - fRmax1 * (fRmax1 + kRadTolerance * secRMax);
1450  }
1451  else
1452  {
1453  deltaRoi2 = 1.0;
1454  }
1455 
1456  if (nt1 && (deltaRoi2 > 0.0))
1457  {
1458  // Equation quadratic => 2 roots : second root must be leaving
1459 
1460  b = nt2 / nt1;
1461  c = nt3 / nt1;
1462  d = b * b - c;
1463 
1464  if (d >= 0)
1465  {
1466  // Check if on outer cone & heading outwards
1467  // NOTE: Should use rho-rout>-kRadTolerance*0.5
1468 
1469  if (nt3 > -halfRadTolerance && nt2 >= 0)
1470  {
1471  risec = std::sqrt(t3) * secRMax;
1472  aConvex = true;
1473  aNormalVector = UVector3(p.x / risec, p.y / risec, -tanRMax / secRMax);
1474  return snxt = 0;
1475  }
1476  else
1477  {
1478  sider = kRMax ;
1479  if (b > 0)
1480  {
1481  srd = -b - std::sqrt(d);
1482  }
1483  else
1484  {
1485  srd = c / (-b + std::sqrt(d));
1486  }
1487 
1488  zi = p.z + srd * v.z;
1489  ri = tanRMax * zi + rMaxAv;
1490 
1491  if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1492  {
1493  // An intersection within the tolerance
1494  // we will Store it in case it is good -
1495  //
1496  slentol = srd;
1497  sidetol = kRMax;
1498  }
1499  if ((ri < 0) || (srd < halfRadTolerance))
1500  {
1501  // Safety: if both roots -ve ensure that srd cannot `win'
1502  // distance to out
1503 
1504  if (b > 0)
1505  {
1506  sr2 = c / (-b - std::sqrt(d));
1507  }
1508  else
1509  {
1510  sr2 = -b + std::sqrt(d);
1511  }
1512  zi = p.z + sr2 * v.z;
1513  ri = tanRMax * zi + rMaxAv;
1514 
1515  if ((ri >= 0) && (sr2 > halfRadTolerance))
1516  {
1517  srd = sr2;
1518  }
1519  else
1520  {
1521  srd = UUtils::kInfinity;
1522 
1523  if ((-halfRadTolerance <= sr2) && (sr2 <= halfRadTolerance))
1524  {
1525  // An intersection within the tolerance.
1526  // Storing it in case it is good.
1527 
1528  slentol = sr2;
1529  sidetol = kRMax;
1530  }
1531  }
1532  }
1533  }
1534  }
1535  else
1536  {
1537  // No intersection with outer cone & not parallel
1538  // -> already outside, no intersection
1539 
1540  risec = std::sqrt(t3) * secRMax;
1541  aConvex = true;
1542  aNormalVector = UVector3(p.x / risec, p.y / risec, -tanRMax / secRMax);
1543  return snxt = 0.0;
1544  }
1545  }
1546  else if (nt2 && (deltaRoi2 > 0.0))
1547  {
1548  // Linear case (only one intersection) => point outside outer cone
1549 
1550  risec = std::sqrt(t3) * secRMax;
1551  aConvex = true;
1552  aNormalVector = UVector3(p.x / risec, p.y / risec, -tanRMax / secRMax);
1553  return snxt = 0.0;
1554  }
1555  else
1556  {
1557  // No intersection -> parallel to outer cone
1558  // => Z or inner cone intersection
1559 
1560  srd = UUtils::kInfinity;
1561  }
1562 
1563  // Check possible intersection within tolerance
1564 
1565  if (slentol <= halfCarTolerance)
1566  {
1567  // An intersection within the tolerance was found.
1568  // We must accept it only if the momentum points outwards.
1569  //
1570  // UVector3 ptTol; // The point of the intersection
1571  // ptTol= p + slentol*v;
1572  // ri=tanRMax*zi+rMaxAv;
1573  //
1574  // Calculate a normal vector, as below
1575 
1576  xi = p.x + slentol * v.x;
1577  yi = p.y + slentol * v.y;
1578  risec = std::sqrt(xi * xi + yi * yi) * secRMax;
1579  UVector3 norm = UVector3(xi / risec, yi / risec, -tanRMax / secRMax);
1580 
1581  if (norm.Dot(v) > 0) // We will leave the Cone immediatelly
1582  {
1583  aNormalVector = norm.Unit();
1584  aConvex = true;
1585  return snxt = 0.0;
1586  }
1587  else // On the surface, but not heading out so we ignore this intersection
1588  {
1589  // (as it is within tolerance).
1590  slentol = UUtils::kInfinity;
1591  }
1592  }
1593 
1594  // Inner Cone intersection
1595 
1596  if (fRmin1 || fRmin2)
1597  {
1598  nt1 = t1 - (tanRMin * v.z) * (tanRMin * v.z);
1599 
1600  if (nt1)
1601  {
1602  rMinAv = (fRmin1 + fRmin2) * 0.5;
1603  rin = tanRMin * p.z + rMinAv;
1604  nt2 = t2 - tanRMin * v.z * rin;
1605  nt3 = t3 - rin * rin;
1606 
1607  // Equation quadratic => 2 roots : first root must be leaving
1608 
1609  b = nt2 / nt1;
1610  c = nt3 / nt1;
1611  d = b * b - c;
1612 
1613  if (d >= 0.0)
1614  {
1615  // NOTE: should be rho-rin<kRadTolerance*0.5,
1616  // but using squared versions for efficiency
1617 
1618  if (nt3 < kRadTolerance * (rin + kRadTolerance * 0.25))
1619  {
1620  if (nt2 < 0.0)
1621  {
1622  aConvex = false;
1623  risec = std::sqrt(p.x * p.x + p.y * p.y) * secRMin;
1624  aNormalVector = UVector3(-p.x / risec, -p.y / risec, tanRMin / secRMin);
1625  return snxt = 0.0;
1626  }
1627  }
1628  else
1629  {
1630  if (b > 0)
1631  {
1632  sr2 = -b - std::sqrt(d);
1633  }
1634  else
1635  {
1636  sr2 = c / (-b + std::sqrt(d));
1637  }
1638  zi = p.z + sr2 * v.z;
1639  ri = tanRMin * zi + rMinAv;
1640 
1641  if ((ri >= 0.0) && (-halfRadTolerance <= sr2) && (sr2 <= halfRadTolerance))
1642  {
1643  // An intersection within the tolerance
1644  // storing it in case it is good.
1645 
1646  slentol = sr2;
1647  sidetol = kRMax;
1648  }
1649  if ((ri < 0) || (sr2 < halfRadTolerance))
1650  {
1651  if (b > 0)
1652  {
1653  sr3 = c / (-b - std::sqrt(d));
1654  }
1655  else
1656  {
1657  sr3 = -b + std::sqrt(d);
1658  }
1659 
1660  // Safety: if both roots -ve ensure that srd cannot `win'
1661  // distancetoout
1662 
1663  if (sr3 > halfRadTolerance)
1664  {
1665  if (sr3 < srd)
1666  {
1667  zi = p.z + sr3 * v.z;
1668  ri = tanRMin * zi + rMinAv;
1669 
1670  if (ri >= 0.0)
1671  {
1672  srd = sr3;
1673  sider = kRMin;
1674  }
1675  }
1676  }
1677  else if (sr3 > -halfRadTolerance)
1678  {
1679  // Intersection in tolerance. Store to check if it's good
1680 
1681  slentol = sr3;
1682  sidetol = kRMin;
1683  }
1684  }
1685  else if ((sr2 < srd) && (sr2 > halfCarTolerance))
1686  {
1687  srd = sr2;
1688  sider = kRMin;
1689  }
1690  else if (sr2 > -halfCarTolerance)
1691  {
1692  // Intersection in tolerance. Store to check if it's good
1693 
1694  slentol = sr2;
1695  sidetol = kRMin;
1696  }
1697  if (slentol <= halfCarTolerance)
1698  {
1699  // An intersection within the tolerance was found.
1700  // We must accept it only if the momentum points outwards.
1701 
1702  UVector3 norm;
1703 
1704  // Calculate a normal vector, as below
1705 
1706  xi = p.x + slentol * v.x;
1707  yi = p.y + slentol * v.y;
1708  if (sidetol == kRMax)
1709  {
1710  risec = std::sqrt(xi * xi + yi * yi) * secRMax;
1711  norm = UVector3(xi / risec, yi / risec, -tanRMax / secRMax);
1712  }
1713  else
1714  {
1715  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
1716  norm = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
1717  }
1718  if (norm.Dot(v) > 0)
1719  {
1720  // We will leave the cone immediately
1721 
1722  aNormalVector = norm.Unit();
1723  aConvex = true;
1724  return snxt = 0.0;
1725  }
1726  else
1727  {
1728  // On the surface, but not heading out so we ignore this
1729  // intersection (as it is within tolerance).
1730 
1731  slentol = UUtils::kInfinity;
1732  }
1733  }
1734  }
1735  }
1736  }
1737  }
1738 
1739  // Linear case => point outside inner cone ---> outer cone intersect
1740  //
1741  // Phi Intersection
1742 
1743  if (!fPhiFullCone)
1744  {
1745  // add angle calculation with correction
1746  // of the difference in domain of atan2 and Sphi
1747 
1748  vphi = std::atan2(v.y, v.x);
1749 
1750  if (vphi < fSPhi - halfAngTolerance)
1751  {
1752  vphi += 2 * UUtils::kPi;
1753  }
1754  else if (vphi > fSPhi + fDPhi + halfAngTolerance)
1755  {
1756  vphi -= 2 * UUtils::kPi;
1757  }
1758 
1759  if (p.x || p.y) // Check if on z axis (rho not needed later)
1760  {
1761  // pDist -ve when inside
1762 
1763  pDistS = p.x * sinSPhi - p.y * cosSPhi;
1764  pDistE = -p.x * sinEPhi + p.y * cosEPhi;
1765 
1766  // Comp -ve when in direction of outwards normal
1767 
1768  compS = -sinSPhi * v.x + cosSPhi * v.y;
1769  compE = sinEPhi * v.x - cosEPhi * v.y;
1770 
1771  sidephi = kNull;
1772 
1773  if (((fDPhi <= UUtils::kPi) && ((pDistS <= halfCarTolerance)
1774  && (pDistE <= halfCarTolerance)))
1775  || ((fDPhi > UUtils::kPi) && !((pDistS > halfCarTolerance)
1776  && (pDistE > halfCarTolerance))))
1777  {
1778  // Inside both phi *full* planes
1779  if (compS < 0)
1780  {
1781  sphi = pDistS / compS;
1782  if (sphi >= -halfCarTolerance)
1783  {
1784  xi = p.x + sphi * v.x;
1785  yi = p.y + sphi * v.y;
1786 
1787  // Check intersecting with correct half-plane
1788  // (if not -> no intersect)
1789  //
1790  if ((std::fabs(xi) <= VUSolid::Tolerance())
1791  && (std::fabs(yi) <= VUSolid::Tolerance()))
1792  {
1793  sidephi = kSPhi;
1794  if ((fSPhi - halfAngTolerance <= vphi)
1795  && (fSPhi + fDPhi + halfAngTolerance >= vphi))
1796  {
1797  sphi = UUtils::kInfinity;
1798  }
1799  }
1800  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
1801  {
1802  sphi = UUtils::kInfinity;
1803  }
1804  else
1805  {
1806  sidephi = kSPhi;
1807  if (pDistS > -halfCarTolerance)
1808  {
1809  sphi = 0.0; // Leave by sphi immediately
1810  }
1811  }
1812  }
1813  else
1814  {
1815  sphi = UUtils::kInfinity;
1816  }
1817  }
1818  else
1819  {
1820  sphi = UUtils::kInfinity;
1821  }
1822 
1823  if (compE < 0)
1824  {
1825  sphi2 = pDistE / compE;
1826 
1827  // Only check further if < starting phi intersection
1828  //
1829  if ((sphi2 > -halfCarTolerance) && (sphi2 < sphi))
1830  {
1831  xi = p.x + sphi2 * v.x;
1832  yi = p.y + sphi2 * v.y;
1833 
1834  // Check intersecting with correct half-plane
1835 
1836  if ((std::fabs(xi) <= VUSolid::Tolerance())
1837  && (std::fabs(yi) <= VUSolid::Tolerance()))
1838  {
1839  // Leaving via ending phi
1840 
1841  if (!((fSPhi - halfAngTolerance <= vphi)
1842  && (fSPhi + fDPhi + halfAngTolerance >= vphi)))
1843  {
1844  sidephi = kEPhi;
1845  if (pDistE <= -halfCarTolerance)
1846  {
1847  sphi = sphi2;
1848  }
1849  else
1850  {
1851  sphi = 0.0;
1852  }
1853  }
1854  }
1855  else // Check intersecting with correct half-plane
1856  if (yi * cosCPhi - xi * sinCPhi >= 0)
1857  {
1858  // Leaving via ending phi
1859 
1860  sidephi = kEPhi;
1861  if (pDistE <= -halfCarTolerance)
1862  {
1863  sphi = sphi2;
1864  }
1865  else
1866  {
1867  sphi = 0.0;
1868  }
1869  }
1870  }
1871  }
1872  }
1873  else
1874  {
1875  sphi = UUtils::kInfinity;
1876  }
1877  }
1878  else
1879  {
1880  // On z axis + travel not || to z axis -> if phi of vector direction
1881  // within phi of shape, Step limited by rmax, else Step =0
1882 
1883  if ((fSPhi - halfAngTolerance <= vphi)
1884  && (vphi <= fSPhi + fDPhi + halfAngTolerance))
1885  {
1886  sphi = UUtils::kInfinity;
1887  }
1888  else
1889  {
1890  sidephi = kSPhi ; // arbitrary
1891  sphi = 0.0;
1892  }
1893  }
1894  if (sphi < snxt) // Order intersecttions
1895  {
1896  snxt = sphi;
1897  side = sidephi;
1898  }
1899  }
1900  if (srd < snxt) // Order intersections
1901  {
1902  snxt = srd ;
1903  side = sider;
1904  }
1905 
1906  switch (side)
1907  {
1908  // Note: returned vector not normalised
1909  case kRMax: // (divide by frmax for Unit vector)
1910  xi = p.x + snxt * v.x;
1911  yi = p.y + snxt * v.y;
1912  risec = std::sqrt(xi * xi + yi * yi) * secRMax;
1913  aNormalVector = UVector3(xi / risec, yi / risec, -tanRMax / secRMax);
1914  aConvex = true;
1915  break;
1916  case kRMin:
1917  xi = p.x + snxt * v.x;
1918  yi = p.y + snxt * v.y;
1919  risec = std::sqrt(xi * xi + yi * yi) * secRMin;
1920  aNormalVector = UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
1921  aConvex = false; // Rmin is inconvex
1922  break;
1923  case kSPhi:
1924  if (fDPhi <= UUtils::kPi)
1925  {
1926  aNormalVector = UVector3(sinSPhi, -cosSPhi, 0);
1927  aConvex = true;
1928  }
1929  else
1930  {
1931  aNormalVector = UVector3(sinSPhi, -cosSPhi, 0);
1932  aConvex = false;
1933  }
1934  break;
1935  case kEPhi:
1936  if (fDPhi <= UUtils::kPi)
1937  {
1938  aNormalVector = UVector3(-sinEPhi, cosEPhi, 0);
1939  aConvex = true;
1940  }
1941  else
1942  {
1943  aNormalVector = UVector3(-sinEPhi, cosEPhi, 0);
1944  aConvex = false;
1945  }
1946  break;
1947  case kPZ:
1948  aNormalVector = UVector3(0, 0, 1);
1949  aConvex = true;
1950  break;
1951  case kMZ:
1952  aNormalVector = UVector3(0, 0, -1);
1953  aConvex = true;
1954  break;
1955  default:
1956  cout << std::endl;
1957  // DumpInfo();
1958  std::ostringstream message;
1959  int oldprc = message.precision(16);
1960  message << "Undefined side for valid surface normal to solid."
1961  << std::endl
1962  << "Position:" << std::endl << std::endl
1963  << "p.x = " << p.x << " mm" << std::endl
1964  << "p.y = " << p.y << " mm" << std::endl
1965  << "p.z = " << p.z << " mm" << std::endl << std::endl
1966  << "pho at z = " << std::sqrt(p.x * p.x + p.y * p.y)
1967  << " mm" << std::endl << std::endl;
1968  if (p.x != 0. || p.y != 0.)
1969  {
1970  message << "point phi = " << std::atan2(p.y, p.x) / (UUtils::kPi / 180.0)
1971  << " degree" << std::endl << std::endl;
1972  }
1973  message << "Direction:" << std::endl << std::endl
1974  << "v.x = " << v.x << std::endl
1975  << "v.y = " << v.y << std::endl
1976  << "v.z = " << v.z << std::endl << std::endl
1977  << "Proposed distance :" << std::endl << std::endl
1978  << "snxt = " << snxt << " mm" << std::endl;
1979  message.precision(oldprc);
1980  UUtils::Exception("UCons::DistanceToOut()", "UGeomSolids", Warning, 1, message.str().c_str());
1981  break;
1982  }
1983  if (snxt < halfCarTolerance)
1984  {
1985  snxt = 0.;
1986  }
1987 
1988  return snxt;
1989 }
1990 
1992 //
1993 // Calculate distance (<=actual) to closest surface of shape from inside
1994 
1995 double UCons::SafetyFromInside(const UVector3& p, bool) const
1996 {
1997  double safe = 0.0, rho, safeZ;
1998 
1999 #ifdef UCSGDEBUG
2000  if (Inside(p) == eOutside)
2001  {
2002  int oldprc = cout.precision(16);
2003  cout << std::endl;
2004  DumpInfo();
2005  cout << "Position:" << std::endl << std::endl;
2006  cout << "p.x = " << p.x << " mm" << std::endl;
2007  cout << "p.y = " << p.y << " mm" << std::endl;
2008  cout << "p.z = " << p.z << " mm" << std::endl << std::endl;
2009  cout << "pho at z = " << std::sqrt(p.x * p.x + p.y * p.y)
2010  << " mm" << std::endl << std::endl;
2011  if ((p.x != 0.) || (p.x != 0.))
2012  {
2013  cout << "point phi = " << std::atan2(p.y, p.x) / degree
2014  << " degree" << std::endl << std::endl;
2015  }
2016  cout.precision(oldprc);
2017  UUtils::Exception("UCons::UCons()", "UGeomSolids", Warning, 1, message.str().c_str());
2018 
2019  }
2020 #endif
2021 
2022  rho = std::sqrt(p.x * p.x + p.y * p.y);
2023 
2024  safe = SafetyFromInsideR(p, rho, false);
2025  safeZ = fDz - std::fabs(p.z);
2026 
2027  if (safeZ < safe)
2028  {
2029  safe = safeZ;
2030  }
2031  if (safe < 0)
2032  {
2033  safe = 0;
2034  }
2035 
2036  return safe;
2037 }
2038 
2039 
2041 //
2042 // GetEntityType
2043 
2045 {
2046  return std::string("Cons");
2047 }
2048 
2050 //
2051 // Make a clone of the object
2052 //
2054 {
2055 
2056  return new UCons(*this);
2057 }
2058 
2060 //
2061 // Stream object contents to an output stream
2062 
2063 std::ostream& UCons::StreamInfo(std::ostream& os) const
2064 {
2065  int oldprc = os.precision(16);
2066  os << "-----------------------------------------------------------\n"
2067  << " *** Dump for solid - " << GetName() << " ***\n"
2068  << " ===================================================\n"
2069  << " Solid type: UCons\n"
2070  << " Parameters: \n"
2071  << " inside -fDz radius : " << fRmin1 << " mm \n"
2072  << " outside -fDz radius: " << fRmax1 << " mm \n"
2073  << " inside +fDz radius : " << fRmin2 << " mm \n"
2074  << " outside +fDz radius: " << fRmax2 << " mm \n"
2075  << " half length in Z : " << fDz << " mm \n"
2076  << " starting angle of segment: " << fSPhi / (UUtils::kPi / 180.0) << " degrees \n"
2077  << " delta angle of segment : " << fDPhi / (UUtils::kPi / 180.0) << " degrees \n"
2078  << "-----------------------------------------------------------\n";
2079  os.precision(oldprc);
2080 
2081  return os;
2082 }
2083 
2084 
2085 
2087 //
2088 // GetPointOnSurface
2089 
2091 {
2092  // declare working variables
2093  //
2094  double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2095  double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2096  rone = (fRmax1 - fRmax2) / (2.*fDz);
2097  rtwo = (fRmin1 - fRmin2) / (2.*fDz);
2098  qone = 0.;
2099  qtwo = 0.;
2100  if (fRmax1 != fRmax2)
2101  {
2102  qone = fDz * (fRmax1 + fRmax2) / (fRmax1 - fRmax2);
2103  }
2104  if (fRmin1 != fRmin2)
2105  {
2106  qtwo = fDz * (fRmin1 + fRmin2) / (fRmin1 - fRmin2);
2107  }
2108  slin = std::sqrt(UUtils::sqr(fRmin1 - fRmin2) + UUtils::sqr(2.*fDz));
2109  slout = std::sqrt(UUtils::sqr(fRmax1 - fRmax2) + UUtils::sqr(2.*fDz));
2110  Aone = 0.5 * fDPhi * (fRmax2 + fRmax1) * slout;
2111  Atwo = 0.5 * fDPhi * (fRmin2 + fRmin1) * slin;
2112  Athree = 0.5 * fDPhi * (fRmax1 * fRmax1 - fRmin1 * fRmin1);
2113  Afour = 0.5 * fDPhi * (fRmax2 * fRmax2 - fRmin2 * fRmin2);
2114  Afive = fDz * (fRmax1 - fRmin1 + fRmax2 - fRmin2);
2115 
2116  phi = UUtils::Random(fSPhi, fSPhi + fDPhi);
2117  cosu = std::cos(phi);
2118  sinu = std::sin(phi);
2121 
2122  if ((fSPhi == 0.) && fPhiFullCone)
2123  {
2124  Afive = 0.;
2125  }
2126  chose = UUtils::Random(0., Aone + Atwo + Athree + Afour + 2.*Afive);
2127 
2128  if ((chose >= 0.) && (chose < Aone))
2129  {
2130  if (fRmin1 != fRmin2)
2131  {
2132  zRand = UUtils::Random(-1.*fDz, fDz);
2133  return UVector3(rtwo * cosu * (qtwo - zRand),
2134  rtwo * sinu * (qtwo - zRand), zRand);
2135  }
2136  else
2137  {
2138  return UVector3(fRmin1 * cosu, fRmin2 * sinu,
2139  UUtils::Random(-1.*fDz, fDz));
2140  }
2141  }
2142  else if ((chose >= Aone) && (chose <= Aone + Atwo))
2143  {
2144  if (fRmax1 != fRmax2)
2145  {
2146  zRand = UUtils::Random(-1.*fDz, fDz);
2147  return UVector3(rone * cosu * (qone - zRand),
2148  rone * sinu * (qone - zRand), zRand);
2149  }
2150  else
2151  {
2152  return UVector3(fRmax1 * cosu, fRmax2 * sinu,
2153  UUtils::Random(-1.*fDz, fDz));
2154  }
2155  }
2156  else if ((chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree))
2157  {
2158  return UVector3(rRand1 * cosu, rRand1 * sinu, -1 * fDz);
2159  }
2160  else if ((chose >= Aone + Atwo + Athree)
2161  && (chose < Aone + Atwo + Athree + Afour))
2162  {
2163  return UVector3(rRand2 * cosu, rRand2 * sinu, fDz);
2164  }
2165  else if ((chose >= Aone + Atwo + Athree + Afour)
2166  && (chose < Aone + Atwo + Athree + Afour + Afive))
2167  {
2168  zRand = UUtils::Random(-1.*fDz, fDz);
2169  rRand1 = UUtils::Random(fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2),
2170  fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2));
2171  return UVector3(rRand1 * std::cos(fSPhi),
2172  rRand1 * std::sin(fSPhi), zRand);
2173  }
2174  else
2175  {
2176  zRand = UUtils::Random(-1.*fDz, fDz);
2177  rRand1 = UUtils::Random(fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2),
2178  fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2));
2179  return UVector3(rRand1 * std::cos(fSPhi + fDPhi),
2180  rRand1 * std::sin(fSPhi + fDPhi), zRand);
2181  }
2182 }
2183 
2184 
2185 void UCons::Extent(UVector3& aMin, UVector3& aMax) const
2186 {
2187  double max = fRmax1 > fRmax2 ? fRmax1 : fRmax2;
2188  aMin = UVector3(-max, -max, -fDz);
2189  aMax = UVector3(max, max, fDz);
2190 }
2191 
2192 void UCons::GetParametersList(int, double* aArray) const
2193 {
2194  aArray[0] = GetInnerRadiusMinusZ();
2195  aArray[1] = GetOuterRadiusMinusZ();
2196  aArray[2] = GetInnerRadiusPlusZ();
2197  aArray[3] = GetOuterRadiusPlusZ();
2198  aArray[4] = GetZHalfLength();
2199  aArray[5] = GetStartPhiAngle();
2200  aArray[6] = GetDeltaPhiAngle();
2201 
2202 }
double SafetyToPhi(const UVector3 &p, const double rho, bool &outside) const
double sinSPhi
Definition: UCons.hh:201
double SafetyFromOutside(const UVector3 &p, bool precise) const
Definition: UCons.cc:1259
bool fPhiFullCone
Definition: UCons.hh:206
static double frTolerance
Definition: VUSolid.hh:31
Definition: UCons.cc:29
VUSolid::EnumInside Inside(const UVector3 &p) const
void Initialize()
UGeometryType GetEntityType() const
Definition: UCons.cc:2044
Definition: UCons.cc:33
double fSPhi
Definition: UCons.hh:197
std::ostream & StreamInfo(std::ostream &os) const
Definition: UCons.cc:2063
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UCons.cc:440
UCons & operator=(const UCons &rhs)
Definition: UCons.cc:130
const std::string & GetName() const
Definition: VUSolid.hh:103
double cosHDPhiOT
Definition: UCons.hh:201
Definition: UCons.cc:29
~UCons()
Definition: UCons.cc:106
double secRMax
Definition: UCons.hh:210
double cosCPhi
Definition: UCons.hh:201
double SafetyFromInsideR(const UVector3 &p, const double rho, bool) const
double SafetyFromInside(const UVector3 &p, bool precise) const
Definition: UCons.cc:1995
double GetOuterRadiusPlusZ() const
VUSolid * Clone() const
Definition: UCons.cc:2053
static double Tolerance()
Definition: VUSolid.hh:127
double x
Definition: UVector3.hh:136
ENorm
Definition: G4Cons.cc:75
Definition: UCons.hh:49
double tanRMin
Definition: UCons.hh:210
ESide
Definition: UCons.hh:187
double tanRMax
Definition: UCons.hh:210
static const double kInfinity
Definition: UUtils.hh:53
double kAngTolerance
Definition: UCons.hh:193
double GetStartPhiAngle() const
bool Normal(const UVector3 &p, UVector3 &n) const
Definition: UCons.cc:174
double GetZHalfLength() const
Definition: UCons.cc:33
double sinCPhi
Definition: UCons.hh:201
static double faTolerance
Definition: VUSolid.hh:32
double fDPhi
Definition: UCons.hh:197
Definition: UCons.cc:29
ENorm
Definition: UCons.hh:191
Definition: UCons.cc:29
double fRmin2
Definition: UCons.hh:197
const G4int n
double sinEPhi
Definition: UCons.hh:201
double Dot(const UVector3 &) const
Definition: UVector3.hh:257
Definition: UCons.cc:33
Definition: UCons.cc:29
double cosSPhi
Definition: UCons.hh:201
Definition: UCons.cc:29
double GetInnerRadiusPlusZ() const
T sqr(const T &x)
Definition: UUtils.hh:137
double GetRadiusInRing(double rmin, double rmax)
Definition: UUtils.hh:155
double fRmax1
Definition: UCons.hh:197
double GetOuterRadiusMinusZ() const
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UCons.cc:2185
double kRadTolerance
Definition: UCons.hh:193
virtual void GetParametersList(int, double *) const
Definition: UCons.cc:2192
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double fDz
Definition: UCons.hh:197
static const double kPi
Definition: UUtils.hh:48
UVector3 Unit() const
Definition: UVector3.cc:80
void CheckPhiAngles(double sPhi, double dPhi)
double fRmin1
Definition: UCons.hh:197
Definition: UCons.cc:29
double GetDeltaPhiAngle() const
double secRMin
Definition: UCons.hh:210
static const double degree
Definition: G4SIunits.hh:125
std::string UGeometryType
Definition: UTypes.hh:39
double GetInnerRadiusMinusZ() const
Definition: UCons.cc:33
double cosHDPhiIT
Definition: UCons.hh:201
UCons()
Definition: UCons.cc:92
UVector3 GetPointOnSurface() const
Definition: UCons.cc:2090
ESide
Definition: G4Cons.cc:71
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 e3
UVector3 ApproxSurfaceNormal(const UVector3 &p) const
Definition: UCons.cc:298
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double y
Definition: UVector3.hh:137
double fRmax2
Definition: UCons.hh:197
double cosEPhi
Definition: UCons.hh:201
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UCons.cc:1343
Definition: UCons.cc:33