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