Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
UCons Class Reference

#include <UCons.hh>

Inheritance diagram for UCons:
VUSolid

Public Member Functions

 UCons (const std::string &pName, double pRmin1, double pRmax1, double pRmin2, double pRmax2, double pDz, double pSPhi, double pDPhi)
 
 ~UCons ()
 
double GetInnerRadiusMinusZ () const
 
double GetOuterRadiusMinusZ () const
 
double GetInnerRadiusPlusZ () const
 
double GetOuterRadiusPlusZ () const
 
double GetZHalfLength () const
 
double GetStartPhiAngle () const
 
double GetDeltaPhiAngle () const
 
void SetInnerRadiusMinusZ (double Rmin1)
 
void SetOuterRadiusMinusZ (double Rmax1)
 
void SetInnerRadiusPlusZ (double Rmin2)
 
void SetOuterRadiusPlusZ (double Rmax2)
 
void SetZHalfLength (double newDz)
 
void SetStartPhiAngle (double newSPhi, bool trig=true)
 
void SetDeltaPhiAngle (double newDPhi)
 
double GetCubicVolume ()
 
double GetSurfaceArea ()
 
bool Normal (const UVector3 &p, UVector3 &n) const
 
double DistanceToIn (const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
 
double SafetyFromOutside (const UVector3 &p, bool precise) const
 
double DistanceToOut (const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
 
double SafetyFromInside (const UVector3 &p, bool precise) const
 
UGeometryType GetEntityType () const
 
UVector3 GetPointOnSurface () const
 
VUSolid * Clone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
void Extent (UVector3 &aMin, UVector3 &aMax) const
 
virtual void GetParametersList (int, double *) const
 
virtual void ComputeBBox (UBBox *, bool)
 
VUSolid::EnumInside Inside (const UVector3 &p) const
 
 UCons ()
 
 UCons (const UCons &rhs)
 
UCons & operator= (const UCons &rhs)
 
double GetRmin1 () const
 
double GetRmax1 () const
 
double GetRmin2 () const
 
double GetRmax2 () const
 
double GetDz () const
 
double GetSPhi () const
 
double GetDPhi () const
 
- Public Member Functions inherited from VUSolid
 VUSolid ()
 
 VUSolid (const std::string &name)
 
virtual ~VUSolid ()
 
double GetCarTolerance () const
 
double GetRadTolerance () const
 
double GetAngTolerance () const
 
void SetCarTolerance (double eps)
 
void SetRadTolerance (double eps)
 
void SetAngTolerance (double eps)
 
virtual void ExtentAxis (EAxisType aAxis, double &aMin, double &aMax) const
 
const std::string & GetName () const
 
void SetName (const std::string &aName)
 
virtual void SamplePointsInside (int, UVector3 *) const
 
virtual void SamplePointsOnSurface (int, UVector3 *) const
 
virtual void SamplePointsOnEdge (int, UVector3 *) const
 
double EstimateCubicVolume (int nStat, double epsilon) const
 
double EstimateSurfaceArea (int nStat, double ell) const
 

Additional Inherited Members

- Public Types inherited from VUSolid
enum  EnumInside { eInside =0, eSurface =1, eOutside =2 }
 
enum  EAxisType { eXaxis =0, eYaxis =1, eZaxis =2 }
 
- Static Public Member Functions inherited from VUSolid
static double Tolerance ()
 
- Static Protected Attributes inherited from VUSolid
static double fgTolerance = 1.0E-9
 
static double frTolerance = 1.0E-9
 
static double faTolerance = 1.0E-9
 

Detailed Description

Definition at line 49 of file UCons.hh.

Constructor & Destructor Documentation

UCons::UCons ( const std::string &  pName,
double  pRmin1,
double  pRmax1,
double  pRmin2,
double  pRmax2,
double  pDz,
double  pSPhi,
double  pDPhi 
)

Definition at line 40 of file UCons.cc.

References UUtils::Exception(), FatalErrorInArguments, VUSolid::faTolerance, VUSolid::frTolerance, and VUSolid::GetName().

45  : VUSolid(pName.c_str()), fRmin1(pRmin1), fRmin2(pRmin2),
46  fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
47 {
48  kRadTolerance = frTolerance;
49  kAngTolerance = faTolerance;
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 }
VUSolid::VUSolid
VUSolid()
Definition: VUSolid.cc:18
VUSolid::frTolerance
static double frTolerance
Definition: VUSolid.hh:31
FatalErrorInArguments
Definition: UUtils.hh:36
VUSolid::GetName
const std::string & GetName() const
Definition: VUSolid.hh:103
VUSolid::faTolerance
static double faTolerance
Definition: VUSolid.hh:32
UUtils::Exception
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:177
UCons::~UCons ( )

Definition at line 107 of file UCons.cc.

108 {
109 }
UCons::UCons ( )

Definition at line 93 of file UCons.cc.

Referenced by Clone().

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 }
VUSolid::VUSolid
VUSolid()
Definition: VUSolid.cc:18
UCons::UCons ( const UCons &  rhs)

Definition at line 115 of file UCons.cc.

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 }
VUSolid::VUSolid
VUSolid()
Definition: VUSolid.cc:18

