PatternTools/src/TwoTrackMinimumDistanceHelixLine.cc

0001 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistanceHelixLine.h"
0002 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0003 #include "MagneticField/Engine/interface/MagneticField.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include <iostream>
0006 #include <iomanip>
0007
0008 using namespace std;
0009
0010 bool TwoTrackMinimumDistanceHelixLine::updateCoeffs() {
0011   bool isFirstALine = firstGTP->charge() == 0. || firstGTP->magneticField().inTesla(firstGTP->position()).z() == 0.;
0012   bool isSecondALine = secondGTP->charge() == 0. || secondGTP->magneticField().inTesla(secondGTP->position()).z() == 0.;
0013   if (isFirstALine && !isSecondALine) {
0014     theL = firstGTP;
0015     theH = secondGTP;
0016   } else if (!isFirstALine && isSecondALine) {
0017     theH = firstGTP;
0018     theL = secondGTP;
0019   } else {
0020     edm::LogWarning("TwoTrackMinimumDistanceHelixLine")
0021         << "Error in track charge: "
0022         << "One of the tracks has to be charged, and the other not." << endl
0023         << "Track Charges: " << firstGTP->charge() << " and " << secondGTP->charge();
0024     return true;
0025   }
0026
0027   Hn = theH->momentum().mag();
0028   Ln = theL->momentum().mag();
0029
0030   if (Hn == 0. || Ln == 0.) {
0031     edm::LogWarning("TwoTrackMinimumDistanceHelixLine") << "Momentum of input trajectory is zero.";
0032     return true;
0033   };
0034
0035   GlobalPoint lOrig = theL->position();
0036   GlobalPoint hOrig = theH->position();
0037   posDiff = GlobalVector((lOrig - hOrig).basicVector());
0038   X = posDiff.x();
0039   Y = posDiff.y();
0040   Z = posDiff.z();
0041   theLp = theL->momentum();
0042   px = theLp.x();
0043   px2 = px * px;
0044   py = theLp.y();
0045   py2 = py * py;
0046   pz = theLp.z();
0047   pz2 = pz * pz;
0048
0049   const double Bc2kH = theH->magneticField().inTesla(hOrig).z() * 2.99792458e-3;
0050   //   MagneticField::inInverseGeV ( hOrig ).z();
0051
0052   if (Bc2kH == 0.) {
0053     edm::LogWarning("TwoTrackMinimumDistanceHelixLine") << "Magnetic field at point " << hOrig << " is zero.";
0054     return true;
0055   };
0056
0057   theh = -Hn / (theH->charge() * Bc2kH) * sqrt(1 - (((theH->momentum().z() * theH->momentum().z()) / (Hn * Hn))));
0058
0059   thetanlambdaH = -theH->momentum().z() / (theH->charge() * Bc2kH * theh);
0060
0061   thePhiH0 = theH->momentum().phi();
0062   thesinPhiH0 = sin(thePhiH0);
0063   thecosPhiH0 = cos(thePhiH0);
0064
0065   aa = (X + theh * thesinPhiH0) * (py2 + pz2) - px * (py * Y + pz * Z);
0066   bb = (Y - theh * thecosPhiH0) * (px2 + pz2) - py * (px * X + pz * Z);
0067   cc = pz * theh * thetanlambdaH;
0068   dd = theh * px * py;
0069   ee = theh * (px2 - py2);
0070   ff = (px2 + py2) * theh * thetanlambdaH * thetanlambdaH;
0071
0072   baseFct = thetanlambdaH * (Z * (px2 + py2) - pz * (px * X + py * Y));
0073   baseDer = -ff;
0074   return false;
0075 }
0076
0077 bool TwoTrackMinimumDistanceHelixLine::oneIteration(double& thePhiH, double& fct, double& derivative) const {
0078   double thesinPhiH = sin(thePhiH);
0079   double thecosPhiH = cos(thePhiH);
0080
0081   // Fonction of which the root is to be found:
0082
0083   fct = baseFct;
0084   fct -= ff * (thePhiH - thePhiH0);
0085   fct += thecosPhiH * aa;
0086   fct += thesinPhiH * bb;
0087   fct += cc * (thePhiH - thePhiH0) * (px * thecosPhiH + py * thesinPhiH);
0088   fct += cc * (px * (thesinPhiH - thesinPhiH0) - py * (thecosPhiH - thecosPhiH0));
0089   fct += dd * (thesinPhiH * (thesinPhiH - thesinPhiH0) - thecosPhiH * (thecosPhiH - thecosPhiH0));
0090   fct += ee * thecosPhiH * thesinPhiH;
0091
0092   // Its derivative:
0093
0094   derivative = baseDer;
0095   derivative += -thesinPhiH * aa;
0096   derivative += thecosPhiH * bb;
0097   derivative += cc * (thePhiH - thePhiH0) * (py * thecosPhiH - px * thesinPhiH);
0098   derivative += 2 * cc * (px * thecosPhiH + py * thesinPhiH);
0099   derivative += dd * (4 * thecosPhiH * thesinPhiH - thecosPhiH * thesinPhiH0 - thesinPhiH * thecosPhiH0);
0100   derivative += ee * (thecosPhiH * thecosPhiH - thesinPhiH * thesinPhiH);
0101
0102   return false;
0103 }
0104
0105 bool TwoTrackMinimumDistanceHelixLine::calculate(const GlobalTrajectoryParameters& theFirstGTP,
0106                                                  const GlobalTrajectoryParameters& theSecondGTP,
0107                                                  const float qual) {
0108   pointsUpdated = false;
0109   firstGTP = &theFirstGTP;
0110   secondGTP = &theSecondGTP;
0111
0112   if (updateCoeffs()) {
0113     finalPoints();
0114     return true;
0115   };
0116
0117   double fctVal, derVal, dPhiH;
0118   thePhiH = thePhiH0;
0119
0120   double x1 = thePhiH0 - M_PI, x2 = thePhiH0 + M_PI;
0121   for (int j = 1; j <= themaxiter; ++j) {
0122     oneIteration(thePhiH, fctVal, derVal);
0123     dPhiH = fctVal / derVal;
0124     thePhiH -= dPhiH;
0125     if ((x1 - thePhiH) * (thePhiH - x2) < 0.0) {
0126       LogDebug("TwoTrackMinimumDistanceHelixLine") << "Jumped out of brackets in root finding. Will be moved closer.";
0127       thePhiH += (dPhiH * 0.8);
0128     }
0129     if (fabs(dPhiH) < qual) {
0130       finalPoints();
0131       return false;
0132     }
0133   }
0134   LogDebug("TwoTrackMinimumDistanceHelixLine") << "Number of steps exceeded. Has not converged.";
0135   finalPoints();
0136   return true;
0137 }
0138
0139 double TwoTrackMinimumDistanceHelixLine::firstAngle() const {
0140   if (firstGTP == theL)
0141     return theL->momentum().phi();
0142   else
0143     return thePhiH;
0144 }
0145
0146 double TwoTrackMinimumDistanceHelixLine::secondAngle() const {
0147   if (secondGTP == theL)
0148     return theL->momentum().phi();
0149   else
0150     return thePhiH;
0151 }
0152
0153 pair<GlobalPoint, GlobalPoint> TwoTrackMinimumDistanceHelixLine::points() const {
0154   if (firstGTP == theL)
0155     return pair<GlobalPoint, GlobalPoint>(linePoint, helixPoint);
0156   else
0157     return pair<GlobalPoint, GlobalPoint>(helixPoint, linePoint);
0158 }
0159
0160 pair<double, double> TwoTrackMinimumDistanceHelixLine::pathLength() const {
0161   if (firstGTP == theL)
0162     return pair<double, double>(linePath, helixPath);
0163   else
0164     return pair<double, double>(helixPath, linePath);
0165 }
0166
0167 void TwoTrackMinimumDistanceHelixLine::finalPoints() {
0168   if (pointsUpdated)
0169     return;
0170   helixPoint = GlobalPoint(theH->position().x() + theh * (sin(thePhiH) - thesinPhiH0),
0171                            theH->position().y() + theh * (-cos(thePhiH) + thecosPhiH0),
0172                            theH->position().z() + theh * (thetanlambdaH * (thePhiH - thePhiH0)));
0173   helixPath = (thePhiH - thePhiH0) * (theh * sqrt(1 + thetanlambdaH * thetanlambdaH));
0174
0175   GlobalVector diff((theL->position() - helixPoint).basicVector());
0176   tL = (-diff.dot(theLp)) / (Ln * Ln);
0177   linePoint =
0178       GlobalPoint(theL->position().x() + tL * px, theL->position().y() + tL * py, theL->position().z() + tL * pz);
0179   linePath = tL * theLp.mag();
0180   pointsUpdated = true;
0181 }