MuonAssociators/interface/L1MuonMatcherAlgo.h

0001 #ifndef MuonAnalysis_MuonAssociators_interface_L1MuonMatcherAlgo_h
0002 #define MuonAnalysis_MuonAssociators_interface_L1MuonMatcherAlgo_h
0003 /**
0004   \class    L1MuonMatcherAlgo L1MuonMatcherAlgo.h "MuonAnalysis/MuonAssociators/interface/L1MuonMatcherAlgo.h"
0005   \brief    Matcher of reconstructed objects to L1 Muons
0006
0007   \author   Giovanni Petrucciani
0008   \version  $Id: L1MuonMatcherAlgo.h,v 1.8 2010/07/01 07:40:52 gpetrucc Exp $
0009 */
0010 #include <cmath>
0011 #include <type_traits>
0012 #include "DataFormats/Math/interface/deltaR.h"
0013 #include "DataFormats/Math/interface/deltaPhi.h"
0014 #include "DataFormats/Candidate/interface/Candidate.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0017 #include "DataFormats/L1Trigger/interface/Muon.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0023 #include "CommonTools/Utils/interface/AnySelector.h"
0024 #include "MuonAnalysis/MuonAssociators/interface/PropagateToMuonSetup.h"
0025
0026 template <edm::Transition Tr>
0027 class L1MuonMatcherAlgoT {
0028 public:
0029   explicit L1MuonMatcherAlgoT(const edm::ParameterSet &iConfig, edm::ConsumesCollector);
0030   ~L1MuonMatcherAlgoT() = default;
0031
0032   /// Call this method at the beginning of each run, to initialize geometry, magnetic field and propagators
0033   void init(const edm::EventSetup &iSetup);
0034
0035   /// Extrapolate reco::Track to the muon station 2, return an invalid TSOS if it fails
0036   TrajectoryStateOnSurface extrapolate(const reco::Track &tk) const { return prop_.extrapolate(tk); }
0037
0038   /// Extrapolate reco::Candidate to the muon station 2, return an invalid TSOS if it fails
0039   TrajectoryStateOnSurface extrapolate(const reco::Candidate &tk) const { return prop_.extrapolate(tk); }
0040
0041   /// Extrapolate a SimTrack to the muon station 2, return an invalid TSOS if it fails. Requires SimVertices to know where to start from.
0042   TrajectoryStateOnSurface extrapolate(const SimTrack &tk, const edm::SimVertexContainer &vtx) const {
0043     return prop_.extrapolate(tk, vtx);
0044   }
0045
0046   /// Extrapolate a FreeTrajectoryState to the muon station 2, return an invalid TSOS if it fails
0047   TrajectoryStateOnSurface extrapolate(const FreeTrajectoryState &state) const { return prop_.extrapolate(state); }
0048
0049   /// Return the propagator to second muon station (in case it's needed)
0050   PropagateToMuon &propagatorToMuon() { return prop_; }
0051   /// Return the propagator to second muon station (in case it's needed)
0052   const PropagateToMuon &propagatorToMuon() const { return prop_; }
0053
0054   /// Try to match one track to one L1. Return true if succeeded (and update deltaR, deltaPhi and propagated TSOS accordingly)
0055   /// The preselection cut on L1, if specified in the config, is applied before the match
0056   bool match(const reco::Track &tk,
0057              const l1extra::L1MuonParticle &l1,
0058              float &deltaR,
0059              float &deltaPhi,
0060              TrajectoryStateOnSurface &propagated) const {
0061     propagated = extrapolate(tk);
0062     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : false;
0063   }
0064
0065   /// Try to match one track to one L1. Return true if succeeded (and update deltaR, deltaPhi and propagated TSOS accordingly)
0066   /// The preselection cut on L1, if specified in the config, is applied before the match
0067   bool match(const reco::Candidate &c,
0068              const l1extra::L1MuonParticle &l1,
0069              float &deltaR,
0070              float &deltaPhi,
0071              TrajectoryStateOnSurface &propagated) const {
0072     propagated = extrapolate(c);
0073     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : false;
0074   }
0075
0076   /// Try to match one simtrack to one L1. Return true if succeeded (and update deltaR, deltaPhi and propagated TSOS accordingly)
0077   /// The preselection cut on L1, if specified in the config, is applied before the match
0078   bool match(const SimTrack &tk,
0079              const edm::SimVertexContainer &vtxs,
0080              const l1extra::L1MuonParticle &l1,
0081              float &deltaR,
0082              float &deltaPhi,
0083              TrajectoryStateOnSurface &propagated) const {
0084     propagated = extrapolate(tk, vtxs);
0085     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : false;
0086   }
0087
0088   /// Try to match one track to one L1. Return true if succeeded (and update deltaR, deltaPhi accordingly)
0089   /// The preselection cut on L1, if specified in the config, is applied before the match
0090   bool match(TrajectoryStateOnSurface &propagated,
0091              const l1extra::L1MuonParticle &l1,
0092              float &deltaR,
0093              float &deltaPhi) const;
0094
0095   // Methods to match with vectors of legacy L1 muon trigger objects
0096
0097   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0098   /// Returns -1 if the match fails
0099   /// The preselection cut on L1, if specified in the config, is applied before the match
0100   int match(const reco::Track &tk,
0101             const std::vector<l1extra::L1MuonParticle> &l1,
0102             float &deltaR,
0103             float &deltaPhi,
0104             TrajectoryStateOnSurface &propagated) const {
0105     propagated = extrapolate(tk);
0106     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
0107   }
0108
0109   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0110   /// Returns -1 if the match fails
0111   /// The preselection cut on L1, if specified in the config, is applied before the match
0112   int match(const reco::Candidate &c,
0113             const std::vector<l1extra::L1MuonParticle> &l1,
0114             float &deltaR,
0115             float &deltaPhi,
0116             TrajectoryStateOnSurface &propagated) const {
0117     propagated = extrapolate(c);
0118     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
0119   }
0120
0121   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0122   /// Returns -1 if the match fails
0123   /// The preselection cut on L1, if specified in the config, is applied before the match
0124   int match(const SimTrack &tk,
0125             const edm::SimVertexContainer &vtxs,
0126             const std::vector<l1extra::L1MuonParticle> &l1,
0127             float &deltaR,
0128             float &deltaPhi,
0129             TrajectoryStateOnSurface &propagated) const {
0130     propagated = extrapolate(tk, vtxs);
0131     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
0132   }
0133
0134   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi accordingly)
0135   /// Returns -1 if the match fails
0136   /// The preselection cut on L1, if specified in the config, is applied before the match
0137   int match(TrajectoryStateOnSurface &propagated,
0138             const std::vector<l1extra::L1MuonParticle> &l1,
0139             float &deltaR,
0140             float &deltaPhi) const;
0141
0142   // Methods to match with vectors of stage2 L1 muon trigger objects
0143
0144   /// Find the best match to stage2 L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0145   /// Returns -1 if the match fails
0146   /// The preselection cut on stage2 L1, if specified in the config, is applied before the match
0147   int match(const reco::Track &tk,
0148             const std::vector<l1t::Muon> &l1,
0149             float &deltaR,
0150             float &deltaPhi,
0151             TrajectoryStateOnSurface &propagated) const {
0152     propagated = extrapolate(tk);
0153     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
0154   }
0155
0156   /// Find the best match to stage2 L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0157   /// Returns -1 if the match fails
0158   /// The preselection cut on stage2 L1, if specified in the config, is applied before the match
0159   int match(const reco::Candidate &c,
0160             const std::vector<l1t::Muon> &l1,
0161             float &deltaR,
0162             float &deltaPhi,
0163             TrajectoryStateOnSurface &propagated) const {
0164     propagated = extrapolate(c);
0165     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
0166   }
0167
0168   /// Find the best match to stage2 L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0169   /// Returns -1 if the match fails
0170   /// The preselection cut on stage 2 L1, if specified in the config, is applied before the match
0171   int match(const SimTrack &tk,
0172             const edm::SimVertexContainer &vtxs,
0173             const std::vector<l1t::Muon> &l1,
0174             float &deltaR,
0175             float &deltaPhi,
0176             TrajectoryStateOnSurface &propagated) const {
0177     propagated = extrapolate(tk, vtxs);
0178     return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
0179   }
0180
0181   /// Find the best match to stage 2 L1, and return its index in the vector (and update deltaR, deltaPhi accordingly)
0182   /// Returns -1 if the match fails
0183   /// The preselection cut on stage 2 L1, if specified in the config, is applied before the match
0184   int match(TrajectoryStateOnSurface &propagated,
0185             const std::vector<l1t::Muon> &l1,
0186             float &deltaR,
0187             float &deltaPhi) const;
0188
0189   // Generic matching
0190
0191   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0192   /// Returns -1 if the match fails
0193   /// Only the objects passing the selector will be allowed for the match.
0194   /// If you don't need a selector, just use an AnySelector (CommonTools/Utils) which accepts everything
0195   template <typename Collection, typename Selector>
0196   int matchGeneric(const reco::Track &tk,
0197                    const Collection &l1,
0198                    const Selector &sel,
0199                    float &deltaR,
0200                    float &deltaPhi,
0201                    TrajectoryStateOnSurface &propagated) const {
0202     propagated = extrapolate(tk);
0203     return propagated.isValid() ? matchGeneric(propagated, l1, sel, deltaR, deltaPhi) : -1;
0204   }
0205
0206   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi and propagated TSOS accordingly)
0207   /// Returns -1 if the match fails
0208   /// Only the objects passing the selector will be allowed for the match.
0209   /// If you don't need a selector, just use an AnySelector (CommonTools/Utils) which accepts everything
0210   template <typename Collection, typename Selector>
0211   int matchGeneric(const reco::Candidate &c,
0212                    const Collection &l1,
0213                    const Selector &sel,
0214                    float &deltaR,
0215                    float &deltaPhi,
0216                    TrajectoryStateOnSurface &propagated) const {
0217     propagated = extrapolate(c);
0218     return propagated.isValid() ? matchGeneric(propagated, l1, sel, deltaR, deltaPhi) : -1;
0219   }
0220
0221   /// Find the best match to L1, and return its index in the vector (and update deltaR, deltaPhi accordingly)
0222   /// Returns -1 if the match fails
0223   /// Only the objects passing the selector will be allowed for the match.
0224   /// The selector defaults to an AnySelector (CommonTools/Utils) which just accepts everything
0225   template <typename Collection, typename Selector>
0226   int matchGeneric(TrajectoryStateOnSurface &propagated,
0227                    const Collection &l1,
0228                    const Selector &sel,
0229                    float &deltaR,
0230                    float &deltaPhi) const;
0231
0232   /// Add this offset to the L1 phi before doing the match, to correct for different scales in L1 vs offline
0233   void setL1PhiOffset(double l1PhiOffset) { l1PhiOffset_ = l1PhiOffset; }
0234
0235 private:
0236   template <class T>
0237   int genericQuality(T const &cand) const {
0238     return 0;
0239   }
0240
0241   int genericQuality(l1extra::L1MuonParticle const &cand) const { return cand.gmtMuonCand().rank(); }
0242   int genericQuality(l1t::Muon const &cand) const { return cand.hwQual(); }
0243
0244   PropagateToMuonSetupT<Tr> propSetup_;
0245   PropagateToMuon prop_;
0246
0247   bool useStage2L1_;
0248
0249   typedef StringCutObjectSelector<reco::Candidate, true> L1Selector;
0250   /// Preselection cut to apply to L1 candidates before matching
0251   L1Selector preselectionCut_;
0252
0253   /// Matching cuts
0254   double deltaR2_, deltaPhi_, deltaEta_;
0255
0256   /// Sort by deltaPhi or deltaEta instead of deltaR
0257   enum SortBy { SortByDeltaR = 0, SortByDeltaPhi, SortByDeltaEta, SortByPt, SortByQual };
0258   SortBy sortBy_;
0259
0260   /// offset to be added to the L1 phi before the match
0261   double l1PhiOffset_;
0262 };
0263
0264 template <edm::Transition Tr>
0265 template <typename Collection, typename Selector>
0266 int L1MuonMatcherAlgoT<Tr>::matchGeneric(TrajectoryStateOnSurface &propagated,
0267                                          const Collection &l1s,
0268                                          const Selector &sel,
0269                                          float &deltaR,
0270                                          float &deltaPhi) const {
0271   typedef typename Collection::value_type obj;
0272
0273   int match = -1;
0274   double minDeltaPhi = deltaPhi_;
0275   double minDeltaEta = deltaEta_;
0276   double minDeltaR2 = deltaR2_;
0277   double maxPt = -1.0;
0278   int maxQual = 0;
0279   GlobalPoint pos = propagated.globalPosition();
0280   for (int i = 0, n = l1s.size(); i < n; ++i) {
0281     const obj &l1 = l1s[i];
0282     if (sel(l1)) {
0283       double thisDeltaPhi = ::deltaPhi(double(pos.phi()), l1.phi() + l1PhiOffset_);
0284       double thisDeltaEta = pos.eta() - l1.eta();
0285       double thisDeltaR2 = ::deltaR2(double(pos.eta()), double(pos.phi()), l1.eta(), l1.phi() + l1PhiOffset_);
0286       int thisQual = genericQuality<obj>(l1);
0287       double thisPt = l1.pt();
0288       if ((fabs(thisDeltaPhi) < deltaPhi_) && (fabs(thisDeltaEta) < deltaEta_) &&
0289           (thisDeltaR2 < deltaR2_)) {  // check all
0290         bool betterMatch = (match == -1);
0291         switch (sortBy_) {
0292           case SortByDeltaR:
0293             betterMatch = (thisDeltaR2 < minDeltaR2);
0294             break;
0295           case SortByDeltaEta:
0296             betterMatch = (fabs(thisDeltaEta) < fabs(minDeltaEta));
0297             break;
0298           case SortByDeltaPhi:
0299             betterMatch = (fabs(thisDeltaPhi) < fabs(minDeltaPhi));
0300             break;
0301           case SortByPt:
0302             betterMatch = (thisPt > maxPt);
0303             break;
0304           // Quality is an int, adding sorting by pT in case of identical qualities
0305           case SortByQual:
0306             betterMatch = (thisQual > maxQual || ((thisQual == maxQual) && (thisPt > maxPt)));
0307             break;
0308         }
0309         if (betterMatch) {  // sort on one
0310           match = i;
0311           deltaR = std::sqrt(thisDeltaR2);
0312           deltaPhi = thisDeltaPhi;
0313           minDeltaR2 = thisDeltaR2;
0314           minDeltaEta = thisDeltaEta;
0315           minDeltaPhi = thisDeltaPhi;
0316           maxQual = thisQual;
0317           maxPt = thisPt;
0318         }
0319       }
0320     }
0321   }
0322   return match;
0323 }
0324
0325 template <edm::Transition Tr>
0326 L1MuonMatcherAlgoT<Tr>::L1MuonMatcherAlgoT(const edm::ParameterSet &iConfig, edm::ConsumesCollector iCollector)
0327     : propSetup_(iConfig, iCollector),
0328       useStage2L1_(iConfig.existsAs<bool>("useStage2L1") ? iConfig.getParameter<bool>("useStage2L1") : false),
0329       preselectionCut_(
0330           (iConfig.existsAs<std::string>("preselection")) ? iConfig.getParameter<std::string>("preselection") : ""),
0331       deltaR2_(std::pow(iConfig.getParameter<double>("maxDeltaR"), 2)),
0332       deltaPhi_(iConfig.existsAs<double>("maxDeltaPhi") ? iConfig.getParameter<double>("maxDeltaPhi") : 10),
0333       deltaEta_(iConfig.existsAs<double>("maxDeltaEta") ? iConfig.getParameter<double>("maxDeltaEta") : 10),
0334       l1PhiOffset_(iConfig.existsAs<double>("l1PhiOffset") ? iConfig.getParameter<double>("l1PhiOffset") : 0) {
0335   bool reqQual = iConfig.existsAs<bool>("sortByQual") && iConfig.getParameter<bool>("sortByQual");
0336   bool reqPhi = iConfig.existsAs<bool>("sortByDeltaPhi") && iConfig.getParameter<bool>("sortByDeltaPhi");
0337   bool reqEta = iConfig.existsAs<bool>("sortByDeltaEta") && iConfig.getParameter<bool>("sortByDeltaEta");
0338   bool reqPt = iConfig.existsAs<bool>("sortByPt") && iConfig.getParameter<bool>("sortByPt");
0339   std::string sortBy = iConfig.existsAs<std::string>("sortBy") ? iConfig.getParameter<std::string>("sortBy") : "";
0340   if (reqPhi + reqEta + reqPt + reqQual > 1)
0341     throw cms::Exception("Configuration")
0342         << "L1MuonMatcherAlgoT: Can't set more than one 'sortBy<XXX>' parameter to True.\n";
0343   if (sortBy == "deltaPhi") {
0344     if (reqEta || reqPt || reqQual)
0345       throw cms::Exception("Configuration") << "L1MuonMatcherAlgoT: Can't set sortBy = 'deltaPhi' and set also another "
0346                                                "'sortBy<XXX>' parameter to True.\n";
0347     else
0348       reqPhi = true;
0349   }
0350   if (sortBy == "deltaEta") {
0351     if (reqPhi || reqPt || reqQual)
0352       throw cms::Exception("Configuration") << "L1MuonMatcherAlgoT: Can't set sortBy = 'deltaEta' and set also another "
0353                                                "'sortBy<XXX>' parameter to True.\n";
0354     else
0355       reqEta = true;
0356   }
0357   if (sortBy == "pt") {
0358     if (reqPhi || reqEta || reqQual)
0359       throw cms::Exception("Configuration")
0360           << "L1MuonMatcherAlgoT: Can't set sortBy = 'pt' and set also another 'sortBy<XXX>' parameter to True.\n";
0361     else
0362       reqPt = true;
0363   }
0364   if (sortBy == "quality") {
0365     if (reqPhi || reqEta || reqPt)
0366       throw cms::Exception("Configuration")
0367           << "L1MuonMatcherAlgoT: Can't set sortBy = 'quality' and set also another 'sortBy<XXX>' parameter to True.\n";
0368     else
0369       reqQual = true;
0370   }
0371   if (sortBy == "deltaR") {
0372     if (reqPhi || reqEta || reqPt || reqQual)
0373       throw cms::Exception("Configuration")
0374           << "L1MuonMatcherAlgoT: Can't set sortBy = 'deltaR' and set also another 'sortBy<XXX>' parameter to True.\n";
0375   }
0376   // so, if we're here there's no ambiguity in what the user may want. either everything is false, or exactly one req is true.
0377   if (reqEta)
0378     sortBy_ = SortByDeltaEta;
0379   else if (reqPhi)
0380     sortBy_ = SortByDeltaPhi;
0381   else if (reqQual)
0382     sortBy_ = SortByQual;
0383   else if (reqPt)
0384     sortBy_ = SortByPt;
0385   else
0386     sortBy_ = SortByDeltaR;
0387 }
0388
0389 template <edm::Transition Tr>
0390 void L1MuonMatcherAlgoT<Tr>::init(const edm::EventSetup &iSetup) {
0391   prop_ = propSetup_.init(iSetup);
0392 }
0393
0394 template <edm::Transition Tr>
0395 bool L1MuonMatcherAlgoT<Tr>::match(TrajectoryStateOnSurface &propagated,
0396                                    const l1extra::L1MuonParticle &l1,
0397                                    float &deltaR,
0398                                    float &deltaPhi) const {
0399   if (preselectionCut_(l1)) {
0400     GlobalPoint pos = propagated.globalPosition();
0401     double thisDeltaPhi = ::deltaPhi(double(pos.phi()), l1.phi() + l1PhiOffset_);
0402     double thisDeltaEta = pos.eta() - l1.eta();
0403     double thisDeltaR2 = ::deltaR2(double(pos.eta()), double(pos.phi()), l1.eta(), l1.phi() + l1PhiOffset_);
0404     if ((fabs(thisDeltaPhi) < deltaPhi_) && (fabs(thisDeltaEta) < deltaEta_) && (thisDeltaR2 < deltaR2_)) {
0405       deltaR = std::sqrt(thisDeltaR2);
0406       deltaPhi = thisDeltaPhi;
0407       return true;
0408     }
0409   }
0410   return false;
0411 }
0412
0413 template <edm::Transition Tr>
0414 int L1MuonMatcherAlgoT<Tr>::match(TrajectoryStateOnSurface &propagated,
0415                                   const std::vector<l1extra::L1MuonParticle> &l1s,
0416                                   float &deltaR,
0417                                   float &deltaPhi) const {
0418   return matchGeneric(propagated, l1s, preselectionCut_, deltaR, deltaPhi);
0419 }
0420
0421 template <edm::Transition Tr>
0422 int L1MuonMatcherAlgoT<Tr>::match(TrajectoryStateOnSurface &propagated,
0423                                   const std::vector<l1t::Muon> &l1s,
0424                                   float &deltaR,
0425                                   float &deltaPhi) const {
0426   return matchGeneric(propagated, l1s, preselectionCut_, deltaR, deltaPhi);
0427 }
0428
0429 using L1MuonMatcherAlgo = L1MuonMatcherAlgoT<edm::Transition::Event>;
0430
0431 #endif