Member Function Documentation

VUSolid * UCons::Clone ( ) const
virtual

Implements VUSolid.

Definition at line 2098 of file UCons.cc.

References UCons().

2099 {
2100 
2101  return new UCons(*this);
2102 }
UCons::UCons
UCons()
Definition: UCons.cc:93
virtual void UCons::ComputeBBox ( UBBox *  ,
bool   
)
inlinevirtual

Implements VUSolid.

Definition at line 122 of file UCons.hh.

122 {}
double UCons::DistanceToIn ( const UVector3 &  p,
const UVector3 &  v,
double  aPstep = UUtils::kInfinity 
) const
virtual

Implements VUSolid.

Definition at line 441 of file UCons.cc.

References test::b, test::c, UVector3::Dot(), G4INCL::Math::max(), plottest35::t1, VUSolid::Tolerance(), UVector3::x, UVector3::y, and UVector3::z.

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 }
UCons::DistanceToIn
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UCons.cc:441
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
UVector3::Dot
double Dot(const UVector3 &) const
Definition: UVector3.hh:257
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
plottest35.t1
tuple t1
Definition: plottest35.py:33
UVector3::z
double z
Definition: UVector3.hh:138
UVector3::y
double y
Definition: UVector3.hh:137
double UCons::DistanceToOut ( const UVector3 &  aPoint,
const UVector3 &  aDirection,
UVector3 &  aNormalVector,
bool &  aConvex,
double  aPstep = UUtils::kInfinity 
) const
virtual

Implements VUSolid.

Definition at line 1355 of file UCons.cc.

References test::b, test::c, UVector3::Dot(), UUtils::Exception(), plottest35::t1, VUSolid::Tolerance(), UVector3::Unit(), Warning, UVector3::x, UVector3::y, and UVector3::z.

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 }
p
const char * p
Definition: xmltok.h:285
Warning
Definition: UUtils.hh:36
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
UVector3::Dot
double Dot(const UVector3 &) const
Definition: UVector3.hh:257
UVector3::Unit
UVector3 Unit() const
Definition: UVector3.cc:80
plottest35.t1
tuple t1
Definition: plottest35.py:33
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
void UCons::Extent ( UVector3 &  aMin,
UVector3 &  aMax 
) const
virtual

Implements VUSolid.

Definition at line 2230 of file UCons.cc.

References G4INCL::Math::max().

2231 {
2232  double max = fRmax1 > fRmax2 ? fRmax1 : fRmax2;
2233  aMin = UVector3(-max, -max, -fDz);
2234  aMax = UVector3(max, max, fDz);
2235 }
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
double UCons::GetCubicVolume ( )
inline
double UCons::GetDeltaPhiAngle ( ) const
inline

Referenced by G4UCons::GetDeltaPhiAngle(), and GetParametersList().

double UCons::GetDPhi ( ) const
inline
double UCons::GetDz ( ) const
inline

Referenced by UPolycone::InsideSection().

UGeometryType UCons::GetEntityType ( ) const
virtual

Implements VUSolid.

Definition at line 2089 of file UCons.cc.

2090 {
2091  return std::string("Cons");
2092 }
double UCons::GetInnerRadiusMinusZ ( ) const
inline

Referenced by G4UCons::GetInnerRadiusMinusZ(), and GetParametersList().

double UCons::GetInnerRadiusPlusZ ( ) const
inline

Referenced by G4UCons::GetInnerRadiusPlusZ(), and GetParametersList().

double UCons::GetOuterRadiusMinusZ ( ) const
inline

Referenced by G4UCons::GetOuterRadiusMinusZ(), and GetParametersList().

double UCons::GetOuterRadiusPlusZ ( ) const
inline

Referenced by G4UCons::GetOuterRadiusPlusZ(), and GetParametersList().

void UCons::GetParametersList ( int  ,
double *  aArray 
) const
virtual

