L1MuonMatcherAlgo.h

Go to the documentation of this file.

#ifndef MuonAnalysis_MuonAssociators_interface_L1MuonMatcherAlgo_h
#define MuonAnalysis_MuonAssociators_interface_L1MuonMatcherAlgo_h

#include <cmath>
#include <type_traits>
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/Math/interface/deltaPhi.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
#include "DataFormats/L1Trigger/interface/Muon.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
#include "CommonTools/Utils/interface/AnySelector.h"
#include "MuonAnalysis/MuonAssociators/interface/PropagateToMuonSetup.h"

template <edm::Transition Tr>
class L1MuonMatcherAlgoT {
public:
  explicit L1MuonMatcherAlgoT(const edm::ParameterSet &iConfig, edm::ConsumesCollector);
  ~L1MuonMatcherAlgoT() = default;

  void init(const edm::EventSetup &iSetup);

  TrajectoryStateOnSurface extrapolate(const reco::Track &tk) const { return prop_.extrapolate(tk); }

  TrajectoryStateOnSurface extrapolate(const reco::Candidate &tk) const { return prop_.extrapolate(tk); }

  TrajectoryStateOnSurface extrapolate(const SimTrack &tk, const edm::SimVertexContainer &vtx) const {
    return prop_.extrapolate(tk, vtx);
  }

  TrajectoryStateOnSurface extrapolate(const FreeTrajectoryState &state) const { return prop_.extrapolate(state); }

  PropagateToMuon &propagatorToMuon() { return prop_; }
  const PropagateToMuon &propagatorToMuon() const { return prop_; }

  bool match(const reco::Track &tk,
             const l1extra::L1MuonParticle &l1,
             float &deltaR,
             float &deltaPhi,
             TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : false;
  }

  bool match(const reco::Candidate &c,
             const l1extra::L1MuonParticle &l1,
             float &deltaR,
             float &deltaPhi,
             TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(c);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : false;
  }

  bool match(const SimTrack &tk,
             const edm::SimVertexContainer &vtxs,
             const l1extra::L1MuonParticle &l1,
             float &deltaR,
             float &deltaPhi,
             TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk, vtxs);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : false;
  }

  bool match(TrajectoryStateOnSurface &propagated,
             const l1extra::L1MuonParticle &l1,
             float &deltaR,
             float &deltaPhi) const;

  // Methods to match with vectors of legacy L1 muon trigger objects

  int match(const reco::Track &tk,
            const std::vector<l1extra::L1MuonParticle> &l1,
            float &deltaR,
            float &deltaPhi,
            TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
  }

  int match(const reco::Candidate &c,
            const std::vector<l1extra::L1MuonParticle> &l1,
            float &deltaR,
            float &deltaPhi,
            TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(c);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
  }

  int match(const SimTrack &tk,
            const edm::SimVertexContainer &vtxs,
            const std::vector<l1extra::L1MuonParticle> &l1,
            float &deltaR,
            float &deltaPhi,
            TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk, vtxs);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
  }

  int match(TrajectoryStateOnSurface &propagated,
            const std::vector<l1extra::L1MuonParticle> &l1,
            float &deltaR,
            float &deltaPhi) const;

  // Methods to match with vectors of stage2 L1 muon trigger objects

  int match(const reco::Track &tk,
            const std::vector<l1t::Muon> &l1,
            float &deltaR,
            float &deltaPhi,
            TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
  }

  int match(const reco::Candidate &c,
            const std::vector<l1t::Muon> &l1,
            float &deltaR,
            float &deltaPhi,
            TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(c);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
  }

  int match(const SimTrack &tk,
            const edm::SimVertexContainer &vtxs,
            const std::vector<l1t::Muon> &l1,
            float &deltaR,
            float &deltaPhi,
            TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk, vtxs);
    return propagated.isValid() ? match(propagated, l1, deltaR, deltaPhi) : -1;
  }

  int match(TrajectoryStateOnSurface &propagated,
            const std::vector<l1t::Muon> &l1,
            float &deltaR,
            float &deltaPhi) const;

  // Generic matching

  template <typename Collection, typename Selector>
  int matchGeneric(const reco::Track &tk,
                   const Collection &l1,
                   const Selector &sel,
                   float &deltaR,
                   float &deltaPhi,
                   TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(tk);
    return propagated.isValid() ? matchGeneric(propagated, l1, sel, deltaR, deltaPhi) : -1;
  }

  template <typename Collection, typename Selector>
  int matchGeneric(const reco::Candidate &c,
                   const Collection &l1,
                   const Selector &sel,
                   float &deltaR,
                   float &deltaPhi,
                   TrajectoryStateOnSurface &propagated) const {
    propagated = extrapolate(c);
    return propagated.isValid() ? matchGeneric(propagated, l1, sel, deltaR, deltaPhi) : -1;
  }

  template <typename Collection, typename Selector>
  int matchGeneric(TrajectoryStateOnSurface &propagated,
                   const Collection &l1,
                   const Selector &sel,
                   float &deltaR,
                   float &deltaPhi) const;

  void setL1PhiOffset(double l1PhiOffset) { l1PhiOffset_ = l1PhiOffset; }

