#include <IterativeHelixExtrapolatorToLine.h>

Inheritance diagram for IterativeHelixExtrapolatorToLine:

Public Member Functions
virtual DirectionType	direction (double s) const
DirectionTypeDouble	directionInDouble (double s) const
	IterativeHelixExtrapolatorToLine (const PositionType &point, const DirectionType &direction, const float curvature, const PropagationDirection propDir=anyDirection)
virtual std::pair< bool, double >	pathLength (const Line &line) const
virtual std::pair< bool, double >	pathLength (const GlobalPoint &point) const
virtual PositionType	position (double s) const
PositionTypeDouble	positionInDouble (double s) const
virtual	~IterativeHelixExtrapolatorToLine ()
Private Member Functions
template<class T >
std::pair< bool, double >	genericPathLength (const T &object) const
	common functionality for extrapolation to line or point
Private Attributes
double	theCachedCDPhi
double	theCachedDPhi
double	theCachedS
double	theCachedSDPhi
double	theCosPhi0
double	theCosTheta
const PropagationDirection	thePropDir
HelixExtrapolatorToLine2Order	theQuadraticSolutionFromStart
const double	theRho
double	theSinPhi0
double	theSinTheta
const double	theX0
const double	theY0
const double	theZ0

Detailed Description

Calculates closest approach of a helix to a line or a point by iterative use of a 2nd order expansion of the helix.

Definition at line 12 of file IterativeHelixExtrapolatorToLine.h.

Constructor & Destructor Documentation

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

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

Definition at line 5 of file IterativeHelixExtrapolatorToLine.cc.

References L1TEmulatorMonitor_cff::p, mathSSE::sqrt(), Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

                                     :
  theX0(point.x()),
  theY0(point.y()),
  theZ0(point.z()),
  theRho(curvature),
  theQuadraticSolutionFromStart(point,direction,curvature,propDir),
  thePropDir(propDir),
  theCachedS(0),
  theCachedDPhi(0.),
  theCachedSDPhi(0.),
  theCachedCDPhi(1.)
{
  //
  // 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);
  theCosPhi0 = px/pt;
  theSinPhi0 = py/pt;
  theCosTheta = pz/p;
  theSinTheta = pt/p;
}

virtual IterativeHelixExtrapolatorToLine::~IterativeHelixExtrapolatorToLine ( ) [inline, virtual]

Definition at line 21 of file IterativeHelixExtrapolatorToLine.h.

{}

Member Function Documentation

HelixLineExtrapolation::DirectionType IterativeHelixExtrapolatorToLine::direction ( double s ) const [virtual]

Direction at pathlength s from the starting point.

Implements HelixLineExtrapolation.

Definition at line 163 of file IterativeHelixExtrapolatorToLine.cc.

References directionInDouble().

Referenced by AnalyticalTrajectoryExtrapolatorToLine::propagateWithHelix(), and AnalyticalImpactPointExtrapolator::propagateWithHelix().

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

HelixLineExtrapolation::DirectionTypeDouble IterativeHelixExtrapolatorToLine::directionInDouble ( double s ) const

Direction at pathlength s from the starting point.

Definition at line 174 of file IterativeHelixExtrapolatorToLine.cc.

References funct::cos(), HelixExtrapolatorToLine2Order::directionInDouble(), asciidump::s, funct::sin(), theCachedCDPhi, theCachedDPhi, theCachedS, theCachedSDPhi, theCosPhi0, theCosTheta, theQuadraticSolutionFromStart, theRho, theSinPhi0, and theSinTheta.

Referenced by direction(), and genericPathLength().

                                                                   {
  //
  // Calculate delta phi (if not already available)
  //
  if ( s!=theCachedS ) {
    theCachedS = s;
    theCachedDPhi = theCachedS*theRho*theSinTheta;
    theCachedSDPhi = sin(theCachedDPhi);
    theCachedCDPhi = cos(theCachedDPhi);
  }

  if ( fabs(theCachedDPhi)>1.e-4 ) {
    // full helix formula
    return DirectionTypeDouble(theCosPhi0*theCachedCDPhi-theSinPhi0*theCachedSDPhi,
                               theSinPhi0*theCachedCDPhi+theCosPhi0*theCachedSDPhi,
                               theCosTheta/theSinTheta);
  }
  else {
    // 1st order
    return theQuadraticSolutionFromStart.directionInDouble(theCachedS);
  }
}

template<class T >

std::pair< bool, double > IterativeHelixExtrapolatorToLine::genericPathLength ( const T & object ) const [private]

common functionality for extrapolation to line or point

Definition at line 38 of file IterativeHelixExtrapolatorToLine.cc.

