#include <UTubs.hh>

Inheritance diagram for UTubs:

Public Member Functions
	UTubs (const std::string &pName, double pRMin, double pRMax, double pDz, double pSPhi, double pDPhi)

virtual	~UTubs ()

double	GetInnerRadius () const

double	GetOuterRadius () const

double	GetZHalfLength () const

double	GetStartPhiAngle () const

double	GetDeltaPhiAngle () const

void	SetInnerRadius (double newRMin)

void	SetOuterRadius (double newRMax)

void	SetZHalfLength (double newDz)

void	SetStartPhiAngle (double newSPhi, bool trig=true)

void	SetDeltaPhiAngle (double newDPhi)

double	Capacity ()

double	SurfaceArea ()

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

bool	Normal (const UVector3 &p, UVector3 &normal) const

double	DistanceToIn (const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const

double	SafetyFromInside (const UVector3 &p, bool precise=false) const

double	DistanceToOut (const UVector3 &p, const UVector3 &v, UVector3 &n, bool &validNorm, double aPstep=UUtils::kInfinity) const

double	SafetyFromOutside (const UVector3 &p, bool precise=false) 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)

	UTubs ()

	UTubs (const UTubs &rhs)

UTubs &	operator= (const UTubs &rhs)

double	GetRMin () const

double	GetRMax () 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

Protected Types
enum	ESide { kNull, kRMin, kRMax, kSPhi, kEPhi, kPZ, kMZ }

enum	ENorm { kNRMin, kNRMax, kNSPhi, kNEPhi, kNZ }

Protected Member Functions
void	Initialize ()

void	CheckSPhiAngle (double sPhi)

void	CheckDPhiAngle (double dPhi)

void	CheckPhiAngles (double sPhi, double dPhi)

void	InitializeTrigonometry ()

virtual UVector3	ApproxSurfaceNormal (const UVector3 &p) const

Protected Attributes
double	fCubicVolume

double	fSurfaceArea

double	kRadTolerance

double	kAngTolerance

double	fRMin

double	fRMax

double	fDz

double	fSPhi

double	fDPhi

double	fSinCPhi

double	fCosCPhi

double	fCosHDPhiOT

double	fCosHDPhiIT

double	fSinSPhi

double	fCosSPhi

double	fSinEPhi

double	fCosEPhi

double	fSinSPhiDPhi

double	fCosSPhiDPhi

bool	fPhiFullTube

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 47 of file UTubs.hh.

Member Enumeration Documentation

enum UTubs::ENorm

protected

Enumerator
kNRMin
kNRMax
kNSPhi
kNEPhi
kNZ

Definition at line 164 of file UTubs.hh.

164 {kNRMin, kNRMax, kNSPhi, kNEPhi, kNZ};

UTubs::kNRMax

Definition: UTubs.hh:164

UTubs::kNRMin

Definition: UTubs.hh:164

UTubs::kNEPhi

Definition: UTubs.hh:164

UTubs::kNZ

Definition: UTubs.hh:164

UTubs::kNSPhi

Definition: UTubs.hh:164

enum UTubs::ESide

protected

Enumerator
kNull
kRMin
kRMax
kSPhi
kEPhi
kPZ
kMZ

Definition at line 160 of file UTubs.hh.

160 {kNull, kRMin, kRMax, kSPhi, kEPhi, kPZ, kMZ};

UTubs::kRMax

Definition: UTubs.hh:160

UTubs::kEPhi

Definition: UTubs.hh:160

UTubs::kMZ

Definition: UTubs.hh:160

UTubs::kSPhi

Definition: UTubs.hh:160

UTubs::kPZ

Definition: UTubs.hh:160

UTubs::kNull

Definition: UTubs.hh:160

UTubs::kRMin

Definition: UTubs.hh:160

Constructor & Destructor Documentation

UTubs::UTubs	(	const std::string &	pName,
		double	pRMin,
		double	pRMax,
		double	pDz,
		double	pSPhi,
		double	pDPhi
	)

Definition at line 30 of file UTubs.cc.

References CheckPhiAngles(), UUtils::Exception(), FatalErrorInArguments, VUSolid::faTolerance, fCubicVolume, VUSolid::frTolerance, fSurfaceArea, VUSolid::GetName(), kAngTolerance, and kRadTolerance.

   : VUSolid(pName.c_str()), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
 {
   fCubicVolume=0.;
   fSurfaceArea=0.;
   kRadTolerance = frTolerance;
   kAngTolerance = faTolerance;
 
   if (pDz <= 0) // Check z-len
   {
     std::ostringstream message;
     message << "Negative Z half-length (" << pDz << ") in solid: " << GetName();
     UUtils::Exception("UTubs::UTubs()", "GeomSolids0002", FatalErrorInArguments, 1, message.str().c_str());
   }
   if ((pRMin >= pRMax) || (pRMin < 0))   // Check radii
   {
     std::ostringstream message;
     message << "Invalid values for radii in solid: " << GetName()
             << std::endl
             << "                pRMin = " << pRMin << ", pRMax = " << pRMax;
     UUtils::Exception("UTubs::UTubs()", "GeomSolids0002", FatalErrorInArguments, 1, message.str().c_str());
   }
 
   // Check angles
   //
   CheckPhiAngles(pSPhi, pDPhi);
 }

UTubs::~UTubs ( )

virtual

Definition at line 81 of file UTubs.cc.

82 {

83 }

UTubs::UTubs ( )

Definition at line 66 of file UTubs.cc.

Referenced by Clone().

   : VUSolid(""), fCubicVolume(0.), fSurfaceArea(0.),
     kRadTolerance(0.), kAngTolerance(0.),
     fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
     fSinCPhi(0.), fCosCPhi(0.), fCosHDPhiOT(0.), fCosHDPhiIT(0.),
     fSinSPhi(0.), fCosSPhi(0.), fSinEPhi(0.), fCosEPhi(0.),
     fSinSPhiDPhi(0.), fCosSPhiDPhi(0.),
     fPhiFullTube(false)
 {
 }

UTubs::UTubs ( const UTubs & rhs )

Definition at line 89 of file UTubs.cc.

   : VUSolid(rhs),
     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
     kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
     fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
     fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
     fSinCPhi(rhs.fSinCPhi), fCosCPhi(rhs.fSinCPhi),
     fCosHDPhiOT(rhs.fCosHDPhiOT), fCosHDPhiIT(rhs.fCosHDPhiOT),
     fSinSPhi(rhs.fSinSPhi), fCosSPhi(rhs.fCosSPhi),
     fSinEPhi(rhs.fSinEPhi), fCosEPhi(rhs.fCosEPhi),
     fSinSPhiDPhi(rhs.fSinSPhiDPhi), fCosSPhiDPhi(rhs.fCosSPhiDPhi),
     fPhiFullTube(rhs.fPhiFullTube)
 {
 }

Member Function Documentation

UVector3 UTubs::ApproxSurfaceNormal ( const UVector3 & p ) const

protectedvirtual

Definition at line 473 of file UTubs.cc.

References UUtils::Exception(), fCosSPhi, fCosSPhiDPhi, fDPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinSPhi, fSinSPhiDPhi, fSPhi, kNEPhi, kNRMax, kNRMin, kNSPhi, kNZ, Warning, UVector3::x, UVector3::y, and UVector3::z.

