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

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

const 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

const GlobalTrajectoryParameters *	secondGTP

double	thecosPhiH0

const GlobalTrajectoryParameters *	theH

double	theh

const 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::TwoTrackMinimumDistanceHelixLine ( )

inline

Definition at line 19 of file TwoTrackMinimumDistanceHelixLine.h.

19 : theH(nullptr), theL(nullptr), themaxiter(12), pointsUpdated(false) {}

TwoTrackMinimumDistanceHelixLine::theH

const GlobalTrajectoryParameters * theH

Definition: TwoTrackMinimumDistanceHelixLine.h:44

TwoTrackMinimumDistanceHelixLine::pointsUpdated

bool pointsUpdated

Definition: TwoTrackMinimumDistanceHelixLine.h:59

TwoTrackMinimumDistanceHelixLine::theL

const GlobalTrajectoryParameters * theL

Definition: TwoTrackMinimumDistanceHelixLine.h:44

TwoTrackMinimumDistanceHelixLine::themaxiter

int themaxiter

Definition: TwoTrackMinimumDistanceHelixLine.h:53

◆ ~TwoTrackMinimumDistanceHelixLine()

TwoTrackMinimumDistanceHelixLine::~TwoTrackMinimumDistanceHelixLine ( )

inline

Definition at line 20 of file TwoTrackMinimumDistanceHelixLine.h.

20 {}

Member Function Documentation

◆ calculate()

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 105 of file TwoTrackMinimumDistanceHelixLine.cc.

References dqmiolumiharvest::j, LogDebug, M_PI, testProducerWithPsetDescEmpty_cfi::x1, and testProducerWithPsetDescEmpty_cfi::x2.

                                                                    {
   pointsUpdated = false;
   firstGTP = &theFirstGTP;
   secondGTP = &theSecondGTP;
 
   if (updateCoeffs()) {
     finalPoints();
     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) {
       finalPoints();
       return false;
     }
   }
   LogDebug("TwoTrackMinimumDistanceHelixLine") << "Number of steps exceeded. Has not converged.";
   finalPoints();
   return true;
 }

◆ finalPoints()

void TwoTrackMinimumDistanceHelixLine::finalPoints ( )

private

Definition at line 167 of file TwoTrackMinimumDistanceHelixLine.cc.

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

                                                    {
   if (pointsUpdated)
     return;
   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;
 }

◆ firstAngle()

double TwoTrackMinimumDistanceHelixLine::firstAngle ( ) const

Definition at line 139 of file TwoTrackMinimumDistanceHelixLine.cc.

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

◆ oneIteration()

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

private

Definition at line 77 of file TwoTrackMinimumDistanceHelixLine.cc.

References funct::cos(), createTree::dd, funct::derivative(), alignCSCRings::ff, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, 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;
 }

◆ pathLength()

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

Definition at line 160 of file TwoTrackMinimumDistanceHelixLine.cc.

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

◆ points()

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 153 of file TwoTrackMinimumDistanceHelixLine.cc.

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

◆ secondAngle()

double TwoTrackMinimumDistanceHelixLine::secondAngle ( ) const

Definition at line 146 of file TwoTrackMinimumDistanceHelixLine.cc.

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

◆ updateCoeffs()

bool TwoTrackMinimumDistanceHelixLine::updateCoeffs ( )

private

Definition at line 10 of file TwoTrackMinimumDistanceHelixLine.cc.

References funct::cos(), createTree::dd, alignCSCRings::ff, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, funct::sin(), mathSSE::sqrt(), X, BeamSpotPI::Y, and BeamSpotPI::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;
 }