#include <HelixArbitraryPlaneCrossing.h>

Inheritance diagram for HelixArbitraryPlaneCrossing:

Public Types
typedef Basic3DVector< double >	DirectionTypeDouble
typedef Basic3DVector< double >	PositionTypeDouble
Public Member Functions
virtual DirectionType	direction (double s) const
DirectionTypeDouble	directionInDouble (double s) const
	HelixArbitraryPlaneCrossing (const PositionType &point, const DirectionType &direction, const float curvature, const PropagationDirection propDir=alongMomentum)
virtual std::pair< bool, double >	pathLength (const Plane &plane)
virtual PositionType	position (double s) const
PositionTypeDouble	positionInDouble (double s) const
virtual	~HelixArbitraryPlaneCrossing ()
Private Member Functions
bool	notAtSurface (const Plane &, const PositionTypeDouble &, const float) const dso_internal
Private Attributes
double	theCachedCDPhi
double	theCachedDPhi
double	theCachedS
double	theCachedSDPhi
double	theCosPhi0
double	theCosTheta
const PropagationDirection	thePropDir
HelixArbitraryPlaneCrossing2Order	theQuadraticCrossingFromStart
const double	theRho
double	theSinPhi0
double	theSinTheta
const double	theX0
const double	theY0
const double	theZ0
Static Private Attributes
static const float	theMaxDistToPlane = 1.e-4f
static const float	theNumericalPrecision = 5.e-7f

Detailed Description

Calculates intersections of a helix with planes of any orientation.

Definition at line 10 of file HelixArbitraryPlaneCrossing.h.

Member Typedef Documentation

typedef Basic3DVector<double> HelixArbitraryPlaneCrossing::DirectionTypeDouble

Definition at line 38 of file HelixArbitraryPlaneCrossing.h.

typedef Basic3DVector<double> HelixArbitraryPlaneCrossing::PositionTypeDouble

Definition at line 37 of file HelixArbitraryPlaneCrossing.h.

Constructor & Destructor Documentation

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

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

Definition at line 9 of file HelixArbitraryPlaneCrossing.cc.

References p2, mathSSE::sqrt(), theCosPhi0, theCosTheta, theSinPhi0, theSinTheta, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

                                                                                             :
  theQuadraticCrossingFromStart(point,direction,curvature,propDir),
  theX0(point.x()),
  theY0(point.y()),
  theZ0(point.z()),
  theRho(curvature),
  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 pt2 = px*px+py*py;
  double p2 = pt2+pz*pz;
  double pI = 1./sqrt(p2);
  double ptI = 1./sqrt(pt2);
  theCosPhi0 = px*ptI;
  theSinPhi0 = py*ptI;
  theCosTheta = pz*pI;
  theSinTheta = pt2*ptI*pI;
}

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

Definition at line 19 of file HelixArbitraryPlaneCrossing.h.

{}

Member Function Documentation

HelixPlaneCrossing::DirectionType HelixArbitraryPlaneCrossing::direction ( double s ) const [virtual]

Direction at pathlength s from the starting point.

Implements HelixPlaneCrossing.

Definition at line 169 of file HelixArbitraryPlaneCrossing.cc.

References dir, directionInDouble(), Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by AnalyticalPropagator::propagateParametersOnPlane().

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

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

Direction at pathlength s from the starting point.

Definition at line 178 of file HelixArbitraryPlaneCrossing.cc.

References abs, funct::cos(), HelixArbitraryPlaneCrossing2Order::directionInDouble(), alignCSCRings::e, alignCSCRings::s, funct::sin(), theCachedCDPhi, theCachedDPhi, theCachedS, theCachedSDPhi, theCosPhi0, theCosTheta, theQuadraticCrossingFromStart, theRho, theSinPhi0, theSinTheta, and unlikely.

Referenced by direction(), and pathLength().

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

  if ( std::abs(theCachedDPhi)>1.e-4 ) {
    // full helix formula
    return DirectionTypeDouble(theCosPhi0*theCachedCDPhi-theSinPhi0*theCachedSDPhi,
                               theSinPhi0*theCachedCDPhi+theCosPhi0*theCachedSDPhi,
                               theCosTheta/theSinTheta);
  }
  else {
    // 2nd order
    return theQuadraticCrossingFromStart.directionInDouble(theCachedS);
  }
}

bool HelixArbitraryPlaneCrossing::notAtSurface	(	const Plane &	plane,
		const PositionTypeDouble &	point,
		const float	maxDist
	)		const `[inline, private]`

Iteration control: check for significant distance to plane.

Definition at line 202 of file HelixArbitraryPlaneCrossing.cc.

References abs, Plane::localZ(), insertMaterial::maxDist, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by pathLength().

                                                                           {
  float dz = plane.localZ(Surface::GlobalPoint(point.x(),point.y(),point.z()));
  return std::abs(dz)>maxDist;
}

std::pair< bool, double > HelixArbitraryPlaneCrossing::pathLength ( const Plane & plane ) [virtual]

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