Referenced by Normal().

 {
   ENorm side;
   UVector3 norm;
   double rho, phi;
   double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin;
 
   rho = std::sqrt(p.x * p.x + p.y * p.y);
 
   distRMin = std::fabs(rho - fRMin);
   distRMax = std::fabs(rho - fRMax);
   distZ   = std::fabs(std::fabs(p.z) - fDz);
 
   if (distRMin < distRMax) // First minimum
   {
     if (distZ < distRMin)
     {
       distMin = distZ;
       side    = kNZ;
     }
     else
     {
       distMin = distRMin;
       side    = kNRMin  ;
     }
   }
   else
   {
     if (distZ < distRMax)
     {
       distMin = distZ;
       side    = kNZ ;
     }
     else
     {
       distMin = distRMax;
       side    = kNRMax  ;
     }
   }
   if (!fPhiFullTube &&  rho)  // Protected against (0,0,z)
   {
     phi = std::atan2(p.y, p.x);
 
     if (phi < 0)
     {
       phi += 2 * UUtils::kPi;
     }
 
     if (fSPhi < 0)
     {
       distSPhi = std::fabs(phi - (fSPhi + 2 * UUtils::kPi)) * rho;
     }
     else
     {
       distSPhi = std::fabs(phi - fSPhi) * rho;
     }
     distEPhi = std::fabs(phi - fSPhi - fDPhi) * rho;
 
     if (distSPhi < distEPhi) // Find new minimum
     {
       if (distSPhi < distMin)
       {
         side = kNSPhi;
       }
     }
     else
     {
       if (distEPhi < distMin)
       {
         side = kNEPhi;
       }
     }
   }
   switch (side)
   {
     case kNRMin : // Inner radius
     {
       norm = UVector3(-p.x / rho, -p.y / rho, 0);
       break;
     }
     case kNRMax : // Outer radius
     {
       norm = UVector3(p.x / rho, p.y / rho, 0);
       break;
     }
     case kNZ :    // + or - dz
     {
       if (p.z > 0)
       {
         norm = UVector3(0, 0, 1);
       }
       else
       {
         norm = UVector3(0, 0, -1);
       }
       break;
     }
     case kNSPhi:
     {
       norm = UVector3(fSinSPhi, -fCosSPhi, 0);
       break;
     }
     case kNEPhi:
     {
       norm = UVector3(-fSinSPhiDPhi, fCosSPhiDPhi, 0);
       break;
     }
     default:      // Should never reach this case ...
     {
       // DumpInfo();
       UUtils::Exception("UTubs::ApproxSurfaceNormal()",
                         "GeomSolids1002", Warning, 1,
                         "Undefined side for valid surface normal to solid.");
       break;
     }
   }
   return norm;
 }

double UTubs::Capacity ( )

inlinevirtual

Implements VUSolid.

void UTubs::CheckDPhiAngle ( double dPhi )

inlineprotected

void UTubs::CheckPhiAngles	(	double	sPhi,
		double	dPhi
	)

inlineprotected

Referenced by UTubs().

void UTubs::CheckSPhiAngle ( double sPhi )

inlineprotected

VUSolid * UTubs::Clone ( ) const

virtual

Implements VUSolid.

Definition at line 1647 of file UTubs.cc.

References UTubs().

 {
   return new UTubs(*this);
 }

virtual void UTubs::ComputeBBox	(	UBBox *	,
		bool
	)

inlinevirtual

Implements VUSolid.

Definition at line 106 of file UTubs.hh.

106 {}

double UTubs::DistanceToIn	(	const UVector3 &	p,
		const UVector3 &	v,
		double	aPstep = `UUtils::kInfinity`
	)		const

virtual

Implements VUSolid.

Definition at line 614 of file UTubs.cc.

References test::b, test::c, fCosCPhi, fCosEPhi, fCosHDPhiIT, fCosSPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinCPhi, fSinEPhi, fSinSPhi, kRadTolerance, plottest35::t1, VUSolid::Tolerance(), UVector3::x, UVector3::y, and UVector3::z.

