#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.

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 j, LogDebug, and M_PI.

 {
   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.

 {
   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.

 {
   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.

 {
   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.

 {
   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, and Gflash::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;
 }