#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
double	dv
double	overDu
double	pv
double	su
Range	theIpRangeMinus
Range	theIpRangePlus
Rotation	theRotation
double	theTolerance
double	u1u2

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 93 of file ThirdHitPredictionFromInvParabola.h.

References overDu, pv, and u1u2.

Referenced by findPointAtCurve(), and predV().

{
  return -charge*pv*overDu - u1u2*impactParameter;
}

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

Definition at line 99 of file ThirdHitPredictionFromInvParabola.h.

References dv, overDu, and su.

Referenced by findPointAtCurve(), and predV().

{
  return charge*dv*overDu - su*impactParameter;
}

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

Definition at line 67 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, dv, PixelRecoRange< T >::intersection(), ipFromCurvature(), overDu, p1, p2, pv, PixelRecoRange< T >::sort(), su, theIpRangeMinus, theIpRangePlus, theRotation, ThirdHitPredictionFromInvParabola::MappedPoint< T >::u(), u1u2, Vector3DBase< T, FrameTag >::unit(), ThirdHitPredictionFromInvParabola::MappedPoint< T >::v(), 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); 

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

  u1u2 = p1.u()*p2.u();
  overDu = 1./(p2.u() - p1.u());
  pv = p1.v()*p2.u() - p2.v()*p1.u();
  dv = p2.v() - p1.v();
  su = p2.u() + p1.u();

  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 `[inline, private]`

Definition at line 105 of file ThirdHitPredictionFromInvParabola.h.

References overDu, pv, and u1u2.

Referenced by init().

{
  double overU1u2 = 1./u1u2;
  double inInf = -charge*pv*overDu*overU1u2;
  return inInf-curvature*overU1u2*0.5;
}

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 171 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 91 of file ThirdHitPredictionFromInvParabola.cc.

References PixelRecoRange< T >::empty(), findPointAtCurve(), PixelRecoRange< T >::intersection(), M_PI, PixelRecoRange< T >::max(), PixelRecoRange< T >::min(), 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 119 of file ThirdHitPredictionFromInvParabola.cc.

References DeDxDiscriminatorTools::charge(), PixelRecoRange< T >::empty(), i, PixelRecoRange< T >::intersection(), M_PI, PixelRecoRange< T >::max(), PixelRecoRange< T >::min(), 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;

}