PatternTools/src/TwoTrackMinimumDistanceLineLine.cc

0001 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistanceLineLine.h"
0002 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0003 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "MagneticField/Engine/interface/MagneticField.h"
0006
0007 bool TwoTrackMinimumDistanceLineLine::calculate(const GlobalTrajectoryParameters& theG,
0008                                                 const GlobalTrajectoryParameters& theH) {
0009   GlobalPoint gOrig = theG.position();
0010   GlobalPoint hOrig = theH.position();
0011   bool isLineG = theG.magneticField().inTesla(gOrig).z() == 0 || theG.charge() == 0;
0012   bool isLineH = theH.magneticField().inTesla(hOrig).z() == 0 || theH.charge() == 0;
0013   if (!(isLineG && isLineH)) {
0014     edm::LogWarning("TwoTrackMinimumDistanceLineLine")
0015         << "charge of input track is not zero or field non zero"
0016         << "\n positions: " << gOrig << " , " << hOrig << "\n Bz fields: " << theG.magneticField().inTesla(gOrig).z()
0017         << " , " << theH.magneticField().inTesla(hOrig).z() << "\n charges: " << theG.charge() << " , "
0018         << theH.charge();
0019     return true;
0020   };
0021
0022   GlobalVector gVec = theG.momentum();
0023   const double gMag = theG.momentum().mag();
0024   double gMag2 = gMag * gMag;
0025   GlobalVector hVec = theH.momentum();
0026   const double hMag = theH.momentum().mag();
0027   double hMag2 = hMag * hMag;
0028
0029   if (gMag == 0. || hMag == 0.) {
0030     edm::LogWarning("TwoTrackMinimumDistanceLineLine") << "momentum of input trajectory is zero.";
0031     return true;
0032   };
0033
0034   phiG = theG.momentum().phi();
0035   phiH = theH.momentum().phi();
0036
0037   double gVec_Dot_hVec = gVec.dot(hVec);
0038   double norm = gVec_Dot_hVec * gVec_Dot_hVec - gMag2 * hMag2;
0039
0040   if (norm == 0) {
0041     edm::LogWarning("TwoTrackMinimumDistanceLineLine") << "Tracks are parallel.";
0042     return true;
0043   }
0044
0045   GlobalVector posDiff = gOrig - hOrig;
0046
0047   double tG = (posDiff.dot(gVec) * hMag2 - gVec_Dot_hVec * posDiff.dot(hVec)) / norm;
0048   double tH = (gVec_Dot_hVec * posDiff.dot(gVec) - posDiff.dot(hVec) * gMag2) / norm;
0049
0050   gPos = GlobalPoint(gOrig + tG * gVec);
0051   hPos = GlobalPoint(hOrig + tH * hVec);
0052   //   cout << "TwoTrackMinimumDistanceLineLine "<<gPos<<hPos<<endl;
0053
0054   pathG = tG * gMag;
0055   pathH = tH * hMag;
0056   return false;
0057 }
0058
0059 std::pair<GlobalPoint, GlobalPoint> TwoTrackMinimumDistanceLineLine::points() const {
0060   return std::pair<GlobalPoint, GlobalPoint>(gPos, hPos);
0061 }
0062
0063 std::pair<double, double> TwoTrackMinimumDistanceLineLine::pathLength() const {
0064   return std::pair<double, double>(pathG, pathH);
0065 }