private:
  template <class T>
  int genericQuality(T const &cand) const {
    return 0;
  }

  int genericQuality(l1extra::L1MuonParticle const &cand) const { return cand.gmtMuonCand().rank(); }
  int genericQuality(l1t::Muon const &cand) const { return cand.hwQual(); }

  PropagateToMuonSetupT<Tr> propSetup_;
  PropagateToMuon prop_;

  bool useStage2L1_;

  typedef StringCutObjectSelector<reco::Candidate, true> L1Selector;
  L1Selector preselectionCut_;

  double deltaR2_, deltaPhi_, deltaEta_;

  enum SortBy { SortByDeltaR = 0, SortByDeltaPhi, SortByDeltaEta, SortByPt, SortByQual };
  SortBy sortBy_;

  double l1PhiOffset_;
};

template <edm::Transition Tr>
template <typename Collection, typename Selector>
int L1MuonMatcherAlgoT<Tr>::matchGeneric(TrajectoryStateOnSurface &propagated,
                                         const Collection &l1s,
                                         const Selector &sel,
                                         float &deltaR,
                                         float &deltaPhi) const {
  typedef typename Collection::value_type obj;

  int match = -1;
  double minDeltaPhi = deltaPhi_;
  double minDeltaEta = deltaEta_;
  double minDeltaR2 = deltaR2_;
  double maxPt = -1.0;
  int maxQual = 0;
  GlobalPoint pos = propagated.globalPosition();
  for (int i = 0, n = l1s.size(); i < n; ++i) {
    const obj &l1 = l1s[i];
    if (sel(l1)) {
      double thisDeltaPhi = ::deltaPhi(double(pos.phi()), l1.phi() + l1PhiOffset_);
      double thisDeltaEta = pos.eta() - l1.eta();
      double thisDeltaR2 = ::deltaR2(double(pos.eta()), double(pos.phi()), l1.eta(), l1.phi() + l1PhiOffset_);
      int thisQual = genericQuality<obj>(l1);
      double thisPt = l1.pt();
      if ((fabs(thisDeltaPhi) < deltaPhi_) && (fabs(thisDeltaEta) < deltaEta_) &&
          (thisDeltaR2 < deltaR2_)) {  // check all
        bool betterMatch = (match == -1);
        switch (sortBy_) {
          case SortByDeltaR:
            betterMatch = (thisDeltaR2 < minDeltaR2);
            break;
          case SortByDeltaEta:
            betterMatch = (fabs(thisDeltaEta) < fabs(minDeltaEta));
            break;
          case SortByDeltaPhi:
            betterMatch = (fabs(thisDeltaPhi) < fabs(minDeltaPhi));
            break;
          case SortByPt:
            betterMatch = (thisPt > maxPt);
            break;
          // Quality is an int, adding sorting by pT in case of identical qualities
          case SortByQual:
            betterMatch = (thisQual > maxQual || ((thisQual == maxQual) && (thisPt > maxPt)));
            break;
        }
        if (betterMatch) {  // sort on one
          match = i;
          deltaR = std::sqrt(thisDeltaR2);
          deltaPhi = thisDeltaPhi;
          minDeltaR2 = thisDeltaR2;
          minDeltaEta = thisDeltaEta;
          minDeltaPhi = thisDeltaPhi;
          maxQual = thisQual;
          maxPt = thisPt;
        }
      }
    }
  }
  return match;
}