Referenced by UEnclosingCylinder::DistanceTo().

 {
   double snxt = UUtils::kInfinity;      // snxt = default return value
   double tolORMin2, tolIRMax2;  // 'generous' radii squared
   double tolORMax2, tolIRMin2, tolODz, tolIDz;
   const double dRmax = 100.*fRMax;
 
   static const double halfCarTolerance = 0.5 * VUSolid::Tolerance();
   static const double halfRadTolerance = 0.5 * kRadTolerance;
 
   // Intersection point variables
   //
   double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp;
   double t1, t2, t3, b, c, d;    // Quadratic solver variables
 
   // Calculate tolerant rmin and rmax
 
   if (fRMin > kRadTolerance)
   {
     tolORMin2 = (fRMin - halfRadTolerance) * (fRMin - halfRadTolerance);
     tolIRMin2 = (fRMin + halfRadTolerance) * (fRMin + halfRadTolerance);
   }
   else
   {
     tolORMin2 = 0.0;
     tolIRMin2 = 0.0;
   }
   tolORMax2 = (fRMax + halfRadTolerance) * (fRMax + halfRadTolerance);
   tolIRMax2 = (fRMax - halfRadTolerance) * (fRMax - halfRadTolerance);
 
   // Intersection with Z surfaces
 
   tolIDz = fDz - halfCarTolerance;
   tolODz = fDz + halfCarTolerance;
 
   if (std::fabs(p.z) >= tolIDz)
   {
     if (p.z * v.z < 0)   // at +Z going in -Z or visa versa
     {
       sd = (std::fabs(p.z) - fDz) / std::fabs(v.z); // Z intersect distance
 
       if (sd < 0.0)
       {
         sd = 0.0;
       }
 
       xi   = p.x + sd * v.x;              // Intersection coords
       yi   = p.y + sd * v.y;
       rho2 = xi * xi + yi * yi;
 
       // Check validity of intersection
 
       if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
       {
         if (!fPhiFullTube && rho2)
         {
           // Psi = angle made with central (average) phi of shape
           //
           inum   = xi * fCosCPhi + yi * fSinCPhi;
           iden   = std::sqrt(rho2);
           cosPsi = inum / iden;
           if (cosPsi >= fCosHDPhiIT)
           {
             return sd;
           }
         }
         else
         {
           return sd;
         }
       }
     }
     else
     {
       if (snxt < halfCarTolerance)
       {
         snxt = 0;
       }
       return snxt;  // On/outside extent, and heading away
       // -> cannot intersect
     }
   }
 
   // -> Can not intersect z surfaces
   //
   // Intersection with rmax (possible return) and rmin (must also check phi)
   //
   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
   //
   // Intersects with x^2+y^2=R^2
   //
   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
   //            t1                t2                t3
 
   t1 = 1.0 - v.z * v.z;
   t2 = p.x * v.x + p.y * v.y;
   t3 = p.x * p.x + p.y * p.y;
 
   if (t1 > 0)          // Check not || to z axis
   {
     b = t2 / t1;
     c = t3 - fRMax * fRMax;
     if ((t3 >= tolORMax2) && (t2 < 0)) // This also handles the tangent case
     {
       // Try outer cylinder intersection
       //          c=(t3-fRMax*fRMax)/t1;
 
       c /= t1;
       d = b * b - c;
 
       if (d >= 0) // If real root
       {
         sd = c / (-b + std::sqrt(d));
         if (sd >= 0)  // If 'forwards'
         {
           if (sd > dRmax) // Avoid rounding errors due to precision issues on
           {
             // 64 bits systems. Split long distances and recompute
             double fTerm = sd - std::fmod(sd, dRmax);
             sd = fTerm + DistanceToIn(p + fTerm * v, v);
           }
           // Check z intersection
           //
           zi = p.z + sd * v.z;
           if (std::fabs(zi) <= tolODz)
           {
             // Z ok. Check phi intersection if reqd
             //
             if (fPhiFullTube)
             {
               return sd;
             }
             else
             {
               xi     = p.x + sd * v.x;
               yi     = p.y + sd * v.y;
               cosPsi = (xi * fCosCPhi + yi * fSinCPhi) / fRMax;
               if (cosPsi >= fCosHDPhiIT)
               {
                 return sd;
               }
             }
           } //  end if std::fabs(zi)
         }   //  end if (sd>=0)
       }     //  end if (d>=0)
     }       //  end if (r>=fRMax)
     else
     {
       // Inside outer radius :
       // check not inside, and heading through tubs (-> 0 to in)
 
       if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z) <= tolIDz))
       {
         // Inside both radii, delta r -ve, inside z extent
 
         if (!fPhiFullTube)
         {
           inum   = p.x * fCosCPhi + p.y * fSinCPhi;
           iden   = std::sqrt(t3);
           cosPsi = inum / iden;
           if (cosPsi >= fCosHDPhiIT)
           {
             // In the old version, the small negative tangent for the point
             // on surface was not taken in account, and returning 0.0 ...
             // New version: check the tangent for the point on surface and
             // if no intersection, return UUtils::kInfinity, if intersection instead
             // return sd.
             //
             c = t3 - fRMax * fRMax;
             if (c <= 0.0)
             {
               return 0.0;
             }
             else
             {
               c = c / t1;
               d = b * b - c;
               if (d >= 0.0)
               {
                 snxt = c / (-b + std::sqrt(d)); // using safe solution
                 // for quadratic equation
                 if (snxt < halfCarTolerance)
                 {
                   snxt = 0;
                 }
                 return snxt;
               }
               else
               {
                 return UUtils::kInfinity;
               }
             }
           }
         }
         else
         {
           // In the old version, the small negative tangent for the point
           // on surface was not taken in account, and returning 0.0 ...
           // New version: check the tangent for the point on surface and
           // if no intersection, return UUtils::kInfinity, if intersection instead
           // return sd.
           //
           c = t3 - fRMax * fRMax;
           if (c <= 0.0)
           {
             return 0.0;
           }
           else
           {
             c = c / t1;
             d = b * b - c;
             if (d >= 0.0)
             {
               snxt = c / (-b + std::sqrt(d)); // using safe solution
               // for quadratic equation
               if (snxt < halfCarTolerance)
               {
                 snxt = 0;
               }
               return snxt;
             }
             else
             {
               return UUtils::kInfinity;
             }
           }
         } // end if  (!fPhiFullTube)
       }  // end if   (t3>tolIRMin2)
     }    // end if   (Inside Outer Radius)
     if (fRMin)     // Try inner cylinder intersection
     {
       c = (t3 - fRMin * fRMin) / t1;
       d = b * b - c;
       if (d >= 0.0)    // If real root
       {
         // Always want 2nd root - we are outside and know rmax Hit was bad
         // - If on surface of rmin also need farthest root
 
         sd = (b > 0.) ? c / (-b - std::sqrt(d)) : (-b + std::sqrt(d));
         if (sd >= -halfCarTolerance)  // check forwards
         {
           // Check z intersection
           //
           if (sd < 0.0)
           {
             sd = 0.0;
           }
           if (sd > dRmax) // Avoid rounding errors due to precision issues seen
           {
             // 64 bits systems. Split long distances and recompute
             double fTerm = sd - std::fmod(sd, dRmax);
             sd = fTerm + DistanceToIn(p + fTerm * v, v);
           }
           zi = p.z + sd * v.z;
           if (std::fabs(zi) <= tolODz)
           {
             // Z ok. Check phi
             //
             if (fPhiFullTube)
             {
               return sd;
             }
             else
             {
               xi     = p.x + sd * v.x;
               yi     = p.y + sd * v.y;
               cosPsi = (xi * fCosCPhi + yi * fSinCPhi) / fRMin;
               if (cosPsi >= fCosHDPhiIT)
               {
                 // Good inner radius isect
                 // - but earlier phi isect still possible
 
                 snxt = sd;
               }
             }
           }       //    end if std::fabs(zi)
         }         //    end if (sd>=0)
       }           //    end if (d>=0)
     }             //    end if (fRMin)
   }
 
   // Phi segment intersection
   //
   // o Tolerant of points inside phi planes by up to VUSolid::Tolerance()*0.5
   //
   // o NOTE: Large duplication of code between sphi & ephi checks
   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
   //            intersection check <=0 -> >=0
   //         -> use some form of loop Construct ?
   //
   if (!fPhiFullTube)
   {
     // First phi surface (Starting phi)
     //
     Comp    = v.x * fSinSPhi - v.y * fCosSPhi;
 
     if (Comp < 0)    // Component in outwards normal dirn
     {
       Dist = (p.y * fCosSPhi - p.x * fSinSPhi);
 
       if (Dist < halfCarTolerance)
       {
         sd = Dist / Comp;
 
         if (sd < snxt)
         {
           if (sd < 0)
           {
             sd = 0.0;
           }
           zi = p.z + sd * v.z;
           if (std::fabs(zi) <= tolODz)
           {
             xi   = p.x + sd * v.x;
             yi   = p.y + sd * v.y;
             rho2 = xi * xi + yi * yi;
 
             if (((rho2 >= tolIRMin2) && (rho2 <= tolIRMax2))
                 || ((rho2 > tolORMin2) && (rho2 < tolIRMin2)
                     && (v.y * fCosSPhi - v.x * fSinSPhi > 0)
                     && (v.x * fCosSPhi + v.y * fSinSPhi >= 0))
                 || ((rho2 > tolIRMax2) && (rho2 < tolORMax2)
                     && (v.y * fCosSPhi - v.x * fSinSPhi > 0)
                     && (v.x * fCosSPhi + v.y * fSinSPhi < 0)))
             {
               // z and r intersections good
               // - check intersecting with correct half-plane
               //
               if ((yi * fCosCPhi - xi * fSinCPhi) <= halfCarTolerance)
               {
                 snxt = sd;
               }
             }
           }
         }
       }
     }
 
     // Second phi surface (Ending phi)
 
     Comp    = -(v.x * fSinEPhi - v.y * fCosEPhi);
 
     if (Comp < 0)   // Component in outwards normal dirn
     {
       Dist = -(p.y * fCosEPhi - p.x * fSinEPhi);
 
       if (Dist < halfCarTolerance)
       {
         sd = Dist / Comp;
 
         if (sd < snxt)
         {
           if (sd < 0)
           {
             sd = 0;
           }
           zi = p.z + sd * v.z;
           if (std::fabs(zi) <= tolODz)
           {
             xi   = p.x + sd * v.x;
             yi   = p.y + sd * v.y;
             rho2 = xi * xi + yi * yi;
             if (((rho2 >= tolIRMin2) && (rho2 <= tolIRMax2))
                 || ((rho2 > tolORMin2)  && (rho2 < tolIRMin2)
                     && (v.x * fSinEPhi - v.y * fCosEPhi > 0)
                     && (v.x * fCosEPhi + v.y * fSinEPhi >= 0))
                 || ((rho2 > tolIRMax2) && (rho2 < tolORMax2)
                     && (v.x * fSinEPhi - v.y * fCosEPhi > 0)
                     && (v.x * fCosEPhi + v.y * fSinEPhi < 0)))
             {
               // z and r intersections good
               // - check intersecting with correct half-plane
               //
               if ((yi * fCosCPhi - xi * fSinCPhi) >= 0)
               {
                 snxt = sd;
               }
             }                        //?? >=-halfCarTolerance
           }
         }
       }
     }        // Comp < 0
   }          // !fPhiFullTube
   if (snxt < halfCarTolerance)
   {
     snxt = 0;
   }
   return snxt;
 }

