#include <ThirdHitPredictionFromInvParabola.h>

Classes
class	MappedPoint
Public Types
typedef PixelRecoRange< float >	Range
typedef TkRotation< double >	Rotation
Public Member Functions
void	init (const GlobalPoint &P1, const GlobalPoint &P2, double ip, double curv)
Range	operator() (double radius, int charge) const
Range	rangeRPhi (double radius, int charge) const
Range	rangeRPhiSlow (double radius, int charge, int nIter=5) const
	ThirdHitPredictionFromInvParabola (const GlobalPoint &P1, const GlobalPoint &P2, double ip, double curv, double tolerance)
Private Types
typedef MappedPoint< double >	PointUV
Private Member Functions
double	coeffA (double impactParameter, int charge) const
double	coeffB (double impactParameter, int charge) const
PointUV	findPointAtCurve (double radius, int charge, double ip) const
double	ipFromCurvature (double curvature, int charge) const
double	predV (double u, double ip, int charge) const
Private Attributes
PointUV	p1
PointUV	p2
Range	theIpRangeMinus
Range	theIpRangePlus
Rotation	theRotation
double	theTolerance

Detailed Description

Definition at line 24 of file ThirdHitPredictionFromInvParabola.h.

Member Typedef Documentation

typedef MappedPoint<double> ThirdHitPredictionFromInvParabola::PointUV [private]

Definition at line 79 of file ThirdHitPredictionFromInvParabola.h.

typedef PixelRecoRange<float> ThirdHitPredictionFromInvParabola::Range

Definition at line 29 of file ThirdHitPredictionFromInvParabola.h.

typedef TkRotation<double> ThirdHitPredictionFromInvParabola::Rotation

Definition at line 28 of file ThirdHitPredictionFromInvParabola.h.

Constructor & Destructor Documentation

ThirdHitPredictionFromInvParabola::ThirdHitPredictionFromInvParabola	(	const GlobalPoint &	P1,
		const GlobalPoint &	P2,
		double	ip,
		double	curv,
		double	tolerance
	)

Definition at line 25 of file ThirdHitPredictionFromInvParabola.cc.

References init().

  : theTolerance(torlerance)
{
  init(P1,P2,ip,fabs(curv));
}

Member Function Documentation

double ThirdHitPredictionFromInvParabola::coeffA	(	double	impactParameter,
		int	charge
	)		const `[inline, private]`

Definition at line 91 of file ThirdHitPredictionFromInvParabola.h.

References p1, p2, ThirdHitPredictionFromInvParabola::MappedPoint< T >::u(), and ThirdHitPredictionFromInvParabola::MappedPoint< T >::v().

Referenced by findPointAtCurve(), and predV().

{
  double u1u2 = p1.u()*p2.u();
  double du = p2.u() - p1.u();
  double pv = p1.v()*p2.u() - p2.v()*p1.u();
  return -charge*pv/du - u1u2*impactParameter;
}

double ThirdHitPredictionFromInvParabola::coeffB	(	double	impactParameter,
		int	charge
	)		const `[inline, private]`

Definition at line 100 of file ThirdHitPredictionFromInvParabola.h.

References p1, p2, ThirdHitPredictionFromInvParabola::MappedPoint< T >::u(), and ThirdHitPredictionFromInvParabola::MappedPoint< T >::v().

Referenced by findPointAtCurve(), and predV().

{
  double dv = p2.v() - p1.v();
  double du = p2.u() - p1.u();
  double su = p2.u() + p1.u();
  return charge*dv/du - su*impactParameter;
}

ThirdHitPredictionFromInvParabola::PointUV ThirdHitPredictionFromInvParabola::findPointAtCurve	(	double	radius,
		int	charge,
		double	ip
	)		const `[private]`

Definition at line 61 of file ThirdHitPredictionFromInvParabola.cc.

References funct::A, alpha, coeffA(), coeffB(), delta, csvReporter::r, mathSSE::sqrt(), theRotation, and v.

Referenced by rangeRPhi().

