File indexing completed on 2023-05-31 00:38:06
0001 #include "MuonAnalysis/MuonAssociators/interface/MatcherUsingTracksAlgorithm.h"
0002
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005
0006 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008
0009 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0011 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0013 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0014 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
0015
0016 #include "FWCore/Utilities/interface/Exception.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 MatcherUsingTracksAlgorithm::MatcherUsingTracksAlgorithm(const edm::ParameterSet &iConfig, edm::ConsumesCollector iC)
0020 : whichTrack1_(None),
0021 whichTrack2_(None),
0022 whichState1_(AtVertex),
0023 whichState2_(AtVertex),
0024 srcCut_(iConfig.existsAs<std::string>("srcPreselection") ? iConfig.getParameter<std::string>("srcPreselection")
0025 : ""),
0026 matchedCut_(iConfig.existsAs<std::string>("matchedPreselection")
0027 ? iConfig.getParameter<std::string>("matchedPreselection")
0028 : ""),
0029 requireSameCharge_(iConfig.existsAs<bool>("requireSameCharge") ? iConfig.getParameter<bool>("requireSameCharge")
0030 : false),
0031 magfieldToken_(iC.esConsumes()),
0032 propagatorToken_(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0033 geometryToken_(iC.esConsumes()) {
0034 std::string algo = iConfig.getParameter<std::string>("algorithm");
0035 if (algo == "byTrackRef") {
0036 algo_ = ByTrackRef;
0037 } else if (algo == "byPropagatingSrc") {
0038 algo_ = ByPropagatingSrc;
0039 } else if (algo == "byPropagatingMatched") {
0040 algo_ = ByPropagatingMatched;
0041 } else if (algo == "byDirectComparison") {
0042 algo_ = ByDirectComparison;
0043 } else
0044 throw cms::Exception("Configuration") << "Value '" << algo << "' for algorithm not yet implemented.\n";
0045
0046 getConf(iConfig, "src", whichTrack1_, whichState1_);
0047 getConf(iConfig, "matched", whichTrack2_, whichState2_);
0048
0049 if (algo_ == ByTrackRef) {
0050
0051 if (whichTrack1_ == None || whichTrack2_ == None)
0052 throw cms::Exception("Configuration") << "Algorithm 'byTrackRef' needs tracks not to be 'none'.\n";
0053 } else if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingMatched || algo_ == ByDirectComparison) {
0054
0055 maxLocalPosDiff_ = iConfig.getParameter<double>("maxDeltaLocalPos");
0056 maxGlobalMomDeltaR_ = iConfig.getParameter<double>("maxDeltaR");
0057 maxGlobalMomDeltaEta_ =
0058 iConfig.existsAs<double>("maxDeltaEta") ? iConfig.getParameter<double>("maxDeltaEta") : maxGlobalMomDeltaR_;
0059 maxGlobalMomDeltaPhi_ =
0060 iConfig.existsAs<double>("maxDeltaPhi") ? iConfig.getParameter<double>("maxDeltaPhi") : maxGlobalMomDeltaR_;
0061 maxGlobalDPtRel_ = iConfig.getParameter<double>("maxDeltaPtRel");
0062
0063
0064 std::string sortBy = iConfig.getParameter<std::string>("sortBy");
0065 if (sortBy == "deltaLocalPos")
0066 sortBy_ = LocalPosDiff;
0067 else if (sortBy == "deltaPtRel")
0068 sortBy_ = GlobalDPtRel;
0069 else if (sortBy == "deltaR")
0070 sortBy_ = GlobalMomDeltaR;
0071 else if (sortBy == "deltaEta")
0072 sortBy_ = GlobalMomDeltaEta;
0073 else if (sortBy == "deltaPhi")
0074 sortBy_ = GlobalMomDeltaPhi;
0075 else if (sortBy == "chi2")
0076 sortBy_ = Chi2;
0077 else
0078 throw cms::Exception("Configuration")
0079 << "Parameter 'sortBy' must be one of: deltaLocalPos, deltaPtRel, deltaR, chi2.\n";
0080
0081 if (algo_ == ByPropagatingSrc) {
0082 if (whichTrack2_ == None || whichState2_ == AtVertex) {
0083 algo_ = ByPropagatingSrcTSCP;
0084
0085 }
0086 } else if (algo_ == ByPropagatingMatched) {
0087 if (whichTrack1_ == None || whichState1_ == AtVertex) {
0088 algo_ = ByPropagatingMatchedTSCP;
0089
0090 }
0091 } else if (algo_ == ByDirectComparison) {
0092 bool firstAtVertex = (whichTrack1_ == None || whichState1_ == AtVertex);
0093 bool secAtVertex = (whichTrack2_ == None || whichState2_ == AtVertex);
0094 if (firstAtVertex) {
0095 if (!secAtVertex)
0096 throw cms::Exception("Configuration")
0097 << "When using 'byDirectComparison' with 'src' at vertex (or None), 'matched' must be at vertex too.\n";
0098 } else {
0099 if (secAtVertex)
0100 throw cms::Exception("Configuration")
0101 << "When using 'byDirectComparison' with 'src' not at vertex, 'matched' can't be at vertex or None.\n";
0102 if (whichState1_ != whichState2_)
0103 throw cms::Exception("Configuration") << "You can't use 'byDirectComparison' with non-matching states.\n";
0104 }
0105 }
0106
0107 useChi2_ = iConfig.existsAs<bool>("computeChi2") ? iConfig.getParameter<bool>("computeChi2") : false;
0108 if (useChi2_) {
0109 if (whichTrack1_ == None && whichTrack2_ == None)
0110 throw cms::Exception("Configuration") << "Can't compute chi2s if both tracks are set to 'none'.\n";
0111 maxChi2_ = iConfig.getParameter<double>("maxChi2");
0112 chi2DiagonalOnly_ = iConfig.getParameter<bool>("chi2DiagonalOnly");
0113 if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingMatched) {
0114 chi2UseVertex_ = iConfig.getParameter<bool>("chi2UsePosition");
0115 } else {
0116 chi2UseVertex_ = iConfig.getParameter<bool>("chi2UseVertex");
0117 if (algo_ == ByDirectComparison) {
0118 std::string choice = iConfig.getParameter<std::string>("chi2MomentumForDxy");
0119 if (choice == "src")
0120 chi2FirstMomentum_ = true;
0121 else if (choice != "matched")
0122 throw cms::Exception("Configuration") << "chi2MomentumForDxy must be 'src' or 'matched'\n";
0123 }
0124 }
0125 } else
0126 maxChi2_ = 1;
0127
0128 if (sortBy_ == Chi2 && !useChi2_)
0129 throw cms::Exception("Configuration") << "Can't sort by chi2s if 'computeChi2s' is not set to true.\n";
0130 }
0131 }
0132
0133 void MatcherUsingTracksAlgorithm::getConf(const edm::ParameterSet &iConfig,
0134 const std::string &whatFor,
0135 WhichTrack &whichTrack,
0136 WhichState &whichState) {
0137 std::string s_whichTrack = iConfig.getParameter<std::string>(whatFor + "Track");
0138 if (s_whichTrack == "none") {
0139 whichTrack = None;
0140 } else if (s_whichTrack == "tracker") {
0141 whichTrack = TrackerTk;
0142 } else if (s_whichTrack == "muon") {
0143 whichTrack = MuonTk;
0144 } else if (s_whichTrack == "global") {
0145 whichTrack = GlobalTk;
0146 } else
0147 throw cms::Exception("Configuration") << "Parameter 'useTrack' must be 'none', 'tracker', 'muon', 'global'\n";
0148 if ((whichTrack != None) && (algo_ != ByTrackRef)) {
0149 std::string s_whichState = iConfig.getParameter<std::string>(whatFor + "State");
0150 if (s_whichState == "atVertex") {
0151 whichState = AtVertex;
0152 } else if (s_whichState == "innermost") {
0153 whichState = Innermost;
0154 } else if (s_whichState == "outermost") {
0155 whichState = Outermost;
0156 } else
0157 throw cms::Exception("Configuration") << "Parameter 'useState' must be 'atVertex', 'innermost', 'outermost'\n";
0158 }
0159 }
0160
0161
0162
0163 bool MatcherUsingTracksAlgorithm::match(const reco::Candidate &c1,
0164 const reco::Candidate &c2,
0165 float &deltR,
0166 float &deltEta,
0167 float &deltPhi,
0168 float &deltaLocalPos,
0169 float &deltaPtRel,
0170 float &chi2) const {
0171 if (!(srcCut_(c1) && matchedCut_(c2)))
0172 return false;
0173 if (requireSameCharge_ && (c1.charge() != c2.charge()))
0174 return false;
0175 switch (algo_) {
0176 case ByTrackRef: {
0177 reco::TrackRef t1 = getTrack(c1, whichTrack1_);
0178 reco::TrackRef t2 = getTrack(c2, whichTrack2_);
0179 if (t1.isNonnull()) {
0180 if (t1 == t2)
0181 return true;
0182 if (t1.id() != t2.id()) {
0183 edm::LogWarning("MatcherUsingTracksAlgorithm")
0184 << "Trying to match by reference tracks coming from different collections.\n";
0185 }
0186 }
0187 }
0188 [[fallthrough]];
0189 case ByPropagatingSrc: {
0190 FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
0191 TrajectoryStateOnSurface target = targetState(c2, whichTrack2_, whichState2_);
0192 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
0193 }
0194 case ByPropagatingMatched: {
0195 FreeTrajectoryState start = startingState(c2, whichTrack2_, whichState2_);
0196 TrajectoryStateOnSurface target = targetState(c1, whichTrack1_, whichState1_);
0197 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
0198 }
0199 case ByPropagatingSrcTSCP: {
0200 FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
0201 FreeTrajectoryState target = startingState(c2, whichTrack2_, whichState2_);
0202 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
0203 }
0204 case ByPropagatingMatchedTSCP: {
0205 FreeTrajectoryState start = startingState(c2, whichTrack2_, whichState2_);
0206 FreeTrajectoryState target = startingState(c1, whichTrack1_, whichState1_);
0207 return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
0208 }
0209
0210 case ByDirectComparison: {
0211 FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
0212 FreeTrajectoryState otherstart = startingState(c2, whichTrack2_, whichState2_);
0213 return matchByDirectComparison(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
0214 }
0215 }
0216 return false;
0217 }
0218
0219
0220
0221
0222 int MatcherUsingTracksAlgorithm::match(const reco::Candidate &c1,
0223 const edm::View<reco::Candidate> &c2s,
0224 float &deltR,
0225 float &deltEta,
0226 float &deltPhi,
0227 float &deltaLocalPos,
0228 float &deltaPtRel,
0229 float &chi2) const {
0230 if (!srcCut_(c1))
0231 return -1;
0232
0233
0234 FreeTrajectoryState start;
0235 TrajectoryStateOnSurface target;
0236 int match = -1;
0237
0238
0239 if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingSrcTSCP || algo_ == ByPropagatingMatchedTSCP ||
0240 algo_ == ByDirectComparison) {
0241 start = startingState(c1, whichTrack1_, whichState1_);
0242 } else if (algo_ == ByPropagatingMatched)
0243 target = targetState(c1, whichTrack1_, whichState1_);
0244
0245
0246 edm::View<reco::Candidate>::const_iterator it, ed;
0247 int i;
0248 for (it = c2s.begin(), ed = c2s.end(), i = 0; it != ed; ++it, ++i) {
0249 if (!matchedCut_(*it))
0250 continue;
0251 if (requireSameCharge_ && (c1.charge() != it->charge()))
0252 continue;
0253 bool exit = false;
0254 switch (algo_) {
0255 case ByTrackRef: {
0256 reco::TrackRef t1 = getTrack(c1, whichTrack1_);
0257 reco::TrackRef t2 = getTrack(*it, whichTrack2_);
0258 if (t1.isNonnull()) {
0259 if (t1 == t2) {
0260 match = i;
0261 exit = true;
0262 }
0263 if (t1.id() != t2.id()) {
0264 edm::LogWarning("MatcherUsingTracksAlgorithm")
0265 << "Trying to match by reference tracks coming from different collections.\n";
0266 }
0267 }
0268 } break;
0269 case ByPropagatingSrc:
0270 case ByPropagatingMatched: {
0271 if (algo_ == ByPropagatingMatched)
0272 start = startingState(*it, whichTrack2_, whichState2_);
0273 else if (algo_ == ByPropagatingSrc)
0274 target = targetState(*it, whichTrack2_, whichState2_);
0275 if (matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
0276 match = i;
0277 }
0278 } break;
0279 case ByDirectComparison: {
0280 FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
0281 if (matchByDirectComparison(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
0282 match = i;
0283 }
0284 } break;
0285 case ByPropagatingSrcTSCP: {
0286 FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
0287 if (matchWithPropagation(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
0288 match = i;
0289 }
0290 } break;
0291 case ByPropagatingMatchedTSCP: {
0292 FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
0293 if (matchWithPropagation(otherstart, start, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
0294 match = i;
0295 }
0296 } break;
0297 }
0298 if (exit)
0299 break;
0300 }
0301
0302 return match;
0303 }
0304
0305 void MatcherUsingTracksAlgorithm::init(const edm::EventSetup &iSetup) {
0306 magfield_ = iSetup.getHandle(magfieldToken_);
0307 propagator_ = iSetup.getHandle(propagatorToken_);
0308 geometry_ = iSetup.getHandle(geometryToken_);
0309 }
0310
0311 reco::TrackRef MatcherUsingTracksAlgorithm::getTrack(const reco::Candidate &reco, WhichTrack whichTrack) const {
0312 reco::TrackRef tk;
0313 const reco::RecoCandidate *rc = dynamic_cast<const reco::RecoCandidate *>(&reco);
0314 if (rc == nullptr)
0315 throw cms::Exception("Invalid Data") << "Input object is not a RecoCandidate.\n";
0316 switch (whichTrack) {
0317 case TrackerTk:
0318 tk = rc->track();
0319 break;
0320 case MuonTk:
0321 tk = rc->standAloneMuon();
0322 break;
0323 case GlobalTk:
0324 tk = rc->combinedMuon();
0325 break;
0326 default:
0327 break;
0328 }
0329 return tk;
0330 }
0331
0332 FreeTrajectoryState MatcherUsingTracksAlgorithm::startingState(const reco::Candidate &reco,
0333 WhichTrack whichTrack,
0334 WhichState whichState) const {
0335 FreeTrajectoryState ret;
0336 if (whichTrack != None) {
0337 reco::TrackRef tk = getTrack(reco, whichTrack);
0338 if (tk.isNull()) {
0339 ret = FreeTrajectoryState();
0340 } else {
0341 switch (whichState) {
0342 case AtVertex:
0343 ret = trajectoryStateTransform::initialFreeState(*tk, magfield_.product());
0344 break;
0345 case Innermost:
0346 ret = trajectoryStateTransform::innerFreeState(*tk, magfield_.product());
0347 break;
0348 case Outermost:
0349 ret = trajectoryStateTransform::outerFreeState(*tk, magfield_.product());
0350 break;
0351 }
0352 }
0353 } else {
0354 ret = FreeTrajectoryState(GlobalPoint(reco.vx(), reco.vy(), reco.vz()),
0355 GlobalVector(reco.px(), reco.py(), reco.pz()),
0356 reco.charge(),
0357 magfield_.product());
0358 }
0359 return ret;
0360 }
0361
0362 TrajectoryStateOnSurface MatcherUsingTracksAlgorithm::targetState(const reco::Candidate &reco,
0363 WhichTrack whichTrack,
0364 WhichState whichState) const {
0365 TrajectoryStateOnSurface ret;
0366 reco::TrackRef tk = getTrack(reco, whichTrack);
0367 if (tk.isNonnull()) {
0368 switch (whichState) {
0369 case Innermost:
0370 ret = trajectoryStateTransform::innerStateOnSurface(*tk, *geometry_, magfield_.product());
0371 break;
0372 case Outermost:
0373 ret = trajectoryStateTransform::outerStateOnSurface(*tk, *geometry_, magfield_.product());
0374 break;
0375 default:
0376 break;
0377 }
0378 }
0379 return ret;
0380 }
0381
0382 bool MatcherUsingTracksAlgorithm::matchWithPropagation(const FreeTrajectoryState &start,
0383 const TrajectoryStateOnSurface &target,
0384 float &lastDeltaR,
0385 float &lastDeltaEta,
0386 float &lastDeltaPhi,
0387 float &lastDeltaLocalPos,
0388 float &lastGlobalDPtRel,
0389 float &lastChi2) const {
0390 if ((start.momentum().mag() == 0) || !target.isValid())
0391 return false;
0392
0393 TrajectoryStateOnSurface tsos = propagator_->propagate(start, target.surface());
0394
0395 bool isBest = false;
0396 if (tsos.isValid()) {
0397 float thisLocalPosDiff = (tsos.localPosition() - target.localPosition()).mag();
0398 float thisGlobalMomDeltaR = deltaR(tsos.globalMomentum(), target.globalMomentum());
0399 float thisGlobalMomDeltaPhi = fabs(deltaPhi(tsos.globalMomentum().barePhi(), target.globalMomentum().barePhi()));
0400 float thisGlobalMomDeltaEta = fabs(tsos.globalMomentum().eta() - target.globalMomentum().eta());
0401 float thisGlobalDPtRel =
0402 (tsos.globalMomentum().perp() - target.globalMomentum().perp()) / target.globalMomentum().perp();
0403
0404 if ((thisLocalPosDiff < maxLocalPosDiff_) && (thisGlobalMomDeltaR < maxGlobalMomDeltaR_) &&
0405 (thisGlobalMomDeltaEta < maxGlobalMomDeltaEta_) && (thisGlobalMomDeltaPhi < maxGlobalMomDeltaPhi_) &&
0406 (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
0407 float thisChi2 = useChi2_ ? getChi2(target, tsos, chi2DiagonalOnly_, chi2UseVertex_) : 0;
0408 if (thisChi2 >= maxChi2_)
0409 return false;
0410 switch (sortBy_) {
0411 case LocalPosDiff:
0412 isBest = (thisLocalPosDiff < lastDeltaLocalPos);
0413 break;
0414 case GlobalMomDeltaR:
0415 isBest = (thisGlobalMomDeltaR < lastDeltaR);
0416 break;
0417 case GlobalMomDeltaEta:
0418 isBest = (thisGlobalMomDeltaEta < lastDeltaEta);
0419 break;
0420 case GlobalMomDeltaPhi:
0421 isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi);
0422 break;
0423 case GlobalDPtRel:
0424 isBest = (thisGlobalDPtRel < lastGlobalDPtRel);
0425 break;
0426 case Chi2:
0427 isBest = (thisChi2 < lastChi2);
0428 break;
0429 }
0430 if (isBest) {
0431 lastDeltaLocalPos = thisLocalPosDiff;
0432 lastDeltaR = thisGlobalMomDeltaR;
0433 lastDeltaEta = thisGlobalMomDeltaEta;
0434 lastDeltaPhi = thisGlobalMomDeltaPhi;
0435 lastGlobalDPtRel = thisGlobalDPtRel;
0436 lastChi2 = thisChi2;
0437 }
0438 }
0439 }
0440
0441 return isBest;
0442 }
0443
0444 bool MatcherUsingTracksAlgorithm::matchWithPropagation(const FreeTrajectoryState &start,
0445 const FreeTrajectoryState &target,
0446 float &lastDeltaR,
0447 float &lastDeltaEta,
0448 float &lastDeltaPhi,
0449 float &lastDeltaLocalPos,
0450 float &lastGlobalDPtRel,
0451 float &lastChi2) const {
0452 if ((start.momentum().mag() == 0) || (target.momentum().mag() == 0))
0453 return false;
0454 TSCPBuilderNoMaterial propagator;
0455 try {
0456 TrajectoryStateClosestToPoint tscp = propagator(start, target.position());
0457
0458
0459 bool isBest = false;
0460 float thisLocalPosDiff = (tscp.position() - target.position()).mag();
0461 float thisGlobalMomDeltaR = deltaR(tscp.momentum(), target.momentum());
0462 float thisGlobalMomDeltaPhi = fabs(deltaPhi(tscp.momentum().barePhi(), target.momentum().barePhi()));
0463 float thisGlobalMomDeltaEta = fabs(tscp.momentum().eta() - target.momentum().eta());
0464 float thisGlobalDPtRel = (tscp.momentum().perp() - target.momentum().perp()) / target.momentum().perp();
0465
0466 if ((thisLocalPosDiff < maxLocalPosDiff_) && (thisGlobalMomDeltaR < maxGlobalMomDeltaR_) &&
0467 (thisGlobalMomDeltaEta < maxGlobalMomDeltaEta_) && (thisGlobalMomDeltaPhi < maxGlobalMomDeltaPhi_) &&
0468 (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
0469 float thisChi2 = useChi2_ ? getChi2(target, tscp, chi2DiagonalOnly_, chi2UseVertex_) : 0;
0470 if (thisChi2 >= maxChi2_)
0471 return false;
0472 switch (sortBy_) {
0473 case LocalPosDiff:
0474 isBest = (thisLocalPosDiff < lastDeltaLocalPos);
0475 break;
0476 case GlobalMomDeltaR:
0477 isBest = (thisGlobalMomDeltaR < lastDeltaR);
0478 break;
0479 case GlobalMomDeltaEta:
0480 isBest = (thisGlobalMomDeltaEta < lastDeltaEta);
0481 break;
0482 case GlobalMomDeltaPhi:
0483 isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi);
0484 break;
0485 case GlobalDPtRel:
0486 isBest = (thisGlobalDPtRel < lastGlobalDPtRel);
0487 break;
0488 case Chi2:
0489 isBest = (thisChi2 < lastChi2);
0490 break;
0491 }
0492 if (isBest) {
0493 lastDeltaLocalPos = thisLocalPosDiff;
0494 lastDeltaR = thisGlobalMomDeltaR;
0495 lastDeltaEta = thisGlobalMomDeltaEta;
0496 lastDeltaPhi = thisGlobalMomDeltaPhi;
0497 lastGlobalDPtRel = thisGlobalDPtRel;
0498 lastChi2 = thisChi2;
0499 }
0500 }
0501
0502 return isBest;
0503 } catch (const TrajectoryStateException &err) { return false; }
0504 }
0505
0506 bool MatcherUsingTracksAlgorithm::matchByDirectComparison(const FreeTrajectoryState &start,
0507 const FreeTrajectoryState &target,
0508 float &lastDeltaR,
0509 float &lastDeltaEta,
0510 float &lastDeltaPhi,
0511 float &lastDeltaLocalPos,
0512 float &lastGlobalDPtRel,
0513 float &lastChi2) const {
0514 if ((start.momentum().mag() == 0) || target.momentum().mag() == 0)
0515 return false;
0516
0517 bool isBest = false;
0518 float thisLocalPosDiff = (start.position() - target.position()).mag();
0519 float thisGlobalMomDeltaR = deltaR(start.momentum(), target.momentum());
0520 float thisGlobalMomDeltaPhi = fabs(deltaPhi(start.momentum().barePhi(), target.momentum().barePhi()));
0521 float thisGlobalMomDeltaEta = fabs(start.momentum().eta() - target.momentum().eta());
0522 float thisGlobalDPtRel = (start.momentum().perp() - target.momentum().perp()) / target.momentum().perp();
0523
0524 if ((thisLocalPosDiff < maxLocalPosDiff_) && (thisGlobalMomDeltaR < maxGlobalMomDeltaR_) &&
0525 (thisGlobalMomDeltaEta < maxGlobalMomDeltaEta_) && (thisGlobalMomDeltaPhi < maxGlobalMomDeltaPhi_) &&
0526 (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
0527 float thisChi2 = useChi2_ ? getChi2(start, target, chi2DiagonalOnly_, chi2UseVertex_, chi2FirstMomentum_) : 0;
0528 if (thisChi2 >= maxChi2_)
0529 return false;
0530 switch (sortBy_) {
0531 case LocalPosDiff:
0532 isBest = (thisLocalPosDiff < lastDeltaLocalPos);
0533 break;
0534 case GlobalMomDeltaR:
0535 isBest = (thisGlobalMomDeltaR < lastDeltaR);
0536 break;
0537 case GlobalMomDeltaEta:
0538 isBest = (thisGlobalMomDeltaEta < lastDeltaEta);
0539 break;
0540 case GlobalMomDeltaPhi:
0541 isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi);
0542 break;
0543 case GlobalDPtRel:
0544 isBest = (thisGlobalDPtRel < lastGlobalDPtRel);
0545 break;
0546 case Chi2:
0547 isBest = (thisChi2 < lastChi2);
0548 break;
0549 }
0550 if (isBest) {
0551 lastDeltaLocalPos = thisLocalPosDiff;
0552 lastDeltaR = thisGlobalMomDeltaR;
0553 lastDeltaEta = thisGlobalMomDeltaEta;
0554 lastDeltaPhi = thisGlobalMomDeltaPhi;
0555 lastGlobalDPtRel = thisGlobalDPtRel;
0556 lastChi2 = thisChi2;
0557 }
0558 }
0559
0560 return isBest;
0561 }
0562
0563 double MatcherUsingTracksAlgorithm::getChi2(const FreeTrajectoryState &start,
0564 const FreeTrajectoryState &other,
0565 bool diagonalOnly,
0566 bool useVertex,
0567 bool useFirstMomentum) {
0568 if (!start.hasError() && !other.hasError())
0569 throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
0570 AlgebraicSymMatrix55 cov;
0571 if (start.hasError())
0572 cov += start.curvilinearError().matrix();
0573 if (other.hasError())
0574 cov += other.curvilinearError().matrix();
0575 cropAndInvert(cov, diagonalOnly, !useVertex);
0576 GlobalVector p1 = start.momentum(), p2 = other.momentum();
0577 GlobalPoint x1 = start.position(), x2 = other.position();
0578 GlobalVector p = useFirstMomentum ? p1 : p2;
0579 double pt = p.perp(), pm = p.mag();
0580 double dsz =
0581 (x1.z() - x2.z()) * pt / pm - ((x1.x() - x2.x()) * p.x() + (x1.y() - x2.y()) * p.y()) / pt * p.z() / pm;
0582 double dxy = (-(x1.x() - x2.x()) * p.y() + (x1.y() - x2.y()) * p.x()) / pt;
0583 AlgebraicVector5 diff(start.charge() / p1.mag() - other.charge() / p2.mag(),
0584 p1.theta() - p2.theta(),
0585 (p1.phi() - p2.phi()).value(),
0586 dxy,
0587 dsz);
0588 return ROOT::Math::Similarity(diff, cov);
0589 }
0590
0591 double MatcherUsingTracksAlgorithm::getChi2(const FreeTrajectoryState &start,
0592 const TrajectoryStateClosestToPoint &other,
0593 bool diagonalOnly,
0594 bool useVertex) {
0595 if (!start.hasError() && !other.hasError())
0596 throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
0597 double pt;
0598 AlgebraicSymMatrix55 cov;
0599 if (start.hasError())
0600 cov += PerigeeConversions::ftsToPerigeeError(start).covarianceMatrix();
0601 if (other.hasError())
0602 cov += other.perigeeError().covarianceMatrix();
0603 cropAndInvert(cov, diagonalOnly, !useVertex);
0604 AlgebraicVector5 pgpar1 = PerigeeConversions::ftsToPerigeeParameters(start, other.referencePoint(), pt).vector();
0605 AlgebraicVector5 pgpar2 = other.perigeeParameters().vector();
0606 AlgebraicVector5 diff(pgpar1 - pgpar2);
0607 return ROOT::Math::Similarity(diff, cov);
0608 }
0609
0610 double MatcherUsingTracksAlgorithm::getChi2(const TrajectoryStateOnSurface &start,
0611 const TrajectoryStateOnSurface &other,
0612 bool diagonalOnly,
0613 bool usePosition) {
0614 if (!start.hasError() && !other.hasError())
0615 throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
0616 AlgebraicSymMatrix55 cov;
0617 if (start.hasError())
0618 cov += start.localError().matrix();
0619 if (other.hasError())
0620 cov += other.localError().matrix();
0621 cropAndInvert(cov, diagonalOnly, !usePosition);
0622 AlgebraicVector5 diff(start.localParameters().mixedFormatVector() - other.localParameters().mixedFormatVector());
0623 return ROOT::Math::Similarity(diff, cov);
0624 }
0625
0626 void MatcherUsingTracksAlgorithm::cropAndInvert(AlgebraicSymMatrix55 &cov, bool diagonalOnly, bool top3by3only) {
0627 if (!top3by3only) {
0628 if (diagonalOnly) {
0629 for (size_t i = 0; i < 5; ++i) {
0630 for (size_t j = i + 1; j < 5; ++j) {
0631 cov(i, j) = 0;
0632 }
0633 }
0634 }
0635 cov.Invert();
0636 } else {
0637
0638 AlgebraicSymMatrix33 momCov = cov.Sub<AlgebraicSymMatrix33>(0, 0);
0639 if (diagonalOnly) {
0640 momCov(0, 1) = 0;
0641 momCov(0, 2) = 0;
0642 momCov(1, 2) = 0;
0643 }
0644
0645 momCov.Invert();
0646
0647 cov.Place_at(momCov, 0, 0);
0648
0649 for (size_t i = 3; i < 5; ++i) {
0650 for (size_t j = i; j < 5; ++j) {
0651 cov(i, j) = 0;
0652 }
0653 }
0654 }
0655 }