#include <HelixExtrapolatorToLine2Order.h>

Inheritance diagram for HelixExtrapolatorToLine2Order:

Public Member Functions
virtual DirectionType	direction (double s) const
	Direction at pathlength s from the starting point. More...

DirectionTypeDouble	directionInDouble (double s) const
	Direction at pathlength s from the starting point in double precision. More...

	HelixExtrapolatorToLine2Order (const PositionType &point, const DirectionType &direction, const float curvature, const PropagationDirection propDir=alongMomentum)
	Constructor using point, direction and (transverse!) curvature. More...

	HelixExtrapolatorToLine2Order (const double &x0, const double &y0, const double &z0, const double &cosPhi0, const double &sinPhi0, const double &cosTheta, const double &sinTheta, const double &rho, const PropagationDirection propDir=alongMomentum)
	Fast constructor (for use by IterativeHelixExtrapolatorToLine). More...

virtual std::pair< bool, double >	pathLength (const GlobalPoint &point) const

virtual std::pair< bool, double >	pathLength (const Line &line) const

virtual PositionType	position (double s) const
	Position at pathlength s from the starting point. More...

PositionTypeDouble	positionInDouble (double s) const
	Position at pathlength s from the starting point in double precision. More...

virtual	~HelixExtrapolatorToLine2Order ()

Public Member Functions inherited from HelixLineExtrapolation
virtual	~HelixLineExtrapolation ()=default

Private Member Functions
virtual std::pair< bool, double >	pathLengthFromCoefficients (const double ceq[4]) const
	common part for propagation to point and line More...

int	solve2ndOrder (const double ceq[], double sol[]) const
	Solutions of 2nd order equation. More...

int	solve3rdOrder (const double ceq[], double sol[]) const
	Solutions of 3rd order equation. More...

Private Attributes
DirectionTypeDouble	theDirection

const PositionTypeDouble	thePosition

const PropagationDirection	thePropDir

const double	theRho

double	theSinTheta

Additional Inherited Members
Public Types inherited from HelixLineExtrapolation
typedef Basic3DVector< float >	DirectionType

typedef Basic3DVector< double >	DirectionTypeDouble

typedef Basic3DVector< float >	PositionType

typedef Basic3DVector< double >	PositionTypeDouble

Detailed Description

Calculates intersections of a helix with planes of any orientation using a parabolic approximation.

Definition at line 11 of file HelixExtrapolatorToLine2Order.h.

Constructor & Destructor Documentation

HelixExtrapolatorToLine2Order::HelixExtrapolatorToLine2Order	(	const PositionType &	point,
		const DirectionType &	direction,
		const float	curvature,
		const PropagationDirection	propDir = `alongMomentum`
	)

Constructor using point, direction and (transverse!) curvature.

Definition at line 7 of file HelixExtrapolatorToLine2Order.cc.

References AlCaHLTBitMon_ParallelJobs::p, EnergyCorrector::pt, mathSSE::sqrt(), theDirection, theSinTheta, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

                                                                      :
   thePosition(point),
   theRho(curvature),
   thePropDir(propDir)
 {
   //
   // Components of direction vector (with correct normalisation)
   //
   double px = direction.x();
   double py = direction.y();
   double pz = direction.z();
   double pt = px*px+py*py;
   double p = sqrt(pt+pz*pz);
   pt = sqrt(pt);
   theDirection = DirectionTypeDouble(px/pt,py/pt,pz/pt);
   theSinTheta = pt/p;
 }

HelixExtrapolatorToLine2Order::HelixExtrapolatorToLine2Order	(	const double &	x0,
		const double &	y0,
		const double &	z0,
		const double &	cosPhi0,
		const double &	sinPhi0,
		const double &	cosTheta,
		const double &	sinTheta,
		const double &	rho,
		const PropagationDirection	propDir = `alongMomentum`
	)

inline

Fast constructor (for use by IterativeHelixExtrapolatorToLine).