{
  //
  // assume u=(1-alpha^2/2)/r v=alpha/r
  // solve qudratic equation neglecting aplha^4 term
  //
  double A = coeffA(ip,c);
  double B = coeffB(ip,c);

  double overR = 1./r;
  double ipOverR = ip*overR;

  double delta = 1-4*(0.5*B+ipOverR)*(-B+A*r-ipOverR);
  double sqrtdelta = (delta > 0) ? std::sqrt(delta) : 0.;
  double alpha = (c>0)?  (-c+sqrtdelta)/(B+2*ipOverR) :  (-c-sqrtdelta)/(B+2*ipOverR);

  double v = alpha*overR;
  double d2 = overR*overR - v*v;
  double u = (d2 > 0) ? std::sqrt(d2) : 0.;

  return PointUV(u,v,&theRotation);
}

void ThirdHitPredictionFromInvParabola::init	(	const GlobalPoint &	P1,
		const GlobalPoint &	P2,
		double	ip,
		double	curv
	)

Definition at line 34 of file ThirdHitPredictionFromInvParabola.cc.

References abs, PixelRecoRange< T >::intersection(), ipFromCurvature(), p1, p2, PixelRecoRange< T >::sort(), theIpRangeMinus, theIpRangePlus, theRotation, Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

Referenced by ThirdHitPredictionFromInvParabola().

{
//  GlobalVector aX = GlobalVector( P2.x()-P1.x(), P2.y()-P1.y(), 0.).unit();
  GlobalVector aX = GlobalVector( P1.x(), P1.y(), 0.).unit();
  GlobalVector aY( -aX.y(), aX.x(), 0.); 
  GlobalVector aZ( 0., 0., 1.);
  theRotation = Rotation(aX,aY,aZ); 

  p1 = PointUV(Point2D(P1.x(),P1.y()), &theRotation);
  p2 = PointUV(Point2D(P2.x(),P2.y()), &theRotation);

  Range ipRange(-ip, ip); 
  ipRange.sort();
  
  double ipIntyPlus = ipFromCurvature(0.,1);
  double ipCurvPlus = ipFromCurvature(std::abs(curv), 1);
  double ipCurvMinus = ipFromCurvature(std::abs(curv), -1);

  
  Range ipRangePlus = Range(ipIntyPlus, ipCurvPlus); ipRangePlus.sort();
  Range ipRangeMinus = Range(-ipIntyPlus, ipCurvMinus); ipRangeMinus.sort();

  theIpRangePlus  = ipRangePlus.intersection(ipRange);
  theIpRangeMinus = ipRangeMinus.intersection(ipRange);
}

double ThirdHitPredictionFromInvParabola::ipFromCurvature	(	double	curvature,
		int	charge
	)		const `[private]`

Definition at line 165 of file ThirdHitPredictionFromInvParabola.cc.

References p1, p2, ThirdHitPredictionFromInvParabola::MappedPoint< T >::u(), and ThirdHitPredictionFromInvParabola::MappedPoint< T >::v().

Referenced by init().

{
  double u1u2 = p1.u()*p2.u();
  double du = p2.u() - p1.u();
  double pv = p1.v()*p2.u() - p2.v()*p1.u();

  double inInf = -charge*pv/(du*u1u2);
  return inInf-curvature/(2.*u1u2);
}

Range ThirdHitPredictionFromInvParabola::operator()	(	double	radius,
		int	charge
	)		const `[inline]`

Definition at line 35 of file ThirdHitPredictionFromInvParabola.h.

References rangeRPhi().

{ return rangeRPhi(radius,charge); }

double ThirdHitPredictionFromInvParabola::predV	(	double	u,
		double	ip,
		int	charge
	)		const `[private]`

Definition at line 178 of file ThirdHitPredictionFromInvParabola.cc.

References coeffA(), coeffB(), and funct::sqr().

Referenced by rangeRPhiSlow().

{
  return -charge*( coeffA(ip,charge) - coeffB(ip,charge)*u - ip*sqr(u));
}

ThirdHitPredictionFromInvParabola::Range ThirdHitPredictionFromInvParabola::rangeRPhi	(	double	radius,
		int	charge
	)		const

Definition at line 85 of file ThirdHitPredictionFromInvParabola.cc.

