File indexing completed on 2024-04-06 12:31:35
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 }
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
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
0133 if (!(theTTMDhh.calculate(sta, stb, .0001))) {
0134 points_ = theTTMDhh.points();
0135 return true;
0136 };
0137 };
0138
0139
0140 bool cairStat = theIniAlgo.calculate(sta, stb);
0141
0142 if (!cairStat) {
0143 edm::LogWarning("TwoTrackMinimumDistance") << "Computation HelixHelix::CAIR failed.";
0144 if (theModus == SlowMode) {
0145 if (!(theTTMDhh.calculate(sta, stb, .0001))) {
0146 points_ = theTTMDhh.points();
0147 return true;
0148 }
0149 };
0150
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
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_); }