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