#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 ()

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

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

19 {}

Member Function Documentation

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

virtual

Direction at pathlength s from the starting point.

Implements HelixPlaneCrossing.

Definition at line 200 of file HelixArbitraryPlaneCrossing.cc.

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

Referenced by PixelForwardLayer::computeCrossings(), and PixelForwardLayerPhase1::computeCrossings().

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

References funct::abs(), HelixArbitraryPlaneCrossing2Order::directionInDouble(), alignCSCRings::e, alignCSCRings::s, 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 ) {  // very very unlikely!
     theCachedS = s;
     theCachedDPhi = theCachedS*theRho*theSinTheta;
     vdt::fast_sincos(theCachedDPhi,theCachedSDPhi,theCachedCDPhi);
   }
 
   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 232 of file HelixArbitraryPlaneCrossing.cc.

References funct::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 67 of file HelixArbitraryPlaneCrossing.cc.

References funct::abs(), alongMomentum, anyDirection, constexpr, directionInDouble(), align_cfg::iteration, Plane::localZ(), LogDebug, Basic3DVector< T >::mag(), maxiter, notAtSurface(), oppositeToMomentum, HelixArbitraryPlaneCrossing2Order::pathLength(), GloballyPositioned< T >::position(), positionInDouble(), theCosTheta, theMaxDistToPlane, theNumericalPrecision, thePropDir, theQuadraticCrossingFromStart, theRho, theSinTheta, theX0, theY0, theZ0, and unlikely.

Referenced by PixelForwardLayer::computeCrossings(), PixelForwardLayerPhase1::computeCrossings(), and PathToPlane2Order::operator()().

                                                           {
   //
   // Constants used for control of convergence
   //
   constexpr int maxIterations = 20;
   //
   // 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 
   //
   
   float dz = plane.localZ(Surface::GlobalPoint(theX0,theY0,theZ0));
   if (std::abs(dz)<safeMaxDist) return std::make_pair(true,0.);
 
   bool notFail;
   double dSTotal;
   // Use existing 2nd order object at first pass
   std::tie(notFail,dSTotal) = theQuadraticCrossingFromStart.pathLength(plane);
   if unlikely(!notFail) return std::make_pair(notFail,dSTotal);
   auto  xnew = positionInDouble(dSTotal);
 
   auto propDir = thePropDir;
   auto newDir = dSTotal>=0 ? alongMomentum : oppositeToMomentum;
   if ( propDir == anyDirection ) {
       propDir = newDir;
   }  else {
      if unlikely( newDir!=propDir )  return std::pair<bool,double>(false,0);
   }
  
 
   //
   // Prepare iterations: count and total pathlength
   //
   auto iteration = maxIterations;
   while ( notAtSurface(plane,xnew,safeMaxDist) ) {
     //
     // return empty solution vector if no convergence after maxIterations iterations
     //
     if unlikely( --iteration == 0 ) {
       LogDebug("HelixArbitraryPlaneCrossing") << "pathLength : no convergence";
       return std::pair<bool,double>(false,0);
     }
 
     //
     // create temporary object for subsequent passes.
     auto  pnew = directionInDouble(dSTotal);
     HelixArbitraryPlaneCrossing2Order quadraticCrossing(xnew.x(),xnew.y(),xnew.z(),
                               pnew.x(),pnew.y(),
                               theCosTheta,theSinTheta,
                               theRho,
                               anyDirection);
       
     auto  deltaS2 = quadraticCrossing.pathLength(plane);
    
      
     if unlikely( !deltaS2.first )  return deltaS2;
     //
     // Calculate and sort total pathlength (max. 2 solutions)
     //
     dSTotal += deltaS2.second;
     auto 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);
   }
   //
   // Return result
   //
   maxiter(iteration);
 
   return std::make_pair(true,dSTotal);
 }

HelixPlaneCrossing::PositionType HelixArbitraryPlaneCrossing::position ( double s ) const

virtual

Position at pathlength s from the starting point.

Implements HelixPlaneCrossing.

Definition at line 153 of file HelixArbitraryPlaneCrossing.cc.

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

Referenced by PixelForwardLayer::computeCrossings(), and PixelForwardLayerPhase1::computeCrossings().

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

References funct::abs(), alignCSCRings::e, connectstrParser::o, HelixArbitraryPlaneCrossing2Order::positionInDouble(), alignCSCRings::s, 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;
     vdt::fast_sincos(theCachedDPhi,theCachedSDPhi,theCachedCDPhi);
   }
   //
   // 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);
   }
 }