double UTubs::DistanceToOut	(	const UVector3 &	p,
		const UVector3 &	v,
		UVector3 &	n,
		bool &	validNorm,
		double	aPstep = `UUtils::kInfinity`
	)		const

virtual

Implements VUSolid.

Definition at line 1109 of file UTubs.cc.

References test::b, test::c, UUtils::Exception(), fCosCPhi, fCosEPhi, fCosSPhi, fDPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinCPhi, fSinEPhi, fSinSPhi, fSPhi, kAngTolerance, kEPhi, kMZ, kNull, kPZ, kRadTolerance, kRMax, kRMin, kSPhi, plottest35::t1, VUSolid::Tolerance(), Warning, UVector3::x, UVector3::y, and UVector3::z.

 {
   ESide side = kNull , sider = kNull, sidephi = kNull;
   double snxt, srd = UUtils::kInfinity, sphi = UUtils::kInfinity, pdist;
   double deltaR, t1, t2, t3, b, c, d2, roMin2;
 
   static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
   static const double halfAngTolerance = kAngTolerance * 0.5;
 
   // Vars for phi intersection:
 
   double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2;
 
   // Z plane intersection
 
   if (v.z > 0)
   {
     pdist = fDz - p.z;
     if (pdist > halfCarTolerance)
     {
       snxt = pdist / v.z;
       side = kPZ;
     }
     else
     {
       //      if (calcNorm)
       {
         n        = UVector3(0, 0, 1);
         validNorm = true;
       }
       return snxt = 0;
     }
   }
   else if (v.z < 0)
   {
     pdist = fDz + p.z;
 
     if (pdist > halfCarTolerance)
     {
       snxt = -pdist / v.z;
       side = kMZ;
     }
     else
     {
       //      if (calcNorm)
       {
         n        = UVector3(0, 0, -1);
         validNorm = true;
       }
       return snxt = 0.0;
     }
   }
   else
   {
     snxt = UUtils::kInfinity;   // Travel perpendicular to z axis
     side = kNull;
   }
 
   // Radial Intersections
   //
   // Find intersection with cylinders at rmax/rmin
   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
   //
   // Intersects with x^2+y^2=R^2
   //
   // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
   //
   //            t1                t2                    t3
 
   t1   = 1.0 - v.z * v.z;   // since v normalised
   t2   = p.x * v.x + p.y * v.y;
   t3   = p.x * p.x + p.y * p.y;
 
   if (snxt > 10 * (fDz + fRMax))
   {
     roi2 = 2 * fRMax * fRMax;
   }
   else
   {
     roi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3;  // radius^2 on +-fDz
   }
 
   if (t1 > 0)   // Check not parallel
   {
     // Calculate srd, r exit distance
 
     if ((t2 >= 0.0) && (roi2 > fRMax * (fRMax + kRadTolerance)))
     {
       // Delta r not negative => leaving via rmax
 
       deltaR = t3 - fRMax * fRMax;
 
       // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
       // - avoid sqrt for efficiency
 
       if (deltaR < -kRadTolerance * fRMax)
       {
         b    = t2 / t1;
         c    = deltaR / t1;
         d2    = b * b - c;
         if (d2 >= 0)
         {
           srd = c / (-b - std::sqrt(d2));
         }
         else
         {
           srd = 0.;
         }
         sider = kRMax;
       }
       else
       {
         // On tolerant boundary & heading outwards (or perpendicular to)
         // outer radial surface -> leaving immediately
 
         //        if ( calcNorm )
         {
           n        = UVector3(p.x / fRMax, p.y / fRMax, 0);
           validNorm = true;
         }
         return snxt = 0; // Leaving by rmax immediately
       }
     }
     else if (t2 < 0.)   // i.e.  t2 < 0; Possible rmin intersection
     {
       roMin2 = t3 - t2 * t2 / t1; // min ro2 of the plane of movement
 
       if (fRMin && (roMin2 < fRMin * (fRMin - kRadTolerance)))
       {
         deltaR = t3 - fRMin * fRMin;
         b     = t2 / t1;
         c     = deltaR / t1;
         d2     = b * b - c;
 
         if (d2 >= 0)    // Leaving via rmin
         {
           // NOTE: SHould use rho-rmin>kRadTolerance*0.5
           // - avoid sqrt for efficiency
 
           if (deltaR > kRadTolerance * fRMin)
           {
             srd = c / (-b + std::sqrt(d2));
             sider = kRMin;
           }
           else
           {
             //if ( calcNorm )
             {
               validNorm = false;  // Concave side
             }
             return snxt = 0.0;
           }
         }
         else    // No rmin intersect -> must be rmax intersect
         {
           deltaR = t3 - fRMax * fRMax;
           c    = deltaR / t1;
           d2    = b * b - c;
           if (d2 >= 0.)
           {
             srd    = -b + std::sqrt(d2);
             sider = kRMax;
           }
           else // Case: On the border+t2<kRadTolerance
             //       (v is perpendicular to the surface)
           {
             // if (calcNorm)
             {
               n = UVector3(p.x / fRMax, p.y / fRMax, 0);
               validNorm = true;
             }
             return snxt = 0.0;
           }
         }
       }
       else if (roi2 > fRMax * (fRMax + kRadTolerance))
         // No rmin intersect -> must be rmax intersect
       {
         deltaR = t3 - fRMax * fRMax;
         b     = t2 / t1;
         c     = deltaR / t1;
         d2     = b * b - c;
         if (d2 >= 0)
         {
           srd    = -b + std::sqrt(d2);
           sider = kRMax;
         }
         else // Case: On the border+t2<kRadTolerance
           //       (v is perpendicular to the surface)
         {
           //          if (calcNorm)
           {
             n = UVector3(p.x / fRMax, p.y / fRMax, 0);
             validNorm = true;
           }
           return snxt = 0.0;
         }
       }
     }
 
     // Phi Intersection
 
     if (!fPhiFullTube)
     {
       // add angle calculation with correction
       // of the difference in domain of atan2 and Sphi
       //
       vphi = std::atan2(v.y, v.x);
 
       if (vphi < fSPhi - halfAngTolerance)
       {
         vphi += 2 * UUtils::kPi;
       }
       else if (vphi > fSPhi + fDPhi + halfAngTolerance)
       {
         vphi -= 2 * UUtils::kPi;
       }
 
 
       if (p.x || p.y)    // Check if on z axis (rho not needed later)
       {
         // pDist -ve when inside
 
         pDistS = p.x * fSinSPhi - p.y * fCosSPhi;
         pDistE = -p.x * fSinEPhi + p.y * fCosEPhi;
 
         // Comp -ve when in direction of outwards normal
 
         compS  = -fSinSPhi * v.x + fCosSPhi * v.y;
         compE  =  fSinEPhi * v.x - fCosEPhi * v.y;
 
         sidephi = kNull;
 
         if (((fDPhi <= UUtils::kPi) && ((pDistS <= halfCarTolerance)
                                         && (pDistE <= halfCarTolerance)))
             || ((fDPhi >  UUtils::kPi) && !((pDistS > halfCarTolerance)
                                             && (pDistE >  halfCarTolerance))))
         {
           // Inside both phi *full* planes
 
           if (compS < 0)
           {
             sphi = pDistS / compS;
 
             if (sphi >= -halfCarTolerance)
             {
               xi = p.x + sphi * v.x;
               yi = p.y + sphi * v.y;
 
               // Check intersecting with correct half-plane
               // (if not -> no intersect)
               //
               if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
               {
                 sidephi = kSPhi;
                 if (((fSPhi - halfAngTolerance) <= vphi)
                     && ((fSPhi + fDPhi + halfAngTolerance) >= vphi))
                 {
                   sphi = UUtils::kInfinity;
                 }
               }
               else if (yi * fCosCPhi - xi * fSinCPhi >= 0)
               {
                 sphi = UUtils::kInfinity;
               }
               else
               {
                 sidephi = kSPhi;
                 if (pDistS > -halfCarTolerance)
                 {
                   sphi = 0.0; // Leave by sphi immediately
                 }
               }
             }
             else
             {
               sphi = UUtils::kInfinity;
             }
           }
           else
           {
             sphi = UUtils::kInfinity;
           }
 
           if (compE < 0)
           {
             sphi2 = pDistE / compE;
 
             // Only check further if < starting phi intersection
             //
             if ((sphi2 > -halfCarTolerance) && (sphi2 < sphi))
             {
               xi = p.x + sphi2 * v.x;
               yi = p.y + sphi2 * v.y;
 
               if ((std::fabs(xi) <= VUSolid::Tolerance()) && (std::fabs(yi) <= VUSolid::Tolerance()))
               {
                 // Leaving via ending phi
                 //
                 if (!((fSPhi - halfAngTolerance <= vphi)
                       && (fSPhi + fDPhi + halfAngTolerance >= vphi)))
                 {
                   sidephi = kEPhi;
                   if (pDistE <= -halfCarTolerance)
                   {
                     sphi = sphi2;
                   }
                   else
                   {
                     sphi = 0.0;
                   }
                 }
               }
               else    // Check intersecting with correct half-plane
 
                 if ((yi * fCosCPhi - xi * fSinCPhi) >= 0)
                 {
                   // Leaving via ending phi
                   //
                   sidephi = kEPhi;
                   if (pDistE <= -halfCarTolerance)
                   {
                     sphi = sphi2;
                   }
                   else
                   {
                     sphi = 0.0;
                   }
                 }
             }
           }
         }
         else
         {
           sphi = UUtils::kInfinity;
         }
       }
       else
       {
         // On z axis + travel not || to z axis -> if phi of vector direction
         // within phi of shape, Step limited by rmax, else Step =0
 
         if ((fSPhi - halfAngTolerance <= vphi)
             && (vphi <= fSPhi + fDPhi + halfAngTolerance))
         {
           sphi = UUtils::kInfinity;
         }
         else
         {
           sidephi = kSPhi; // arbitrary
           sphi    = 0.0;
         }
       }
       if (sphi < snxt)  // Order intersecttions
       {
         snxt = sphi;
         side = sidephi;
       }
     }
     if (srd < snxt) // Order intersections
     {
       snxt = srd;
       side = sider;
     }
   }
   //  if (calcNorm)
   {
     switch (side)
     {
       case kRMax:
         // Note: returned vector not normalised
         // (divide by fRMax for Unit vector)
         //
         xi = p.x + snxt * v.x;
         yi = p.y + snxt * v.y;
         n = UVector3(xi / fRMax, yi / fRMax, 0);
         validNorm = true;
         break;
 
       case kRMin:
         validNorm = false;  // Rmin is inconvex
         break;
 
       case kSPhi:
         if (fDPhi <= UUtils::kPi)
         {
           n        = UVector3(fSinSPhi, -fCosSPhi, 0);
           validNorm = true;
         }
         else
         {
           validNorm = false;
         }
         break;
 
       case kEPhi:
         if (fDPhi <= UUtils::kPi)
         {
           n = UVector3(-fSinEPhi, fCosEPhi, 0);
           validNorm = true;
         }
         else
         {
           validNorm = false;
         }
         break;
 
       case kPZ:
         n        = UVector3(0, 0, 1);
         validNorm = true;
         break;
 
       case kMZ:
         n        = UVector3(0, 0, -1);
         validNorm = true;
         break;
 
       default:
         cout << std::endl;
         // DumpInfo();
         std::ostringstream message;
         int oldprc = message.precision(16);
         message << "Undefined side for valid surface normal to solid."
                 << std::endl
                 << "Position:"  << std::endl << std::endl
                 << "p.x = "  << p.x << " mm" << std::endl
                 << "p.y = "  << p.y << " mm" << std::endl
                 << "p.z = "  << p.z << " mm" << std::endl << std::endl
                 << "Direction:" << std::endl << std::endl
                 << "v.x = "  << v.x << std::endl
                 << "v.y = "  << v.y << std::endl
                 << "v.z = "  << v.z << std::endl << std::endl
                 << "Proposed distance :" << std::endl << std::endl
                 << "snxt = "    << snxt << " mm" << std::endl;
         message.precision(oldprc);
         UUtils::Exception("UTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
                           Warning, 1, message.str().c_str());
         break;
     }
   }
   if (snxt < halfCarTolerance)
   {
     snxt = 0;
   }
 
   return snxt;
 }

