#include <HelixBarrelPlaneCrossingByCircle.h>

Inheritance diagram for HelixBarrelPlaneCrossingByCircle:

Public Member Functions
virtual DirectionType	direction (double s) const

	HelixBarrelPlaneCrossingByCircle (const PositionType &pos, const DirectionType &dir, double rho, PropagationDirection propDir=alongMomentum)

	HelixBarrelPlaneCrossingByCircle (const GlobalPoint &pos, const GlobalVector &dir, double rho, PropagationDirection propDir=alongMomentum)

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

virtual PositionType	position (double s) const

Private Types
typedef Basic2DVector< double >	Vector2D

Private Member Functions
bool	chooseSolution (const Vector2D &d1, const Vector2D &d2)

void	init ()

Private Attributes
double	theActualDir

double	theCosTheta

Vector2D	theD

double	theDmag

PropagationDirection	thePropDir

double	theRho

double	theS

double	theSinTheta

DirectionType	theStartingDir

PositionType	theStartingPos

double	theXCenter

double	theYCenter

bool	useStraightLine

Additional Inherited Members
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...

Detailed Description

Computes the crossing of a helix with a barrel plane. Exact if the magnetic field is parallel to the plane.

Definition at line 13 of file HelixBarrelPlaneCrossingByCircle.h.

Member Typedef Documentation

typedef Basic2DVector<double> HelixBarrelPlaneCrossingByCircle::Vector2D

private

Definition at line 34 of file HelixBarrelPlaneCrossingByCircle.h.

Constructor & Destructor Documentation

HelixBarrelPlaneCrossingByCircle::HelixBarrelPlaneCrossingByCircle	(	const PositionType &	pos,
		const DirectionType &	dir,
		double	rho,
		PropagationDirection	propDir = `alongMomentum`
	)

Definition at line 11 of file HelixBarrelPlaneCrossingByCircle.cc.

References init().

                                                             :
   theStartingPos( pos), theStartingDir(dir), 
   theRho(rho), thePropDir(propDir) { init();}

HelixBarrelPlaneCrossingByCircle::HelixBarrelPlaneCrossingByCircle	(	const GlobalPoint &	pos,
		const GlobalVector &	dir,
		double	rho,
		PropagationDirection	propDir = `alongMomentum`
	)

Definition at line 18 of file HelixBarrelPlaneCrossingByCircle.cc.

References init().

                                                             :
   theStartingPos( pos.basicVector()), theStartingDir(dir.basicVector()), 
   theRho(rho), thePropDir(propDir) { init();}

Member Function Documentation

bool HelixBarrelPlaneCrossingByCircle::chooseSolution	(	const Vector2D &	d1,
		const Vector2D &	d2
	)

private

Definition at line 121 of file HelixBarrelPlaneCrossingByCircle.cc.

References alongMomentum, anyDirection, Basic2DVector< T >::mag2(), mathSSE::samesign(), theActualDir, theD, thePropDir, theStartingDir, Basic2DVector< T >::x(), Basic3DVector< T >::x(), Basic2DVector< T >::y(), and Basic3DVector< T >::y().

Referenced by pathLength().

 {
   bool solved;
   double momProj1 = theStartingDir.x()*d1.x() + theStartingDir.y()*d1.y();
   double momProj2 = theStartingDir.x()*d2.x() + theStartingDir.y()*d2.y();
 
   if ( thePropDir == anyDirection ) {
     solved = true;
     if (d1.mag2()<d2.mag2()) {
       theD = d1;
       theActualDir = (momProj1 > 0) ? 1. : -1.;
     }
     else {
       theD = d2;
       theActualDir = (momProj2 > 0) ? 1. : -1.;
     }
 
   }
   else {
     double propSign = thePropDir==alongMomentum ? 1 : -1;
     if (!mathSSE::samesign(momProj1,momProj2)) {
       // if different signs return the positive one
       solved = true;
       theD         = mathSSE::samesign(momProj1,propSign) ? d1 : d2;
       theActualDir = propSign;
     }
     else if (mathSSE::samesign(momProj1,propSign)) {
       // if both positive, return the shortest
       solved = true;
       theD = (d1.mag2()<d2.mag2()) ? d1 : d2;
       theActualDir = propSign;
     }
     else solved = false;
   }
   return solved;
 }

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