Definition at line 20 of file HelixExtrapolatorToLine2Order.h.

                                                                     :
     thePosition(x0,y0,z0),
     theDirection(cosPhi0,sinPhi0,cosTheta/sinTheta),
     theSinTheta(sinTheta),
     theRho(rho), 
     thePropDir(propDir) {}

virtual HelixExtrapolatorToLine2Order::~HelixExtrapolatorToLine2Order ( )

inlinevirtual

Definition at line 32 of file HelixExtrapolatorToLine2Order.h.

References direction(), directionInDouble(), dso_internal, geometryCSVtoXML::line, pathLength(), pathLengthFromCoefficients(), position(), positionInDouble(), alignCSCRings::s, solve2ndOrder(), and solve3rdOrder().

32 {}

Member Function Documentation

HelixLineExtrapolation::DirectionType HelixExtrapolatorToLine2Order::direction ( double s ) const

virtual

Direction at pathlength s from the starting point.

Implements HelixLineExtrapolation.

Definition at line 226 of file HelixExtrapolatorToLine2Order.cc.

References directionInDouble().

Referenced by ~HelixExtrapolatorToLine2Order().

                                                         {
   // use double precision result
 //   DirectionTypeDouble dir = directionInDouble(s);
 //   return DirectionType(dir.x(),dir.y(),dir.z());
    return DirectionType(directionInDouble(s));
 }

HelixExtrapolatorToLine2Order::DirectionTypeDouble HelixExtrapolatorToLine2Order::directionInDouble ( double s ) const

Direction at pathlength s from the starting point in double precision.

Definition at line 236 of file HelixExtrapolatorToLine2Order.cc.

References theDirection, theRho, theSinTheta, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by direction(), IterativeHelixExtrapolatorToLine::directionInDouble(), and ~HelixExtrapolatorToLine2Order().

                                                                 {
   // based on delta phi
   double dph = s*theRho*theSinTheta;
   return DirectionTypeDouble(theDirection.x()-(theDirection.y()+0.5*theDirection.x()*dph)*dph,
                  theDirection.y()+(theDirection.x()-0.5*theDirection.y()*dph)*dph,
                  theDirection.z());
 }

std::pair< bool, double > HelixExtrapolatorToLine2Order::pathLength ( const GlobalPoint & point ) const

virtual

Propagation status (true if valid) and (signed) path length along the helix from the starting point to the closest approach to the point. The starting point is given in the constructor.

Implements HelixLineExtrapolation.

Definition at line 32 of file HelixExtrapolatorToLine2Order.cc.

References Basic3DVector< T >::dot(), Basic3DVector< T >::mag2(), pathLengthFromCoefficients(), position(), theDirection, thePosition, theRho, Basic3DVector< T >::x(), and Basic3DVector< T >::y().

Referenced by IterativeHelixExtrapolatorToLine::genericPathLength(), and ~HelixExtrapolatorToLine2Order().

                                                                          {
   //
   PositionTypeDouble position(point);
   DirectionTypeDouble helix2(-0.5*theRho*theDirection.y(),
                  0.5*theRho*theDirection.x(),
                  0.);
   DirectionTypeDouble deltaPos(thePosition-position);
   //
   // coefficients of 3rd order equation
   //
   double ceq[4];
   ceq[3] = 2*helix2.mag2();
   // ceq[2] = 3*theDirection.dot(helix2) = 0 since they are orthogonal 
   ceq[2] = 0.;
   ceq[1] = theDirection.mag2()+2*deltaPos.dot(helix2);
   ceq[0] = deltaPos.dot(theDirection);
   //
   return pathLengthFromCoefficients(ceq);
 }

std::pair< bool, double > HelixExtrapolatorToLine2Order::pathLength ( const Line & line ) const

virtual

Propagation status (true if valid) and (signed) path length along the helix from the starting point to the closest approach to the line. The starting point is given in the constructor.

Implements HelixLineExtrapolation.

Definition at line 56 of file HelixExtrapolatorToLine2Order.cc.

