#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, 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, FWPFMaths::sgn(), funct::sqr(), mathSSE::sqrt(), 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, 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. * 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, 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);
  }
}