#include <ThirdHitPredictionFromCircle.h>

Classes
class	HelixRZ

Public Types
typedef PixelRecoRange< float >	Range

Public Member Functions
double	angle (double curvature, double radius) const

Range	curvature (double transverseIP) const

double	curvature (const Basic2DVector< double > &thirdPoint) const

Range	operator() (Range curvature, double radius) const

double	phi (double curvature, double radius) const

	ThirdHitPredictionFromCircle (const GlobalPoint &P1, const GlobalPoint &P2, float tolerance)

double	transverseIP (const Basic2DVector< double > &thirdPoint) const

Private Member Functions
double	invCenterOnAxis (const Basic2DVector< double > &thirdPoint) const

Private Attributes
Basic2DVector< double >	axis

Basic2DVector< double >	center

double	delta

double	delta2

Basic2DVector< double >	p1

double	theTolerance

Friends
class	HelixRZ

Detailed Description

Definition at line 11 of file ThirdHitPredictionFromCircle.h.

Member Typedef Documentation

typedef PixelRecoRange<float> ThirdHitPredictionFromCircle::Range

Definition at line 14 of file ThirdHitPredictionFromCircle.h.

Constructor & Destructor Documentation

ThirdHitPredictionFromCircle::ThirdHitPredictionFromCircle	(	const GlobalPoint &	P1,
		const GlobalPoint &	P2,
		float	tolerance
	)

Definition at line 33 of file ThirdHitPredictionFromCircle.cc.

References axis, center, delta, delta2, diffTreeTool::diff, p1, p2, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

   : p1(P1.x(), P1.y()), theTolerance(tolerance)
 {
   Point2D p2(P2.x(), P2.y());
   Point2D diff = 0.5 * (p2 - p1);
   delta2 = diff.mag2();
   delta = std::sqrt(delta2);
   axis = Point2D(-diff.y(), diff.x()) / delta;
   center = p1 + diff;
 }

Member Function Documentation

double ThirdHitPredictionFromCircle::angle	(	double	curvature,
		double	radius
	)		const

Definition at line 68 of file ThirdHitPredictionFromCircle.cc.

References abs, axis, center, funct::cos(), delta2, alignCSCRings::e, Basic2DVector< T >::mag2(), FWPFMaths::sgn(), funct::sin(), funct::sqr(), and unlikely.

Referenced by MatchedHitRZCorrectionFromBending::tibMatchedHitZFixup().

 {
   if (unlikely(std::abs(curvature) < 1.0e-5)) {
     double sin = (center * axis) / radius;
     return sin / clamped_sqrt(1 - sqr(sin));
   } else {
     double radius2 = sqr(1.0 / curvature);
     double orthog = clamped_sqrt(radius2 - delta2);
     Basic2DVector<double> lcenter = center - sgn(curvature) * orthog * axis;
  
     double cos = (radius2 + sqr(radius) - lcenter.mag2()) *
                  curvature / (2. * radius);
     return - cos / clamped_sqrt(1 - sqr(cos));
  }
 }

ThirdHitPredictionFromCircle::Range ThirdHitPredictionFromCircle::curvature ( double transverseIP ) const

Definition at line 96 of file ThirdHitPredictionFromCircle.cc.

References abs, axis, center, delta, delta2, alignCSCRings::e, FWPFMaths::sgn(), funct::sqr(), mathSSE::sqrt(), std::swap(), tmp, transverseIP(), unlikely, Basic2DVector< T >::x(), and Basic2DVector< T >::y().

Referenced by ThirdHitPredictionFromCircle::HelixRZ::rAtZ(), and ThirdHitPredictionFromCircle::HelixRZ::zAtR().

 {
   // this is a mess.  Use a CAS and lots of drawings to verify...
 
   transverseIP = std::abs(transverseIP);
   double transverseIP2 = sqr(transverseIP);
   double tip = axis * center;
   double tip2 = sqr(tip);
   double lip = axis.x() * center.y() - axis.y() * center.x();
   double lip2 = sqr(lip);
 
   double origin = std::sqrt(tip2 + lip2);
   double tmp1 = lip2 + tip2 - transverseIP2;
   double tmp2 = 2. * (tip - transverseIP) * (tip + transverseIP);
   double tmp3 = 2. * delta * origin;
   double tmp4 = tmp1 + delta2;
   double tmp5 = 2. * delta * lip;
 
   // I am probably being overly careful here with border cases
   // but you never know what crap you might get fed
 
   double u1, u2;
   if (unlikely(tmp4 - tmp5 < 1.0e-5)) {
     u1 = -0.;   // yes, I am making use of signed zero
     u2 = +0.;   // -> no -ffast-math please
   } else {
     if (unlikely(std::abs(tmp2) < 1.0e-5)) {
       // the denominator is zero
       // this means that one of the tracks will be straight
       // and the other can be computed from the limit of the equation
       double tmp = lip2 - delta2;
       u1 = INFINITY;    // and require 1 / sqrt(inf^2 + x) = 0 (with x > 0)
       u2 = (sqr(0.5 * tmp) - delta2 * tip2) / (tmp * tip);
       if (tip < 0)
         std::swap(u1, u2);
     } else {
       double tmp6 = (tmp4 - tmp5) * (tmp4 + tmp5);
       if (unlikely(tmp6 < 1.0e-5)) {
         u1 = -0.;
         u2 = +0.;
       } else {
         double tmp7 = tmp6 > 0 ? (transverseIP * std::sqrt(tmp6) / tmp2) : 0.;
         double tmp8 = tip * (tmp1 - delta2) / tmp2;
         // the two quadratic solutions
         u1 = tmp8 + tmp7;
         u2 = tmp8 - tmp7;
       }
     }
 
     if (tmp4 <= std::abs(tmp3)) {
       if ((tmp3 < 0) == (tip < 0))
         u2 = +0.;
       else
         u1 = -0.;
     }
   }
 
   return Range(sgn(u1) / std::sqrt(sqr(u1) + delta2),
                sgn(u2) / std::sqrt(sqr(u2) + delta2));
 }

