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

#include <USphere.hh>

Inheritance diagram for USphere:
VUSolid

Public Member Functions

 USphere (const std::string &pName, double pRmin, double pRmax, double pSPhi, double pDPhi, double pSTheta, double pDTheta)
 
 ~USphere ()
 
double GetInnerRadius () const
 
double GetOuterRadius () const
 
double GetStartPhiAngle () const
 
double GetDeltaPhiAngle () const
 
double GetStartThetaAngle () const
 
double GetDeltaThetaAngle () const
 
void SetInnerRadius (double newRMin)
 
void SetOuterRadius (double newRmax)
 
void SetStartPhiAngle (double newSphi, bool trig=true)
 
void SetDeltaPhiAngle (double newDphi)
 
void SetStartThetaAngle (double newSTheta)
 
void SetDeltaThetaAngle (double newDTheta)
 
double Capacity ()
 
double SurfaceArea ()
 
VUSolid::EnumInside Inside (const UVector3 &p) const
 
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 aAccurate=false) const
 
double DistanceToOut (const UVector3 &p, const UVector3 &v, UVector3 &n, bool &validNorm, double aPstep=UUtils::kInfinity) const
 
double SafetyFromInside (const UVector3 &p, bool aAccurate=false) const
 
UGeometryType GetEntityType () const
 
UVector3 GetPointOnSurface () const
 
VUSolid * Clone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
UVisExtent GetExtent () const
 
void Extent (UVector3 &aMin, UVector3 &aMax) const
 
void GetParametersList (int, double *) const
 
virtual void ComputeBBox (UBBox *, bool)
 
 USphere (const USphere &rhs)
 
USphere & operator= (const USphere &rhs)
 
double GetRmin () const
 
double GetRmax () const
 
double GetSPhi () const
 
double GetDPhi () const
 
double GetSTheta () const
 
double GetDTheta () const
 
double GetInsideRadius () const
 
void SetInsideRadius (double newRmin)
 
- 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 55 of file USphere.hh.

Constructor & Destructor Documentation

USphere::USphere ( const std::string &  pName,
double  pRmin,
double  pRmax,
double  pSPhi,
double  pDPhi,
double  pSTheta,
double  pDTheta 
)

Definition at line 34 of file USphere.cc.

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

Referenced by Clone().

38  : VUSolid(pName), fCubicVolume(0.),
39  fSurfaceArea(0.),fEpsilon(2.e-11),
40  fFullPhiSphere(true), fFullThetaSphere(true)
41 {
42  kAngTolerance = faTolerance;
43 
44  // Check radii and Set radial tolerances
45 
46  kRadTolerance = frTolerance;
47  if ((pRmin >= pRmax) || (pRmax < 1.1 * kRadTolerance) || (pRmin < 0))
48  {
49  std::ostringstream message;
50  message << "Invalid radii for Solid: " << GetName() << std::endl
51  << " pRmin = " << pRmin << ", pRmax = " << pRmax;
52  UUtils::Exception("USphere::USphere()", "GeomSolids0002",
53  FatalErrorInArguments, 1, message.str().c_str());
54  }
55  fRmin = pRmin;
56  fRmax = pRmax;
57  fRminTolerance = (fRmin) ? std::max(kRadTolerance, fEpsilon * fRmin) : 0;
58  kTolerance = std::max(kRadTolerance, fEpsilon * fRmax);
59 
60  // Check angles
61 
62  CheckPhiAngles(pSPhi, pDPhi);
63  CheckThetaAngles(pSTheta, pDTheta);
64 }
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
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
UUtils::Exception
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:177
USphere::~USphere ( )

Definition at line 70 of file USphere.cc.

71 {
72 }
USphere::USphere ( const USphere &  rhs)

Definition at line 78 of file USphere.cc.

79  : VUSolid(rhs), fCubicVolume(0.),fSurfaceArea(0.),
80  fRminTolerance(rhs.fRminTolerance),
81  kTolerance(rhs.kTolerance), kAngTolerance(rhs.kAngTolerance),
82  kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
83  fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
84  fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
85  sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
86  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
87  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
88  sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
89  hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
90  sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
91  sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
92  tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
93  tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
94  fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
95  fFullSphere(rhs.fFullSphere)
96 {
97 }
VUSolid::VUSolid
VUSolid()
Definition: VUSolid.cc:18

Member Function Documentation

double USphere::Capacity ( )
inlinevirtual

Implements VUSolid.

Definition at line 489 of file USphere.hh.

490 {
491  if (fCubicVolume != 0.)
492  {
493  ;
494  }
495  else
496  {
497  fCubicVolume = fDPhi * (std::cos(fSTheta) - std::cos(fSTheta + fDTheta)) *
498  (fRmax * fRmax * fRmax - fRmin * fRmin * fRmin) / 3.;
499  }
500  return fCubicVolume;
501 }
VUSolid * USphere::Clone ( ) const
virtual

Implements VUSolid.

Definition at line 2805 of file USphere.cc.

References USphere().

2806 {
2807  return new USphere(*this);
2808 }
USphere::USphere
USphere(const std::string &pName, double pRmin, double pRmax, double pSPhi, double pDPhi, double pSTheta, double pDTheta)
Definition: USphere.cc:34
virtual void USphere::ComputeBBox ( UBBox *  ,
bool   
)
inlinevirtual

Implements VUSolid.

Definition at line 129 of file USphere.hh.

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

Implements VUSolid.

Definition at line 611 of file USphere.cc.

References test::b, test::c, UUtils::Infinity(), plottest35::t1, VUSolid::Tolerance(), UVector3::x, UVector3::y, and UVector3::z.

