#include <TwoTrackMinimumDistanceHelixLine.h>

Public Member Functions
bool	calculate (const GlobalTrajectoryParameters &, const GlobalTrajectoryParameters &, const float qual=.0001)

double	firstAngle () const

std::pair< double, double >	pathLength () const

std::pair< GlobalPoint, GlobalPoint >	points () const

double	secondAngle () const

	TwoTrackMinimumDistanceHelixLine ()

	~TwoTrackMinimumDistanceHelixLine ()

Private Member Functions
void	finalPoints () const

bool	oneIteration (double &thePhiH, double &fct, double &derivative) const

bool	updateCoeffs ()

Private Attributes
double	aa

double	baseDer

double	baseFct

double	bb

double	cc

double	dd

double	ee

double	ff

GlobalTrajectoryParameters *	firstGTP

double	helixPath

GlobalPoint	helixPoint

double	Hn

double	linePath

GlobalPoint	linePoint

double	Ln

bool	pointsUpdated

GlobalVector	posDiff

double	px

double	px2

double	py

double	py2

double	pz

double	pz2

GlobalTrajectoryParameters *	secondGTP

double	thecosPhiH0

GlobalTrajectoryParameters *	theH

double	theh

GlobalTrajectoryParameters *	theL

GlobalVector	theLp

int	themaxiter

double	thePhiH

double	thePhiH0

double	thesinPhiH0

double	thetanlambdaH

double	tL

double	X

double	Y

double	Z

Detailed Description

This is a helper class for TwoTrackMinimumDistance, for the case where one of the tracks is charged and the other not. No user should need direct access to this class. It implements a Newton method for finding the minimum distance between two tracks.

Definition at line 17 of file TwoTrackMinimumDistanceHelixLine.h.

Constructor & Destructor Documentation

TwoTrackMinimumDistanceHelixLine::TwoTrackMinimumDistanceHelixLine ( )

inline

Definition at line 21 of file TwoTrackMinimumDistanceHelixLine.h.

21 : theH(0), theL(0), themaxiter(12),

22 pointsUpdated(false){}

TwoTrackMinimumDistanceHelixLine::pointsUpdated

bool pointsUpdated

Definition: TwoTrackMinimumDistanceHelixLine.h:63

TwoTrackMinimumDistanceHelixLine::theL

GlobalTrajectoryParameters * theL

Definition: TwoTrackMinimumDistanceHelixLine.h:48

TwoTrackMinimumDistanceHelixLine::theH

GlobalTrajectoryParameters * theH

Definition: TwoTrackMinimumDistanceHelixLine.h:48

TwoTrackMinimumDistanceHelixLine::themaxiter

int themaxiter

Definition: TwoTrackMinimumDistanceHelixLine.h:57

TwoTrackMinimumDistanceHelixLine::~TwoTrackMinimumDistanceHelixLine ( )

inline

Definition at line 23 of file TwoTrackMinimumDistanceHelixLine.h.

References calculate(), firstAngle(), pathLength(), points(), and secondAngle().

23 {}

Member Function Documentation

bool TwoTrackMinimumDistanceHelixLine::calculate	(	const GlobalTrajectoryParameters &	theFirstGTP,
		const GlobalTrajectoryParameters &	theSecondGTP,
		const float	qual = `.0001`
	)

Calculates the PCA between a charged particle (helix) and a neutral particle (line). The order of the trajectories (helix-line or line-helix) is irrelevent, and will be conserved.

Definition at line 113 of file TwoTrackMinimumDistanceHelixLine.cc.

References LogDebug, and M_PI.

Referenced by ~TwoTrackMinimumDistanceHelixLine().

 {
   pointsUpdated = false;
   firstGTP  = (GlobalTrajectoryParameters *) &theFirstGTP;
   secondGTP = (GlobalTrajectoryParameters *) &theSecondGTP;
 
   if ( updateCoeffs () )
   {
     return true;
   };
 
   double fctVal, derVal, dPhiH;
   thePhiH = thePhiH0;
   
   double x1=thePhiH0-M_PI, x2=thePhiH0+M_PI;
   for (int j=1; j<=themaxiter; ++j) { 
     oneIteration(thePhiH, fctVal, derVal);
     dPhiH=fctVal/derVal;
     thePhiH -= dPhiH;
     if ((x1-thePhiH)*(thePhiH-x2) < 0.0) {
       LogDebug ("TwoTrackMinimumDistanceHelixLine")
         << "Jumped out of brackets in root finding. Will be moved closer.";
       thePhiH += (dPhiH*0.8);
     }
     if (fabs(dPhiH) < qual) {return false;}
   }
   LogDebug ("TwoTrackMinimumDistanceHelixLine")
     <<"Number of steps exceeded. Has not converged.";
   return true;
 }