double ThirdHitPredictionFromCircle::curvature ( const Basic2DVector< double > & thirdPoint ) const

Definition at line 167 of file ThirdHitPredictionFromCircle.cc.

References delta2, invCenterOnAxis(), FWPFMaths::sgn(), funct::sqr(), and mathSSE::sqrt().

 {
   double invDist = invCenterOnAxis(p2);
   double invDist2 = sqr(invDist);
   double curv = std::sqrt(invDist2 / (1. + invDist2 * delta2));
   return sgn(invDist) * curv;
 }

double ThirdHitPredictionFromCircle::invCenterOnAxis ( const Basic2DVector< double > & thirdPoint ) const

private

Definition at line 157 of file ThirdHitPredictionFromCircle.cc.

References a, axis, b, center, delta, diffTreeTool::diff, Basic2DVector< T >::mag(), p1, Basic2DVector< T >::x(), and Basic2DVector< T >::y().

Referenced by curvature(), and transverseIP().

 {
   Point2D delta = p2 - p1;
   Point2D axis2 = Point2D(-delta.y(), delta.x()) / delta.mag();
   Point2D diff = p1 + 0.5 * delta - center;
   double a = diff.y() * axis2.x() - diff.x() * axis2.y();
   double b = axis.y() * axis2.x() - axis.x() * axis2.y();
   return b / a;
 }

ThirdHitPredictionFromCircle::Range ThirdHitPredictionFromCircle::operator()	(	Range	curvature,
		double	radius
	)		const

Definition at line 85 of file ThirdHitPredictionFromCircle.cc.

References M_PI, phi(), theTolerance, and unlikely.

 {
   double phi1 = phi(curvature.second, radius);
   double phi2 = phi(curvature.first, radius);
 
   while(unlikely(phi2 <  phi1)) phi2 += 2. * M_PI; 
 
   return Range(phi1 * radius - theTolerance, phi2 * radius + theTolerance);
 }

double ThirdHitPredictionFromCircle::phi	(	double	curvature,
		double	radius
	)		const

Definition at line 45 of file ThirdHitPredictionFromCircle.cc.

References abs, axis, center, funct::cos(), delta2, alignCSCRings::e, M_PI, Basic2DVector< T >::mag2(), Basic2DVector< T >::phi(), CosmicsPD_Skims::radius, FWPFMaths::sgn(), funct::sqr(), mathSSE::sqrt(), and unlikely.

Referenced by operator()(), and ThirdHitPredictionFromCircle::HelixRZ::rAtZ().

 {
   double phi;
   if (unlikely(std::abs(curvature) < 1.0e-5)) {
     double cos = (center * axis) / radius;
     phi = axis.phi() - clamped_acos(cos);
   } else {
     double sign = sgn(curvature);
     double radius2 = sqr(1.0 / curvature);
     double orthog = clamped_sqrt(radius2 - delta2);
     Basic2DVector<double> lcenter = center - sign * orthog * axis;
     double rc2 = lcenter.mag2();
     double cos = (rc2 + sqr(radius) - radius2) /
       (2. *std:: sqrt(rc2) * radius);
     phi = lcenter.phi() + sign * clamped_acos(cos);
  }
 
   while(unlikely(phi >= M_PI)) phi -= 2. * M_PI;
   while(unlikely(phi < -M_PI)) phi += 2. * M_PI;
 
   return phi;
 }

double ThirdHitPredictionFromCircle::transverseIP ( const Basic2DVector< double > & thirdPoint ) const

Definition at line 175 of file ThirdHitPredictionFromCircle.cc.

References abs, axis, center, delta2, alignCSCRings::e, invCenterOnAxis(), mag(), CosmicsPD_Skims::radius, funct::sqr(), mathSSE::sqrt(), and unlikely.

Referenced by curvature().

 {
   double invDist = invCenterOnAxis(p2);
   if (unlikely(std::abs(invDist) < 1.0e-5))
     return std::abs(p2 * axis);
   else {
     double dist = 1.0 / invDist;
     double radius = std::sqrt(sqr(dist) + delta2);
     return std::abs((center + axis * dist).mag() - radius);
   }
 }