Implements HelixPlaneCrossing.

Definition at line 43 of file HelixArbitraryPlaneCrossing.cc.

References alongMomentum, anyDirection, directionInDouble(), first, align_cfg::iteration, Basic3DVector< T >::mag(), notAtSurface(), oppositeToMomentum, HelixArbitraryPlaneCrossing2Order::pathLength(), GloballyPositioned< T >::position(), positionInDouble(), theCosPhi0, theCosTheta, theMaxDistToPlane, theNumericalPrecision, thePropDir, theQuadraticCrossingFromStart, theRho, theSinPhi0, theSinTheta, theX0, theY0, theZ0, unlikely, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by PathToPlane2Order::operator()().

                                                          {
  //
  // Constants used for control of convergence
  //
  const int maxIterations(100);
  //
  // maximum distance to plane (taking into account numerical precision)
  //
  float maxNumDz = theNumericalPrecision*plane.position().mag();
  float safeMaxDist = (theMaxDistToPlane>maxNumDz?theMaxDistToPlane:maxNumDz);
  //
  // 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
  //
  unsigned int iteration(maxIterations+1);
  double dSTotal(0.);
  //
  bool first = true;
  while ( notAtSurface(plane,xnew,safeMaxDist) ) {
    //
    // return empty solution vector if no convergence after maxIterations iterations
    //
    if unlikely( --iteration == 0 ) {
      edm::LogInfo("HelixArbitraryPlaneCrossing") << "pathLength : no convergence";
      return std::pair<bool,double>(false,0);
    }

    //
    // Use existing 2nd order object at first pass, create temporary object
    // for subsequent passes.
    //
    std::pair<bool,double> deltaS2;
    if unlikely( first ) {
      first = false;
      deltaS2 = theQuadraticCrossingFromStart.pathLength(plane);
    }
    else {
      HelixArbitraryPlaneCrossing2Order quadraticCrossing(xnew.x(),xnew.y(),xnew.z(),
                                                          pnew.x(),pnew.y(),
                                                          theCosTheta,theSinTheta,
                                                          theRho,
                                                          anyDirection);
      
      deltaS2 = quadraticCrossing.pathLength(plane);
    }
     
    if unlikely( !deltaS2.first )  return deltaS2;
    //
    // Calculate and sort total pathlength (max. 2 solutions)
    //
    dSTotal += deltaS2.second;
    PropagationDirection newDir = dSTotal>=0 ? alongMomentum : oppositeToMomentum;
    if ( propDir == anyDirection ) {
      propDir = newDir;
    }
    else {
      if unlikely( newDir!=propDir )  return std::pair<bool,double>(false,0);
    }
    //
    // Step forward by dSTotal.
    //
    xnew = positionInDouble(dSTotal);
    pnew = directionInDouble(dSTotal);
  }
  //
  // Return result
  //
  return std::pair<bool,double>(true,dSTotal);
}

HelixPlaneCrossing::PositionType HelixArbitraryPlaneCrossing::position ( double s ) const [virtual]

Position at pathlength s from the starting point.

Implements HelixPlaneCrossing.

Definition at line 121 of file HelixArbitraryPlaneCrossing.cc.

References pos, positionInDouble(), Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by AnalyticalPropagator::propagateParametersOnPlane().

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

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

Position at pathlength s from the starting point.

Definition at line 130 of file HelixArbitraryPlaneCrossing.cc.

References abs, funct::cos(), alignCSCRings::e, python::connectstrParser::o, HelixArbitraryPlaneCrossing2Order::positionInDouble(), alignCSCRings::s, funct::sin(), theCachedCDPhi, theCachedDPhi, theCachedS, theCachedSDPhi, theCosPhi0, theCosTheta, theQuadraticCrossingFromStart, theRho, theSinPhi0, theSinTheta, theX0, theY0, theZ0, and unlikely.

Referenced by pathLength(), and position().

                                                             {
  //
  // Calculate delta phi (if not already available)
  //
  if unlikely( s!=theCachedS ) {
    theCachedS = s;
    theCachedDPhi = theCachedS*theRho*theSinTheta;
    theCachedSDPhi = sin(theCachedDPhi);
    theCachedCDPhi = cos(theCachedDPhi);
  }
  //
  // Calculate with appropriate formulation of full helix formula or with 
  //   2nd order approximation.
  //
//    if ( fabs(theCachedDPhi)>1.e-1 ) {
  if ( std::abs(theCachedDPhi)>1.e-4 ) {
    // "standard" helix formula
    double o = 1./theRho;
    return PositionTypeDouble(theX0+(-theSinPhi0*(1.-theCachedCDPhi)+theCosPhi0*theCachedSDPhi)*o,
                              theY0+( theCosPhi0*(1.-theCachedCDPhi)+theSinPhi0*theCachedSDPhi)*o,
                              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 2nd order.
    return theQuadraticCrossingFromStart.positionInDouble(theCachedS);
  }
}