00001 #include "MuonAnalysis/MuonAssociators/interface/MatcherUsingTracksAlgorithm.h"
00002
00003 #include "DataFormats/Math/interface/deltaR.h"
00004 #include "DataFormats/Math/interface/deltaPhi.h"
00005
00006 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008
00009 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00011 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00013 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00014 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
00015
00016 #include "FWCore/Utilities/interface/Exception.h"
00017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00018
00019
00020 MatcherUsingTracksAlgorithm::MatcherUsingTracksAlgorithm(const edm::ParameterSet & iConfig) :
00021 whichTrack1_(None), whichTrack2_(None),
00022 whichState1_(AtVertex), whichState2_(AtVertex),
00023 srcCut_(iConfig.existsAs<std::string>("srcPreselection") ? iConfig.getParameter<std::string>("srcPreselection") : ""),
00024 matchedCut_(iConfig.existsAs<std::string>("matchedPreselection") ? iConfig.getParameter<std::string>("matchedPreselection") : ""),
00025 requireSameCharge_(iConfig.existsAs<bool>("requireSameCharge") ? iConfig.getParameter<bool>("requireSameCharge") : false)
00026 {
00027 std::string algo = iConfig.getParameter<std::string>("algorithm");
00028 if (algo == "byTrackRef") { algo_ = ByTrackRef; }
00029 else if (algo == "byPropagatingSrc") { algo_ = ByPropagatingSrc; }
00030 else if (algo == "byPropagatingMatched") { algo_ = ByPropagatingMatched; }
00031 else if (algo == "byDirectComparison") { algo_ = ByDirectComparison; }
00032 else throw cms::Exception("Configuration") << "Value '" << algo << "' for algorithm not yet implemented.\n";
00033
00034 getConf(iConfig, "src", whichTrack1_, whichState1_);
00035 getConf(iConfig, "matched", whichTrack2_, whichState2_);
00036
00037 if (algo_ == ByTrackRef) {
00038
00039 if (whichTrack1_ == None || whichTrack2_ == None) throw cms::Exception("Configuration") << "Algorithm 'byTrackRef' needs tracks not to be 'none'.\n";
00040 } else if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingMatched || algo_ == ByDirectComparison) {
00041
00042 maxLocalPosDiff_ = iConfig.getParameter<double>("maxDeltaLocalPos");
00043 maxGlobalMomDeltaR_ = iConfig.getParameter<double>("maxDeltaR");
00044 maxGlobalMomDeltaEta_ = iConfig.existsAs<double>("maxDeltaEta") ? iConfig.getParameter<double>("maxDeltaEta") : maxGlobalMomDeltaR_;
00045 maxGlobalMomDeltaPhi_ = iConfig.existsAs<double>("maxDeltaPhi") ? iConfig.getParameter<double>("maxDeltaPhi") : maxGlobalMomDeltaR_;
00046 maxGlobalDPtRel_ = iConfig.getParameter<double>("maxDeltaPtRel");
00047
00048
00049 std::string sortBy = iConfig.getParameter<std::string>("sortBy");
00050 if (sortBy == "deltaLocalPos") sortBy_ = LocalPosDiff;
00051 else if (sortBy == "deltaPtRel") sortBy_ = GlobalDPtRel;
00052 else if (sortBy == "deltaR") sortBy_ = GlobalMomDeltaR;
00053 else if (sortBy == "deltaEta") sortBy_ = GlobalMomDeltaEta;
00054 else if (sortBy == "deltaPhi") sortBy_ = GlobalMomDeltaPhi;
00055 else if (sortBy == "chi2") sortBy_ = Chi2;
00056 else throw cms::Exception("Configuration") << "Parameter 'sortBy' must be one of: deltaLocalPos, deltaPtRel, deltaR, chi2.\n";
00057
00058 if (algo_ == ByPropagatingSrc) {
00059 if (whichTrack2_ == None || whichState2_ == AtVertex) {
00060 algo_ = ByPropagatingSrcTSCP;
00061
00062 }
00063 } else if (algo_ == ByPropagatingMatched) {
00064 if (whichTrack1_ == None || whichState1_ == AtVertex) {
00065 algo_ = ByPropagatingMatchedTSCP;
00066
00067 }
00068 } else if (algo_ == ByDirectComparison) {
00069 bool firstAtVertex = (whichTrack1_ == None || whichState1_ == AtVertex);
00070 bool secAtVertex = (whichTrack2_ == None || whichState2_ == AtVertex);
00071 if (firstAtVertex) {
00072 if (!secAtVertex) throw cms::Exception("Configuration") << "When using 'byDirectComparison' with 'src' at vertex (or None), 'matched' must be at vertex too.\n";
00073 } else {
00074 if (secAtVertex) throw cms::Exception("Configuration") << "When using 'byDirectComparison' with 'src' not at vertex, 'matched' can't be at vertex or None.\n";
00075 if (whichState1_ != whichState2_) throw cms::Exception("Configuration") << "You can't use 'byDirectComparison' with non-matching states.\n";
00076 }
00077 }
00078
00079 useChi2_ = iConfig.existsAs<bool>("computeChi2") ? iConfig.getParameter<bool>("computeChi2") : false;
00080 if (useChi2_) {
00081 if (whichTrack1_ == None && whichTrack2_ == None) throw cms::Exception("Configuration") << "Can't compute chi2s if both tracks are set to 'none'.\n";
00082 maxChi2_ = iConfig.getParameter<double>("maxChi2");
00083 chi2DiagonalOnly_ = iConfig.getParameter<bool>("chi2DiagonalOnly");
00084 if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingMatched) {
00085 chi2UseVertex_ = iConfig.getParameter<bool>("chi2UsePosition");
00086 } else {
00087 chi2UseVertex_ = iConfig.getParameter<bool>("chi2UseVertex");
00088 if (algo_ == ByDirectComparison) {
00089 std::string choice = iConfig.getParameter<std::string>("chi2MomentumForDxy");
00090 if (choice == "src") chi2FirstMomentum_ = true;
00091 else if (choice != "matched") throw cms::Exception("Configuration") << "chi2MomentumForDxy must be 'src' or 'matched'\n";
00092 }
00093 }
00094 } else maxChi2_ = 1;
00095
00096 if (sortBy_ == Chi2 && !useChi2_) throw cms::Exception("Configuration") << "Can't sort by chi2s if 'computeChi2s' is not set to true.\n";
00097 }
00098 }
00099
00100 void
00101 MatcherUsingTracksAlgorithm::getConf(const edm::ParameterSet & iConfig, const std::string &whatFor, WhichTrack &whichTrack, WhichState &whichState) {
00102 std::string s_whichTrack = iConfig.getParameter<std::string>(whatFor+"Track");
00103 if (s_whichTrack == "none") { whichTrack = None; }
00104 else if (s_whichTrack == "tracker") { whichTrack = TrackerTk; }
00105 else if (s_whichTrack == "muon") { whichTrack = MuonTk; }
00106 else if (s_whichTrack == "global") { whichTrack = GlobalTk; }
00107 else throw cms::Exception("Configuration") << "Parameter 'useTrack' must be 'none', 'tracker', 'muon', 'global'\n";
00108 if ((whichTrack != None) && (algo_ != ByTrackRef)) {
00109 std::string s_whichState = iConfig.getParameter<std::string>(whatFor+"State");
00110 if (s_whichState == "atVertex") { whichState = AtVertex; }
00111 else if (s_whichState == "innermost") { whichState = Innermost; }
00112 else if (s_whichState == "outermost") { whichState = Outermost; }
00113 else throw cms::Exception("Configuration") << "Parameter 'useState' must be 'atVertex', 'innermost', 'outermost'\n";
00114 }
00115 }
00116
00119 bool
00120 MatcherUsingTracksAlgorithm::match(const reco::Candidate &c1, const reco::Candidate &c2, float &deltR, float &deltEta, float &deltPhi, float &deltaLocalPos, float &deltaPtRel, float &chi2) const {
00121 if (!(srcCut_(c1) && matchedCut_(c2))) return false;
00122 if (requireSameCharge_ && (c1.charge() != c2.charge())) return false;
00123 switch (algo_) {
00124 case ByTrackRef: {
00125 reco::TrackRef t1 = getTrack(c1, whichTrack1_);
00126 reco::TrackRef t2 = getTrack(c2, whichTrack2_);
00127 if (t1.isNonnull()) {
00128 if (t1 == t2) return true;
00129 if (t1.id() != t2.id()) { edm::LogWarning("MatcherUsingTracksAlgorithm") << "Trying to match by reference tracks coming from different collections.\n"; }
00130 }
00131 }
00132 case ByPropagatingSrc: {
00133 FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
00134 TrajectoryStateOnSurface target = targetState(c2, whichTrack2_, whichState2_);
00135 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00136 }
00137 case ByPropagatingMatched: {
00138 FreeTrajectoryState start = startingState(c2, whichTrack2_, whichState2_);
00139 TrajectoryStateOnSurface target = targetState(c1, whichTrack1_, whichState1_);
00140 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00141 }
00142 case ByPropagatingSrcTSCP: {
00143 FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
00144 FreeTrajectoryState target = startingState(c2, whichTrack2_, whichState2_);
00145 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00146 }
00147 case ByPropagatingMatchedTSCP: {
00148 FreeTrajectoryState start = startingState(c2, whichTrack2_, whichState2_);
00149 FreeTrajectoryState target = startingState(c1, whichTrack1_, whichState1_);
00150 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00151 }
00152
00153 case ByDirectComparison: {
00154 FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
00155 FreeTrajectoryState otherstart = startingState(c2, whichTrack2_, whichState2_);
00156 return matchByDirectComparison(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00157 }
00158 }
00159 return false;
00160 }
00161
00165 int
00166 MatcherUsingTracksAlgorithm::match(const reco::Candidate &c1, const edm::View<reco::Candidate> &c2s, float &deltR, float &deltEta, float &deltPhi, float &deltaLocalPos, float &deltaPtRel, float &chi2) const {
00167 if (!srcCut_(c1)) return -1;
00168
00169
00170 FreeTrajectoryState start; TrajectoryStateOnSurface target;
00171 int match = -1;
00172
00173
00174 if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingSrcTSCP || algo_ == ByPropagatingMatchedTSCP || algo_ == ByDirectComparison) {
00175 start = startingState(c1, whichTrack1_, whichState1_);
00176 } else if (algo_ == ByPropagatingMatched) target = targetState(c1, whichTrack1_, whichState1_);
00177
00178
00179 edm::View<reco::Candidate>::const_iterator it, ed; int i;
00180 for (it = c2s.begin(), ed = c2s.end(), i = 0; it != ed; ++it, ++i) {
00181 if (!matchedCut_(*it)) continue;
00182 if (requireSameCharge_ && (c1.charge() != it->charge()) ) continue;
00183 bool exit = false;
00184 switch (algo_) {
00185 case ByTrackRef: {
00186 reco::TrackRef t1 = getTrack(c1, whichTrack1_);
00187 reco::TrackRef t2 = getTrack(*it, whichTrack2_);
00188 if (t1.isNonnull()) {
00189 if (t1 == t2) { match = i; exit = true; }
00190 if (t1.id() != t2.id()) { edm::LogWarning("MatcherUsingTracksAlgorithm") << "Trying to match by reference tracks coming from different collections.\n"; }
00191 }
00192 } break;
00193 case ByPropagatingSrc:
00194 case ByPropagatingMatched: {
00195 if (algo_ == ByPropagatingMatched) start = startingState(*it, whichTrack2_, whichState2_);
00196 else if (algo_ == ByPropagatingSrc) target = targetState(*it, whichTrack2_, whichState2_);
00197 if (matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00198 match = i;
00199 }
00200 } break;
00201 case ByDirectComparison: {
00202 FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
00203 if (matchByDirectComparison(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00204 match = i;
00205 }
00206 } break;
00207 case ByPropagatingSrcTSCP: {
00208 FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
00209 if (matchWithPropagation(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00210 match = i;
00211 }
00212 } break;
00213 case ByPropagatingMatchedTSCP: {
00214 FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
00215 if (matchWithPropagation(otherstart, start, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00216 match = i;
00217 }
00218 } break;
00219 }
00220 if (exit) break;
00221 }
00222
00223 return match;
00224 }
00225
00226 void
00227 MatcherUsingTracksAlgorithm::init(const edm::EventSetup & iSetup) {
00228 iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
00229 iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator_);
00230 iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_);
00231 }
00232
00233
00234 reco::TrackRef
00235 MatcherUsingTracksAlgorithm::getTrack(const reco::Candidate &reco, WhichTrack whichTrack) const {
00236 reco::TrackRef tk;
00237 const reco::RecoCandidate *rc = dynamic_cast<const reco::RecoCandidate *>(&reco);
00238 if (rc == 0) throw cms::Exception("Invalid Data") << "Input object is not a RecoCandidate.\n";
00239 switch (whichTrack) {
00240 case TrackerTk: tk = rc->track(); break;
00241 case MuonTk : tk = rc->standAloneMuon(); break;
00242 case GlobalTk : tk = rc->combinedMuon(); break;
00243 default: break;
00244 }
00245 return tk;
00246 }
00247
00248 FreeTrajectoryState
00249 MatcherUsingTracksAlgorithm::startingState(const reco::Candidate &reco, WhichTrack whichTrack, WhichState whichState) const {
00250 FreeTrajectoryState ret;
00251 if (whichTrack != None) {
00252 reco::TrackRef tk = getTrack(reco, whichTrack);
00253 if (tk.isNull()) {
00254 ret = FreeTrajectoryState();
00255 } else {
00256 switch (whichState) {
00257 case AtVertex: ret = trajectoryStateTransform::initialFreeState(*tk, magfield_.product()); break;
00258 case Innermost: ret = trajectoryStateTransform::innerFreeState( *tk, magfield_.product()); break;
00259 case Outermost: ret = trajectoryStateTransform::outerFreeState( *tk, magfield_.product()); break;
00260 }
00261 }
00262 } else {
00263 ret = FreeTrajectoryState( GlobalPoint( reco.vx(), reco.vy(), reco.vz()),
00264 GlobalVector(reco.px(), reco.py(), reco.pz()),
00265 reco.charge(),
00266 magfield_.product());
00267 }
00268 return ret;
00269 }
00270
00271 TrajectoryStateOnSurface
00272 MatcherUsingTracksAlgorithm::targetState(const reco::Candidate &reco, WhichTrack whichTrack, WhichState whichState) const {
00273 TrajectoryStateOnSurface ret;
00274 reco::TrackRef tk = getTrack(reco, whichTrack);
00275 if (tk.isNonnull()) {
00276 switch (whichState) {
00277 case Innermost: ret = trajectoryStateTransform::innerStateOnSurface( *tk, *geometry_, magfield_.product()); break;
00278 case Outermost: ret = trajectoryStateTransform::outerStateOnSurface( *tk, *geometry_, magfield_.product()); break;
00279 default: break;
00280 }
00281 }
00282 return ret;
00283 }
00284
00285 bool
00286 MatcherUsingTracksAlgorithm::matchWithPropagation(const FreeTrajectoryState &start,
00287 const TrajectoryStateOnSurface &target,
00288 float &lastDeltaR,
00289 float &lastDeltaEta,
00290 float &lastDeltaPhi,
00291 float &lastDeltaLocalPos,
00292 float &lastGlobalDPtRel,
00293 float &lastChi2) const
00294 {
00295 if ((start.momentum().mag() == 0) || !target.isValid()) return false;
00296
00297 TrajectoryStateOnSurface tsos = propagator_->propagate(start, target.surface());
00298
00299 bool isBest = false;
00300 if (tsos.isValid()) {
00301 float thisLocalPosDiff = (tsos.localPosition()-target.localPosition()).mag();
00302 float thisGlobalMomDeltaR = deltaR(tsos.globalMomentum(), target.globalMomentum());
00303 float thisGlobalMomDeltaPhi = fabs(deltaPhi(tsos.globalMomentum().phi(), target.globalMomentum().phi()));
00304 float thisGlobalMomDeltaEta = fabs(tsos.globalMomentum().eta() - target.globalMomentum().eta());
00305 float thisGlobalDPtRel = (tsos.globalMomentum().perp() - target.globalMomentum().perp())/target.globalMomentum().perp();
00306
00307 if ((thisLocalPosDiff < maxLocalPosDiff_) &&
00308 (thisGlobalMomDeltaR < maxGlobalMomDeltaR_) &&
00309 (thisGlobalMomDeltaEta < maxGlobalMomDeltaEta_) &&
00310 (thisGlobalMomDeltaPhi < maxGlobalMomDeltaPhi_) &&
00311 (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
00312 float thisChi2 = useChi2_ ? getChi2(target, tsos, chi2DiagonalOnly_, chi2UseVertex_) : 0;
00313 if (thisChi2 >= maxChi2_) return false;
00314 switch (sortBy_) {
00315 case LocalPosDiff: isBest = (thisLocalPosDiff < lastDeltaLocalPos); break;
00316 case GlobalMomDeltaR: isBest = (thisGlobalMomDeltaR < lastDeltaR); break;
00317 case GlobalMomDeltaEta: isBest = (thisGlobalMomDeltaEta < lastDeltaEta); break;
00318 case GlobalMomDeltaPhi: isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi); break;
00319 case GlobalDPtRel: isBest = (thisGlobalDPtRel < lastGlobalDPtRel); break;
00320 case Chi2: isBest = (thisChi2 < lastChi2); break;
00321 }
00322 if (isBest) {
00323 lastDeltaLocalPos = thisLocalPosDiff;
00324 lastDeltaR = thisGlobalMomDeltaR;
00325 lastDeltaEta = thisGlobalMomDeltaEta;
00326 lastDeltaPhi = thisGlobalMomDeltaPhi;
00327 lastGlobalDPtRel = thisGlobalDPtRel;
00328 lastChi2 = thisChi2;
00329 }
00330 }
00331 }
00332
00333 return isBest;
00334 }
00335
00336 bool
00337 MatcherUsingTracksAlgorithm::matchWithPropagation(const FreeTrajectoryState &start,
00338 const FreeTrajectoryState &target,
00339 float &lastDeltaR,
00340 float &lastDeltaEta,
00341 float &lastDeltaPhi,
00342 float &lastDeltaLocalPos,
00343 float &lastGlobalDPtRel,
00344 float &lastChi2) const
00345 {
00346 if ((start.momentum().mag() == 0) || (target.momentum().mag() == 0)) return false;
00347 TSCPBuilderNoMaterial propagator;
00348 try {
00349 TrajectoryStateClosestToPoint tscp = propagator(start, target.position());
00350
00351
00352 bool isBest = false;
00353 float thisLocalPosDiff = (tscp.position()-target.position()).mag();
00354 float thisGlobalMomDeltaR = deltaR(tscp.momentum(), target.momentum());
00355 float thisGlobalMomDeltaPhi = fabs(deltaPhi(tscp.momentum().phi(), target.momentum().phi()));
00356 float thisGlobalMomDeltaEta = fabs(tscp.momentum().eta() - target.momentum().eta());
00357 float thisGlobalDPtRel = (tscp.momentum().perp() - target.momentum().perp())/target.momentum().perp();
00358
00359 if ((thisLocalPosDiff < maxLocalPosDiff_) &&
00360 (thisGlobalMomDeltaR < maxGlobalMomDeltaR_) &&
00361 (thisGlobalMomDeltaEta < maxGlobalMomDeltaEta_) &&
00362 (thisGlobalMomDeltaPhi < maxGlobalMomDeltaPhi_) &&
00363 (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
00364 float thisChi2 = useChi2_ ? getChi2(target, tscp, chi2DiagonalOnly_, chi2UseVertex_) : 0;
00365 if (thisChi2 >= maxChi2_) return false;
00366 switch (sortBy_) {
00367 case LocalPosDiff: isBest = (thisLocalPosDiff < lastDeltaLocalPos); break;
00368 case GlobalMomDeltaR: isBest = (thisGlobalMomDeltaR < lastDeltaR); break;
00369 case GlobalMomDeltaEta: isBest = (thisGlobalMomDeltaEta < lastDeltaEta); break;
00370 case GlobalMomDeltaPhi: isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi); break;
00371 case GlobalDPtRel: isBest = (thisGlobalDPtRel < lastGlobalDPtRel); break;
00372 case Chi2: isBest = (thisChi2 < lastChi2); break;
00373 }
00374 if (isBest) {
00375 lastDeltaLocalPos = thisLocalPosDiff;
00376 lastDeltaR = thisGlobalMomDeltaR;
00377 lastDeltaEta = thisGlobalMomDeltaEta;
00378 lastDeltaPhi = thisGlobalMomDeltaPhi;
00379 lastGlobalDPtRel = thisGlobalDPtRel;
00380 lastChi2 = thisChi2;
00381 }
00382 }
00383
00384 return isBest;
00385 } catch (const TrajectoryStateException &err) { return false; }
00386 }
00387
00388
00389 bool
00390 MatcherUsingTracksAlgorithm::matchByDirectComparison(const FreeTrajectoryState &start,
00391 const FreeTrajectoryState &target,
00392 float &lastDeltaR,
00393 float &lastDeltaEta,
00394 float &lastDeltaPhi,
00395 float &lastDeltaLocalPos,
00396 float &lastGlobalDPtRel,
00397 float &lastChi2) const
00398 {
00399 if ((start.momentum().mag() == 0) || target.momentum().mag() == 0) return false;
00400
00401 bool isBest = false;
00402 float thisLocalPosDiff = (start.position()-target.position()).mag();
00403 float thisGlobalMomDeltaR = deltaR(start.momentum(), target.momentum());
00404 float thisGlobalMomDeltaPhi = fabs(deltaPhi(start.momentum().phi(), target.momentum().phi()));
00405 float thisGlobalMomDeltaEta = fabs(start.momentum().eta() - target.momentum().eta());
00406 float thisGlobalDPtRel = (start.momentum().perp() - target.momentum().perp())/target.momentum().perp();
00407
00408 if ((thisLocalPosDiff < maxLocalPosDiff_) &&
00409 (thisGlobalMomDeltaR < maxGlobalMomDeltaR_) &&
00410 (thisGlobalMomDeltaEta < maxGlobalMomDeltaEta_) &&
00411 (thisGlobalMomDeltaPhi < maxGlobalMomDeltaPhi_) &&
00412 (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
00413 float thisChi2 = useChi2_ ? getChi2(start, target, chi2DiagonalOnly_, chi2UseVertex_, chi2FirstMomentum_) : 0;
00414 if (thisChi2 >= maxChi2_) return false;
00415 switch (sortBy_) {
00416 case LocalPosDiff: isBest = (thisLocalPosDiff < lastDeltaLocalPos); break;
00417 case GlobalMomDeltaR: isBest = (thisGlobalMomDeltaR < lastDeltaR); break;
00418 case GlobalMomDeltaEta: isBest = (thisGlobalMomDeltaEta < lastDeltaEta); break;
00419 case GlobalMomDeltaPhi: isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi); break;
00420 case GlobalDPtRel: isBest = (thisGlobalDPtRel < lastGlobalDPtRel); break;
00421 case Chi2: isBest = (thisChi2 < lastChi2); break;
00422 }
00423 if (isBest) {
00424 lastDeltaLocalPos = thisLocalPosDiff;
00425 lastDeltaR = thisGlobalMomDeltaR;
00426 lastDeltaEta = thisGlobalMomDeltaEta;
00427 lastDeltaPhi = thisGlobalMomDeltaPhi;
00428 lastGlobalDPtRel = thisGlobalDPtRel;
00429 lastChi2 = thisChi2;
00430 }
00431 }
00432
00433 return isBest;
00434 }
00435
00436 double
00437 MatcherUsingTracksAlgorithm::getChi2(const FreeTrajectoryState &start, const FreeTrajectoryState &other, bool diagonalOnly, bool useVertex, bool useFirstMomentum) {
00438 if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
00439 AlgebraicSymMatrix55 cov;
00440 if (start.hasError()) cov += start.curvilinearError().matrix();
00441 if (other.hasError()) cov += other.curvilinearError().matrix();
00442 cropAndInvert(cov, diagonalOnly, !useVertex);
00443 GlobalVector p1 = start.momentum(), p2 = other.momentum();
00444 GlobalPoint x1 = start.position(), x2 = other.position();
00445 GlobalVector p = useFirstMomentum ? p1 : p2; double pt = p.perp(), pm = p.mag();
00446 double dsz = (x1.z()-x2.z())*pt/pm - ( (x1.x()-x2.x())*p.x() + (x1.y()-x2.y())*p.y() ) / pt * p.z()/pm;
00447 double dxy = ( -(x1.x()-x2.x())*p.y() + (x1.y()-x2.y())*p.x() )/pt;
00448 AlgebraicVector5 diff(start.charge()/p1.mag()-other.charge()/p2.mag(),
00449 p1.theta()-p2.theta(),
00450 (p1.phi()-p2.phi()).value(),
00451 dxy,
00452 dsz);
00453 return ROOT::Math::Similarity(diff, cov);
00454 }
00455
00456 double
00457 MatcherUsingTracksAlgorithm::getChi2(const FreeTrajectoryState &start, const TrajectoryStateClosestToPoint &other, bool diagonalOnly, bool useVertex) {
00458 if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
00459 double pt;
00460 AlgebraicSymMatrix55 cov;
00461 if (start.hasError()) cov += PerigeeConversions::ftsToPerigeeError(start).covarianceMatrix();
00462 if (other.hasError()) cov += other.perigeeError().covarianceMatrix();
00463 cropAndInvert(cov, diagonalOnly, !useVertex);
00464 AlgebraicVector5 pgpar1 = PerigeeConversions::ftsToPerigeeParameters(start,other.referencePoint(),pt).vector();
00465 AlgebraicVector5 pgpar2 = other.perigeeParameters().vector();
00466 AlgebraicVector5 diff(pgpar1-pgpar2);
00467 return ROOT::Math::Similarity(diff, cov);
00468 }
00469
00470 double
00471 MatcherUsingTracksAlgorithm::getChi2(const TrajectoryStateOnSurface &start, const TrajectoryStateOnSurface &other, bool diagonalOnly, bool usePosition) {
00472 if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
00473 AlgebraicSymMatrix55 cov;
00474 if (start.hasError()) cov += start.localError().matrix();
00475 if (other.hasError()) cov += other.localError().matrix();
00476 cropAndInvert(cov, diagonalOnly, !usePosition);
00477 AlgebraicVector5 diff(start.localParameters().mixedFormatVector() - other.localParameters().mixedFormatVector());
00478 return ROOT::Math::Similarity(diff, cov);
00479 }
00480
00481 void
00482 MatcherUsingTracksAlgorithm::cropAndInvert(AlgebraicSymMatrix55 &cov, bool diagonalOnly, bool top3by3only) {
00483 if (!top3by3only) {
00484 if (diagonalOnly) {
00485 for (size_t i = 0; i < 5; ++i) { for (size_t j = i+1; j < 5; ++j) {
00486 cov(i,j) = 0;
00487 } }
00488 }
00489 cov.Invert();
00490 } else {
00491
00492 AlgebraicSymMatrix33 momCov = cov.Sub<AlgebraicSymMatrix33>(0,0);
00493 if (diagonalOnly) { momCov(0,1) = 0; momCov(0,2) = 0; momCov(1,2) = 0; }
00494
00495 momCov.Invert();
00496
00497 cov.Place_at(momCov,0,0);
00498
00499 for (size_t i = 3; i < 5; ++i) { for (size_t j = i; j < 5; ++j) {
00500 cov(i,j) = 0;
00501 } }
00502 }
00503 }
00504
00505