virtual

Returns the direction along the helix that corresponds to path length "s" from the starting point. As for position, the direction of the crossing with a plane (if it exists!) is given by direction( pathLength( plane)).

Implements HelixPlaneCrossing.

Definition at line 183 of file HelixBarrelPlaneCrossingByCircle.cc.

References funct::cos(), phi, funct::sin(), mathSSE::sqrt(), theDmag, theRho, theS, theSinTheta, theStartingDir, tmp, useStraightLine, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by TSCPBuilderNoMaterial::createFTSatTransverseImpactPointCharged(), and TrackKinematicStatePropagator::propagateToTheTransversePCACharged().

 {
   if(useStraightLine){return theStartingDir;}
   else{
     double sinPhi, cosPhi;
     if ( s==theS) {
       double tmp = 0.5*theDmag*theRho;
       if (s < 0) tmp = -tmp;
       // protect sqrt
       sinPhi = 1.-(tmp*tmp);
       if ( sinPhi<0 )  sinPhi = 0.;
       sinPhi = 2.*tmp*sqrt(sinPhi);
       cosPhi = 1.-2.*(tmp*tmp);
     }
     else {
       double phi = s*theSinTheta*theRho;
       sinPhi = sin(phi);
       cosPhi = cos(phi);
     }
     return DirectionType(theStartingDir.x()*cosPhi-theStartingDir.y()*sinPhi,
              theStartingDir.x()*sinPhi+theStartingDir.y()*cosPhi,
              theStartingDir.z());
   }
 }

void HelixBarrelPlaneCrossingByCircle::init ( void )

private

Definition at line 24 of file HelixBarrelPlaneCrossingByCircle.cc.

References Basic3DVector< T >::mag(), python.connectstrParser::o, Basic3DVector< T >::perp(), theCosTheta, theRho, theSinTheta, theStartingDir, theStartingPos, theXCenter, theYCenter, useStraightLine, Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by HelixBarrelPlaneCrossingByCircle().

 {
   double pabsI = 1./theStartingDir.mag();
   double pt   = theStartingDir.perp();
   theCosTheta = theStartingDir.z()*pabsI;
   theSinTheta = pt*pabsI;
 
   // protect for zero curvature case
   const double sraightLineCutoff = 1.e-7;
   if (fabs(theRho) < sraightLineCutoff && 
       fabs(theRho)*theStartingPos.perp()  < sraightLineCutoff) {
     useStraightLine = true;
   }else{
     // circle parameters
     // position of center of curvature is on a line perpendicular
     // to startingDir and at a distance 1/curvature.
     double o = 1./(pt*theRho);
     theXCenter = theStartingPos.x() - theStartingDir.y()*o;
     theYCenter = theStartingPos.y() + theStartingDir.x()*o;
     useStraightLine = false;
   }
 }