References Line::direction(), Basic3DVector< T >::dot(), pathLengthFromCoefficients(), Line::position(), theDirection, thePosition, theRho, Basic3DVector< T >::x(), and Basic3DVector< T >::y().

                                                                  {
   //
   // Auxiliary vectors. Assumes that line.direction().mag()=1 !
   //
   PositionTypeDouble linePosition(line.position());
   DirectionTypeDouble lineDirection(line.direction());
   DirectionTypeDouble helix2(-0.5*theRho*theDirection.y(),
                  0.5*theRho*theDirection.x(),
                  0.);
   DirectionTypeDouble deltaPos(thePosition-linePosition);
   DirectionTypeDouble helix1p(theDirection-lineDirection*theDirection.dot(lineDirection));
   DirectionTypeDouble helix2p(helix2-lineDirection*helix2.dot(lineDirection));
   //
   // coefficients of 3rd order equation
   //
   double ceq[4];
   ceq[3] = 2*helix2.dot(helix2p);
   // ceq[2] = 3*helix1.dot(helix1p); 
   // since theDirection.dot(helix2)==0 equivalent to
   ceq[2] = 3*theDirection.dot(lineDirection)*helix2.dot(lineDirection);
   ceq[1] = theDirection.dot(helix1p)+2*deltaPos.dot(helix2p);
   ceq[0] = deltaPos.dot(helix1p);
   //
   return pathLengthFromCoefficients(ceq);
 }

std::pair< bool, double > HelixExtrapolatorToLine2Order::pathLengthFromCoefficients ( const double ceq[4] ) const

privatevirtual

common part for propagation to point and line

Definition at line 86 of file HelixExtrapolatorToLine2Order.cc.

References alongMomentum, anyDirection, mps_fire::i, oppositeToMomentum, StEvtSolProducer_cfi::solutions, solve3rdOrder(), thePropDir, and theSinTheta.

Referenced by pathLength(), and ~HelixExtrapolatorToLine2Order().

 {
   //
   // Solution of 3rd order equation
   //
   double solutions[3];
   unsigned int nRaw = solve3rdOrder(ceq,solutions);
   //
   // check compatibility with propagation direction
   //
   unsigned int nDir(0);
   for ( unsigned int i=0; i<nRaw; i++ ) {
     if ( thePropDir==anyDirection ||
      (solutions[i]>=0&&thePropDir==alongMomentum) ||
      (solutions[i]<=0&&thePropDir==oppositeToMomentum) )
       solutions[nDir++] = solutions[i];
   }
   if ( nDir==0 )  return std::make_pair(false,0.);
   //
   // check 2nd derivative
   //
   unsigned int nMin(0);
   for ( unsigned int i=0; i<nDir; i++ ) {
     double st = solutions[i];
     double deri2 = (3*ceq[3]*st+2*ceq[2])*st+ceq[1];
     if ( deri2>0. )  solutions[nMin++] = st;
   }
   if ( nMin==0 )  return std::make_pair(false,0.);
   //
   // choose smallest path length
   //
   double dSt = solutions[0];
   for ( unsigned int i=1; i<nMin; i++ ) {
     if ( fabs(solutions[i])<fabs(dSt) )  dSt = solutions[i];
   }
 
   return std::make_pair(true,dSt/theSinTheta);
 }

HelixLineExtrapolation::PositionType HelixExtrapolatorToLine2Order::position ( double s ) const

virtual

Position at pathlength s from the starting point.

Implements HelixLineExtrapolation.

Definition at line 205 of file HelixExtrapolatorToLine2Order.cc.

References positionInDouble().

Referenced by pathLength(), and ~HelixExtrapolatorToLine2Order().

                                                        {
   // use double precision result
 //   PositionTypeDouble pos = positionInDouble(s);
 //   return PositionType(pos.x(),pos.y(),pos.z());
   return PositionType(positionInDouble(s));
 }

HelixExtrapolatorToLine2Order::PositionTypeDouble HelixExtrapolatorToLine2Order::positionInDouble ( double s ) const

Position at pathlength s from the starting point in double precision.

Definition at line 215 of file HelixExtrapolatorToLine2Order.cc.

