PatternTools/src/TwoTrackMinimumDistance.cc

0001 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007
0008 using namespace std;
0009
0010 namespace {
0011   inline GlobalPoint mean(pair<GlobalPoint, GlobalPoint> pr) {
0012     return GlobalPoint(0.5 * (pr.first.basicVector() + pr.second.basicVector()));
0013   }
0014
0015   inline double dist(pair<GlobalPoint, GlobalPoint> pr) { return (pr.first - pr.second).mag(); }
0016 }  // namespace
0017
0018 double TwoTrackMinimumDistance::firstAngle() const {
0019   if (!status_)
0020     throw cms::Exception(
0021         "TrackingTools/PatternTools",
0022         "TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
0023   switch (theCharge) {
0024     case (hh):
0025       return theTTMDhh.firstAngle();
0026       break;
0027     case (hl):
0028       return theTTMDhl.firstAngle();
0029       break;
0030     case (ll):
0031       return theTTMDll.firstAngle();
0032       break;
0033   }
0034   return 0;
0035 }
0036
0037 double TwoTrackMinimumDistance::secondAngle() const {
0038   if (!status_)
0039     throw cms::Exception(
0040         "TrackingTools/PatternTools",
0041         "TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
0042   switch (theCharge) {
0043     case (hh):
0044       return theTTMDhh.secondAngle();
0045       break;
0046     case (hl):
0047       return theTTMDhl.secondAngle();
0048       break;
0049     case (ll):
0050       return theTTMDll.secondAngle();
0051       break;
0052   }
0053   return 0;
0054 }
0055
0056 pair<double, double> TwoTrackMinimumDistance::pathLength() const {
0057   if (!status_)
0058     throw cms::Exception(
0059         "TrackingTools/PatternTools",
0060         "TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
0061   switch (theCharge) {
0062     case (hh):
0063       return theTTMDhh.pathLength();
0064       break;
0065     case (hl):
0066       return theTTMDhl.pathLength();
0067       break;
0068     case (ll):
0069       return theTTMDll.pathLength();
0070       break;
0071   }
0072   return std::pair<double, double>(0, 0);
0073 }
0074
0075 pair<GlobalPoint, GlobalPoint> TwoTrackMinimumDistance::points() const {
0076   if (!status_)
0077     throw cms::Exception(
0078         "TrackingTools/PatternTools",
0079         "TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
0080   return points_;
0081 }
0082
0083 bool TwoTrackMinimumDistance::calculate(const TrajectoryStateOnSurface& sta, const TrajectoryStateOnSurface& stb) {
0084   return calculate(sta.globalParameters(), stb.globalParameters());
0085 }
0086
0087 bool TwoTrackMinimumDistance::calculate(const FreeTrajectoryState& sta, const FreeTrajectoryState& stb) {
0088   //  pair<GlobalPoint, GlobalPoint> ret  = theIniAlgo.points ( sta, stb );
0089   return calculate(sta.parameters(), stb.parameters());
0090 }
0091
0092 bool TwoTrackMinimumDistance::calculate(const GlobalTrajectoryParameters& sta, const GlobalTrajectoryParameters& stb) {
0093   bool isHelixA = (sta.magneticField().inTesla(sta.position()).z() != 0.) && sta.charge() != 0.;
0094   bool isHelixB = (stb.magneticField().inTesla(stb.position()).z() != 0.) && stb.charge() != 0.;
0095   if (!isHelixA && !isHelixB) {
0096     status_ = pointsLineLine(sta, stb);
0097   } else if (isHelixA && isHelixB) {
0098     status_ = pointsHelixHelix(sta, stb);
0099   } else {
0100     status_ = pointsHelixLine(sta, stb);
0101   }
0102   return status_;
0103 }
0104
0105 bool TwoTrackMinimumDistance::pointsLineLine(const GlobalTrajectoryParameters& sta,
0106                                              const GlobalTrajectoryParameters& stb) {
0107   theCharge = ll;
0108   if (theTTMDll.calculate(sta, stb))
0109     return false;
0110   points_ = theTTMDll.points();
0111   return true;
0112 }
0113
0114 bool TwoTrackMinimumDistance::pointsHelixLine(const GlobalTrajectoryParameters& sta,
0115                                               const GlobalTrajectoryParameters& stb) {
0116   theCharge = hl;
0117   if (theTTMDhl.calculate(sta, stb, 0.000001))
0118     return false;
0119   points_ = theTTMDhl.points();
0120   return true;
0121 }
0122
0123 bool TwoTrackMinimumDistance::pointsHelixHelix(const GlobalTrajectoryParameters& sta,
0124                                                const GlobalTrajectoryParameters& stb) {
0125   if ((sta.position() - stb.position()).mag2() < 1e-7f && (sta.momentum() - stb.momentum()).mag2() < 1e-7f &&
0126       sta.charge() == stb.charge()) {
0127     edm::LogWarning("TwoTrackMinimumDistance") << "comparing track with itself!";
0128   };
0129
0130   theCharge = hh;
0131   if (theModus == FastMode) {
0132     // first we try directly - in FastMode only ...
0133     if (!(theTTMDhh.calculate(sta, stb, .0001))) {
0134       points_ = theTTMDhh.points();
0135       return true;
0136     };
0137   };
0138
0139   // okay. did not work. so we use CAIR, and then TTMD again.
0140   bool cairStat = theIniAlgo.calculate(sta, stb);
0141
0142   if (!cairStat) {  // yes. CAIR may fail.
0143     edm::LogWarning("TwoTrackMinimumDistance") << "Computation HelixHelix::CAIR failed.";
0144     if (theModus == SlowMode) {  // we can still try ttmd here.
0145       if (!(theTTMDhh.calculate(sta, stb, .0001))) {
0146         points_ = theTTMDhh.points();
0147         return true;
0148       }
0149     };
0150     // we can try with more sloppy settings
0151     if (!(theTTMDhh.calculate(sta, stb, .1))) {
0152       points_ = theTTMDhh.points();
0153       return true;
0154     }
0155     return false;
0156     edm::LogWarning("TwoTrackMinimumDistance") << "TwoTrackMinimumDistanceHelixHelix failed";
0157   };
0158
0159   pair<GlobalTrajectoryParameters, GlobalTrajectoryParameters> ini = theIniAlgo.trajectoryParameters();
0160
0161   pair<GlobalPoint, GlobalPoint> inip(ini.first.position(), ini.second.position());
0162   bool isFirstALine = ini.first.charge() == 0. || ini.first.magneticField().inTesla(ini.first.position()).z() == 0.;
0163   bool isSecondALine = ini.second.charge() == 0. || ini.second.magneticField().inTesla(ini.second.position()).z() == 0.;
0164   bool gotDist = false;
0165   if (!isFirstALine && !isSecondALine)
0166     gotDist = theTTMDhh.calculate(ini.first, ini.second, .0001);
0167   else if (isFirstALine && isSecondALine)
0168     gotDist = theTTMDll.calculate(ini.first, ini.second);
0169   else
0170     gotDist = theTTMDhl.calculate(ini.first, ini.second, .0001);
0171   if (gotDist) {
0172     points_ = inip;
0173   } else {
0174     if (!isFirstALine && !isSecondALine)
0175       points_ = theTTMDhh.points();
0176     else if (isFirstALine && isSecondALine)
0177       points_ = theTTMDll.points();
0178     else
0179       points_ = theTTMDhl.points();
0180     // if we are still worse than CAIR, we use CAIR results.
0181     if (dist(points_) > dist(inip))
0182       points_ = inip;
0183   };
0184   return true;
0185 }
0186
0187 GlobalPoint TwoTrackMinimumDistance::crossingPoint() const { return mean(points_); }
0188
0189 float TwoTrackMinimumDistance::distance() const { return dist(points_); }