References alongMomentum, anyDirection, directionInDouble(), first, align_cfg::iteration, oppositeToMomentum, HelixExtrapolatorToLine2Order::pathLength(), positionInDouble(), theCosPhi0, theCosTheta, thePropDir, theQuadraticSolutionFromStart, theRho, theSinPhi0, theSinTheta, theX0, theY0, theZ0, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by pathLength().

                                                                          {
  //
  // Constants used for control of convergence
  //
  const int maxIterations(100);
  //
  // Prepare internal value of the propagation direction and position / direction vectors for iteration 
  //
  PropagationDirection propDir = thePropDir;
  PositionTypeDouble xnew(theX0,theY0,theZ0);
  DirectionTypeDouble pnew(theCosPhi0,theSinPhi0,theCosTheta/theSinTheta);
  //
  // Prepare iterations: count and total pathlength
  //
  int iteration(maxIterations);
  double dSTotal(0.);
  //
  // Convergence criterion: maximal lateral displacement in a step < 1um
  //
  double maxDeltaS2 = 2*1.e-4/theSinTheta/theSinTheta/fabs(theRho);
  //
  bool first(true);
  while ( true ) {
    //
    // return empty solution vector if no convergence after maxIterations iterations
    //
    if ( --iteration<0 ) {
      return std::pair<bool,double>(false,0);
    }
    //
    // Use existing second order object at first pass, create temporary object
    // for subsequent passes.
    //
    std::pair<bool,double> deltaS1;
    if ( first ) {
      first = false;
      deltaS1 = theQuadraticSolutionFromStart.pathLength(object);
    }
    else {
      HelixExtrapolatorToLine2Order linearCrossing(xnew.x(),xnew.y(),xnew.z(),
                                                   pnew.x(),pnew.y(),
                                                   theCosTheta,theSinTheta,
                                                   theRho,anyDirection);
      deltaS1 = linearCrossing.pathLength(object);
    }
    if ( !deltaS1.first )  return deltaS1;
    //
    // Calculate total pathlength
    //
    dSTotal += deltaS1.second;
    PropagationDirection newDir = dSTotal>=0 ? alongMomentum : oppositeToMomentum;
    if ( propDir == anyDirection ) {
      propDir = newDir;
    }
    else {
      if ( newDir!=propDir )  return std::pair<bool,double>(false,0);
    }
    //
    // Check convergence
    //
    if ( deltaS1.second*deltaS1.second<maxDeltaS2 )  break;
    //
    // Step forward by dSTotal.
    //
    xnew = positionInDouble(dSTotal);
    pnew = directionInDouble(dSTotal);
  }
  //
  // Return result
  //
  return std::pair<bool,double>(true,dSTotal);
}

virtual std::pair<bool,double> IterativeHelixExtrapolatorToLine::pathLength ( const GlobalPoint & point ) const [inline, 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 27 of file IterativeHelixExtrapolatorToLine.h.

References genericPathLength().

Referenced by AnalyticalTrajectoryExtrapolatorToLine::propagateWithHelix(), and AnalyticalImpactPointExtrapolator::propagateWithHelix().

                                                                           {
    return genericPathLength(point);
  }

virtual std::pair<bool,double> IterativeHelixExtrapolatorToLine::pathLength ( const Line & line ) const [inline, 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 35 of file IterativeHelixExtrapolatorToLine.h.

References genericPathLength().

                                                                   {
    return genericPathLength(line);
  }

HelixLineExtrapolation::PositionType IterativeHelixExtrapolatorToLine::position ( double s ) const [virtual]

Position at pathlength s from the starting point.

Implements HelixLineExtrapolation.

Definition at line 115 of file IterativeHelixExtrapolatorToLine.cc.

References positionInDouble().

Referenced by AnalyticalTrajectoryExtrapolatorToLine::propagateWithHelix(), and AnalyticalImpactPointExtrapolator::propagateWithHelix().

                                                          {
  // use result in double precision
  return PositionType(positionInDouble(s));
}

HelixLineExtrapolation::PositionTypeDouble IterativeHelixExtrapolatorToLine::positionInDouble ( double s ) const

Position at pathlength s from the starting point.

Definition at line 124 of file IterativeHelixExtrapolatorToLine.cc.

References funct::cos(), HelixExtrapolatorToLine2Order::positionInDouble(), asciidump::s, funct::sin(), theCachedCDPhi, theCachedDPhi, theCachedS, theCachedSDPhi, theCosPhi0, theCosTheta, theQuadraticSolutionFromStart, theRho, theSinPhi0, theSinTheta, theX0, theY0, and theZ0.

Referenced by genericPathLength(), and position().

                                                                  {
  //
  // Calculate delta phi (if not already available)
  //
  if ( s!=theCachedS ) {
    theCachedS = s;
    theCachedDPhi = theCachedS*theRho*theSinTheta;
    theCachedSDPhi = sin(theCachedDPhi);
    theCachedCDPhi = cos(theCachedDPhi);
  }
  //
  // Calculate with appropriate formulation of full helix formula or with 
  //   1st order approximation.
  //
//    if ( fabs(theCachedDPhi)>1.e-1 ) {
  if ( fabs(theCachedDPhi)>1.e-4 ) {
    // "standard" helix formula
    return PositionTypeDouble(theX0+(-theSinPhi0*(1.-theCachedCDPhi)+theCosPhi0*theCachedSDPhi)/theRho,
                              theY0+(theCosPhi0*(1.-theCachedCDPhi)+theSinPhi0*theCachedSDPhi)/theRho,
                              theZ0+theCachedS*theCosTheta);
    }
//    else if ( fabs(theCachedDPhi)>theNumericalPrecision ) {
//      // full helix formula, but avoiding (1-cos(deltaPhi)) for small angles
//      return PositionTypeDouble(theX0+(-theSinPhi0*theCachedSDPhi*theCachedSDPhi/(1.+theCachedCDPhi)+
//                                   theCosPhi0*theCachedSDPhi)/theRho,
//                            theY0+(theCosPhi0*theCachedSDPhi*theCachedSDPhi/(1.+theCachedCDPhi)+
//                                   theSinPhi0*theCachedSDPhi)/theRho,
//                            theZ0+theCachedS*theCosTheta);
//    }
  else {
    // Use 1st order.
    return theQuadraticSolutionFromStart.positionInDouble(theCachedS);
  }
}