References theDirection, thePosition, theRho, theSinTheta, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by position(), IterativeHelixExtrapolatorToLine::positionInDouble(), and ~HelixExtrapolatorToLine2Order().

                                                                {
   // based on path length in the transverse plane
   double st = s*theSinTheta;
   return PositionTypeDouble(thePosition.x()+(theDirection.x()-st*0.5*theRho*theDirection.y())*st,
                 thePosition.y()+(theDirection.y()+st*0.5*theRho*theDirection.x())*st,
                 thePosition.z()+st*theDirection.z());
 }

int HelixExtrapolatorToLine2Order::solve2ndOrder	(	const double	ceq[],
		double	sol[]
	)		const

private

Solutions of 2nd order equation.

Definition at line 172 of file HelixExtrapolatorToLine2Order.cc.

References MillePedeFileConverter_cfg::e, and mathSSE::sqrt().

Referenced by solve3rdOrder(), and ~HelixExtrapolatorToLine2Order().

 {
   //
   double deq1 = coeff[1]*coeff[1];
   double deq2 = coeff[2]*coeff[0];
   if ( fabs(deq1)<FLT_MIN || fabs(deq2/deq1)>1.e-6 ) {
     //
     // Standard solution for quadratic equations
     //
     double deq = deq1+2*deq2;
     if ( deq<0. )  return 0;
     double ceq = -0.5*(coeff[1]+(coeff[1]>0?1:-1)*sqrt(deq));
     solutions[0] = -2*ceq/coeff[2];
     solutions[1] = coeff[0]/ceq;
     return 2;
   }
   else {
     //
     // Solution by expansion of sqrt(1+deq)
     //
     double ceq = coeff[1]/coeff[2];
     double deq = deq2/deq1;
     deq *= (1-deq/2);
     solutions[0] = -ceq*deq;
     solutions[1] = ceq*(2+deq);
     return 2;
   }
 }

int HelixExtrapolatorToLine2Order::solve3rdOrder	(	const double	ceq[],
		double	sol[]
	)		const

private

Solutions of 3rd order equation.

Definition at line 126 of file HelixExtrapolatorToLine2Order.cc.

References a, b, funct::cos(), mps_fire::i, M_PI, phi, funct::pow(), lumiQueryAPI::q, alignCSCRings::r, mps_fire::result, solve2ndOrder(), and mathSSE::sqrt().

Referenced by pathLengthFromCoefficients(), and ~HelixExtrapolatorToLine2Order().

 {
   //
   // Real 3rd order equation? Follow numerical recipes ..
   //
   if ( fabs(ceq[3])>FLT_MIN ) {
     int result(0);
     double q = (ceq[2]*ceq[2]-3*ceq[3]*ceq[1]) / (ceq[3]*ceq[3]) / 9.;
     double r = (2*ceq[2]*ceq[2]*ceq[2]-9*ceq[3]*ceq[2]*ceq[1]+27*ceq[3]*ceq[3]*ceq[0]) 
       / (ceq[3]*ceq[3]*ceq[3]) / 54.;
     double q3 = q*q*q;
     if ( r*r<q3 ) {
       double phi = acos(r/sqrt(q3))/3.;
       double rootq = sqrt(q);
       for ( int i=0; i<3; i++ ) {
     solutions[i] = -2*rootq*cos(phi) - ceq[2]/ceq[3]/3.;
     phi += 2./3.*M_PI;
       }
       result = 3;
     }
     else {
       double a = pow(fabs(r)+sqrt(r*r-q3),1./3.);
       if ( r>0. ) a *= -1;
       double b = fabs(a)>FLT_MIN ? q/a : 0.;
       solutions[0] = a + b - ceq[2]/ceq[3]/3.;
       result = 1;
     }
     return result;
   }
   //
   // Second order equation
   //
   else if ( fabs(ceq[2])>FLT_MIN ) {
     return solve2ndOrder(ceq,solutions);
   }
   else {
     //
     // Special case: linear equation
     //
     solutions[0] = -ceq[0]/ceq[1];
     return 1;
   }
 }