#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
double	theh
GlobalTrajectoryParameters *	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.

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

TwoTrackMinimumDistanceHelixLine::~TwoTrackMinimumDistanceHelixLine ( ) [inline]

Definition at line 23 of file TwoTrackMinimumDistanceHelixLine.h.

{}

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

References funct::cos(), diffTreeTool::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 145 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 79 of file TwoTrackMinimumDistanceHelixLine.cc.

References funct::cos(), createTree::dd, createTree::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 166 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 157 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 151 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, createTree::ff, funct::sin(), mathSSE::sqrt(), X, and Gflash::Z.

{

  if (firstGTP->charge() == 0. && secondGTP->charge() != 0.) {
    theL= firstGTP;
    theH= secondGTP;
  } else if (firstGTP->charge() != 0. && secondGTP->charge() == 0.) {
    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;
}