std::pair< bool, double > HelixBarrelPlaneCrossingByCircle::pathLength ( const 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 48 of file HelixBarrelPlaneCrossingByCircle.cc.

References funct::A, funct::C, chooseSolution(), RealQuadEquation::first, RealQuadEquation::hasSolution, Plane::localZ(), Basic2DVector< T >::mag(), n, Plane::normalVector(), StraightLinePlaneCrossing::pathLength(), RealQuadEquation::second, theActualDir, theD, theDmag, thePropDir, theRho, theS, theSinTheta, theStartingDir, theStartingPos, theXCenter, theYCenter, useStraightLine, PV3DBase< T, PVType, FrameType >::x(), Basic3DVector< T >::x(), PV3DBase< T, PVType, FrameType >::y(), and Basic3DVector< T >::y().

Referenced by TSCPBuilderNoMaterial::createFTSatTransverseImpactPointCharged(), and TrackKinematicStatePropagator::propagateToTheTransversePCACharged().

 {
   typedef std::pair<bool,double>     ResultType;
   
   if(useStraightLine){
     // switch to straight line case
     StraightLinePlaneCrossing slc( theStartingPos, 
                    theStartingDir, thePropDir);
     return slc.pathLength( plane);
   }
   
 
   // plane parameters
   GlobalVector n = plane.normalVector();
   double distToPlane = -plane.localZ( GlobalPoint( theStartingPos));
   //    double distToPlane = (plane.position().x()-theStartingPos.x()) * n.x() +
   //                         (plane.position().y()-theStartingPos.y()) * n.y();
   double nx = n.x();  // convert to double
   double ny = n.y();  // convert to double
   double distCx = theStartingPos.x() - theXCenter;
   double distCy = theStartingPos.y() - theYCenter;
 
   double nfac, dfac;
   double A, B, C;
   bool solveForX;
   if (fabs(nx) > fabs(ny)) {
     solveForX = false;
     nfac = ny/nx;
     dfac = distToPlane/nx;
     B = distCy - nfac*distCx;  // only part of B, may have large cancelation
     C = (2.*distCx + dfac)*dfac;
   }
   else {
     solveForX = true;
     nfac = nx/ny;
     dfac = distToPlane/ny;
     B = distCx - nfac*distCy; // only part of B, may have large cancelation
     C = (2.*distCy + dfac)*dfac;
   }
   B -= nfac*dfac; B *= 2;  // the rest of B (normally small)
   A = 1.+ nfac*nfac;
 
   double dx1, dx2, dy1, dy2;
   RealQuadEquation equation(A,B,C);
   if (!equation.hasSolution) return ResultType( false, 0.);
   else {
     if (solveForX) {
       dx1 = equation.first;
       dx2 = equation.second;
       dy1 = dfac - nfac*dx1;
       dy2 = dfac - nfac*dx2;
     }
     else {
       dy1 = equation.first;
       dy2 = equation.second;
       dx1 = dfac - nfac*dy1;
       dx2 = dfac - nfac*dy2;
     }
   }
   bool solved = chooseSolution( Vector2D(dx1, dy1), Vector2D(dx2, dy2));
   if (solved) {
     theDmag = theD.mag();
     // protect asin (taking some safety margin)
     double sinAlpha = 0.5*theDmag*theRho;
     if ( sinAlpha>(1.-10*DBL_EPSILON) )  sinAlpha = 1.-10*DBL_EPSILON;
     else if ( sinAlpha<-(1.-10*DBL_EPSILON) )  sinAlpha = -(1.-10*DBL_EPSILON);
     theS = theActualDir*2./(theRho*theSinTheta) * asin( sinAlpha);
     return ResultType( true, theS);
   }
   else return ResultType( false, 0.);
 }

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

virtual

Returns the position along the helix that corresponds to path length "s" from the starting point. If s is obtained from the pathLength method the position is the destination point, i.e. the position of the crossing with a plane (if it exists!) is given by position( pathLength( plane)).

Implements HelixPlaneCrossing.

Definition at line 160 of file HelixBarrelPlaneCrossingByCircle.cc.

References funct::cos(), phi, funct::sin(), theCosTheta, theD, theRho, theS, theSinTheta, theStartingDir, theStartingPos, theXCenter, theYCenter, Basic3DVector< T >::unit(), useStraightLine, Basic2DVector< T >::x(), Basic3DVector< T >::x(), Basic2DVector< T >::y(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by TSCPBuilderNoMaterial::createFTSatTransverseImpactPointCharged(), and TrackKinematicStatePropagator::propagateToTheTransversePCACharged().

 {
   if(useStraightLine){
     return PositionType(theStartingPos+s*theStartingDir.unit());
   }else{
     if ( s==theS) {
       return PositionType( theStartingPos.x() + theD.x(),
                theStartingPos.y() + theD.y(), 
                theStartingPos.z() + s*theCosTheta);
     }
     else {
       double phi = s*theSinTheta*theRho;
       double x1Shift = theStartingPos.x() - theXCenter;
       double y1Shift = theStartingPos.y() - theYCenter;
       
       return PositionType(x1Shift*cos(phi)-y1Shift*sin(phi) + theXCenter,
               x1Shift*sin(phi)+y1Shift*cos(phi) + theYCenter,
               theStartingPos.z() + s*theCosTheta);
     }
   }
 }