Back to home page

Project CMSSW displayed by LXR
 
 

    

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     // validate the config
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     // read matching cuts
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     // choice of sorting variable
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     // validate the config
0081     if (algo_ == ByPropagatingSrc) {
0082       if (whichTrack2_ == None || whichState2_ == AtVertex) {
0083         algo_ = ByPropagatingSrcTSCP;
0084         //throw cms::Exception("Configuration") << "Destination track must be non-null, and state must not be 'AtVertex' (not yet).\n";
0085       }
0086     } else if (algo_ == ByPropagatingMatched) {
0087       if (whichTrack1_ == None || whichState1_ == AtVertex) {
0088         algo_ = ByPropagatingMatchedTSCP;
0089         //throw cms::Exception("Configuration") << "Destination track must be non-null, and state must not be 'AtVertex' (not yet).\n";
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 /// Try to match one track to another one. Return true if succeeded.
0162 /// For matches not by ref, it will update deltaR, deltaLocalPos and deltaPtRel if the match suceeded
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 /// Find the best match to another candidate, and return its index in the vector
0220 /// For matches not by ref, it will update deltaR, deltaLocalPos and deltaPtRel if the match suceeded
0221 /// Returns -1 if the match fails
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   // working and output variables
0234   FreeTrajectoryState start;
0235   TrajectoryStateOnSurface target;
0236   int match = -1;
0237 
0238   // pre-fetch some states if needed
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   // loop on the  collection
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;  // just to make gcc happy
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;  // just to make gcc happy
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     }  // if match
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   /*2.2.X*/ try {
0456     TrajectoryStateClosestToPoint tscp = propagator(start, target.position());
0457     // if (!tscp.isValid()) return false;  // in 3.1.X
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     }  // if match
0501 
0502     return isBest;
0503     /*2.2.X*/ } 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     }  // if match
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;  // needed by pgconvert
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       // get 3x3 covariance
0638       AlgebraicSymMatrix33 momCov = cov.Sub<AlgebraicSymMatrix33>(0, 0);  // get 3x3 matrix
0639       if (diagonalOnly) {
0640         momCov(0, 1) = 0;
0641         momCov(0, 2) = 0;
0642         momCov(1, 2) = 0;
0643       }
0644       // invert
0645       momCov.Invert();
0646       // place it
0647       cov.Place_at(momCov, 0, 0);
0648       // zero the rest
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!