Implements VUSolid.

Definition at line 2237 of file UCons.cc.

References GetDeltaPhiAngle(), GetInnerRadiusMinusZ(), GetInnerRadiusPlusZ(), GetOuterRadiusMinusZ(), GetOuterRadiusPlusZ(), GetStartPhiAngle(), and GetZHalfLength().

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 }
UCons::GetOuterRadiusPlusZ
double GetOuterRadiusPlusZ() const
UCons::GetStartPhiAngle
double GetStartPhiAngle() const
UCons::GetZHalfLength
double GetZHalfLength() const
UCons::GetInnerRadiusPlusZ
double GetInnerRadiusPlusZ() const
UCons::GetOuterRadiusMinusZ
double GetOuterRadiusMinusZ() const
UCons::GetDeltaPhiAngle
double GetDeltaPhiAngle() const
UCons::GetInnerRadiusMinusZ
double GetInnerRadiusMinusZ() const
UVector3 UCons::GetPointOnSurface ( ) const
virtual

Implements VUSolid.

Definition at line 2135 of file UCons.cc.

References UUtils::GetRadiusInRing(), UUtils::Random(), and UUtils::sqr().

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);
2164  rRand1 = UUtils::GetRadiusInRing(fRmin1, fRmin2);
2165  rRand2 = UUtils::GetRadiusInRing(fRmax1, fRmax2);
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 }
UUtils::sqr
T sqr(const T &x)
Definition: UUtils.hh:142
UUtils::GetRadiusInRing
double GetRadiusInRing(double rmin, double rmax)
Definition: UUtils.hh:160
UUtils::Random
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double UCons::GetRmax1 ( ) const
inline

Referenced by UPolycone::InsideSection().

double UCons::GetRmax2 ( ) const
inline

Referenced by UPolycone::InsideSection().

double UCons::GetRmin1 ( ) const
inline

Referenced by UPolycone::InsideSection().

double UCons::GetRmin2 ( ) const
inline

Referenced by UPolycone::InsideSection().

double UCons::GetSPhi ( ) const
inline
double UCons::GetStartPhiAngle ( ) const
inline

Referenced by GetParametersList(), and G4UCons::GetStartPhiAngle().

double UCons::GetSurfaceArea ( )
inline
double UCons::GetZHalfLength ( ) const
inline

Referenced by GetParametersList(), and G4UCons::GetZHalfLength().

VUSolid::EnumInside UCons::Inside ( const UVector3 &  p) const
inlinevirtual

Implements VUSolid.

Definition at line 127 of file UCons.hh.

References VUSolid::eInside, VUSolid::eOutside, VUSolid::eSurface, VUSolid::Tolerance(), UVector3::x, UVector3::y, and UVector3::z.

Referenced by SafetyFromInside().

128  {
129  double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2;
130  VUSolid::EnumInside in;
131  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
132  static const double halfRadTolerance = kRadTolerance * 0.5;
133  static const double halfAngTolerance = kAngTolerance * 0.5;
134 
135  if (std::fabs(p.z) > fDz + halfCarTolerance)
136  {
137  return in = eOutside;
138  }
139  else if (std::fabs(p.z) >= fDz - halfCarTolerance)
140  {
141  in = eSurface;
142  }
143  else
144  {
145  in = eInside;
146  }
147 
148  r2 = p.x * p.x + p.y * p.y;
149  rl = 0.5 * (fRmin2 * (p.z + fDz) + fRmin1 * (fDz - p.z)) / fDz;
150  rh = 0.5 * (fRmax2 * (p.z + fDz) + fRmax1 * (fDz - p.z)) / fDz;
151 
152  // rh2 = rh*rh;
153 
154  tolRMin = rl - halfRadTolerance;
155  if (tolRMin < 0)
156  {
157  tolRMin = 0;
158  }
159  tolRMax = rh + halfRadTolerance;
160 
161  if ((r2 < tolRMin * tolRMin) || (r2 > tolRMax * tolRMax))
162  {
163  return in = eOutside;
164  }
165 
166  if (rl)
167  {
168  tolRMin = rl + halfRadTolerance;
169  }
170  else
171  {
172  tolRMin = 0.0;
173  }
174  tolRMax = rh - halfRadTolerance;
175 
176  if (in == eInside) // else it's eSurface already
177  {
178  if ((r2 < tolRMin * tolRMin) || (r2 >= tolRMax * tolRMax))
179  {
180  in = eSurface;
181  }
182  }
183  if (!fPhiFullCone && ((p.x != 0.0) || (p.y != 0.0)))
184  {
185  pPhi = std::atan2(p.y, p.x);
186 
187  if (pPhi < fSPhi - halfAngTolerance)
188  {
189  pPhi += 2 * UUtils::kPi;
190  }
191  else if (pPhi > fSPhi + fDPhi + halfAngTolerance)
192  {
193  pPhi -= 2 * UUtils::kPi;
194  }
195 
196  if ((pPhi < fSPhi - halfAngTolerance) ||
197  (pPhi > fSPhi + fDPhi + halfAngTolerance))
198  {
199  return in = eOutside;
200  }
201 
202  else if (in == eInside) // else it's eSurface anyway already
203  {
204  if ((pPhi < fSPhi + halfAngTolerance) ||
205  (pPhi > fSPhi + fDPhi - halfAngTolerance))
206  {
207  in = eSurface;
208  }
209  }
210  }
211  else if (!fPhiFullCone)
212  {
213  in = eSurface;
214  }
215 
216  return in;
217  }
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
VUSolid::EnumInside
EnumInside
Definition: VUSolid.hh:23
UVector3::z
double z
Definition: UVector3.hh:138
UVector3::y
double y
Definition: UVector3.hh:137
bool UCons::Normal ( const UVector3 &  p,
UVector3 &  n 
) const
virtual

