#include <HelixArbitraryPlaneCrossing.h>

Inheritance diagram for HelixArbitraryPlaneCrossing:

Public Types
typedef Basic3DVector< double >	DirectionTypeDouble

typedef Basic3DVector< double >	PositionTypeDouble

Public Types inherited from HelixPlaneCrossing
typedef Basic3DVector< float >	DirectionType

typedef Basic3DVector< float >	PositionType
	the helix is passed to the constructor and does not appear in the interface More...

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

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 9 of file HelixArbitraryPlaneCrossing.h.

Member Typedef Documentation

typedef Basic3DVector<double> HelixArbitraryPlaneCrossing::DirectionTypeDouble

Definition at line 37 of file HelixArbitraryPlaneCrossing.h.

typedef Basic3DVector<double> HelixArbitraryPlaneCrossing::PositionTypeDouble

Definition at line 36 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 8 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 ( )

inlinevirtual

Definition at line 18 of file HelixArbitraryPlaneCrossing.h.

18 {}

Member Function Documentation

HelixPlaneCrossing::DirectionType HelixArbitraryPlaneCrossing::direction ( double s ) const

virtual

Direction at pathlength s from the starting point.

Implements HelixPlaneCrossing.

Definition at line 168 of file HelixArbitraryPlaneCrossing.cc.

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

                                                       {
   // 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 177 of file HelixArbitraryPlaneCrossing.cc.

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

Referenced by direction(), and pathLength().

                                                               {
   //
   // Calculate delta phi (if not already available)
   //
   if ( 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

inlineprivate

Iteration control: check for significant distance to plane.

Definition at line 201 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 42 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, 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
   //
   int iteration(maxIterations);
   double dSTotal(0.);
   //
   bool first = true;
   while ( notAtSurface(plane,xnew,safeMaxDist) ) {
     //
     // return empty solution vector if no convergence after maxIterations iterations
     //
     if ( --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 ( 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 ( !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 ( 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 120 of file HelixArbitraryPlaneCrossing.cc.

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

                                                      {
   // 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 129 of file HelixArbitraryPlaneCrossing.cc.

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

Referenced by pathLength(), 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 
   //   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);
   }
 }