void UTubs::Extent	(	UVector3 &	aMin,
		UVector3 &	aMax
	)		const

virtual

Implements VUSolid.

Definition at line 1748 of file UTubs.cc.

References fDz, and fRMax.

 {
   aMin = UVector3(-fRMax, -fRMax, -fDz);
   aMax = UVector3(fRMax, fRMax, fDz);
 }

double UTubs::GetDeltaPhiAngle ( ) const

inline

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

double UTubs::GetDPhi ( ) const

inline

double UTubs::GetDz ( ) const

inline

Referenced by UPolycone::InsideSection().

UGeometryType UTubs::GetEntityType ( ) const

virtual

Implements VUSolid.

Definition at line 1638 of file UTubs.cc.

 {
   return std::string("Tubs");
 }

double UTubs::GetInnerRadius ( ) const

inline

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

double UTubs::GetOuterRadius ( ) const

inline

Referenced by G4UTubs::GetOuterRadius(), and GetParametersList().

void UTubs::GetParametersList	(	int	,
		double *	aArray
	)		const

virtual

Implements VUSolid.

Definition at line 1754 of file UTubs.cc.

References GetDeltaPhiAngle(), GetInnerRadius(), GetOuterRadius(), GetStartPhiAngle(), and GetZHalfLength().

 {
   aArray[0] = GetInnerRadius();
   aArray[1] = GetOuterRadius();
   aArray[2] = GetZHalfLength();
   aArray[3] = GetStartPhiAngle();
   aArray[4] = GetDeltaPhiAngle();
 }