612 {
613  double snxt = UUtils::Infinity(); // snxt = default return value
614  double rho2, rad2, pDotV2d, pDotV3d, pTheta;
615  double tolSTheta = 0., tolETheta = 0.;
616  const double dRmax = 100.*fRmax;
617 
618  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
619  static const double halfAngTolerance = kAngTolerance * 0.5;
620  const double halfTolerance = kTolerance * 0.5;
621  const double halfRminTolerance = fRminTolerance * 0.5;
622  const double tolORMin2 = (fRmin > halfRminTolerance)
623  ? (fRmin - halfRminTolerance) * (fRmin - halfRminTolerance) : 0;
624  const double tolIRMin2 =
625  (fRmin + halfRminTolerance) * (fRmin + halfRminTolerance);
626  const double tolORMax2 =
627  (fRmax + halfTolerance) * (fRmax + halfTolerance);
628  const double tolIRMax2 =
629  (fRmax - halfTolerance) * (fRmax - halfTolerance);
630 
631  // Intersection point
632  //
633  double xi, yi, zi, rhoi, rhoi2, radi2, iTheta;
634 
635  // Phi intersection
636  //
637  double Comp;
638 
639  // Phi precalcs
640  //
641  double Dist, cosPsi;
642 
643  // Theta precalcs
644  //
645  double dist2STheta, dist2ETheta;
646  double t1, t2, b, c, d2, d, sd = UUtils::Infinity();
647 
648  // General Precalcs
649  //
650  rho2 = p.x * p.x + p.y * p.y;
651  rad2 = rho2 + p.z * p.z;
652  pTheta = std::atan2(std::sqrt(rho2), p.z);
653 
654  pDotV2d = p.x * v.x + p.y * v.y;
655  pDotV3d = pDotV2d + p.z * v.z;
656 
657  // Theta precalcs
658  //
659  if (!fFullThetaSphere)
660  {
661  tolSTheta = fSTheta - halfAngTolerance;
662  tolETheta = eTheta + halfAngTolerance;
663  }
664 
665  // Outer spherical shell intersection
666  // - Only if outside tolerant fRmax
667  // - Check for if inside and outer USphere heading through solid (-> 0)
668  // - No intersect -> no intersection with USphere
669  //
670  // Shell eqn: x^2+y^2+z^2=RSPH^2
671  //
672  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
673  //
674  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
675  // => rad2 +2sd(pDotV3d) +sd^2 =R^2
676  //
677  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
678 
679  c = rad2 - fRmax * fRmax;
680 
681  if (c > kTolerance * fRmax)
682  {
683  // If outside tolerant boundary of outer USphere
684  // [should be std::sqrt(rad2)-fRmax > halfTolerance]
685 
686  d2 = pDotV3d * pDotV3d - c;
687 
688  if (d2 >= 0)
689  {
690  sd = -pDotV3d - std::sqrt(d2);
691 
692  if (sd >= 0)
693  {
694  if (sd > dRmax) // Avoid rounding errors due to precision issues seen on
695  {
696  // 64 bits systems. Split long distances and recompute
697  double fTerm = sd - std::fmod(sd, dRmax);
698  sd = fTerm + DistanceToIn(p + fTerm * v, v);
699  }
700  xi = p.x + sd * v.x;
701  yi = p.y + sd * v.y;
702  rhoi = std::sqrt(xi * xi + yi * yi);
703 
704  if (!fFullPhiSphere && rhoi) // Check phi intersection
705  {
706  cosPsi = (xi * cosCPhi + yi * sinCPhi) / rhoi;
707 
708  if (cosPsi >= cosHDPhiOT)
709  {
710  if (!fFullThetaSphere) // Check theta intersection
711  {
712  zi = p.z + sd * v.z;
713 
714  // rhoi & zi can never both be 0
715  // (=>intersect at origin =>fRmax=0)
716  //
717  iTheta = std::atan2(rhoi, zi);
718  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
719  {
720  return snxt = sd;
721  }
722  }
723  else
724  {
725  return snxt = sd;
726  }
727  }
728  }
729  else
730  {
731  if (!fFullThetaSphere) // Check theta intersection
732  {
733  zi = p.z + sd * v.z;
734 
735  // rhoi & zi can never both be 0
736  // (=>intersect at origin => fRmax=0 !)
737  //
738  iTheta = std::atan2(rhoi, zi);
739  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
740  {
741  return snxt = sd;
742  }
743  }
744  else
745  {
746  return snxt = sd;
747  }
748  }
749  }
750  }
751  else // No intersection with USphere
752  {
753  return snxt = UUtils::Infinity();
754  }
755  }
756  else
757  {
758  // Inside outer radius
759  // check not inside, and heading through USphere (-> 0 to in)
760 
761  d2 = pDotV3d * pDotV3d - c;
762 
763  if ((rad2 > tolIRMax2)
764  && ((d2 >= kTolerance * fRmax) && (pDotV3d < 0)))
765  {
766  if (!fFullPhiSphere)
767  {
768  // Use inner phi tolerant boundary -> if on tolerant
769  // phi boundaries, phi intersect code handles leaving/entering checks
770 
771  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
772 
773  if (cosPsi >= cosHDPhiIT)
774  {
775  // inside radii, delta r -ve, inside phi
776 
777  if (!fFullThetaSphere)
778  {
779  if ((pTheta >= tolSTheta + kAngTolerance)
780  && (pTheta <= tolETheta - kAngTolerance))
781  {
782  return snxt = 0;
783  }
784  }
785  else // strictly inside Theta in both cases
786  {
787  return snxt = 0;
788  }
789  }
790  }
791  else
792  {
793  if (!fFullThetaSphere)
794  {
795  if ((pTheta >= tolSTheta + kAngTolerance)
796  && (pTheta <= tolETheta - kAngTolerance))
797  {
798  return snxt = 0;
799  }
800  }
801  else // strictly inside Theta in both cases
802  {
803  return snxt = 0;
804  }
805  }
806  }
807  }
808 
809  // Inner spherical shell intersection
810  // - Always farthest root, because would have passed through outer
811  // surface first.
812  // - Tolerant check if travelling through solid
813 
814  if (fRmin)
815  {
816  c = rad2 - fRmin * fRmin;
817  d2 = pDotV3d * pDotV3d - c;
818 
819  // Within tolerance inner radius of inner USphere
820  // Check for immediate entry/already inside and travelling outwards
821 
822  if ((c > -halfRminTolerance) && (rad2 < tolIRMin2)
823  && ((d2 < fRmin * VUSolid::Tolerance()) || (pDotV3d >= 0)))
824  {
825  if (!fFullPhiSphere)
826  {
827  // Use inner phi tolerant boundary -> if on tolerant
828  // phi boundaries, phi intersect code handles leaving/entering checks
829 
830  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
831  if (cosPsi >= cosHDPhiIT)
832  {
833  // inside radii, delta r -ve, inside phi
834  //
835  if (!fFullThetaSphere)
836  {
837  if ((pTheta >= tolSTheta + kAngTolerance)
838  && (pTheta <= tolETheta - kAngTolerance))
839  {
840  return snxt = 0;
841  }
842  }
843  else
844  {
845  return snxt = 0;
846  }
847  }
848  }
849  else
850  {
851  if (!fFullThetaSphere)
852  {
853  if ((pTheta >= tolSTheta + kAngTolerance)
854  && (pTheta <= tolETheta - kAngTolerance))
855  {
856  return snxt = 0;
857  }
858  }
859  else
860  {
861  return snxt = 0;
862  }
863  }
864  }
865  else // Not special tolerant case
866  {
867  if (d2 >= 0)
868  {
869  sd = -pDotV3d + std::sqrt(d2);
870  if (sd >= halfRminTolerance) // It was >= 0 ??
871  {
872  xi = p.x + sd * v.x;
873  yi = p.y + sd * v.y;
874  rhoi = std::sqrt(xi * xi + yi * yi);
875 
876  if (!fFullPhiSphere && rhoi) // Check phi intersection
877  {
878  cosPsi = (xi * cosCPhi + yi * sinCPhi) / rhoi;
879 
880  if (cosPsi >= cosHDPhiOT)
881  {
882  if (!fFullThetaSphere) // Check theta intersection
883  {
884  zi = p.z + sd * v.z;
885 
886  // rhoi & zi can never both be 0
887  // (=>intersect at origin =>fRmax=0)
888  //
889  iTheta = std::atan2(rhoi, zi);
890  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
891  {
892  snxt = sd;
893  }
894  }
895  else
896  {
897  snxt = sd;
898  }
899  }
900  }
901  else
902  {
903  if (!fFullThetaSphere) // Check theta intersection
904  {
905  zi = p.z + sd * v.z;
906 
907  // rhoi & zi can never both be 0
908  // (=>intersect at origin => fRmax=0 !)
909  //
910  iTheta = std::atan2(rhoi, zi);
911  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
912  {
913  snxt = sd;
914  }
915  }
916  else
917  {
918  snxt = sd;
919  }
920  }
921  }
922  }
923  }
924  }
925 
926  // Phi segment intersection
927  //
928  // o Tolerant of points inside phi planes by up to VUSolid::Tolerance()*0.5
929  //
930  // o NOTE: Large duplication of code between sphi & ephi checks
931  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
932  // intersection check <=0 -> >=0
933  // -> Should use some form of loop Construct
934  //
935  if (!fFullPhiSphere)
936  {
937  // First phi surface ('S'tarting phi)
938  // Comp = Component in outwards normal dirn
939  //
940  Comp = v.x * sinSPhi - v.y * cosSPhi;
941 
942  if (Comp < 0)
943  {
944  Dist = p.y * cosSPhi - p.x * sinSPhi;
945 
946  if (Dist < halfCarTolerance)
947  {
948  sd = Dist / Comp;
949 
950  if (sd < snxt)
951  {
952  if (sd > 0)
953  {
954  xi = p.x + sd * v.x;
955  yi = p.y + sd * v.y;
956  zi = p.z + sd * v.z;
957  rhoi2 = xi * xi + yi * yi ;
958  radi2 = rhoi2 + zi * zi ;
959  }
960  else
961  {
962  sd = 0;
963  xi = p.x;
964  yi = p.y;
965  zi = p.z;
966  rhoi2 = rho2;
967  radi2 = rad2;
968  }
969  if ((radi2 <= tolORMax2)
970  && (radi2 >= tolORMin2)
971  && ((yi * cosCPhi - xi * sinCPhi) <= 0))
972  {
973  // Check theta intersection
974  // rhoi & zi can never both be 0
975  // (=>intersect at origin =>fRmax=0)
976  //
977  if (!fFullThetaSphere)
978  {
979  iTheta = std::atan2(std::sqrt(rhoi2), zi);
980  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
981  {
982  // r and theta intersections good
983  // - check intersecting with correct half-plane
984 
985  if ((yi * cosCPhi - xi * sinCPhi) <= 0)
986  {
987  snxt = sd;
988  }
989  }
990  }
991  else
992  {
993  snxt = sd;
994  }
995  }
996  }
997  }
998  }
999 
1000  // Second phi surface ('E'nding phi)
1001  // Component in outwards normal dirn
1002 
1003  Comp = -(v.x * sinEPhi - v.y * cosEPhi);
1004 
1005  if (Comp < 0)
1006  {
1007  Dist = -(p.y * cosEPhi - p.x * sinEPhi);
1008  if (Dist < halfCarTolerance)
1009  {
1010  sd = Dist / Comp;
1011 
1012  if (sd < snxt)
1013  {
1014  if (sd > 0)
1015  {
1016  xi = p.x + sd * v.x;
1017  yi = p.y + sd * v.y;
1018  zi = p.z + sd * v.z;
1019  rhoi2 = xi * xi + yi * yi ;
1020  radi2 = rhoi2 + zi * zi ;
1021  }
1022  else
1023  {
1024  sd = 0 ;
1025  xi = p.x;
1026  yi = p.y;
1027  zi = p.z;
1028  rhoi2 = rho2 ;
1029  radi2 = rad2 ;
1030  }
1031  if ((radi2 <= tolORMax2)
1032  && (radi2 >= tolORMin2)
1033  && ((yi * cosCPhi - xi * sinCPhi) >= 0))
1034  {
1035  // Check theta intersection
1036  // rhoi & zi can never both be 0
1037  // (=>intersect at origin =>fRmax=0)
1038  //
1039  if (!fFullThetaSphere)
1040  {
1041  iTheta = std::atan2(std::sqrt(rhoi2), zi);
1042  if ((iTheta >= tolSTheta) && (iTheta <= tolETheta))
1043  {
1044  // r and theta intersections good
1045  // - check intersecting with correct half-plane
1046 
1047  if ((yi * cosCPhi - xi * sinCPhi) >= 0)
1048  {
1049  snxt = sd;
1050  }
1051  }
1052  }
1053  else
1054  {
1055  snxt = sd;
1056  }
1057  }
1058  }
1059  }
1060  }
1061  }
1062 
1063  // Theta segment intersection
1064 
1065  if (!fFullThetaSphere)
1066  {
1067 
1068  // Intersection with theta surfaces
1069  // Known failure cases:
1070  // o Inside tolerance of stheta surface, skim
1071  // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1072  //
1073  // To solve: Check 2nd root of etheta surface in addition to stheta
1074  //
1075  // o start/end theta is exactly UUtils::kPi/2
1076  // Intersections with cones
1077  //
1078  // Cone equation: x^2+y^2=z^2tan^2(t)
1079  //
1080  // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1081  //
1082  // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1083  // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1084  //
1085  // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1086 
1087  if (fSTheta)
1088  {
1089  dist2STheta = rho2 - p.z * p.z * tanSTheta2;
1090  }
1091  else
1092  {
1093  dist2STheta = UUtils::Infinity();
1094  }
1095  if (eTheta < UUtils::kPi)
1096  {
1097  dist2ETheta = rho2 - p.z * p.z * tanETheta2;
1098  }
1099  else
1100  {
1101  dist2ETheta = UUtils::Infinity();
1102  }
1103  if (pTheta < tolSTheta)
1104  {
1105  // Inside (theta<stheta-tol) stheta cone
1106  // First root of stheta cone, second if first root -ve
1107 
1108  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1109  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1110  if (t1)
1111  {
1112  b = t2 / t1;
1113  c = dist2STheta / t1;
1114  d2 = b * b - c;
1115 
1116  if (d2 >= 0)
1117  {
1118  d = std::sqrt(d2);
1119  sd = -b - d; // First root
1120  zi = p.z + sd * v.z;
1121 
1122  if ((sd < 0) || (zi * (fSTheta - UUtils::kPi / 2) > 0))
1123  {
1124  sd = -b + d; // Second root
1125  }
1126  if ((sd >= 0) && (sd < snxt))
1127  {
1128  xi = p.x + sd * v.x;
1129  yi = p.y + sd * v.y;
1130  zi = p.z + sd * v.z;
1131  rhoi2 = xi * xi + yi * yi;
1132  radi2 = rhoi2 + zi * zi;
1133  if ((radi2 <= tolORMax2)
1134  && (radi2 >= tolORMin2)
1135  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1136  {
1137  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1138  {
1139  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1140  if (cosPsi >= cosHDPhiOT)
1141  {
1142  snxt = sd;
1143  }
1144  }
1145  else
1146  {
1147  snxt = sd;
1148  }
1149  }
1150  }
1151  }
1152  }
1153 
1154  // Possible intersection with ETheta cone.
1155  // Second >= 0 root should be considered
1156 
1157  if (eTheta < UUtils::kPi)
1158  {
1159  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1160  t2 = pDotV2d - p.z * v.z * tanETheta2;
1161  if (t1)
1162  {
1163  b = t2 / t1;
1164  c = dist2ETheta / t1;
1165  d2 = b * b - c;
1166 
1167  if (d2 >= 0)
1168  {
1169  d = std::sqrt(d2);
1170  sd = -b + d; // Second root
1171 
1172  if ((sd >= 0) && (sd < snxt))
1173  {
1174  xi = p.x + sd * v.x;
1175  yi = p.y + sd * v.y;
1176  zi = p.z + sd * v.z;
1177  rhoi2 = xi * xi + yi * yi;
1178  radi2 = rhoi2 + zi * zi;
1179 
1180  if ((radi2 <= tolORMax2)
1181  && (radi2 >= tolORMin2)
1182  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1183  {
1184  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1185  {
1186  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1187  if (cosPsi >= cosHDPhiOT)
1188  {
1189  snxt = sd;
1190  }
1191  }
1192  else
1193  {
1194  snxt = sd;
1195  }
1196  }
1197  }
1198  }
1199  }
1200  }
1201  }
1202  else if (pTheta > tolETheta)
1203  {
1204  // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
1205  // Inside (theta > etheta+tol) e-theta cone
1206  // First root of etheta cone, second if first root 'imaginary'
1207 
1208  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1209  t2 = pDotV2d - p.z * v.z * tanETheta2;
1210  if (t1)
1211  {
1212  b = t2 / t1;
1213  c = dist2ETheta / t1;
1214  d2 = b * b - c;
1215 
1216  if (d2 >= 0)
1217  {
1218  d = std::sqrt(d2);
1219  sd = -b - d; // First root
1220  zi = p.z + sd * v.z;
1221 
1222  if ((sd < 0) || (zi * (eTheta - UUtils::kPi / 2) > 0))
1223  {
1224  sd = -b + d; // second root
1225  }
1226  if ((sd >= 0) && (sd < snxt))
1227  {
1228  xi = p.x + sd * v.x;
1229  yi = p.y + sd * v.y;
1230  zi = p.z + sd * v.z;
1231  rhoi2 = xi * xi + yi * yi;
1232  radi2 = rhoi2 + zi * zi;
1233 
1234  if ((radi2 <= tolORMax2)
1235  && (radi2 >= tolORMin2)
1236  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1237  {
1238  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1239  {
1240  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1241  if (cosPsi >= cosHDPhiOT)
1242  {
1243  snxt = sd;
1244  }
1245  }
1246  else
1247  {
1248  snxt = sd;
1249  }
1250  }
1251  }
1252  }
1253  }
1254 
1255  // Possible intersection with STheta cone.
1256  // Second >= 0 root should be considered
1257 
1258  if (fSTheta)
1259  {
1260  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1261  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1262  if (t1)
1263  {
1264  b = t2 / t1;
1265  c = dist2STheta / t1;
1266  d2 = b * b - c;
1267 
1268  if (d2 >= 0)
1269  {
1270  d = std::sqrt(d2);
1271  sd = -b + d; // Second root
1272 
1273  if ((sd >= 0) && (sd < snxt))
1274  {
1275  xi = p.x + sd * v.x;
1276  yi = p.y + sd * v.y;
1277  zi = p.z + sd * v.z;
1278  rhoi2 = xi * xi + yi * yi;
1279  radi2 = rhoi2 + zi * zi;
1280 
1281  if ((radi2 <= tolORMax2)
1282  && (radi2 >= tolORMin2)
1283  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1284  {
1285  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1286  {
1287  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1288  if (cosPsi >= cosHDPhiOT)
1289  {
1290  snxt = sd;
1291  }
1292  }
1293  else
1294  {
1295  snxt = sd;
1296  }
1297  }
1298  }
1299  }
1300  }
1301  }
1302  }
1303  else if ((pTheta < tolSTheta + kAngTolerance)
1304  && (fSTheta > halfAngTolerance))
1305  {
1306  // In tolerance of stheta
1307  // If entering through solid [r,phi] => 0 to in
1308  // else try 2nd root
1309 
1310  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1311  if ((t2 >= 0 && tolIRMin2 < rad2 && rad2 < tolIRMax2 && fSTheta < UUtils::kPi / 2)
1312  || (t2 < 0 && tolIRMin2 < rad2 && rad2 < tolIRMax2 && fSTheta > UUtils::kPi / 2)
1313  || (v.z < 0 && tolIRMin2 < rad2 && rad2 < tolIRMax2 && fSTheta == UUtils::kPi / 2))
1314  {
1315  if (!fFullPhiSphere && rho2) // Check phi intersection
1316  {
1317  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
1318  if (cosPsi >= cosHDPhiIT)
1319  {
1320  return 0;
1321  }
1322  }
1323  else
1324  {
1325  return 0;
1326  }
1327  }
1328 
1329  // Not entering immediately/travelling through
1330 
1331  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1332  if (t1)
1333  {
1334  b = t2 / t1;
1335  c = dist2STheta / t1;
1336  d2 = b * b - c;
1337 
1338  if (d2 >= 0)
1339  {
1340  d = std::sqrt(d2);
1341  sd = -b + d;
1342  if ((sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < UUtils::kPi / 2))
1343  {
1344  // ^^^^^^^^^^^^^^^^^^^^^ shouldn't it be >=0 instead ?
1345  xi = p.x + sd * v.x;
1346  yi = p.y + sd * v.y;
1347  zi = p.z + sd * v.z;
1348  rhoi2 = xi * xi + yi * yi;
1349  radi2 = rhoi2 + zi * zi;
1350 
1351  if ((radi2 <= tolORMax2)
1352  && (radi2 >= tolORMin2)
1353  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1354  {
1355  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1356  {
1357  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1358  if (cosPsi >= cosHDPhiOT)
1359  {
1360  snxt = sd;
1361  }
1362  }
1363  else
1364  {
1365  snxt = sd;
1366  }
1367  }
1368  }
1369  }
1370  }
1371  }
1372  else if ((pTheta > tolETheta - kAngTolerance) && (eTheta < UUtils::kPi - kAngTolerance))
1373  {
1374 
1375  // In tolerance of etheta
1376  // If entering through solid [r,phi] => 0 to in
1377  // else try 2nd root
1378 
1379  t2 = pDotV2d - p.z * v.z * tanETheta2;
1380 
1381  if (((t2 < 0) && (eTheta < UUtils::kPi / 2)
1382  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1383  || ((t2 >= 0) && (eTheta > UUtils::kPi / 2)
1384  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1385  || ((v.z > 0) && (eTheta == UUtils::kPi / 2)
1386  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)))
1387  {
1388  if (!fFullPhiSphere && rho2) // Check phi intersection
1389  {
1390  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / std::sqrt(rho2);
1391  if (cosPsi >= cosHDPhiIT)
1392  {
1393  return 0;
1394  }
1395  }
1396  else
1397  {
1398  return 0;
1399  }
1400  }
1401 
1402  // Not entering immediately/travelling through
1403 
1404  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1405  if (t1)
1406  {
1407  b = t2 / t1;
1408  c = dist2ETheta / t1;
1409  d2 = b * b - c;
1410 
1411  if (d2 >= 0)
1412  {
1413  d = std::sqrt(d2);
1414  sd = -b + d;
1415 
1416  if ((sd >= halfCarTolerance)
1417  && (sd < snxt) && (eTheta > UUtils::kPi / 2))
1418  {
1419  xi = p.x + sd * v.x;
1420  yi = p.y + sd * v.y;
1421  zi = p.z + sd * v.z;
1422  rhoi2 = xi * xi + yi * yi;
1423  radi2 = rhoi2 + zi * zi;
1424 
1425  if ((radi2 <= tolORMax2)
1426  && (radi2 >= tolORMin2)
1427  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1428  {
1429  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1430  {
1431  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1432  if (cosPsi >= cosHDPhiOT)
1433  {
1434  snxt = sd;
1435  }
1436  }
1437  else
1438  {
1439  snxt = sd;
1440  }
1441  }
1442  }
1443  }
1444  }
1445  }
1446  else
1447  {
1448  // stheta+tol<theta<etheta-tol
1449  // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1450 
1451  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1452  t2 = pDotV2d - p.z * v.z * tanSTheta2;
1453  if (t1)
1454  {
1455  b = t2 / t1;
1456  c = dist2STheta / t1;
1457  d2 = b * b - c;
1458 
1459  if (d2 >= 0)
1460  {
1461  d = std::sqrt(d2);
1462  sd = -b + d; // second root
1463 
1464  if ((sd >= 0) && (sd < snxt))
1465  {
1466  xi = p.x + sd * v.x;
1467  yi = p.y + sd * v.y;
1468  zi = p.z + sd * v.z;
1469  rhoi2 = xi * xi + yi * yi;
1470  radi2 = rhoi2 + zi * zi;
1471 
1472  if ((radi2 <= tolORMax2)
1473  && (radi2 >= tolORMin2)
1474  && (zi * (fSTheta - UUtils::kPi / 2) <= 0))
1475  {
1476  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1477  {
1478  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1479  if (cosPsi >= cosHDPhiOT)
1480  {
1481  snxt = sd;
1482  }
1483  }
1484  else
1485  {
1486  snxt = sd;
1487  }
1488  }
1489  }
1490  }
1491  }
1492  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1493  t2 = pDotV2d - p.z * v.z * tanETheta2;
1494  if (t1)
1495  {
1496  b = t2 / t1;
1497  c = dist2ETheta / t1;
1498  d2 = b * b - c;
1499 
1500  if (d2 >= 0)
1501  {
1502  d = std::sqrt(d2);
1503  sd = -b + d; // second root
1504 
1505  if ((sd >= 0) && (sd < snxt))
1506  {
1507  xi = p.x + sd * v.x;
1508  yi = p.y + sd * v.y;
1509  zi = p.z + sd * v.z;
1510  rhoi2 = xi * xi + yi * yi;
1511  radi2 = rhoi2 + zi * zi;
1512 
1513  if ((radi2 <= tolORMax2)
1514  && (radi2 >= tolORMin2)
1515  && (zi * (eTheta - UUtils::kPi / 2) <= 0))
1516  {
1517  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1518  {
1519  cosPsi = (xi * cosCPhi + yi * sinCPhi) / std::sqrt(rhoi2);
1520  if (cosPsi >= cosHDPhiOT)
1521  {
1522  snxt = sd;
1523  }
1524  }
1525  else
1526  {
1527  snxt = sd;
1528  }
1529  }
1530  }
1531  }
1532  }
1533  }
1534  }
1535  return snxt;
1536 }
UUtils::Infinity
double Infinity()
Definition: UUtils.hh:177
USphere::DistanceToIn
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: USphere.cc:611
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
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 USphere::DistanceToOut ( const UVector3 &  p,
const UVector3 &  v,
UVector3 &  n,
bool &  validNorm,
double  aPstep = UUtils::kInfinity 
) const
virtual

Implements VUSolid.

Definition at line 1651 of file USphere.cc.

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

1652 {
1653  double snxt = UUtils::Infinity(); // snxt is default return value
1654  double sphi = UUtils::Infinity(), stheta = UUtils::Infinity();
1655  ESide side = kNull, sidephi = kNull, sidetheta = kNull;
1656 
1657  static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
1658  static const double halfAngTolerance = kAngTolerance * 0.5;
1659  const double halfTolerance = kTolerance * 0.5;
1660  const double halfRminTolerance = fRminTolerance * 0.5;
1661  const double rMaxPlus = fRmax + halfTolerance;
1662  const double Rmin_minus = (fRmin) ? fRmin - halfRminTolerance : 0;
1663  double t1, t2;
1664  double b, c, d;
1665 
1666  // Variables for phi intersection:
1667 
1668  double pDistS, compS, pDistE, compE, sphi2, vphi;
1669 
1670  double rho2, rad2, pDotV2d, pDotV3d;
1671 
1672  double xi, yi, zi; // Intersection point
1673 
1674  // Theta precals
1675  //
1676  double rhoSecTheta;
1677  double dist2STheta, dist2ETheta, distTheta;
1678  double d2, sd;
1679 
1680  // General Precalcs
1681  //
1682  rho2 = p.x * p.x + p.y * p.y;
1683  rad2 = rho2 + p.z * p.z;
1684 
1685  pDotV2d = p.x * v.x + p.y * v.y;
1686  pDotV3d = pDotV2d + p.z * v.z;
1687 
1688  // Radial Intersections from USphere::DistanceToIn
1689  //
1690  // Outer spherical shell intersection
1691  // - Only if outside tolerant fRmax
1692  // - Check for if inside and outer USphere heading through solid (-> 0)
1693  // - No intersect -> no intersection with USphere
1694  //
1695  // Shell eqn: x^2+y^2+z^2=RSPH^2
1696  //
1697  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
1698  //
1699  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
1700  // => rad2 +2sd(pDotV3d) +sd^2 =R^2
1701  //
1702  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
1703 
1704  if ((rad2 <= rMaxPlus * rMaxPlus) && (rad2 >= Rmin_minus * Rmin_minus))
1705  {
1706  c = rad2 - fRmax * fRmax;
1707 
1708  if (c < kTolerance * fRmax)
1709  {
1710  // Within tolerant Outer radius
1711  //
1712  // The test is
1713  // rad - fRmax < 0.5*kRadTolerance
1714  // => rad < fRmax + 0.5*kRadTol
1715  // => rad2 < (fRmax + 0.5*kRadTol)^2
1716  // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
1717  // => rad2 - fRmax^2 <~ fRmax*kRadTol
1718 
1719  d2 = pDotV3d * pDotV3d - c;
1720 
1721  if ((c > - kTolerance * fRmax) // on tolerant surface
1722  && ((pDotV3d >= 0) || (d2 < 0))) // leaving outside from Rmax
1723  // not re-entering
1724  {
1725  validNorm = true;
1726  n = UVector3(p.x / fRmax, p.y / fRmax, p.z / fRmax);
1727  return snxt = 0;
1728  }
1729  else
1730  {
1731  snxt = -pDotV3d + std::sqrt(d2); // second root since inside Rmax
1732  side = kRMax;
1733  }
1734  }
1735 
1736  // Inner spherical shell intersection:
1737  // Always first >=0 root, because would have passed
1738  // from outside of Rmin surface .
1739 
1740  if (fRmin)
1741  {
1742  c = rad2 - fRmin * fRmin;
1743  d2 = pDotV3d * pDotV3d - c;
1744 
1745  if (c > - fRminTolerance * fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
1746  {
1747  if ((c < fRminTolerance * fRmin) // leaving from Rmin
1748  && (d2 >= fRminTolerance * fRmin) && (pDotV3d < 0))
1749  {
1750  validNorm = false; // Rmin surface is concave
1751  return snxt = 0;
1752  }
1753  else
1754  {
1755  if (d2 >= 0.)
1756  {
1757  sd = -pDotV3d - std::sqrt(d2);
1758 
1759  if (sd >= 0.) // Always intersect Rmin first
1760  {
1761  snxt = sd;
1762  side = kRMin;
1763  }
1764  }
1765  }
1766  }
1767  }
1768  }
1769 
1770  // Theta segment intersection
1771 
1772  if (!fFullThetaSphere)
1773  {
1774  // Intersection with theta surfaces
1775  //
1776  // Known failure cases:
1777  // o Inside tolerance of stheta surface, skim
1778  // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1779  //
1780  // To solve: Check 2nd root of etheta surface in addition to stheta
1781  //
1782  // o start/end theta is exactly UUtils::kPi/2
1783  //
1784  // Intersections with cones
1785  //
1786  // Cone equation: x^2+y^2=z^2tan^2(t)
1787  //
1788  // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1789  //
1790  // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1791  // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1792  //
1793  // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1794  //
1795 
1796  if (fSTheta) // intersection with first cons
1797  {
1798  if (std::fabs(tanSTheta) > 5. / kAngTolerance) // kons is plane z=0
1799  {
1800  if (v.z > 0.)
1801  {
1802  if (std::fabs(p.z) <= halfTolerance)
1803  {
1804  validNorm = true;
1805  n = UVector3(0., 0., 1.);
1806  return snxt = 0;
1807  }
1808  stheta = -p.z / v.z;
1809  sidetheta = kSTheta;
1810  }
1811  }
1812  else // kons is not plane
1813  {
1814  t1 = 1 - v.z * v.z * (1 + tanSTheta2);
1815  t2 = pDotV2d - p.z * v.z * tanSTheta2; // ~vDotN if p on cons
1816  dist2STheta = rho2 - p.z * p.z * tanSTheta2; // t3
1817 
1818  distTheta = std::sqrt(rho2) - p.z * tanSTheta;
1819 
1820  if (std::fabs(t1) < halfAngTolerance) // 1st order equation,
1821  {
1822  // v parallel to kons
1823  if (v.z > 0.)
1824  {
1825  if (std::fabs(distTheta) < halfTolerance) // p on surface
1826  {
1827  if ((fSTheta < UUtils::kPi / 2) && (p.z > 0.))
1828  {
1829  validNorm = false;
1830  return snxt = 0.;
1831  }
1832  else if ((fSTheta > UUtils::kPi / 2) && (p.z <= 0))
1833  {
1834  validNorm = true;
1835  if (rho2)
1836  {
1837  rhoSecTheta = std::sqrt(rho2 * (1 + tanSTheta2));
1838 
1839  n = UVector3(p.x / rhoSecTheta,
1840  p.y / rhoSecTheta,
1841  std::sin(fSTheta));
1842  }
1843  else n = UVector3(0., 0., 1.);
1844  return snxt = 0.;
1845  }
1846  }
1847  stheta = -0.5 * dist2STheta / t2;
1848  sidetheta = kSTheta;
1849  }
1850  } // 2nd order equation, 1st root of fSTheta cone,
1851  else // 2nd if 1st root -ve
1852  {
1853  if (std::fabs(distTheta) < halfTolerance)
1854  {
1855  if ((fSTheta > UUtils::kPi / 2) && (t2 >= 0.)) // leave
1856  {
1857  validNorm = true;
1858  if (rho2)
1859  {
1860  rhoSecTheta = std::sqrt(rho2 * (1 + tanSTheta2));
1861 
1862  n = UVector3(p.x / rhoSecTheta,
1863  p.y / rhoSecTheta,
1864  std::sin(fSTheta));
1865  }
1866  else
1867  {
1868  n = UVector3(0., 0., 1.);
1869  }
1870  return snxt = 0.;
1871  }
1872  else if ((fSTheta < UUtils::kPi / 2) && (t2 < 0.) && (p.z >= 0.)) // leave
1873  {
1874  validNorm = false;
1875  return snxt = 0.;
1876  }
1877  }
1878  b = t2 / t1;
1879  c = dist2STheta / t1;
1880  d2 = b * b - c;
1881 
1882  if (d2 >= 0.)
1883  {
1884  d = std::sqrt(d2);
1885 
1886  if (fSTheta > UUtils::kPi / 2)
1887  {
1888  sd = -b - d; // First root
1889 
1890  if (((std::fabs(sd) < halfTolerance) && (t2 < 0.))
1891  || (sd < 0.) || ((sd > 0.) && (p.z + sd * v.z > 0.)))
1892  {
1893  sd = -b + d; // 2nd root
1894  }
1895  if ((sd > halfTolerance) && (p.z + sd * v.z <= 0.))
1896  {
1897  stheta = sd;
1898  sidetheta = kSTheta;
1899  }
1900  }
1901  else // sTheta < UUtils::kPi/2, concave surface, no normal
1902  {
1903  sd = -b - d; // First root
1904 
1905  if (((std::fabs(sd) < halfTolerance) && (t2 >= 0.))
1906  || (sd < 0.) || ((sd > 0.) && (p.z + sd * v.z < 0.)))
1907  {
1908  sd = -b + d; // 2nd root
1909  }
1910  if ((sd > halfTolerance) && (p.z + sd * v.z >= 0.))
1911  {
1912  stheta = sd;
1913  sidetheta = kSTheta;
1914  }
1915  }
1916  }
1917  }
1918  }
1919  }
1920  if (eTheta < UUtils::kPi) // intersection with second cons
1921  {
1922  if (std::fabs(tanETheta) > 5. / kAngTolerance) // kons is plane z=0
1923  {
1924  if (v.z < 0.)
1925  {
1926  if (std::fabs(p.z) <= halfTolerance)
1927  {
1928  validNorm = true;
1929  n = UVector3(0., 0., -1.);
1930  return snxt = 0;
1931  }
1932  sd = -p.z / v.z;
1933 
1934  if (sd < stheta)
1935  {
1936  stheta = sd;
1937  sidetheta = kETheta;
1938  }
1939  }
1940  }
1941  else // kons is not plane
1942  {
1943  t1 = 1 - v.z * v.z * (1 + tanETheta2);
1944  t2 = pDotV2d - p.z * v.z * tanETheta2; // ~vDotN if p on cons
1945  dist2ETheta = rho2 - p.z * p.z * tanETheta2; // t3
1946 
1947  distTheta = std::sqrt(rho2) - p.z * tanETheta;
1948 
1949  if (std::fabs(t1) < halfAngTolerance) // 1st order equation,
1950  {
1951  // v parallel to kons
1952  if (v.z < 0.)
1953  {
1954  if (std::fabs(distTheta) < halfTolerance) // p on surface
1955  {
1956  if ((eTheta > UUtils::kPi / 2) && (p.z < 0.))
1957  {
1958  validNorm = false;
1959  return snxt = 0.;
1960  }
1961  else if ((eTheta < UUtils::kPi / 2) && (p.z >= 0))
1962  {
1963  validNorm = true;
1964  if (rho2)
1965  {
1966  rhoSecTheta = std::sqrt(rho2 * (1 + tanETheta2));
1967  n = UVector3(p.x / rhoSecTheta,
1968  p.y / rhoSecTheta,
1969  -sinETheta);
1970  }
1971  else
1972  {
1973  n = UVector3(0., 0., -1.);
1974  }
1975  return snxt = 0.;
1976  }
1977  }
1978  sd = -0.5 * dist2ETheta / t2;
1979 
1980  if (sd < stheta)
1981  {
1982  stheta = sd;
1983  sidetheta = kETheta;
1984  }
1985  }
1986  } // 2nd order equation, 1st root of fSTheta cone
1987  else // 2nd if 1st root -ve
1988  {
1989  if (std::fabs(distTheta) < halfTolerance)
1990  {
1991  if ((eTheta < UUtils::kPi / 2) && (t2 >= 0.)) // leave
1992  {
1993  validNorm = true;
1994  if (rho2)
1995  {
1996  rhoSecTheta = std::sqrt(rho2 * (1 + tanETheta2));
1997  n = UVector3(p.x / rhoSecTheta,
1998  p.y / rhoSecTheta,
1999  -sinETheta);
2000  }
2001  else n = UVector3(0., 0., -1.);
2002  return snxt = 0.;
2003  }
2004  else if ((eTheta > UUtils::kPi / 2)
2005  && (t2 < 0.) && (p.z <= 0.)) // leave
2006  {
2007  validNorm = false;
2008  return snxt = 0.;
2009  }
2010  }
2011  b = t2 / t1;
2012  c = dist2ETheta / t1;
2013  d2 = b * b - c;
2014 
2015  if (d2 >= 0.)
2016  {
2017  d = std::sqrt(d2);
2018 
2019  if (eTheta < UUtils::kPi / 2)
2020  {
2021  sd = -b - d; // First root
2022 
2023  if (((std::fabs(sd) < halfTolerance) && (t2 < 0.))
2024  || (sd < 0.))
2025  {
2026  sd = -b + d; // 2nd root
2027  }
2028  if (sd > halfTolerance)
2029  {
2030  if (sd < stheta)
2031  {
2032  stheta = sd;
2033  sidetheta = kETheta;
2034  }
2035  }
2036  }
2037  else // sTheta+fDTheta > UUtils::kPi/2, concave surface, no normal
2038  {
2039  sd = -b - d; // First root
2040 
2041  if (((std::fabs(sd) < halfTolerance) && (t2 >= 0.))
2042  || (sd < 0.) || ((sd > 0.) && (p.z + sd * v.z > 0.)))
2043  {
2044  sd = -b + d; // 2nd root
2045  }
2046  if ((sd > halfTolerance) && (p.z + sd * v.z <= 0.))
2047  {
2048  if (sd < stheta)
2049  {
2050  stheta = sd;
2051  sidetheta = kETheta;
2052  }
2053  }
2054  }
2055  }
2056  }
2057  }
2058  }
2059 
2060  } // end theta intersections
2061 
2062  // Phi Intersection
2063 
2064  if (!fFullPhiSphere)
2065  {
2066  if (p.x || p.y) // Check if on z axis (rho not needed later)
2067  {
2068  // pDist -ve when inside
2069 
2070  pDistS = p.x * sinSPhi - p.y * cosSPhi;
2071  pDistE = -p.x * sinEPhi + p.y * cosEPhi;
2072 
2073  // Comp -ve when in direction of outwards normal
2074 
2075  compS = -sinSPhi * v.x + cosSPhi * v.y;
2076  compE = sinEPhi * v.x - cosEPhi * v.y;
2077  sidephi = kNull;
2078 
2079  if ((pDistS <= 0) && (pDistE <= 0))
2080  {
2081  // Inside both phi *full* planes
2082 
2083  if (compS < 0)
2084  {
2085  sphi = pDistS / compS;
2086  xi = p.x + sphi * v.x;
2087  yi = p.y + sphi * v.y;
2088 
2089  // Check intersection with correct half-plane (if not -> no intersect)
2090  //
2091  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2092  {
2093  vphi = std::atan2(v.y, v.x);
2094  sidephi = kSPhi;
2095  if (((fSPhi - halfAngTolerance) <= vphi)
2096  && ((ePhi + halfAngTolerance) >= vphi))
2097  {
2098  sphi = UUtils::Infinity();
2099  }
2100  }
2101  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
2102  {
2103  sphi = UUtils::Infinity();
2104  }
2105  else
2106  {
2107  sidephi = kSPhi;
2108  if (pDistS > -halfCarTolerance)
2109  {
2110  sphi = 0; // Leave by sphi
2111  }
2112  }
2113  }
2114  else
2115  {
2116  sphi = UUtils::Infinity();
2117  }
2118 
2119  if (compE < 0)
2120  {
2121  sphi2 = pDistE / compE;
2122  if (sphi2 < sphi) // Only check further if < starting phi intersection
2123  {
2124  xi = p.x + sphi2 * v.x;
2125  yi = p.y + sphi2 * v.y;
2126 
2127  // Check intersection with correct half-plane
2128  //
2129  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2130  {
2131  // Leaving via ending phi
2132  //
2133  vphi = std::atan2(v.y, v.x);
2134 
2135  if (!((fSPhi - halfAngTolerance <= vphi)
2136  && (fSPhi + fDPhi + halfAngTolerance >= vphi)))
2137  {
2138  sidephi = kEPhi;
2139  if (pDistE <= -halfCarTolerance)
2140  {
2141  sphi = sphi2;
2142  }
2143  else
2144  {
2145  sphi = 0.0;
2146  }
2147  }
2148  }
2149  else if ((yi * cosCPhi - xi * sinCPhi) >= 0) // Leaving via ending phi
2150  {
2151  sidephi = kEPhi;
2152  if (pDistE <= -halfCarTolerance)
2153  {
2154  sphi = sphi2;
2155  }
2156  else
2157  {
2158  sphi = 0;
2159  }
2160  }
2161  }
2162  }
2163  }
2164  else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
2165  {
2166  if (pDistS <= pDistE)
2167  {
2168  sidephi = kSPhi;
2169  }
2170  else
2171  {
2172  sidephi = kEPhi;
2173  }
2174  if (fDPhi > UUtils::kPi)
2175  {
2176  if ((compS < 0) && (compE < 0))
2177  {
2178  sphi = 0;
2179  }
2180  else
2181  {
2182  sphi = UUtils::Infinity();
2183  }
2184  }
2185  else
2186  {
2187  // if towards both >=0 then once inside (after error)
2188  // will remain inside
2189 
2190  if ((compS >= 0) && (compE >= 0))
2191  {
2192  sphi = UUtils::Infinity();
2193  }
2194  else
2195  {
2196  sphi = 0;
2197  }
2198  }
2199  }
2200  else if ((pDistS > 0) && (pDistE < 0))
2201  {
2202  // Outside full starting plane, inside full ending plane
2203 
2204  if (fDPhi > UUtils::kPi)
2205  {
2206  if (compE < 0)
2207  {
2208  sphi = pDistE / compE;
2209  xi = p.x + sphi * v.x;
2210  yi = p.y + sphi * v.y;
2211 
2212  // Check intersection in correct half-plane
2213  // (if not -> not leaving phi extent)
2214  //
2215  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2216  {
2217  vphi = std::atan2(v.y, v.x);
2218  sidephi = kSPhi;
2219  if (((fSPhi - halfAngTolerance) <= vphi)
2220  && ((ePhi + halfAngTolerance) >= vphi))
2221  {
2222  sphi = UUtils::Infinity();
2223  }
2224  }
2225  else if ((yi * cosCPhi - xi * sinCPhi) <= 0)
2226  {
2227  sphi = UUtils::Infinity();
2228  }
2229  else // Leaving via Ending phi
2230  {
2231  sidephi = kEPhi;
2232  if (pDistE > -halfCarTolerance)
2233  {
2234  sphi = 0.;
2235  }
2236  }
2237  }
2238  else
2239  {
2240  sphi = UUtils::Infinity();
2241  }
2242  }
2243  else
2244  {
2245  if (compS >= 0)
2246  {
2247  if (compE < 0)
2248  {
2249  sphi = pDistE / compE;
2250  xi = p.x + sphi * v.x;
2251  yi = p.y + sphi * v.y;
2252 
2253  // Check intersection in correct half-plane
2254  // (if not -> remain in extent)
2255  //
2256  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2257  {
2258  vphi = std::atan2(v.y, v.x);
2259  sidephi = kSPhi;
2260  if (((fSPhi - halfAngTolerance) <= vphi)
2261  && ((ePhi + halfAngTolerance) >= vphi))
2262  {
2263  sphi = UUtils::Infinity();
2264  }
2265  }
2266  else if ((yi * cosCPhi - xi * sinCPhi) <= 0)
2267  {
2268  sphi = UUtils::Infinity();
2269  }
2270  else // otherwise leaving via Ending phi
2271  {
2272  sidephi = kEPhi;
2273  }
2274  }
2275  else sphi = UUtils::Infinity();
2276  }
2277  else // leaving immediately by starting phi
2278  {
2279  sidephi = kSPhi;
2280  sphi = 0;
2281  }
2282  }
2283  }
2284  else
2285  {
2286  // Must be pDistS < 0 && pDistE > 0
2287  // Inside full starting plane, outside full ending plane
2288 
2289  if (fDPhi > UUtils::kPi)
2290  {
2291  if (compS < 0)
2292  {
2293  sphi = pDistS / compS;
2294  xi = p.x + sphi * v.x;
2295  yi = p.y + sphi * v.y;
2296 
2297  // Check intersection in correct half-plane
2298  // (if not -> not leaving phi extent)
2299  //
2300  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2301  {
2302  vphi = std::atan2(v.y, v.x);
2303  sidephi = kSPhi;
2304  if (((fSPhi - halfAngTolerance) <= vphi)
2305  && ((ePhi + halfAngTolerance) >= vphi))
2306  {
2307  sphi = UUtils::Infinity();
2308  }
2309  }
2310  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
2311  {
2312  sphi = UUtils::Infinity();
2313  }
2314  else // Leaving via Starting phi
2315  {
2316  sidephi = kSPhi;
2317  if (pDistS > -halfCarTolerance)
2318  {
2319  sphi = 0;
2320  }
2321  }
2322  }
2323  else
2324  {
2325  sphi = UUtils::Infinity();
2326  }
2327  }
2328  else
2329  {
2330  if (compE >= 0)
2331  {
2332  if (compS < 0)
2333  {
2334  sphi = pDistS / compS;
2335  xi = p.x + sphi * v.x;
2336  yi = p.y + sphi * v.y;
2337 
2338  // Check intersection in correct half-plane
2339  // (if not -> remain in extent)
2340  //
2341  if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
2342  {
2343  vphi = std::atan2(v.y, v.x);
2344  sidephi = kSPhi;
2345  if (((fSPhi - halfAngTolerance) <= vphi)
2346  && ((ePhi + halfAngTolerance) >= vphi))
2347  {
2348  sphi = UUtils::Infinity();
2349  }
2350  }
2351  else if ((yi * cosCPhi - xi * sinCPhi) >= 0)
2352  {
2353  sphi = UUtils::Infinity();
2354  }
2355  else // otherwise leaving via Starting phi
2356  {
2357  sidephi = kSPhi;
2358  }
2359  }
2360  else
2361  {
2362  sphi = UUtils::Infinity();
2363  }
2364  }
2365  else // leaving immediately by ending
2366  {
2367  sidephi = kEPhi;
2368  sphi = 0 ;
2369  }
2370  }
2371  }
2372  }
2373  else
2374  {
2375  // On z axis + travel not || to z axis -> if phi of vector direction
2376  // within phi of shape, Step limited by rmax, else Step =0
2377 
2378  if (v.x || v.y)
2379  {
2380  vphi = std::atan2(v.y, v.x);
2381  if ((fSPhi - halfAngTolerance < vphi) && (vphi < ePhi + halfAngTolerance))
2382  {
2383  sphi = UUtils::Infinity();
2384  }
2385  else
2386  {
2387  sidephi = kSPhi; // arbitrary
2388  sphi = 0;
2389  }
2390  }
2391  else // travel along z - no phi intersection
2392  {
2393  sphi = UUtils::Infinity();
2394  }
2395  }
2396  if (sphi < snxt) // Order intersecttions
2397  {
2398  snxt = sphi;
2399  side = sidephi;
2400  }
2401  }
2402  if (stheta < snxt) // Order intersections
2403  {
2404  snxt = stheta;
2405  side = sidetheta;
2406  }
2407 
2408  switch (side)
2409  {
2410  case kRMax:
2411  xi = p.x + snxt * v.x;
2412  yi = p.y + snxt * v.y;
2413  zi = p.z + snxt * v.z;
2414  n = UVector3(xi / fRmax, yi / fRmax, zi / fRmax);
2415  validNorm = true;
2416  break;
2417 
2418  case kRMin:
2419  validNorm = false; // Rmin is concave
2420  break;
2421 
2422  case kSPhi:
2423  if (fDPhi <= UUtils::kPi) // Normal to Phi-
2424  {
2425  n = UVector3(sinSPhi, -cosSPhi, 0);
2426  validNorm = true;
2427  }
2428  else
2429  {
2430  validNorm = false;
2431  }
2432  break;
2433 
2434  case kEPhi:
2435  if (fDPhi <= UUtils::kPi) // Normal to Phi+
2436  {
2437  n = UVector3(-sinEPhi, cosEPhi, 0);
2438  validNorm = true;
2439  }
2440  else
2441  {
2442  validNorm = false;
2443  }
2444  break;
2445 
2446  case kSTheta:
2447  if (fSTheta == UUtils::kPi / 2)
2448  {
2449  n = UVector3(0., 0., 1.);
2450  validNorm = true;
2451  }
2452  else if (fSTheta > UUtils::kPi / 2)
2453  {
2454  xi = p.x + snxt * v.x;
2455  yi = p.y + snxt * v.y;
2456  rho2 = xi * xi + yi * yi;
2457  if (rho2)
2458  {
2459  rhoSecTheta = std::sqrt(rho2 * (1 + tanSTheta2));
2460  n = UVector3(xi / rhoSecTheta, yi / rhoSecTheta,
2461  -tanSTheta / std::sqrt(1 + tanSTheta2));
2462  }
2463  else
2464  {
2465  n = UVector3(0., 0., 1.);
2466  }
2467  validNorm = true;
2468  }
2469  else
2470  {
2471  validNorm = false; // Concave STheta cone
2472  }
2473  break;
2474 
2475  case kETheta:
2476  if (eTheta == UUtils::kPi / 2)
2477  {
2478  n = UVector3(0., 0., -1.);
2479  validNorm = true;
2480  }
2481  else if (eTheta < UUtils::kPi / 2)
2482  {
2483  xi = p.x + snxt * v.x;
2484  yi = p.y + snxt * v.y;
2485  rho2 = xi * xi + yi * yi;
2486  if (rho2)
2487  {
2488  rhoSecTheta = std::sqrt(rho2 * (1 + tanETheta2));
2489  n = UVector3(xi / rhoSecTheta, yi / rhoSecTheta,
2490  -tanETheta / std::sqrt(1 + tanETheta2));
2491  }
2492  else
2493  {
2494  n = UVector3(0., 0., -1.);
2495  }
2496  validNorm = true;
2497  }
2498  else
2499  {
2500  validNorm = false; // Concave ETheta cone
2501  }
2502  break;
2503 
2504  default:
2505  cout << std::endl;
2506 
2507  std::ostringstream message;
2508  int oldprc = message.precision(16);
2509  message << "Undefined side for valid surface normal to solid."
2510  << std::endl
2511  << "Position:" << std::endl << std::endl
2512  << "p.x = " << p.x << " mm" << std::endl
2513  << "p.y = " << p.y << " mm" << std::endl
2514  << "p.z = " << p.z << " mm" << std::endl << std::endl
2515  << "Direction:" << std::endl << std::endl
2516  << "v.x = " << v.x << std::endl
2517  << "v.y = " << v.y << std::endl
2518  << "v.z = " << v.z << std::endl << std::endl
2519  << "Proposed distance :" << std::endl << std::endl
2520  << "snxt = " << snxt << " mm" << std::endl;
2521  message.precision(oldprc);
2522  UUtils::Exception("USphere::DistanceToOut(p,v,..)",
2523  "GeomSolids1002", Warning, 1, message.str().c_str());
2524  break;
2525  }
2526  if (snxt == UUtils::Infinity())
2527  {
2528  cout << std::endl;
2529 
2530  std::ostringstream message;
2531  int oldprc = message.precision(16);
2532  message << "Logic error: snxt = UUtils::Infinity() ???" << std::endl
2533  << "Position:" << std::endl << std::endl
2534  << "p.x = " << p.x << " mm" << std::endl
2535  << "p.y = " << p.y << " mm" << std::endl
2536  << "p.z = " << p.z << " mm" << std::endl << std::endl
2537  << "Rp = " << std::sqrt(p.x * p.x + p.y * p.y + p.z * p.z)
2538  << " mm" << std::endl << std::endl
2539  << "Direction:" << std::endl << std::endl
2540  << "v.x = " << v.x << std::endl
2541  << "v.y = " << v.y << std::endl
2542  << "v.z = " << v.z << std::endl << std::endl
2543  << "Proposed distance :" << std::endl << std::endl
2544  << "snxt = " << snxt << " mm" << std::endl;
2545  message.precision(oldprc);
2546  UUtils::Exception("USphere::DistanceToOut(p,v,..)",
2547  "GeomSolids1002", Warning, 1, message.str().c_str());
2548  }
2549 
2550  return snxt;
2551 }
Warning
Definition: UUtils.hh:36
UUtils::Infinity
double Infinity()
Definition: UUtils.hh:177
VUSolid::Tolerance
static double Tolerance()
Definition: VUSolid.hh:127
UVector3::x
double x
Definition: UVector3.hh:136
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
UVector3::z
double z
Definition: UVector3.hh:138
UVector3::y
double y
Definition: UVector3.hh:137
void USphere::Extent ( UVector3 &  aMin,
UVector3 &  aMax 
) const
virtual

Implements VUSolid.

Definition at line 2983 of file USphere.cc.

References UVector3::Set().

2984 {
2985  aMin.Set(-fRmax);
2986  aMax.Set(fRmax);
2987 }
UVector3::Set
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:245
double USphere::GetDeltaPhiAngle ( ) const
inline

Definition at line 230 of file USphere.hh.

Referenced by G4USphere::GetDeltaPhiAngle(), GetDPhi(), and GetParametersList().

231 {
232  return fDPhi;
233 }
double USphere::GetDeltaThetaAngle ( ) const
inline

Definition at line 241 of file USphere.hh.

Referenced by G4USphere::GetDeltaThetaAngle(), GetDTheta(), and GetParametersList().

242 {
243  return fDTheta;
244 }
double USphere::GetDPhi ( ) const
inline

Definition at line 471 of file USphere.hh.

References GetDeltaPhiAngle().

472 {
473  return GetDeltaPhiAngle();
474 }
USphere::GetDeltaPhiAngle
double GetDeltaPhiAngle() const
Definition: USphere.hh:230
double USphere::GetDTheta ( ) const
inline

Definition at line 483 of file USphere.hh.

References GetDeltaThetaAngle().

484 {
485  return GetDeltaThetaAngle();
486 }
USphere::GetDeltaThetaAngle
double GetDeltaThetaAngle() const
Definition: USphere.hh:241
UGeometryType USphere::GetEntityType ( ) const
virtual

Implements VUSolid.

Definition at line 2796 of file USphere.cc.

2797 {
2798  return std::string("Sphere");
2799 }
UVisExtent USphere::GetExtent ( ) const
double USphere::GetInnerRadius ( ) const
inline

Definition at line 212 of file USphere.hh.

Referenced by G4USphere::GetInnerRadius(), and GetParametersList().

213 {
214  return fRmin;
215 }
double USphere::GetInsideRadius ( ) const
inline

Definition at line 206 of file USphere.hh.

Referenced by GetRmin().

207 {
208  return fRmin;
209 }
double USphere::GetOuterRadius ( ) const
inline

Definition at line 218 of file USphere.hh.

Referenced by G4USphere::GetOuterRadius(), GetParametersList(), and GetRmax().

219 {
220  return fRmax;
221 }
void USphere::GetParametersList ( int  ,
double *  aArray 
) const
virtual

Implements VUSolid.

Definition at line 2989 of file USphere.cc.

References GetDeltaPhiAngle(), GetDeltaThetaAngle(), GetInnerRadius(), GetOuterRadius(), GetStartPhiAngle(), and GetStartThetaAngle().

2990 {
2991  aArray[0] = GetInnerRadius();
2992  aArray[1] = GetOuterRadius();
2993  aArray[2] = GetStartPhiAngle();
2994  aArray[3] = GetDeltaPhiAngle();
2995  aArray[4] = GetStartThetaAngle();
2996  aArray[5] = GetDeltaThetaAngle();
2997 }
USphere::GetInnerRadius
double GetInnerRadius() const
Definition: USphere.hh:212
USphere::GetDeltaThetaAngle
double GetDeltaThetaAngle() const
Definition: USphere.hh:241
USphere::GetOuterRadius
double GetOuterRadius() const
Definition: USphere.hh:218
USphere::GetStartThetaAngle
double GetStartThetaAngle() const
Definition: USphere.hh:236
USphere::GetStartPhiAngle
double GetStartPhiAngle() const
Definition: USphere.hh:224
USphere::GetDeltaPhiAngle
double GetDeltaPhiAngle() const
Definition: USphere.hh:230
UVector3 USphere::GetPointOnSurface ( ) const
virtual

Implements VUSolid.

Definition at line 2838 of file USphere.cc.

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

2839 {
2840  double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
2841  double height1, height2, slant1, slant2, costheta, sintheta, rRand;
2842 
2843  height1 = (fRmax - fRmin) * cosSTheta;
2844  height2 = (fRmax - fRmin) * cosETheta;
2845  slant1 = std::sqrt(UUtils::sqr((fRmax - fRmin) * sinSTheta) + height1 * height1);
2846  slant2 = std::sqrt(UUtils::sqr((fRmax - fRmin) * sinETheta) + height2 * height2);
2847  rRand = UUtils::GetRadiusInRing(fRmin, fRmax);
2848 
2849  aOne = fRmax * fRmax * fDPhi * (cosSTheta - cosETheta);
2850  aTwo = fRmin * fRmin * fDPhi * (cosSTheta - cosETheta);
2851  aThr = fDPhi * ((fRmax + fRmin) * sinSTheta) * slant1;
2852  aFou = fDPhi * ((fRmax + fRmin) * sinETheta) * slant2;
2853  aFiv = 0.5 * fDTheta * (fRmax * fRmax - fRmin * fRmin);
2854 
2855  phi = UUtils::Random(fSPhi, ePhi);
2856  cosphi = std::cos(phi);
2857  sinphi = std::sin(phi);
2858  costheta = UUtils::Random(cosETheta, cosSTheta);
2859  sintheta = std::sqrt(1. - UUtils::sqr(costheta));
2860 
2861  if (fFullPhiSphere)
2862  {
2863  aFiv = 0;
2864  }
2865  if (fSTheta == 0)
2866  {
2867  aThr = 0;
2868  }
2869  if (eTheta == UUtils::kPi)
2870  {
2871  aFou = 0;
2872  }
2873  if (fSTheta == UUtils::kPi / 2)
2874  {
2875  aThr = UUtils::kPi * (fRmax * fRmax - fRmin * fRmin);
2876  }
2877  if (eTheta == UUtils::kPi / 2)
2878  {
2879  aFou = UUtils::kPi * (fRmax * fRmax - fRmin * fRmin);
2880  }
2881 
2882  chose = UUtils::Random(0., aOne + aTwo + aThr + aFou + 2.*aFiv);
2883  if ((chose >= 0.) && (chose < aOne))
2884  {
2885  return UVector3(fRmax * sintheta * cosphi,
2886  fRmax * sintheta * sinphi, fRmax * costheta);
2887  }
2888  else if ((chose >= aOne) && (chose < aOne + aTwo))
2889  {
2890  return UVector3(fRmin * sintheta * cosphi,
2891  fRmin * sintheta * sinphi, fRmin * costheta);
2892  }
2893  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr))
2894  {
2895  if (fSTheta != UUtils::kPi / 2)
2896  {
2897  zRand = UUtils::Random(fRmin * cosSTheta, fRmax * cosSTheta);
2898  return UVector3(tanSTheta * zRand * cosphi,
2899  tanSTheta * zRand * sinphi, zRand);
2900  }
2901  else
2902  {
2903  return UVector3(rRand * cosphi, rRand * sinphi, 0.);
2904  }
2905  }
2906  else if ((chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + aThr + aFou))
2907  {
2908  if (eTheta != UUtils::kPi / 2)
2909  {
2910  zRand = UUtils::Random(fRmin * cosETheta, fRmax * cosETheta);
2911  return UVector3(tanETheta * zRand * cosphi,
2912  tanETheta * zRand * sinphi, zRand);
2913  }
2914  else
2915  {
2916  return UVector3(rRand * cosphi, rRand * sinphi, 0.);
2917  }
2918  }
2919  else if ((chose >= aOne + aTwo + aThr + aFou) && (chose < aOne + aTwo + aThr + aFou + aFiv))
2920  {
2921  return UVector3(rRand * sintheta * cosSPhi,
2922  rRand * sintheta * sinSPhi, rRand * costheta);
2923  }
2924  else
2925  {
2926  return UVector3(rRand * sintheta * cosEPhi,
2927  rRand * sintheta * sinEPhi, rRand * costheta);
2928  }
2929 }
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 USphere::GetRmax ( ) const
inline

Definition at line 459 of file USphere.hh.

References GetOuterRadius().

460 {
461  return GetOuterRadius();
462 }
USphere::GetOuterRadius
double GetOuterRadius() const
Definition: USphere.hh:218
double USphere::GetRmin ( ) const
inline

Definition at line 453 of file USphere.hh.

References GetInsideRadius().

454 {
455  return GetInsideRadius();
456 }
USphere::GetInsideRadius
double GetInsideRadius() const
Definition: USphere.hh:206
double USphere::GetSPhi ( ) const
inline

Definition at line 465 of file USphere.hh.

References GetStartPhiAngle().

466 {
467  return GetStartPhiAngle();
468 }
USphere::GetStartPhiAngle
double GetStartPhiAngle() const
Definition: USphere.hh:224
double USphere::GetStartPhiAngle ( ) const
inline

Definition at line 224 of file USphere.hh.

Referenced by GetParametersList(), GetSPhi(), and G4USphere::GetStartPhiAngle().

225 {
226  return fSPhi;
227 }
double USphere::GetStartThetaAngle ( ) const
inline

Definition at line 236 of file USphere.hh.

Referenced by GetParametersList(), G4USphere::GetStartThetaAngle(), and GetSTheta().

237 {
238  return fSTheta;
239 }
double USphere::GetSTheta ( ) const
inline

Definition at line 477 of file USphere.hh.

References GetStartThetaAngle().

478 {
479  return GetStartThetaAngle();
480 }
USphere::GetStartThetaAngle
double GetStartThetaAngle() const
Definition: USphere.hh:236
VUSolid::EnumInside USphere::Inside ( const UVector3 &  p) const
virtual

Implements VUSolid.

Definition at line 164 of file USphere.cc.

References VUSolid::eInside, VUSolid::eOutside, VUSolid::eSurface, G4INCL::Math::max(), UVector3::x, UVector3::y, and UVector3::z.

Referenced by SafetyFromInside().

165 {
166  double rho, rho2, rad2, tolRMin, tolRMax;
167  double pPhi, pTheta;
168  VUSolid::EnumInside in = eOutside;
169  static const double halfAngTolerance = kAngTolerance * 0.5;
170  const double halfTolerance = kTolerance * 0.5;
171  const double halfRminTolerance = fRminTolerance * 0.5;
172  const double rMaxMinus = fRmax - halfTolerance;
173  const double rMinPlus = (fRmin > 0) ? fRmin + halfRminTolerance : 0;
174 
175  rho2 = p.x * p.x + p.y * p.y;
176  rad2 = rho2 + p.z * p.z;
177 
178  // Check radial surfaces. Sets 'in'
179 
180  tolRMin = rMinPlus;
181  tolRMax = rMaxMinus;
182 
183  if ((rad2 <= rMaxMinus * rMaxMinus) && (rad2 >= rMinPlus * rMinPlus))
184  {
185  in = eInside;
186  }
187  else
188  {
189  tolRMax = fRmax + halfTolerance; // outside case
190  tolRMin = std::max(fRmin - halfRminTolerance, 0.); // outside case
191  if ((rad2 <= tolRMax * tolRMax) && (rad2 >= tolRMin * tolRMin))
192  {
193  in = eSurface;
194  }
195  else
196  {
197  return in = eOutside;
198  }
199  }
200 
201  // Phi boundaries : Do not check if it has no phi boundary!
202 
203  if (!fFullPhiSphere && rho2) // [fDPhi < 2*UUtils::kPi] and [p.x or p.y]
204  {
205  pPhi = std::atan2(p.y, p.x);
206 
207  if (pPhi < fSPhi - halfAngTolerance)
208  {
209  pPhi += 2 * UUtils::kPi;
210  }
211  else if (pPhi > ePhi + halfAngTolerance)
212  {
213  pPhi -= 2 * UUtils::kPi;
214  }
215 
216  if ((pPhi < fSPhi - halfAngTolerance)
217  || (pPhi > ePhi + halfAngTolerance))
218  {
219  return in = eOutside;
220  }
221 
222  else if (in == eInside) // else it's eSurface anyway already
223  {
224  if ((pPhi < fSPhi + halfAngTolerance)
225  || (pPhi > ePhi - halfAngTolerance))
226  {
227  in = eSurface;
228  }
229  }
230  }
231 
232  // Theta bondaries
233 
234  if ((rho2 || p.z) && (!fFullThetaSphere))
235  {
236  rho = std::sqrt(rho2);
237  pTheta = std::atan2(rho, p.z);
238 
239  if (in == eInside)
240  {
241  if ((pTheta < fSTheta + halfAngTolerance)
242  || (pTheta > eTheta - halfAngTolerance))
243  {
244  if ((pTheta >= fSTheta - halfAngTolerance)
245  && (pTheta <= eTheta + halfAngTolerance))
246  {
247  in = eSurface;
248  }
249  else
250  {
251  in = eOutside;
252  }
253  }
254  }
255  else
256  {
257  if ((pTheta < fSTheta - halfAngTolerance)
258  || (pTheta > eTheta + halfAngTolerance))
259  {
260  in = eOutside;
261  }
262  }
263  }
264  return in;
265 }
UVector3::x
double x
Definition: UVector3.hh:136
VUSolid::EnumInside
EnumInside
Definition: VUSolid.hh:23
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
UVector3::z
double z
Definition: UVector3.hh:138
UVector3::y
double y
Definition: UVector3.hh:137
bool USphere::Normal ( const UVector3 &  p,
UVector3 &  n 
) const
virtual

Implements VUSolid.

Definition at line 273 of file USphere.cc.

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

274 {
275  int noSurfaces = 0;
276  double rho, rho2, radius, pTheta, pPhi = 0.;
277  double distRMin = UUtils::Infinity();
278  double distSPhi = UUtils::Infinity(), distEPhi = UUtils::Infinity();
279  double distSTheta = UUtils::Infinity(), distETheta = UUtils::Infinity();
280  UVector3 nR, nPs, nPe, nTs, nTe, nZ(0., 0., 1.);
281  UVector3 norm, sumnorm(0., 0., 0.);
282 
283  static const double halfCarTolerance = 0.5 * VUSolid::Tolerance();
284  static const double halfAngTolerance = 0.5 * kAngTolerance;
285 
286  rho2 = p.x * p.x + p.y * p.y;
287  radius = std::sqrt(rho2 + p.z * p.z);
288  rho = std::sqrt(rho2);
289 
290  double distRMax = std::fabs(radius - fRmax);
291  if (fRmin) distRMin = std::fabs(radius - fRmin);
292 
293  if (rho && !fFullSphere)
294  {
295  pPhi = std::atan2(p.y, p.x);
296 
297  if (pPhi < fSPhi - halfAngTolerance)
298  {
299  pPhi += 2 * UUtils::kPi;
300  }
301  else if (pPhi > ePhi + halfAngTolerance)
302  {
303  pPhi -= 2 * UUtils::kPi;
304  }
305  }
306  if (!fFullPhiSphere)
307  {
308  if (rho)
309  {
310  distSPhi = std::fabs(pPhi - fSPhi);
311  distEPhi = std::fabs(pPhi - ePhi);
312  }
313  else if (!fRmin)
314  {
315  distSPhi = 0.;
316  distEPhi = 0.;
317  }
318  nPs = UVector3(sinSPhi, -cosSPhi, 0);
319  nPe = UVector3(-sinEPhi, cosEPhi, 0);
320  }
321  if (!fFullThetaSphere)
322  {
323  if (rho)
324  {
325  pTheta = std::atan2(rho, p.z);
326  distSTheta = std::fabs(pTheta - fSTheta);
327  distETheta = std::fabs(pTheta - eTheta);
328 
329  nTs = UVector3(-cosSTheta * p.x / rho,
330  -cosSTheta * p.y / rho,
331  sinSTheta);
332 
333  nTe = UVector3(cosETheta * p.x / rho,
334  cosETheta * p.y / rho,
335  -sinETheta);
336  }
337  else if (!fRmin)
338  {
339  if (fSTheta)
340  {
341  distSTheta = 0.;
342  nTs = UVector3(0., 0., -1.);
343  }
344  if (eTheta < UUtils::kPi)
345  {
346  distETheta = 0.;
347  nTe = UVector3(0., 0., 1.);
348  }
349  }
350  }
351  if (radius)
352  {
353  nR = UVector3(p.x / radius, p.y / radius, p.z / radius);
354  }
355 
356  if (distRMax <= halfCarTolerance)
357  {
358  noSurfaces ++;
359  sumnorm += nR;
360  }
361  if (fRmin && (distRMin <= halfCarTolerance))
362  {
363  noSurfaces ++;
364  sumnorm -= nR;
365  }
366  if (!fFullPhiSphere)
367  {
368  if (distSPhi <= halfAngTolerance)
369  {
370  noSurfaces ++;
371  sumnorm += nPs;
372  }
373  if (distEPhi <= halfAngTolerance)
374  {
375  noSurfaces ++;
376  sumnorm += nPe;
377  }
378  }
379  if (!fFullThetaSphere)
380  {
381  if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
382  {
383  noSurfaces ++;
384  if ((radius <= halfCarTolerance) && fFullPhiSphere)
385  {
386  sumnorm += nZ;
387  }
388  else
389  {
390  sumnorm += nTs;
391  }
392  }
393  if ((distETheta <= halfAngTolerance) && (eTheta < UUtils::kPi))
394  {
395  noSurfaces ++;
396  if ((radius <= halfCarTolerance) && fFullPhiSphere)
397  {
398  sumnorm -= nZ;
399  }
400  else
401  {
402  sumnorm += nTe;
403  }
404  if (sumnorm.z == 0.)
405  {
406  sumnorm += nZ;
407  }
408  }
409  }
410  if (noSurfaces == 0)
411  {
412 #ifdef UDEBUG
413  UUtils::Exception("USphere::SurfaceNormal(p)", "GeomSolids1002",
414  Warning, 1, "Point p is not on surface !?");
415 #endif
416  norm = ApproxSurfaceNormal(p);
417  }
418  else if (noSurfaces == 1)
419  {
420  norm = sumnorm;
421  }
422  else
423  {
424  norm = sumnorm.Unit();
425  }
426  n = norm;
427  return (noSurfaces > 0);
428 }
Warning
Definition: UUtils.hh:36
UUtils::Infinity
double Infinity()
Definition: UUtils.hh:177
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
USphere & USphere::operator= ( const USphere &  rhs)

Definition at line 103 of file USphere.cc.

104 {
105  // Check assignment to self
106  //
107  if (this == &rhs)
108  {
109  return *this;
110  }
111 
112  // Copy base class data
113  //
114  VUSolid::operator=(rhs);
115 
116  // Copy data
117  //
118  fCubicVolume = rhs.fCubicVolume;
119  fSurfaceArea = rhs.fSurfaceArea;
120  fRminTolerance = rhs.fRminTolerance;
121  kTolerance = rhs.kTolerance;
122  kAngTolerance = rhs.kAngTolerance;
123  kRadTolerance = rhs.kRadTolerance;
124  fEpsilon = rhs.fEpsilon;
125  fRmin = rhs.fRmin;
126  fRmax = rhs.fRmax;
127  fSPhi = rhs.fSPhi;
128  fDPhi = rhs.fDPhi;
129  fSTheta = rhs.fSTheta;
130  fDTheta = rhs.fDTheta;
131  sinCPhi = rhs.sinCPhi;
132  cosCPhi = rhs.cosCPhi;
133  cosHDPhiOT = rhs.cosHDPhiOT;
134  cosHDPhiIT = rhs.cosHDPhiIT;
135  sinSPhi = rhs.sinSPhi;
136  cosSPhi = rhs.cosSPhi;
137  sinEPhi = rhs.sinEPhi;
138  cosEPhi = rhs.cosEPhi;
139  hDPhi = rhs.hDPhi;
140  cPhi = rhs.cPhi;
141  ePhi = rhs.ePhi;
142  sinSTheta = rhs.sinSTheta;
143  cosSTheta = rhs.cosSTheta;
144  sinETheta = rhs.sinETheta;
145  cosETheta = rhs.cosETheta;
146  tanSTheta = rhs.tanSTheta;
147  tanSTheta2 = rhs.tanSTheta2;
148  tanETheta = rhs.tanETheta;
149  tanETheta2 = rhs.tanETheta2;
150  eTheta = rhs.eTheta;
151  fFullPhiSphere = rhs.fFullPhiSphere;
152  fFullThetaSphere = rhs.fFullThetaSphere;
153  fFullSphere = rhs.fFullSphere;
154 
155  return *this;
156 }
double USphere::SafetyFromInside ( const UVector3 &  p,
bool  aAccurate = false 
) const
virtual

Implements VUSolid.

Definition at line 2557 of file USphere.cc.

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

2558 {
2559  double safe = 0.0, safeRMin, safeRMax, safePhi, safeTheta;
2560  double rho2, rds, rho;
2561  double pTheta, dTheta1, dTheta2;
2562  rho2 = p.x * p.x + p.y * p.y;
2563  rds = std::sqrt(rho2 + p.z * p.z);
2564  rho = std::sqrt(rho2);
2565 
2566 #ifdef UDEBUG
2567  if (Inside(p) == eOutside)
2568  {
2569  int old_prc = cout.precision(16);
2570  cout << std::endl;
2571 
2572  cout << "Position:" << std::endl << std::endl;
2573  cout << "p.x = " << p.x << " mm" << std::endl;
2574  cout << "p.y = " << p.y << " mm" << std::endl;
2575  cout << "p.z = " << p.z << " mm" << std::endl << std::endl;
2576  cout.precision(old_prc);
2577  UUtils::Exception("USphere::DistanceToOut(p)",
2578  "GeomSolids1002", Warning, 1, "Point p is outside !?");
2579  }
2580 #endif
2581 
2582  //
2583  // Distance to r shells
2584  //
2585  if (fRmin)
2586  {
2587  safeRMin = rds - fRmin;
2588  safeRMax = fRmax - rds;
2589  if (safeRMin < safeRMax)
2590  {
2591  safe = safeRMin;
2592  }
2593  else
2594  {
2595  safe = safeRMax;
2596  }
2597  }
2598  else
2599  {
2600  safe = fRmax - rds;
2601  }
2602 
2603  //
2604  // Distance to phi extent
2605  //
2606  if (!fFullPhiSphere && rho)
2607  {
2608  if ((p.y * cosCPhi - p.x * sinCPhi) <= 0)
2609  {
2610  safePhi = -(p.x * sinSPhi - p.y * cosSPhi);
2611  }
2612  else
2613  {
2614  safePhi = (p.x * sinEPhi - p.y * cosEPhi);
2615  }
2616  if (safePhi < safe)
2617  {
2618  safe = safePhi;
2619  }
2620  }
2621 
2622  //
2623  // Distance to Theta extent
2624  //
2625  if (rds)
2626  {
2627  pTheta = std::acos(p.z / rds);
2628  if (pTheta < 0)
2629  {
2630  pTheta += UUtils::kPi;
2631  }
2632  dTheta1 = pTheta - fSTheta;
2633  dTheta2 = eTheta - pTheta;
2634  if (dTheta1 < dTheta2)
2635  {
2636  safeTheta = rds * std::sin(dTheta1);
2637  }
2638  else
2639  {
2640  safeTheta = rds * std::sin(dTheta2);
2641  }
2642  if (safe > safeTheta)
2643  {
2644  safe = safeTheta;
2645  }
2646  }
2647 
2648  if (safe < 0)
2649  {
2650  safe = 0;
2651  }
2652  return safe;
2653 }
USphere::Inside
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: USphere.cc:164
Warning
Definition: UUtils.hh:36
UVector3::x
double x
Definition: UVector3.hh:136
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 USphere::SafetyFromOutside ( const UVector3 &  p,
bool  aAccurate = false 
) const
virtual

Implements VUSolid.

Definition at line 1546 of file USphere.cc.

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

1547 {
1548  double safe = 0.0, safeRMin, safeRMax, safePhi, safeTheta;
1549  double rho2, rds, rho;
1550  double cosPsi;
1551  double pTheta, dTheta1, dTheta2;
1552  rho2 = p.x * p.x + p.y * p.y;
1553  rds = std::sqrt(rho2 + p.z * p.z);
1554  rho = std::sqrt(rho2);
1555 
1556  //
1557  // Distance to r shells
1558  //
1559  if (fRmin)
1560  {
1561  safeRMin = fRmin - rds;
1562  safeRMax = rds - fRmax;
1563  if (safeRMin > safeRMax)
1564  {
1565  safe = safeRMin;
1566  }
1567  else
1568  {
1569  safe = safeRMax;
1570  }
1571  }
1572  else
1573  {
1574  safe = rds - fRmax;
1575  }
1576 
1577  //
1578  // Distance to phi extent
1579  //
1580  if (!fFullPhiSphere && rho)
1581  {
1582  // Psi=angle from central phi to point
1583  //
1584  cosPsi = (p.x * cosCPhi + p.y * sinCPhi) / rho;
1585  if (cosPsi < std::cos(hDPhi))
1586  {
1587  // Point lies outside phi range
1588  //
1589  if ((p.y * cosCPhi - p.x * sinCPhi) <= 0)
1590  {
1591  safePhi = std::fabs(p.x * sinSPhi - p.y * cosSPhi);
1592  }
1593  else
1594  {
1595  safePhi = std::fabs(p.x * sinEPhi - p.y * cosEPhi);
1596  }
1597  if (safePhi > safe)
1598  {
1599  safe = safePhi;
1600  }
1601  }
1602  }
1603  //
1604  // Distance to Theta extent
1605  //
1606  if ((rds != 0.0) && (!fFullThetaSphere))
1607  {
1608  pTheta = std::acos(p.z / rds);
1609  if (pTheta < 0)
1610  {
1611  pTheta += UUtils::kPi;
1612  }
1613  dTheta1 = fSTheta - pTheta;
1614  dTheta2 = pTheta - eTheta;
1615  if (dTheta1 > dTheta2)
1616  {
1617  if (dTheta1 >= 0) // WHY ???????????
1618  {
1619  safeTheta = rds * std::sin(dTheta1);
1620  if (safe <= safeTheta)
1621  {
1622  safe = safeTheta;
1623  }
1624  }
1625  }
1626  else
1627  {
1628  if (dTheta2 >= 0)
1629  {
1630  safeTheta = rds * std::sin(dTheta2);
1631  if (safe <= safeTheta)
1632  {
1633  safe = safeTheta;
1634  }
1635  }
1636  }
1637  }
1638 
1639  if (safe < 0)
1640  {
1641  safe = 0;
1642  }
1643  return safe;
1644 }
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 USphere::SetDeltaPhiAngle ( double  newDphi)
inline

Definition at line 430 of file USphere.hh.

Referenced by G4USphere::SetDeltaPhiAngle().

431 {
432  CheckPhiAngles(fSPhi, newDPhi);
433  Initialize();
434 }
void USphere::SetDeltaThetaAngle ( double  newDTheta)
inline

Definition at line 444 of file USphere.hh.

Referenced by G4USphere::SetDeltaThetaAngle().

445 {
446  CheckThetaAngles(fSTheta, newDTheta);
447  Initialize();
448 }
void USphere::SetInnerRadius ( double  newRMin)
inline

Definition at line 401 of file USphere.hh.

References SetInsideRadius().

Referenced by G4USphere::SetInnerRadius().

402 {
403  SetInsideRadius(newRmin);
404 }
USphere::SetInsideRadius
void SetInsideRadius(double newRmin)
Definition: USphere.hh:393
void USphere::SetInsideRadius ( double  newRmin)
inline

Definition at line 393 of file USphere.hh.

References G4INCL::Math::max().

Referenced by SetInnerRadius().

394 {
395  fRmin = newRmin;
396  fRminTolerance = (fRmin) ? std::max(kRadTolerance, fEpsilon * fRmin) : 0;
397  Initialize();
398 }
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
void USphere::SetOuterRadius ( double  newRmax)
inline

Definition at line 407 of file USphere.hh.

References G4INCL::Math::max().

Referenced by G4USphere::SetOuterRadius().

408 {
409  fRmax = newRmax;
410  kTolerance = std::max(kRadTolerance, fEpsilon * fRmax);
411  Initialize();
412 }
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
void USphere::SetStartPhiAngle ( double  newSphi,
bool  trig = true 
)
inline

Definition at line 415 of file USphere.hh.

Referenced by G4USphere::SetStartPhiAngle().

416 {
417  // Flag 'compute' can be used to explicitely avoid recomputation of
418  // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
419 
420  CheckSPhiAngle(newSPhi);
421  fFullPhiSphere = false;
422  if (compute)
423  {
424  InitializePhiTrigonometry();
425  }
426  Initialize();
427 }
void USphere::SetStartThetaAngle ( double  newSTheta)
inline

Definition at line 437 of file USphere.hh.

Referenced by G4USphere::SetStartThetaAngle().

438 {
439  CheckThetaAngles(newSTheta, fDTheta);
440  Initialize();
441 }
std::ostream & USphere::StreamInfo ( std::ostream &  os) const
virtual

Implements VUSolid.

Definition at line 2814 of file USphere.cc.

References VUSolid::GetName().

2815 {
2816  int oldprc = os.precision(16);
2817  os << "-----------------------------------------------------------\n"
2818  << " *** Dump for solid - " << GetName() << " ***\n"
2819  << " ===================================================\n"
2820  << " Solid type: USphere\n"
2821  << " Parameters: \n"
2822  << " inner radius: " << fRmin << " mm \n"
2823  << " outer radius: " << fRmax << " mm \n"
2824  << " starting phi of segment : " << fSPhi / (UUtils::kPi / 180.0) << " degrees \n"
2825  << " delta phi of segment : " << fDPhi / (UUtils::kPi / 180.0) << " degrees \n"
2826  << " starting theta of segment: " << fSTheta / (UUtils::kPi / 180.0) << " degrees \n"
2827  << " delta theta of segment : " << fDTheta / (UUtils::kPi / 180.0) << " degrees \n"
2828  << "-----------------------------------------------------------\n";
2829  os.precision(oldprc);
2830 
2831  return os;
2832 }
VUSolid::GetName
const std::string & GetName() const
Definition: VUSolid.hh:103
double USphere::SurfaceArea ( )
virtual

Implements VUSolid.

Definition at line 2935 of file USphere.cc.

2936 {
2937  if (fSurfaceArea != 0.)
2938  {
2939  ;
2940  }
2941  else
2942  {
2943  double Rsq = fRmax * fRmax;
2944  double rsq = fRmin * fRmin;
2945 
2946  fSurfaceArea = fDPhi * (rsq + Rsq) * (cosSTheta - cosETheta);
2947  if (!fFullPhiSphere)
2948  {
2949  fSurfaceArea = fSurfaceArea + fDTheta * (Rsq - rsq);
2950  }
2951  if (fSTheta > 0)
2952  {
2953  double acos1 = std::acos(std::pow(sinSTheta, 2) * std::cos(fDPhi)
2954  + std::pow(cosSTheta, 2));
2955  if (fDPhi > UUtils::kPi)
2956  {
2957  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * (2 * UUtils::kPi - acos1);
2958  }
2959  else
2960  {
2961  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * acos1;
2962  }
2963  }
2964  if (eTheta < UUtils::kPi)
2965  {
2966  double acos2 = std::acos(std::pow(sinETheta, 2) * std::cos(fDPhi)
2967  + std::pow(cosETheta, 2));
2968  if (fDPhi > UUtils::kPi)
2969  {
2970  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * (2 * UUtils::kPi - acos2);
2971  }
2972  else
2973  {
2974  fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * acos2;
2975  }
2976  }
2977  }
2978  return fSurfaceArea;
2979 }
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:57:32 for Geant4.10 by   doxygen 1.8.7