void TwoTrackMinimumDistanceHelixLine::finalPoints ( ) const

private

Definition at line 175 of file TwoTrackMinimumDistanceHelixLine.cc.

References funct::cos(), mps_update::diff, PV3DBase< T, PVType, FrameType >::mag(), funct::sin(), and mathSSE::sqrt().

 {
   helixPoint = GlobalPoint (
       theH->position().x() + theh * ( sin ( thePhiH) - thesinPhiH0 ),
       theH->position().y() + theh * ( - cos ( thePhiH) + thecosPhiH0 ),
       theH->position().z() + theh * ( thetanlambdaH * ( thePhiH- thePhiH0 ))
   );
   helixPath = ( thePhiH- thePhiH0 ) * (theh*sqrt(1+thetanlambdaH*thetanlambdaH)) ;
 
   GlobalVector diff((theL->position() -helixPoint).basicVector());
   tL = ( - diff.dot(theLp)) / (Ln*Ln);
   linePoint = GlobalPoint (
       theL->position().x() + tL * px,
       theL->position().y() + tL * py,
       theL->position().z() + tL * pz );
   linePath = tL * theLp.mag();
   pointsUpdated = true;
 }

double TwoTrackMinimumDistanceHelixLine::firstAngle ( ) const

Definition at line 146 of file TwoTrackMinimumDistanceHelixLine.cc.

Referenced by ~TwoTrackMinimumDistanceHelixLine().

 {
   if (firstGTP==theL) return theL->momentum().phi();
   else return thePhiH;
 }

bool TwoTrackMinimumDistanceHelixLine::oneIteration	(	double &	thePhiH,
		double &	fct,
		double &	derivative
	)		const

private

Definition at line 80 of file TwoTrackMinimumDistanceHelixLine.cc.

References funct::cos(), createTree::dd, alignCSCRings::ff, and funct::sin().

 {
 
   double thesinPhiH = sin(thePhiH);
   double thecosPhiH = cos(thePhiH);
 
     // Fonction of which the root is to be found:
 
   fct = baseFct;
   fct -= ff*(thePhiH - thePhiH0);
   fct += thecosPhiH * aa;
   fct += thesinPhiH * bb;
   fct += cc*(thePhiH - thePhiH0)*(px * thecosPhiH + py * thesinPhiH);
   fct += cc * (px * (thesinPhiH - thesinPhiH0) - py * (thecosPhiH - thecosPhiH0));
   fct += dd * (thesinPhiH* (thesinPhiH - thesinPhiH0) - 
         thecosPhiH*(thecosPhiH - thecosPhiH0));
   fct += ee * thecosPhiH * thesinPhiH;
 
       // Its derivative:
 
   derivative = baseDer;
   derivative += - thesinPhiH * aa;
   derivative += thecosPhiH * bb;
   derivative += cc*(thePhiH - thePhiH0)*(py * thecosPhiH - px * thesinPhiH);
   derivative += 2* cc*(px * thecosPhiH + py * thesinPhiH);
   derivative += dd *(4 * thecosPhiH * thesinPhiH - thecosPhiH * thesinPhiH0 - 
             thesinPhiH * thecosPhiH0);
   derivative += ee * (thecosPhiH*thecosPhiH-thesinPhiH*thesinPhiH);
 
   return false;
 }

pair< double, double > TwoTrackMinimumDistanceHelixLine::pathLength ( ) const

Definition at line 167 of file TwoTrackMinimumDistanceHelixLine.cc.

Referenced by ~TwoTrackMinimumDistanceHelixLine().

 {
   if (!pointsUpdated)finalPoints();
   if (firstGTP==theL) 
     return pair<double, double> (linePath, helixPath);
   else return pair<double, double> (helixPath, linePath);
 }