UVector3 UTubs::GetPointOnSurface ( ) const

virtual

Implements VUSolid.

Definition at line 1679 of file UTubs.cc.

References fCosSPhi, fCosSPhiDPhi, fDPhi, fDz, fRMax, fRMin, fSinSPhi, fSinSPhiDPhi, fSPhi, UUtils::GetRadiusInRing(), and UUtils::Random().

 {
   double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
          aOne, aTwo, aThr, aFou;
   double rRand;
 
   aOne = 2.*fDz * fDPhi * fRMax;
   aTwo = 2.*fDz * fDPhi * fRMin;
   aThr = 0.5 * fDPhi * (fRMax * fRMax - fRMin * fRMin);
   aFou = 2.*fDz * (fRMax - fRMin);
 
   phi   = UUtils::Random(fSPhi, fSPhi + fDPhi);
   cosphi = std::cos(phi);
   sinphi = std::sin(phi);
 
   rRand = UUtils::GetRadiusInRing(fRMin, fRMax);
 
   if ((fSPhi == 0) && (fDPhi == 2 * UUtils::kPi))
   {
     aFou = 0;
   }
 
   chose = UUtils::Random(0., aOne + aTwo + 2.*aThr + 2.*aFou);
 
   if ((chose >= 0) && (chose < aOne))
   {
     xRand = fRMax * cosphi;
     yRand = fRMax * sinphi;
     zRand = UUtils::Random(-1.*fDz, fDz);
     return UVector3(xRand, yRand, zRand);
   }
   else if ((chose >= aOne) && (chose < aOne + aTwo))
   {
     xRand = fRMin * cosphi;
     yRand = fRMin * sinphi;
     zRand = UUtils::Random(-1.*fDz, fDz);
     return UVector3(xRand, yRand, zRand);
   }
   else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr))
   {
     xRand = rRand * cosphi;
     yRand = rRand * sinphi;
     zRand = fDz;
     return UVector3(xRand, yRand, zRand);
   }
   else if ((chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr))
   {
     xRand = rRand * cosphi;
     yRand = rRand * sinphi;
     zRand = -1.*fDz;
     return UVector3(xRand, yRand, zRand);
   }
   else if ((chose >= aOne + aTwo + 2.*aThr)
            && (chose < aOne + aTwo + 2.*aThr + aFou))
   {
     xRand = rRand * fCosSPhi;
     yRand = rRand * fSinSPhi;
     zRand = UUtils::Random(-1.*fDz, fDz);
     return UVector3(xRand, yRand, zRand);
   }
   else
   {
     xRand = rRand * fCosSPhiDPhi;
     yRand = rRand * fSinSPhiDPhi;
     zRand = UUtils::Random(-1.*fDz, fDz);
     return UVector3(xRand, yRand, zRand);
   }
 }

double UTubs::GetRMax ( ) const

inline

Referenced by UPolycone::InsideSection().

double UTubs::GetRMin ( ) const

inline

Referenced by UPolycone::InsideSection().

double UTubs::GetSPhi ( ) const

inline

double UTubs::GetStartPhiAngle ( ) const

inline

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

double UTubs::GetZHalfLength ( ) const

inline

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

void UTubs::Initialize ( )

inlineprotected

void UTubs::InitializeTrigonometry ( )

inlineprotected

VUSolid::EnumInside UTubs::Inside ( const UVector3 & p ) const

inlinevirtual

Implements VUSolid.

Definition at line 153 of file UTubs.cc.

References VUSolid::eInside, VUSolid::eOutside, VUSolid::eSurface, fDPhi, fDz, fPhiFullTube, fRMax, fRMin, fSPhi, kAngTolerance, kRadTolerance, VUSolid::Tolerance(), UVector3::x, UVector3::y, and UVector3::z.

Referenced by SafetyFromInside().

 {
   double r2, pPhi, tolRMin, tolRMax;
   VUSolid::EnumInside in = eOutside;
   static const double halfCarTolerance = VUSolid::Tolerance() * 0.5;
   static const double halfRadTolerance = kRadTolerance * 0.5;
   static const double halfAngTolerance = kAngTolerance * 0.5;
 
   if (std::fabs(p.z) <= fDz - halfCarTolerance)
   {
     r2 = p.x * p.x + p.y * p.y;
 
     if (fRMin)
     {
       tolRMin = fRMin + halfRadTolerance;
     }
     else
     {
       tolRMin = 0;
     }
 
     tolRMax = fRMax - halfRadTolerance;
 
     if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
     {
       if (fPhiFullTube)
       {
         in = eInside;
       }
       else
       {
         // Try inner tolerant phi boundaries (=>inside)
         // if not inside, try outer tolerant phi boundaries
 
         if ((tolRMin == 0) && (std::fabs(p.x) <= halfCarTolerance)
             && (std::fabs(p.y) <= halfCarTolerance))
         {
           in = eSurface;
         }
         else
         {
           pPhi = std::atan2(p.y, p.x);
           if (pPhi < -halfAngTolerance)
           {
             pPhi += 2 * UUtils::kPi;  // 0<=pPhi<2UUtils::kPi
           }
 
           if (fSPhi >= 0)
           {
             if ((std::fabs(pPhi) < halfAngTolerance)
                 && (std::fabs(fSPhi + fDPhi - 2 * UUtils::kPi) < halfAngTolerance))
             {
               pPhi += 2 * UUtils::kPi; // 0 <= pPhi < 2UUtils::kPi
             }
             if ((pPhi >= fSPhi + halfAngTolerance)
                 && (pPhi <= fSPhi + fDPhi - halfAngTolerance))
             {
               in = eInside;
             }
             else if ((pPhi >= fSPhi - halfAngTolerance)
                      && (pPhi <= fSPhi + fDPhi + halfAngTolerance))
             {
               in = eSurface;
             }
           }
           else  // fSPhi < 0
           {
             if ((pPhi <= fSPhi + 2 * UUtils::kPi - halfAngTolerance)
                 && (pPhi >= fSPhi + fDPhi + halfAngTolerance))
             {
               ; //eOutside
             }
             else if ((pPhi <= fSPhi + 2 * UUtils::kPi + halfAngTolerance)
                      && (pPhi >= fSPhi + fDPhi - halfAngTolerance))
             {
               in = eSurface;
             }
             else
             {
               in = eInside;
             }
           }
         }
       }
     }
     else  // Try generous boundaries
     {
       tolRMin = fRMin - halfRadTolerance;
       tolRMax = fRMax + halfRadTolerance;
 
       if (tolRMin < 0)
       {
         tolRMin = 0;
       }
 
       if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
       {
         if (fPhiFullTube || (r2 <= halfRadTolerance * halfRadTolerance))
         {
           // Continuous in phi or on z-axis
           in = eSurface;
         }
         else // Try outer tolerant phi boundaries only
         {
           pPhi = std::atan2(p.y, p.x);
 
           if (pPhi < -halfAngTolerance)
           {
             pPhi += 2 * UUtils::kPi;  // 0<=pPhi<2UUtils::kPi
           }
           if (fSPhi >= 0)
           {
             if ((std::fabs(pPhi) < halfAngTolerance)
                 && (std::fabs(fSPhi + fDPhi - 2 * UUtils::kPi) < halfAngTolerance))
             {
               pPhi += 2 * UUtils::kPi; // 0 <= pPhi < 2UUtils::kPi
             }
             if ((pPhi >= fSPhi - halfAngTolerance)
                 && (pPhi <= fSPhi + fDPhi + halfAngTolerance))
             {
               in = eSurface;
             }
           }
           else  // fSPhi < 0
           {
             if ((pPhi <= fSPhi + 2 * UUtils::kPi - halfAngTolerance)
                 && (pPhi >= fSPhi + fDPhi + halfAngTolerance))
             {
               ; // eOutside
             }
             else
             {
               in = eSurface;
             }
           }
         }
       }
     }
   }
   else if (std::fabs(p.z) <= fDz + halfCarTolerance)
   {
     // Check within tolerant r limits
     r2 = p.x * p.x + p.y * p.y;
     tolRMin = fRMin - halfRadTolerance;
     tolRMax = fRMax + halfRadTolerance;
 
     if (tolRMin < 0)
     {
       tolRMin = 0;
     }
 
     if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
     {
       if (fPhiFullTube || (r2 <= halfRadTolerance * halfRadTolerance))
       {
         // Continuous in phi or on z-axis
         in = eSurface;
       }
       else // Try outer tolerant phi boundaries
       {
         pPhi = std::atan2(p.y, p.x);
 
         if (pPhi < -halfAngTolerance)
         {
           pPhi += 2 * UUtils::kPi;  // 0<=pPhi<2UUtils::kPi
         }
         if (fSPhi >= 0)
         {
           if ((std::fabs(pPhi) < halfAngTolerance)
               && (std::fabs(fSPhi + fDPhi - 2 * UUtils::kPi) < halfAngTolerance))
           {
             pPhi += 2 * UUtils::kPi; // 0 <= pPhi < 2UUtils::kPi
           }
           if ((pPhi >= fSPhi - halfAngTolerance)
               && (pPhi <= fSPhi + fDPhi + halfAngTolerance))
           {
             in = eSurface;
           }
         }
         else  // fSPhi < 0
         {
           if ((pPhi <= fSPhi + 2 * UUtils::kPi - halfAngTolerance)
               && (pPhi >= fSPhi + fDPhi + halfAngTolerance))
           {
             ;
           }
           else
           {
             in = eSurface;
           }
         }
       }
     }
   }
   return in;
 }