References PixelRecoRange< T >::empty(), findPointAtCurve(), PixelRecoRange< T >::intersection(), M_PI, PixelRecoRange< T >::max(), PixelRecoRange< T >::min(), Basic2DVector< T >::phi(), csvReporter::r, PixelRecoRange< T >::sort(), theIpRangeMinus, theIpRangePlus, theTolerance, ThirdHitPredictionFromInvParabola::MappedPoint< T >::unmap(), and ThirdHitPredictionFromInvParabola::MappedPoint< T >::v().

Referenced by operator()().

{
  Range predRPhi(1.,-1.);
  Range ip = (charge > 0) ? theIpRangePlus : theIpRangeMinus;

  PointUV pred_tmp1 = findPointAtCurve(radius,charge,ip.min());
  PointUV pred_tmp2 = findPointAtCurve(radius,charge,ip.max());

  double phi1 = pred_tmp1.unmap().phi();
  while ( phi1 >= M_PI) phi1 -= 2*M_PI;
  while ( phi1 < -M_PI) phi1 += 2*M_PI;
  double phi2 = phi1+radius*(pred_tmp2.v()-pred_tmp1.v()); 
  
  if (ip.empty()) {
    Range r1(phi1*radius-theTolerance, phi1*radius+theTolerance); 
    Range r2(phi2*radius-theTolerance, phi2*radius+theTolerance); 
    predRPhi = r1.intersection(r2);
  } else {
    Range r(phi1, phi2); 
    r.sort();
    predRPhi= Range(radius*r.min()-theTolerance, radius*r.max()+theTolerance);
  }

  return predRPhi;
}

ThirdHitPredictionFromInvParabola::Range ThirdHitPredictionFromInvParabola::rangeRPhiSlow	(	double	radius,
		int	charge,
		int	nIter = `5`
	)		const

Definition at line 113 of file ThirdHitPredictionFromInvParabola.cc.

References DeDxDiscriminatorTools::charge(), PixelRecoRange< T >::empty(), i, PixelRecoRange< T >::intersection(), M_PI, PixelRecoRange< T >::max(), PixelRecoRange< T >::min(), Basic2DVector< T >::phi(), predV(), csvReporter::r, CosmicsPD_Skims::radius, PixelRecoRange< T >::sort(), funct::sqr(), mathSSE::sqrt(), theIpRangeMinus, theIpRangePlus, theRotation, theTolerance, ThirdHitPredictionFromInvParabola::MappedPoint< T >::unmap(), and v.

{
  Range predRPhi(1.,-1.);

  double invr2 = 1/radius/radius;
  double u = sqrt(invr2);
  double v = 0.;

  Range ip = (charge > 0) ? theIpRangePlus : theIpRangeMinus;

  for (int i=0; i < nIter; ++i) {
    v = predV(u, ip.min(), charge); 
    double d2 = invr2-sqr(v);
    u = (d2 > 0) ? sqrt(d2) : 0.;
  }
  PointUV  pred_tmp1(u, v,  &theRotation);
  double phi1 = pred_tmp1.unmap().phi(); 
  while ( phi1 >= M_PI) phi1 -= 2*M_PI;
  while ( phi1 < -M_PI) phi1 += 2*M_PI;


  u = sqrt(invr2); 
  v=0;
  for (int i=0; i < nIter; ++i) {
    v = predV(u, ip.max(), charge); 
    double d2 = invr2-sqr(v);
    u = (d2 > 0) ? sqrt(d2) : 0.;
  }
  PointUV  pred_tmp2(u, v,  &theRotation);
  double phi2 = pred_tmp2.unmap().phi(); 
  while ( phi2-phi1 >= M_PI) phi2 -= 2*M_PI;
  while ( phi2-phi1 < -M_PI) phi2 += 2*M_PI;

// check faster alternative, without while(..) it is enough to:
//  phi2 = phi1+radius*(pred_tmp2.v()-pred_tmp1.v()); 

  if (ip.empty()) {
    Range r1(phi1*radius-theTolerance, phi1*radius+theTolerance); 
    Range r2(phi2*radius-theTolerance, phi2*radius+theTolerance); 
    predRPhi = r1.intersection(r2);
  } else {
    Range r(phi1, phi2); 
    r.sort();
    predRPhi= Range(radius*r.min()-theTolerance, radius*r.max()+theTolerance);
  }
  return predRPhi;

}