TwoTrackMinimumDistanceHelixLine.cc

Go to the documentation of this file.

#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistanceHelixLine.h"
#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
#include <iomanip>

using namespace std;

bool TwoTrackMinimumDistanceHelixLine::updateCoeffs() {
  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;
}

bool TwoTrackMinimumDistanceHelixLine::oneIteration(double& thePhiH, double& fct, double& derivative) const {
  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;
}

bool TwoTrackMinimumDistanceHelixLine::calculate(const GlobalTrajectoryParameters& theFirstGTP,
                                                 const GlobalTrajectoryParameters& theSecondGTP,
                                                 const float qual) {
  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;
}

double TwoTrackMinimumDistanceHelixLine::firstAngle() const {
  if (firstGTP == theL)
    return theL->momentum().phi();
  else
    return thePhiH;
}

double TwoTrackMinimumDistanceHelixLine::secondAngle() const {
  if (secondGTP == theL)
    return theL->momentum().phi();
  else
    return thePhiH;
}

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

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

void TwoTrackMinimumDistanceHelixLine::finalPoints() {
  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;
}