bool UTubs::Normal	(	const UVector3 &	p,
		UVector3 &	normal
	)		const

virtual

Implements VUSolid.

Definition at line 356 of file UTubs.cc.

References ApproxSurfaceNormal(), UUtils::Exception(), fCosSPhi, fCosSPhiDPhi, fDPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinSPhi, fSinSPhiDPhi, fSPhi, kAngTolerance, VUSolid::Tolerance(), UVector3::Unit(), Warning, UVector3::x, UVector3::y, and UVector3::z.

 {
   int noSurfaces = 0;
   double rho, pPhi;
   double distZ, distRMin, distRMax;
   double distSPhi = UUtils::kInfinity, distEPhi = UUtils::kInfinity;
 
   static const double halfCarTolerance = 0.5 * VUSolid::Tolerance();
   static const double halfAngTolerance = 0.5 * kAngTolerance;
 
   UVector3 norm, sumnorm(0., 0., 0.);
   UVector3 nZ = UVector3(0, 0, 1.0);
   UVector3 nR, nPs, nPe;
 
   rho = std::sqrt(p.x * p.x + p.y * p.y);
 
   distRMin = std::fabs(rho - fRMin);
   distRMax = std::fabs(rho - fRMax);
   distZ   = std::fabs(std::fabs(p.z) - fDz);
 
   if (!fPhiFullTube)    // Protected against (0,0,z)
   {
     if (rho > halfCarTolerance)
     {
       pPhi = std::atan2(p.y, p.x);
 
       if (pPhi < fSPhi - halfCarTolerance)
       {
         pPhi += 2 * UUtils::kPi;
       }
       else if (pPhi > fSPhi + fDPhi + halfCarTolerance)
       {
         pPhi -= 2 * UUtils::kPi;
       }
 
       distSPhi = std::fabs(pPhi - fSPhi);
       distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
     }
     else if (!fRMin)
     {
       distSPhi = 0.;
       distEPhi = 0.;
     }
     nPs = UVector3(fSinSPhi, -fCosSPhi, 0);
     nPe = UVector3(-fSinSPhiDPhi /* std::sin(fSPhi+fDPhi)*/, fCosSPhiDPhi, 0);
   }
   if (rho > halfCarTolerance)
   {
     nR = UVector3(p.x / rho, p.y / rho, 0);
   }
 
   if (distRMax <= halfCarTolerance)
   {
     noSurfaces ++;
     sumnorm += nR;
   }
   if (fRMin && (distRMin <= halfCarTolerance))
   {
     noSurfaces ++;
     sumnorm -= nR;
   }
   if (fDPhi < 2 * UUtils::kPi)
   {
     if (distSPhi <= halfAngTolerance)
     {
       noSurfaces ++;
       sumnorm += nPs;
     }
     if (distEPhi <= halfAngTolerance)
     {
       noSurfaces ++;
       sumnorm += nPe;
     }
   }
   if (distZ <= halfCarTolerance)
   {
     noSurfaces ++;
     if (p.z >= 0.)
     {
       sumnorm += nZ;
     }
     else
     {
       sumnorm -= nZ;
     }
   }
   if (noSurfaces == 0)
   {
 #ifdef UDEBUG
     UUtils::Exception("UTubs::SurfaceNormal(p)", "GeomSolids1002",
                       Warning, 1, "Point p is not on surface !?");
     int oldprc = cout.precision(20);
     cout << "UTubs::SN ( " << p.x << ", " << p.y << ", " << p.z << " ); "
          << std::endl << std::endl;
     cout.precision(oldprc);
 #endif
     norm = ApproxSurfaceNormal(p);
   }
   else if (noSurfaces == 1)
   {
     norm = sumnorm;
   }
   else
   {
     norm = sumnorm.Unit();
   }
 
   n = norm;
 
   return noSurfaces; // TODO: return true or false on validity
 }

UTubs & UTubs::operator= ( const UTubs & rhs )

Definition at line 108 of file UTubs.cc.