Implements VUSolid.

Definition at line 175 of file UCons.cc.

References UUtils::Exception(), VUSolid::Tolerance(), UVector3::Unit(), Warning, UVector3::x, UVector3::y, and UVector3::z.

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 }
Warning
Definition: UUtils.hh:36
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
UVector3::Unit
UVector3 Unit() const
Definition: UVector3.cc:80
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
UVector3::y
double y
Definition: UVector3.hh:137
UCons & UCons::operator= ( const UCons &  rhs)

Definition at line 131 of file UCons.cc.

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  //
146  kRadTolerance = rhs.kRadTolerance;
147  kAngTolerance = rhs.kAngTolerance;
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;
163  fPhiFullCone = rhs.fPhiFullCone;
164 
165  Initialize();
166  return *this;
167 }
double UCons::SafetyFromInside ( const UVector3 &  p,
bool  precise 
) const
virtual

Implements VUSolid.

Definition at line 1999 of file UCons.cc.

References python.hepunit::degree, VUSolid::eOutside, UUtils::Exception(), Inside(), Warning, UVector3::x, UVector3::y, and UVector3::z.

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 }
UCons::Inside
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: UCons.hh:127
Warning
Definition: UUtils.hh:36
UVector3::x
double x
Definition: UVector3.hh:136
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
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
UVector3::y
double y
Definition: UVector3.hh:137
double UCons::SafetyFromOutside ( const UVector3 &  p,
bool  precise 
) const
virtual

Implements VUSolid.

Definition at line 1260 of file UCons.cc.

References UVector3::x, UVector3::y, and UVector3::z.

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 }
UVector3::x
double x
Definition: UVector3.hh:136
UVector3::z
double z
Definition: UVector3.hh:138
UVector3::y
double y
Definition: UVector3.hh:137
void UCons::SetDeltaPhiAngle ( double  newDPhi)
inline

Referenced by G4UCons::SetDeltaPhiAngle().

void UCons::SetInnerRadiusMinusZ ( double  Rmin1)
inline

Referenced by G4UCons::SetInnerRadiusMinusZ().

void UCons::SetInnerRadiusPlusZ ( double  Rmin2)
inline

Referenced by G4UCons::SetInnerRadiusPlusZ().

void UCons::SetOuterRadiusMinusZ ( double  Rmax1)
inline

Referenced by G4UCons::SetOuterRadiusMinusZ().

void UCons::SetOuterRadiusPlusZ ( double  Rmax2)
inline

Referenced by G4UCons::SetOuterRadiusPlusZ().

void UCons::SetStartPhiAngle ( double  newSPhi,
bool  trig = true 
)
inline

Referenced by G4UCons::SetStartPhiAngle().

void UCons::SetZHalfLength ( double  newDz)
inline

Referenced by G4UCons::SetZHalfLength().

std::ostream & UCons::StreamInfo ( std::ostream &  os) const
virtual

Implements VUSolid.

Definition at line 2108 of file UCons.cc.

References VUSolid::GetName().

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 }
VUSolid::GetName
const std::string & GetName() const
Definition: VUSolid.hh:103
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:57:31 for Geant4.10 by   doxygen 1.8.7