template <edm::Transition Tr>
L1MuonMatcherAlgoT<Tr>::L1MuonMatcherAlgoT(const edm::ParameterSet &iConfig, edm::ConsumesCollector iCollector)
    : propSetup_(iConfig, iCollector),
      useStage2L1_(iConfig.existsAs<bool>("useStage2L1") ? iConfig.getParameter<bool>("useStage2L1") : false),
      preselectionCut_(
          (iConfig.existsAs<std::string>("preselection")) ? iConfig.getParameter<std::string>("preselection") : ""),
      deltaR2_(std::pow(iConfig.getParameter<double>("maxDeltaR"), 2)),
      deltaPhi_(iConfig.existsAs<double>("maxDeltaPhi") ? iConfig.getParameter<double>("maxDeltaPhi") : 10),
      deltaEta_(iConfig.existsAs<double>("maxDeltaEta") ? iConfig.getParameter<double>("maxDeltaEta") : 10),
      l1PhiOffset_(iConfig.existsAs<double>("l1PhiOffset") ? iConfig.getParameter<double>("l1PhiOffset") : 0) {
  bool reqQual = iConfig.existsAs<bool>("sortByQual") && iConfig.getParameter<bool>("sortByQual");
  bool reqPhi = iConfig.existsAs<bool>("sortByDeltaPhi") && iConfig.getParameter<bool>("sortByDeltaPhi");
  bool reqEta = iConfig.existsAs<bool>("sortByDeltaEta") && iConfig.getParameter<bool>("sortByDeltaEta");
  bool reqPt = iConfig.existsAs<bool>("sortByPt") && iConfig.getParameter<bool>("sortByPt");
  std::string sortBy = iConfig.existsAs<std::string>("sortBy") ? iConfig.getParameter<std::string>("sortBy") : "";
  if (reqPhi + reqEta + reqPt + reqQual > 1)
    throw cms::Exception("Configuration")
        << "L1MuonMatcherAlgoT: Can't set more than one 'sortBy<XXX>' parameter to True.\n";
  if (sortBy == "deltaPhi") {
    if (reqEta || reqPt || reqQual)
      throw cms::Exception("Configuration") << "L1MuonMatcherAlgoT: Can't set sortBy = 'deltaPhi' and set also another "
                                               "'sortBy<XXX>' parameter to True.\n";
    else
      reqPhi = true;
  }
  if (sortBy == "deltaEta") {
    if (reqPhi || reqPt || reqQual)
      throw cms::Exception("Configuration") << "L1MuonMatcherAlgoT: Can't set sortBy = 'deltaEta' and set also another "
                                               "'sortBy<XXX>' parameter to True.\n";
    else
      reqEta = true;
  }
  if (sortBy == "pt") {
    if (reqPhi || reqEta || reqQual)
      throw cms::Exception("Configuration")
          << "L1MuonMatcherAlgoT: Can't set sortBy = 'pt' and set also another 'sortBy<XXX>' parameter to True.\n";
    else
      reqPt = true;
  }
  if (sortBy == "quality") {
    if (reqPhi || reqEta || reqPt)
      throw cms::Exception("Configuration")
          << "L1MuonMatcherAlgoT: Can't set sortBy = 'quality' and set also another 'sortBy<XXX>' parameter to True.\n";
    else
      reqQual = true;
  }
  if (sortBy == "deltaR") {
    if (reqPhi || reqEta || reqPt || reqQual)
      throw cms::Exception("Configuration")
          << "L1MuonMatcherAlgoT: Can't set sortBy = 'deltaR' and set also another 'sortBy<XXX>' parameter to True.\n";
  }
  // 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.
  if (reqEta)
    sortBy_ = SortByDeltaEta;
  else if (reqPhi)
    sortBy_ = SortByDeltaPhi;
  else if (reqQual)
    sortBy_ = SortByQual;
  else if (reqPt)
    sortBy_ = SortByPt;
  else
    sortBy_ = SortByDeltaR;
}

template <edm::Transition Tr>
void L1MuonMatcherAlgoT<Tr>::init(const edm::EventSetup &iSetup) {
  prop_ = propSetup_.init(iSetup);
}

template <edm::Transition Tr>
bool L1MuonMatcherAlgoT<Tr>::match(TrajectoryStateOnSurface &propagated,
                                   const l1extra::L1MuonParticle &l1,
                                   float &deltaR,
                                   float &deltaPhi) const {
  if (preselectionCut_(l1)) {
    GlobalPoint pos = propagated.globalPosition();
    double thisDeltaPhi = ::deltaPhi(double(pos.phi()), l1.phi() + l1PhiOffset_);
    double thisDeltaEta = pos.eta() - l1.eta();
    double thisDeltaR2 = ::deltaR2(double(pos.eta()), double(pos.phi()), l1.eta(), l1.phi() + l1PhiOffset_);
    if ((fabs(thisDeltaPhi) < deltaPhi_) && (fabs(thisDeltaEta) < deltaEta_) && (thisDeltaR2 < deltaR2_)) {
      deltaR = std::sqrt(thisDeltaR2);
      deltaPhi = thisDeltaPhi;
      return true;
    }
  }
  return false;
}

template <edm::Transition Tr>
int L1MuonMatcherAlgoT<Tr>::match(TrajectoryStateOnSurface &propagated,
                                  const std::vector<l1extra::L1MuonParticle> &l1s,
                                  float &deltaR,
                                  float &deltaPhi) const {
  return matchGeneric(propagated, l1s, preselectionCut_, deltaR, deltaPhi);
}

template <edm::Transition Tr>
int L1MuonMatcherAlgoT<Tr>::match(TrajectoryStateOnSurface &propagated,
                                  const std::vector<l1t::Muon> &l1s,
                                  float &deltaR,
                                  float &deltaPhi) const {
  return matchGeneric(propagated, l1s, preselectionCut_, deltaR, deltaPhi);
}

using L1MuonMatcherAlgo = L1MuonMatcherAlgoT<edm::Transition::Event>;

#endif