pair< GlobalPoint, GlobalPoint > TwoTrackMinimumDistanceHelixLine::points ( ) const

Returns the PCA's on the two trajectories. The first point lies on the first trajectory, the second point on the second trajectory.

Definition at line 158 of file TwoTrackMinimumDistanceHelixLine.cc.

Referenced by ~TwoTrackMinimumDistanceHelixLine().

 {
   if (!pointsUpdated)finalPoints();
   if (firstGTP==theL) 
     return pair<GlobalPoint, GlobalPoint> (linePoint, helixPoint);
   else return pair<GlobalPoint, GlobalPoint> (helixPoint, linePoint);
 }

double TwoTrackMinimumDistanceHelixLine::secondAngle ( ) const

Definition at line 152 of file TwoTrackMinimumDistanceHelixLine.cc.

Referenced by ~TwoTrackMinimumDistanceHelixLine().

 {
   if (secondGTP==theL) return theL->momentum().phi();
   else return thePhiH;
 }

bool TwoTrackMinimumDistanceHelixLine::updateCoeffs ( )

private

Definition at line 10 of file TwoTrackMinimumDistanceHelixLine.cc.

References funct::cos(), createTree::dd, alignCSCRings::ff, funct::sin(), mathSSE::sqrt(), X, DOFs::Y, and DOFs::Z.

 {
   bool isFirstALine = firstGTP->charge() == 0. || firstGTP->magneticField().inTesla(firstGTP->position()).z() == 0.;
   bool isSecondALine = secondGTP->charge() == 0. || secondGTP->magneticField().inTesla(secondGTP->position()).z() == 0.;
   if (isFirstALine && !isSecondALine ) {
     theL= firstGTP;
     theH= secondGTP;
   } else if (!isFirstALine && isSecondALine) {
     theH= firstGTP;
     theL= secondGTP;
   } else {
     edm::LogWarning ("TwoTrackMinimumDistanceHelixLine")
       << "Error in track charge: "
       << "One of the tracks has to be charged, and the other not." << endl
       << "Track Charges: "<<firstGTP->charge() << " and " <<secondGTP->charge();
     return true;
   }
 
   Hn = theH->momentum().mag();
   Ln = theL->momentum().mag();
 
   if ( Hn == 0. || Ln == 0. )
   {
     edm::LogWarning ("TwoTrackMinimumDistanceHelixLine")
       << "Momentum of input trajectory is zero.";
     return true;
   };
 
   GlobalPoint lOrig = theL->position();
   GlobalPoint hOrig = theH->position();
   posDiff = GlobalVector((lOrig - hOrig).basicVector());
   X = posDiff.x();
   Y = posDiff.y();
   Z = posDiff.z();
   theLp = theL->momentum();
   px = theLp.x(); px2 = px*px;
   py = theLp.y(); py2 = py*py;
   pz = theLp.z(); pz2 = pz*pz;
   
   const double Bc2kH = theH->magneticField().inTesla(hOrig).z() * 2.99792458e-3;
 //   MagneticField::inInverseGeV ( hOrig ).z();
 
   if ( Bc2kH == 0. )
   {
     edm::LogWarning ("TwoTrackMinimumDistanceHelixLine")
       << "Magnetic field at point " << hOrig << " is zero.";
     return true;
   };
 
   theh= - Hn / (theH->charge() * Bc2kH ) *
      sqrt( 1 - ( ( (theH->momentum().z()*theH->momentum().z()) / (Hn*Hn) )));
 
   thetanlambdaH = - theH->momentum().z() / ( theH->charge() * Bc2kH * theh);
 
   thePhiH0 = theH->momentum().phi();
   thesinPhiH0= sin(thePhiH0);
   thecosPhiH0= cos(thePhiH0);
 
   aa = (X + theh*thesinPhiH0)*(py2 + pz2) - px*(py*Y + pz*Z);
   bb = (Y - theh*thecosPhiH0)*(px2 + pz2) - py*(px*X + pz*Z);
   cc = pz*theh*thetanlambdaH;
   dd = theh* px *py;
   ee = theh*(px2 - py2);
   ff = (px2 + py2)*theh*thetanlambdaH*thetanlambdaH;
 
   baseFct = thetanlambdaH * (Z*(px2+py2) - pz*(px*X + py*Y));
   baseDer = - ff;
   return false;
 }