References fCosCPhi, fCosEPhi, fCosHDPhiIT, fCosHDPhiOT, fCosSPhi, fCosSPhiDPhi, fCubicVolume, fDPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinCPhi, fSinEPhi, fSinSPhi, fSinSPhiDPhi, fSPhi, fSurfaceArea, kAngTolerance, and kRadTolerance.

 {
   // Check assignment to self
   //
   if (this == &rhs)
   {
     return *this;
   }
 
   // Copy base class data
   //
   VUSolid::operator=(rhs);
 
   // Copy data
   //
   fCubicVolume = rhs.fCubicVolume;
   fSurfaceArea = rhs.fSurfaceArea;
   kRadTolerance = rhs.kRadTolerance;
   kAngTolerance = rhs.kAngTolerance;
   fRMin = rhs.fRMin;
   fRMax = rhs.fRMax;
   fDz = rhs.fDz;
   fSPhi = rhs.fSPhi;
   fDPhi = rhs.fDPhi;
   fSinCPhi = rhs.fSinCPhi;
   fCosCPhi = rhs.fSinCPhi;
   fCosHDPhiOT = rhs.fCosHDPhiOT;
   fCosHDPhiIT = rhs.fCosHDPhiOT;
   fSinSPhi = rhs.fSinSPhi;
   fCosSPhi = rhs.fCosSPhi;
   fSinEPhi = rhs.fSinEPhi;
   fCosEPhi = rhs.fCosEPhi;
   fSinSPhiDPhi = rhs.fSinSPhiDPhi;
   fCosSPhiDPhi = rhs.fCosSPhiDPhi;
 
   fPhiFullTube = rhs.fPhiFullTube;
 
   return *this;
 }

double UTubs::SafetyFromInside	(	const UVector3 &	p,
		bool	precise = `false`
	)		const

virtual

Implements VUSolid.

Definition at line 1562 of file UTubs.cc.

References VUSolid::eOutside, UUtils::Exception(), fCosCPhi, fCosEPhi, fCosSPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinCPhi, fSinEPhi, fSinSPhi, Inside(), Warning, UVector3::x, UVector3::y, and UVector3::z.

 {
   double safe = 0.0, rho, safeR1, safeR2, safeZ, safePhi;
   rho = std::sqrt(p.x * p.x + p.y * p.y);
 
 #ifdef UDEBUG
   if (Inside(p) == eOutside)
   {
     int oldprc = cout.precision(16);
     cout << std::endl;
     DumpInfo();
     cout << "Position:" << std::endl << std::endl;
     cout << "p.x = "   << p.x << " mm" << std::endl;
     cout << "p.y = "   << p.y << " mm" << std::endl;
     cout << "p.z = "   << p.z << " mm" << std::endl << std::endl;
     cout.precision(oldprc);
     UUtils::Exception("UTubs::DistanceToOut(p)", "GeomSolids1002",
                       Warning, 1, "Point p is outside !?");
   }
 #endif
 
   if (fRMin)
   {
     safeR1 = rho   - fRMin;
     safeR2 = fRMax - rho;
 
     if (safeR1 < safeR2)
     {
       safe = safeR1;
     }
     else
     {
       safe = safeR2;
     }
   }
   else
   {
     safe = fRMax - rho;
   }
   safeZ = fDz - std::fabs(p.z);
 
   if (safeZ < safe)
   {
     safe = safeZ;
   }
 
   // Check if phi divided, Calc distances closest phi plane
   //
   if (!fPhiFullTube)
   {
     if (p.y * fCosCPhi - p.x * fSinCPhi <= 0)
     {
       safePhi = -(p.x * fSinSPhi - p.y * fCosSPhi);
     }
     else
     {
       safePhi = (p.x * fSinEPhi - p.y * fCosEPhi);
     }
     if (safePhi < safe)
     {
       safe = safePhi;
     }
   }
   if (safe < 0)
   {
     safe = 0;
   }
 
   return safe;
 }

double UTubs::SafetyFromOutside	(	const UVector3 &	p,
		bool	precise = `false`
	)		const

virtual

Implements VUSolid.

Definition at line 1030 of file UTubs.cc.

References fCosCPhi, fCosEPhi, fCosSPhi, fDPhi, fDz, fPhiFullTube, fRMax, fRMin, fSinCPhi, fSinEPhi, fSinSPhi, UVector3::x, UVector3::y, and UVector3::z.

Referenced by UEnclosingCylinder::SafetyFromOutside().

 {
   double safe = 0.0, rho, safe1, safe2, safe3;
   double safePhi, cosPsi;
 
   rho  = std::sqrt(p.x * p.x + p.y * p.y);
   safe1 = fRMin - rho;
   safe2 = rho - fRMax;
   safe3 = std::fabs(p.z) - fDz;
 
   if (safe1 > safe2)
   {
     safe = safe1;
   }
   else
   {
     safe = safe2;
   }
   if (safe3 > safe)
   {
     safe = safe3;
   }
 
   if ((!fPhiFullTube) && (rho))
   {
     // Psi=angle from central phi to point
     //
     cosPsi = (p.x * fCosCPhi + p.y * fSinCPhi) / rho;
 
     if (cosPsi < std::cos(fDPhi * 0.5))
     {
       // Point lies outside phi range
 
       if ((p.y * fCosCPhi - p.x * fSinCPhi) <= 0)
       {
         safePhi = std::fabs(p.x * fSinSPhi - p.y * fCosSPhi);
       }
       else
       {
         safePhi = std::fabs(p.x * fSinEPhi - p.y * fCosEPhi);
       }
       if (safePhi > safe)
       {
         safe = safePhi;
       }
     }
   }
   if (safe < 0)
   {
     safe = 0; return safe; // point is Inside;
   }
   if (!aAccurate) return safe;
   double safsq = 0.0;
   int count = 0;
   if (safe1 > 0)
   {
     safsq += safe1 * safe1;
     count++;
   }
   if (safe2 > 0)
   {
     safsq += safe2 * safe2;
     count++;
   }
   if (safe3 > 0)
   {
     safsq += safe3 * safe3;
     count++;
   }
   if (count == 1) return safe;
   return std::sqrt(safsq);
 }

void UTubs::SetDeltaPhiAngle ( double newDPhi )

inline

Referenced by G4UTubs::SetDeltaPhiAngle().

void UTubs::SetInnerRadius ( double newRMin )

inline

Referenced by G4UTubs::SetInnerRadius().

void UTubs::SetOuterRadius ( double newRMax )

inline

Referenced by G4UTubs::SetOuterRadius().

void UTubs::SetStartPhiAngle	(	double	newSPhi,
		bool	trig = `true`
	)

inline

Referenced by G4UTubs::SetStartPhiAngle().

void UTubs::SetZHalfLength ( double newDz )

inline

Referenced by G4UTubs::SetZHalfLength().

std::ostream & UTubs::StreamInfo ( std::ostream & os ) const

virtual

Implements VUSolid.

Definition at line 1656 of file UTubs.cc.

References fDPhi, fDz, fRMax, fRMin, fSPhi, and VUSolid::GetName().

 {
   int oldprc = os.precision(16);
   os << "-----------------------------------------------------------\n"
      << "       *** Dump for solid - " << GetName() << " ***\n"
      << "       ===================================================\n"
      << " Solid type: UTubs\n"
      << " Parameters: \n"
      << "       inner radius : " << fRMin << " mm \n"
      << "       outer radius : " << fRMax << " mm \n"
      << "       half length Z: " << fDz << " mm \n"
      << "       starting phi : " << fSPhi / (UUtils::kPi / 180.0) << " degrees \n"
      << "       delta phi       : " << fDPhi / (UUtils::kPi / 180.0) << " degrees \n"
      << "-----------------------------------------------------------\n";
   os.precision(oldprc);
 
   return os;
 }