L1TMuonDQMOffline.cc

Go to the documentation of this file.

#include "DQMOffline/L1Trigger/interface/L1TMuonDQMOffline.h"
#include "DataFormats/MuonReco/interface/MuonSelectors.h"

#include <array>

using namespace reco;
using namespace trigger;
using namespace edm;
using namespace std;
using namespace l1t;

//__________RECO-GMT Muon Pair Helper Class____________________________

MuonGmtPair::MuonGmtPair(const reco::Muon* muon,
                         const l1t::Muon* regMu,
                         const PropagateToMuon& propagator,
                         bool useAtVtxCoord)
    : m_muon(muon), m_regMu(regMu) {
  if (m_regMu) {
    if (useAtVtxCoord) {
      m_gmtEta = m_regMu->etaAtVtx();
      m_gmtPhi = m_regMu->phiAtVtx();
    } else {
      m_gmtEta = m_regMu->eta();
      m_gmtPhi = m_regMu->phi();
    }
  } else {
    m_gmtEta = -5.;
    m_gmtPhi = -5.;
  }
  if (m_muon) {
    TrajectoryStateOnSurface trajectory = propagator.extrapolate(*m_muon);
    if (trajectory.isValid()) {
      m_eta = trajectory.globalPosition().eta();
      m_phi = trajectory.globalPosition().phi();
    } else {
      m_eta = 999.;
      m_phi = 999.;
    }
  } else {
    m_eta = 999.;
    m_phi = 999.;
  }
};

MuonGmtPair::MuonGmtPair(const MuonGmtPair& muonGmtPair) {
  m_muon = muonGmtPair.m_muon;
  m_regMu = muonGmtPair.m_regMu;

  m_gmtEta = muonGmtPair.m_gmtEta;
  m_gmtPhi = muonGmtPair.m_gmtPhi;

  m_eta = muonGmtPair.m_eta;
  m_phi = muonGmtPair.m_phi;
}

double MuonGmtPair::dR() {
  float dEta = m_regMu ? (m_gmtEta - m_eta) : 999.;
  float dPhi = m_regMu ? reco::deltaPhi(m_gmtPhi, m_phi) : 999.;
  return sqrt(dEta * dEta + dPhi * dPhi);
}

L1TMuonDQMOffline::EtaRegion MuonGmtPair::etaRegion() const {
  if (std::abs(m_eta) < 0.83)
    return L1TMuonDQMOffline::kEtaRegionBmtf;
  if (std::abs(m_eta) < 1.24)
    return L1TMuonDQMOffline::kEtaRegionOmtf;
  if (std::abs(m_eta) < 2.4)
    return L1TMuonDQMOffline::kEtaRegionEmtf;
  return L1TMuonDQMOffline::kEtaRegionOut;
}

double MuonGmtPair::getDeltaVar(const L1TMuonDQMOffline::ResType type) const {
  if (type == L1TMuonDQMOffline::kResPt)
    return (gmtPt() - pt()) / pt();
  if (type == L1TMuonDQMOffline::kRes1OverPt)
    return (pt() - gmtPt()) / gmtPt();  // (1/gmtPt - 1/pt) / (1/pt)
  if (type == L1TMuonDQMOffline::kResQOverPt)
    return (gmtCharge() * charge() * pt() - gmtPt()) /
           gmtPt();  // (gmtCharge/gmtPt - charge/pt) / (charge/pt) with gmtCharge/charge = gmtCharge*charge
  if (type == L1TMuonDQMOffline::kResPhi)
    return reco::deltaPhi(gmtPhi(), m_phi);
  if (type == L1TMuonDQMOffline::kResEta)
    return gmtEta() - m_eta;
  if (type == L1TMuonDQMOffline::kResCh)
    return gmtCharge() - charge();
  return -999.;
}

double MuonGmtPair::getVar(const L1TMuonDQMOffline::EffType type) const {
  if (type == L1TMuonDQMOffline::kEffPt)
    return pt();
  if (type == L1TMuonDQMOffline::kEffPhi)
    return m_phi;
  if (type == L1TMuonDQMOffline::kEffEta)
    return m_eta;
  return -999.;
}

//__________DQM_base_class_______________________________________________
L1TMuonDQMOffline::L1TMuonDQMOffline(const ParameterSet& ps)
    : m_propagator(ps.getParameter<edm::ParameterSet>("muProp"), consumesCollector()),
      m_effTypes({kEffPt, kEffPhi, kEffEta, kEffVtx}),
      m_resTypes({kResPt, kResQOverPt, kResPhi, kResEta}),
      m_etaRegions({kEtaRegionAll, kEtaRegionBmtf, kEtaRegionOmtf, kEtaRegionEmtf}),
      m_qualLevelsRes({kQualAll}),
      m_effStrings({{kEffPt, "pt"}, {kEffPhi, "phi"}, {kEffEta, "eta"}, {kEffVtx, "vtx"}}),
      m_effLabelStrings({{kEffPt, "p_{T} (GeV)"}, {kEffPhi, "#phi"}, {kEffEta, "#eta"}, {kEffVtx, "# vertices"}}),
      m_resStrings({{kResPt, "pt"},
                    {kRes1OverPt, "1overpt"},
                    {kResQOverPt, "qoverpt"},
                    {kResPhi, "phi"},
                    {kResEta, "eta"},
                    {kResCh, "charge"}}),
      m_resLabelStrings({{kResPt, "(p_{T}^{L1} - p_{T}^{reco})/p_{T}^{reco}"},
                         {kRes1OverPt, "(p_{T}^{reco} - p_{T}^{L1})/p_{T}^{L1}"},
                         {kResQOverPt, "(q^{L1}*q^{reco}*p_{T}^{reco} - p_{T}^{L1})/p_{T}^{L1}"},
                         {kResPhi, "#phi_{L1} - #phi_{reco}"},
                         {kResEta, "#eta_{L1} - #eta_{reco}"},
                         {kResCh, "charge^{L1} - charge^{reco}"}}),
      m_etaStrings({{kEtaRegionAll, "etaMin0_etaMax2p4"},
                    {kEtaRegionBmtf, "etaMin0_etaMax0p83"},
                    {kEtaRegionOmtf, "etaMin0p83_etaMax1p24"},
                    {kEtaRegionEmtf, "etaMin1p24_etaMax2p4"}}),
      m_qualStrings(
          {{kQualAll, "qualAll"}, {kQualOpen, "qualOpen"}, {kQualDouble, "qualDouble"}, {kQualSingle, "qualSingle"}}),
      m_verbose(ps.getUntrackedParameter<bool>("verbose")),
      m_HistFolder(ps.getUntrackedParameter<string>("histFolder")),
      m_TagPtCut(ps.getUntrackedParameter<double>("tagPtCut")),
      m_recoToL1PtCutFactor(ps.getUntrackedParameter<double>("recoToL1PtCutFactor")),
      m_cutsVPSet(ps.getUntrackedParameter<std::vector<edm::ParameterSet>>("cuts")),
      m_MuonInputTag(consumes<reco::MuonCollection>(ps.getUntrackedParameter<InputTag>("muonInputTag"))),
      m_GmtInputTag(consumes<l1t::MuonBxCollection>(ps.getUntrackedParameter<InputTag>("gmtInputTag"))),
      m_VtxInputTag(consumes<VertexCollection>(ps.getUntrackedParameter<InputTag>("vtxInputTag"))),
      m_BsInputTag(consumes<BeamSpot>(ps.getUntrackedParameter<InputTag>("bsInputTag"))),
      m_trigInputTag(consumes<trigger::TriggerEvent>(ps.getUntrackedParameter<InputTag>("trigInputTag"))),
      m_trigProcess(ps.getUntrackedParameter<string>("trigProcess")),
      m_trigProcess_token(consumes<edm::TriggerResults>(ps.getUntrackedParameter<InputTag>("trigProcess_token"))),
      m_trigNames(ps.getUntrackedParameter<vector<string>>("triggerNames")),
      m_effVsPtBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsPtBins")),
      m_effVsPhiBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsPhiBins")),
      m_effVsEtaBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsEtaBins")),
      m_effVsVtxBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsVtxBins")),
      m_useAtVtxCoord(ps.getUntrackedParameter<bool>("useL1AtVtxCoord")),
      m_maxGmtMuonDR(0.3),
      m_minTagProbeDR(0.5),
      m_maxHltMuonDR(0.1) {
  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] ____________ Storage initialization ____________ " << endl;

  for (const auto& cutsPSet : m_cutsVPSet) {
    const auto qCut = cutsPSet.getUntrackedParameter<int>("qualCut");
    QualLevel qLevel = kQualAll;
    if (qCut > 11) {
      qLevel = kQualSingle;
    } else if (qCut > 7) {
      qLevel = kQualDouble;
    } else if (qCut > 3) {
      qLevel = kQualOpen;
    }
    m_cuts.emplace_back(std::make_pair(cutsPSet.getUntrackedParameter<int>("ptCut"), qLevel));
  }
}

//_____________________________________________________________________
L1TMuonDQMOffline::~L1TMuonDQMOffline() {}
//----------------------------------------------------------------------
void L1TMuonDQMOffline::dqmBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] Called beginRun." << endl;
  bool changed = true;
  m_hltConfig.init(run, iSetup, m_trigProcess, changed);
}

//_____________________________________________________________________
void L1TMuonDQMOffline::bookHistograms(DQMStore::IBooker& ibooker, const edm::Run& run, const edm::EventSetup& iSetup) {
  //book histos
  bookControlHistos(ibooker);
  bookEfficiencyHistos(ibooker);
  bookResolutionHistos(ibooker);

  vector<string>::const_iterator trigNamesIt = m_trigNames.begin();
  vector<string>::const_iterator trigNamesEnd = m_trigNames.end();

  for (; trigNamesIt != trigNamesEnd; ++trigNamesIt) {
    TString tNameTmp = TString(*trigNamesIt);  // use TString as it handles regex
    TRegexp tNamePattern = TRegexp(tNameTmp, true);
    int tIndex = -1;

    for (unsigned ipath = 0; ipath < m_hltConfig.size(); ++ipath) {
      TString tmpName = TString(m_hltConfig.triggerName(ipath));
      if (tmpName.Contains(tNamePattern)) {
        tIndex = int(ipath);
        m_trigIndices.push_back(tIndex);
      }
    }
    if (tIndex < 0 && m_verbose)
      cout << "[L1TMuonDQMOffline:] Warning: Could not find trigger " << (*trigNamesIt) << endl;
  }
}

//_____________________________________________________________________
void L1TMuonDQMOffline::analyze(const Event& iEvent, const EventSetup& eventSetup) {
  m_propagator.init(eventSetup);

  Handle<reco::MuonCollection> muons;
  iEvent.getByToken(m_MuonInputTag, muons);
  Handle<BeamSpot> beamSpot;
  iEvent.getByToken(m_BsInputTag, beamSpot);
  Handle<VertexCollection> vertex;
  iEvent.getByToken(m_VtxInputTag, vertex);
  Handle<l1t::MuonBxCollection> gmtCands;
  iEvent.getByToken(m_GmtInputTag, gmtCands);
  Handle<edm::TriggerResults> trigResults;
  iEvent.getByToken(m_trigProcess_token, trigResults);
  edm::Handle<trigger::TriggerEvent> trigEvent;
  iEvent.getByToken(m_trigInputTag, trigEvent);

  const auto nVtx = getNVertices(vertex);
  const Vertex primaryVertex = getPrimaryVertex(vertex, beamSpot);

  getTightMuons(muons, primaryVertex);
  getProbeMuons(trigResults, trigEvent);  // CB add flag to run on orthogonal datasets (no T&P)

  getMuonGmtPairs(gmtCands);

  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] Computing efficiencies" << endl;

  vector<MuonGmtPair>::const_iterator muonGmtPairsIt = m_MuonGmtPairs.begin();
  vector<MuonGmtPair>::const_iterator muonGmtPairsEnd = m_MuonGmtPairs.end();

  // To fill once for global eta and once for TF eta region of the L1T muon.
  // The second entry is a placeholder and will be replaced by the TF eta region of the L1T muon.
  std::array<EtaRegion, 2> regsToFill{{kEtaRegionAll, kEtaRegionAll}};

  for (; muonGmtPairsIt != muonGmtPairsEnd; ++muonGmtPairsIt) {
    // Fill the resolution histograms
    if ((muonGmtPairsIt->etaRegion() != kEtaRegionOut) && (muonGmtPairsIt->gmtPt() > 0)) {
      regsToFill[1] = muonGmtPairsIt->etaRegion();
      m_histoKeyResType histoKeyRes = {kResPt, kEtaRegionAll, kQualAll};
      for (const auto var : m_resTypes) {
        const auto varToFill = muonGmtPairsIt->getDeltaVar(var);
        std::get<0>(histoKeyRes) = var;
        // Fill for the global eta and for TF eta region that the probe muon is in
        for (const auto regToFill : regsToFill) {
          std::get<1>(histoKeyRes) = regToFill;
          for (const auto qualLevel : m_qualLevelsRes) {
            // This assumes that the qualLevel enum has increasing qualities
            // HW quality levels can be 0, 4, 8, or 12
            int qualCut = qualLevel * 4;
            if (muonGmtPairsIt->gmtQual() >= qualCut) {
              std::get<2>(histoKeyRes) = qualLevel;
              m_ResolutionHistos[histoKeyRes]->Fill(varToFill);
            }
          }
        }
      }
    }

    // Fill the efficiency numerator and denominator histograms
    if (muonGmtPairsIt->etaRegion() != kEtaRegionOut) {
      unsigned int cutsCounter = 0;
      for (const auto& cut : m_cuts) {
        const auto gmtPtCut = cut.first;
        const auto qualLevel = cut.second;
        const bool gmtAboveCut = (muonGmtPairsIt->gmtPt() > gmtPtCut);

        // default keys
        m_histoKeyEffDenVarType histoKeyEffDenVar = {kEffPt, gmtPtCut, kEtaRegionAll};
        m_histoKeyEffNumVarType histoKeyEffNumVar = {kEffPt, gmtPtCut, kEtaRegionAll, qualLevel};

        regsToFill[1] = muonGmtPairsIt->etaRegion();
        for (const auto var : m_effTypes) {
          if (var != kEffPt) {
            if (muonGmtPairsIt->pt() < m_recoToL1PtCutFactor * gmtPtCut)
              break;  // efficiency at plateau
          }
          double varToFill;
          if (var == kEffVtx) {
            varToFill = static_cast<double>(nVtx);
          } else {
            varToFill = muonGmtPairsIt->getVar(var);
          }
          // Fill denominators
          if (var == kEffEta) {
            m_EfficiencyDenEtaHistos[gmtPtCut]->Fill(varToFill);
          } else {
            std::get<0>(histoKeyEffDenVar) = var;
            // Fill for the global eta and for TF eta region that the probe muon is in
            for (const auto regToFill : regsToFill) {
              if (var == kEffPt) {
                if (cutsCounter == 0) {
                  m_EfficiencyDenPtHistos[regToFill]->Fill(varToFill);
                }
              } else {
                std::get<2>(histoKeyEffDenVar) = regToFill;
                m_EfficiencyDenVarHistos[histoKeyEffDenVar]->Fill(varToFill);
              }
            }
          }
          // Fill numerators
          std::get<0>(histoKeyEffNumVar) = var;
          // This assumes that the qualLevel enum has increasing qualities
          if (gmtAboveCut && muonGmtPairsIt->gmtQual() >= qualLevel * 4) {
            if (var == kEffEta) {
              m_histoKeyEffNumEtaType histoKeyEffNumEta = {gmtPtCut, qualLevel};
              m_EfficiencyNumEtaHistos[histoKeyEffNumEta]->Fill(varToFill);
            } else {
              std::get<3>(histoKeyEffNumVar) = qualLevel;
              // Fill for the global eta and for TF eta region that the probe muon is in
              for (const auto regToFill : regsToFill) {
                std::get<2>(histoKeyEffNumVar) = regToFill;
                m_EfficiencyNumVarHistos[histoKeyEffNumVar]->Fill(varToFill);
              }
            }
          }
        }
        ++cutsCounter;
      }
    }
  }

  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] Computation finished" << endl;
}

//_____________________________________________________________________
void L1TMuonDQMOffline::bookControlHistos(DQMStore::IBooker& ibooker) {
  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] Booking Control Plot Histos" << endl;

  ibooker.setCurrentFolder(m_HistFolder + "/control_variables");

  m_ControlHistos[kCtrlMuonGmtDeltaR] = ibooker.book1D("MuonGmtDeltaR", "MuonGmtDeltaR; #DeltaR", 50, 0., 0.5);
  m_ControlHistos[kCtrlNTightVsAll] =
      ibooker.book2D("NTightVsAll", "NTightVsAll; # muons; # tight muons", 20, -0.5, 19.5, 16, -0.5, 15.5);
  m_ControlHistos[kCtrlNProbesVsTight] =
      ibooker.book2D("NProbesVsTight", "NProbesVsTight; # tight muons; # probe muons", 8, -0.5, 7.5, 8, -0.5, 7.5);

  m_ControlHistos[kCtrlTagPt] = ibooker.book1D("TagMuonPt", "TagMuonPt; p_{T}", 50, 0., 100.);
  m_ControlHistos[kCtrlTagPhi] = ibooker.book1D("TagMuonPhi", "TagMuonPhi; #phi", 66, -3.3, 3.3);
  m_ControlHistos[kCtrlTagEta] = ibooker.book1D("TagMuonEta", "TagMuonEta; #eta", 50, -2.5, 2.5);

  m_ControlHistos[kCtrlProbePt] = ibooker.book1D("ProbeMuonPt", "ProbeMuonPt; p_{T}", 50, 0., 100.);
  m_ControlHistos[kCtrlProbePhi] = ibooker.book1D("ProbeMuonPhi", "ProbeMuonPhi; #phi", 66, -3.3, 3.3);
  m_ControlHistos[kCtrlProbeEta] = ibooker.book1D("ProbeMuonEta", "ProbeMuonEta; #eta", 50, -2.5, 2.5);

  m_ControlHistos[kCtrlTagProbeDr] =
      ibooker.book1D("TagMuonProbeMuonDeltaR", "TagMuonProbeMuonDeltaR; #DeltaR", 50, 0., 5.0);
  m_ControlHistos[kCtrlTagHltDr] = ibooker.book1D("TagMuonHltDeltaR", "TagMuonHltDeltaR;#DeltaR", 55, 0., 0.11);
}

//_____________________________________________________________________
void L1TMuonDQMOffline::bookEfficiencyHistos(DQMStore::IBooker& ibooker) {
  ibooker.setCurrentFolder(m_HistFolder + "/numerators_and_denominators");

  for (const auto var : m_effTypes) {
    auto histBins = getHistBinsEff(var);
    // histograms for eta variable get a special treatment
    if (var == kEffEta) {
      for (const auto& cut : m_cuts) {
        const auto gmtPtCut = cut.first;
        const auto qualLevel = cut.second;
        std::string name = "effDen_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut);
        m_EfficiencyDenEtaHistos[gmtPtCut] =
            ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
        name = "effNum_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" + m_qualStrings[qualLevel];
        m_histoKeyEffNumEtaType histoKeyEffNumEta = {gmtPtCut, qualLevel};
        m_EfficiencyNumEtaHistos[histoKeyEffNumEta] =
            ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
      }
    } else {
      for (const auto etaReg : m_etaRegions) {
        // denominator histograms for pt variable get a special treatment
        if (var == kEffPt) {
          std::string name = "effDen_" + m_effStrings[var] + "_" + m_etaStrings[etaReg];
          m_EfficiencyDenPtHistos[etaReg] =
              ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
        } else {
          for (const auto& cut : m_cuts) {
            const int gmtPtCut = cut.first;
            std::string name =
                "effDen_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" + m_etaStrings[etaReg];
            m_histoKeyEffDenVarType histoKeyEffDenVar = {var, gmtPtCut, etaReg};
            m_EfficiencyDenVarHistos[histoKeyEffDenVar] =
                ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
          }
        }
        for (const auto& cut : m_cuts) {
          const auto gmtPtCut = cut.first;
          const auto qualLevel = cut.second;
          std::string name = "effNum_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" +
                             m_etaStrings[etaReg] + "_" + m_qualStrings[qualLevel];
          m_histoKeyEffNumVarType histoKeyEffNum = {var, gmtPtCut, etaReg, qualLevel};
          m_EfficiencyNumVarHistos[histoKeyEffNum] =
              ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
        }
      }
    }
  }
}

void L1TMuonDQMOffline::bookResolutionHistos(DQMStore::IBooker& ibooker) {
  if (m_verbose)
    cout << "[L1TMuonOffline:] Booking Resolution Plot Histos" << endl;
  ibooker.setCurrentFolder(m_HistFolder + "/resolution");

  for (const auto var : m_resTypes) {
    auto nbins = std::get<0>(getHistBinsRes(var));
    auto xmin = std::get<1>(getHistBinsRes(var));
    auto xmax = std::get<2>(getHistBinsRes(var));
    for (const auto etaReg : m_etaRegions) {
      for (const auto qualLevel : m_qualLevelsRes) {
        m_histoKeyResType histoKeyRes = {var, etaReg, qualLevel};
        std::string name =
            "resolution_" + m_resStrings[var] + "_" + m_etaStrings[etaReg] + "_" + m_qualStrings[qualLevel];
        m_ResolutionHistos[histoKeyRes] = ibooker.book1D(name, name + ";" + m_resLabelStrings[var], nbins, xmin, xmax);
      }
    }
  }
}

//_____________________________________________________________________
const unsigned int L1TMuonDQMOffline::getNVertices(Handle<VertexCollection>& vertex) {
  unsigned int nVtx = 0;

  if (vertex.isValid()) {
    for (const auto& vertexIt : *vertex) {
      if (vertexIt.isValid() && !vertexIt.isFake()) {
        ++nVtx;
      }
    }
  }
  return nVtx;
}

//_____________________________________________________________________
const reco::Vertex L1TMuonDQMOffline::getPrimaryVertex(Handle<VertexCollection>& vertex, Handle<BeamSpot>& beamSpot) {
  Vertex::Point posVtx;
  Vertex::Error errVtx;

  bool hasPrimaryVertex = false;

  if (vertex.isValid()) {
    vector<Vertex>::const_iterator vertexIt = vertex->begin();
    vector<Vertex>::const_iterator vertexEnd = vertex->end();

    for (; vertexIt != vertexEnd; ++vertexIt) {
      if (vertexIt->isValid() && !vertexIt->isFake()) {
        posVtx = vertexIt->position();
        errVtx = vertexIt->error();
        hasPrimaryVertex = true;
        break;
      }
    }
  }

  if (!hasPrimaryVertex) {
    posVtx = beamSpot->position();
    errVtx(0, 0) = beamSpot->BeamWidthX();
    errVtx(1, 1) = beamSpot->BeamWidthY();
    errVtx(2, 2) = beamSpot->sigmaZ();
  }
  const Vertex primaryVertex(posVtx, errVtx);
  return primaryVertex;
}

//_____________________________________________________________________
void L1TMuonDQMOffline::getTightMuons(edm::Handle<reco::MuonCollection>& muons, const Vertex& vertex) {
  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] Getting tight muons" << endl;
  m_TightMuons.clear();
  MuonCollection::const_iterator muonIt = muons->begin();
  MuonCollection::const_iterator muonEnd = muons->end();

  for (; muonIt != muonEnd; ++muonIt) {
    if (muon::isTightMuon((*muonIt), vertex)) {
      m_TightMuons.push_back(&(*muonIt));
    }
  }
  m_ControlHistos[kCtrlNTightVsAll]->Fill(muons->size(), m_TightMuons.size());
}

//_____________________________________________________________________
void L1TMuonDQMOffline::getProbeMuons(Handle<edm::TriggerResults>& trigResults,
                                      edm::Handle<trigger::TriggerEvent>& trigEvent) {
  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] getting probe muons" << endl;
  m_ProbeMuons.clear();
  std::vector<const reco::Muon*> tagMuonsInHist;

  tagMuonsInHist.clear();

  vector<const reco::Muon*>::const_iterator probeCandIt = m_TightMuons.begin();
  vector<const reco::Muon*>::const_iterator tightMuonsEnd = m_TightMuons.end();

  for (; probeCandIt != tightMuonsEnd; ++probeCandIt) {
    bool isProbe = false;
    vector<const reco::Muon*>::const_iterator tagCandIt = m_TightMuons.begin();
    float deltar = 0.;

    for (; tagCandIt != tightMuonsEnd; ++tagCandIt) {
      bool tagMuonAlreadyInHist = false;
      bool tagHasTrig = false;
      float eta = (*tagCandIt)->eta();
      float phi = (*tagCandIt)->phi();
      float pt = (*tagCandIt)->pt();
      float dEta = eta - (*probeCandIt)->eta();
      float dPhi = phi - (*probeCandIt)->phi();
      deltar = sqrt(dEta * dEta + dPhi * dPhi);

      if ((*tagCandIt) == (*probeCandIt) || deltar < m_minTagProbeDR)
        continue;  // CB has a little bias for closed-by muons
      auto matchHltDeltaR = matchHlt(trigEvent, (*tagCandIt));
      tagHasTrig = (matchHltDeltaR < m_maxHltMuonDR) && (pt > m_TagPtCut);
      isProbe |= tagHasTrig;
      if (tagHasTrig) {
        if (std::distance(m_TightMuons.begin(), m_TightMuons.end()) > 2) {
          for (vector<const reco::Muon*>::const_iterator tagMuonsInHistIt = tagMuonsInHist.begin();
               tagMuonsInHistIt != tagMuonsInHist.end();
               ++tagMuonsInHistIt) {
            if ((*tagCandIt) == (*tagMuonsInHistIt)) {
              tagMuonAlreadyInHist = true;
              break;
            }
          }
          if (tagMuonAlreadyInHist == false)
            tagMuonsInHist.push_back((*tagCandIt));
        }
        if (tagMuonAlreadyInHist == false) {
          m_ControlHistos[kCtrlTagEta]->Fill(eta);
          m_ControlHistos[kCtrlTagPhi]->Fill(phi);
          m_ControlHistos[kCtrlTagPt]->Fill(pt);
          m_ControlHistos[kCtrlTagProbeDr]->Fill(deltar);
          m_ControlHistos[kCtrlTagHltDr]->Fill(matchHltDeltaR);
        }
      }
    }
    if (isProbe)
      m_ProbeMuons.push_back((*probeCandIt));
  }
  m_ControlHistos[kCtrlNProbesVsTight]->Fill(m_TightMuons.size(), m_ProbeMuons.size());
}

//_____________________________________________________________________
void L1TMuonDQMOffline::getMuonGmtPairs(edm::Handle<l1t::MuonBxCollection>& gmtCands) {
  m_MuonGmtPairs.clear();
  if (m_verbose)
    cout << "[L1TMuonDQMOffline:] Getting muon GMT pairs" << endl;

  vector<const reco::Muon*>::const_iterator probeMuIt = m_ProbeMuons.begin();
  vector<const reco::Muon*>::const_iterator probeMuEnd = m_ProbeMuons.end();

  l1t::MuonBxCollection::const_iterator gmtIt;
  l1t::MuonBxCollection::const_iterator gmtEnd = gmtCands->end(0);

  for (; probeMuIt != probeMuEnd; ++probeMuIt) {
    MuonGmtPair pairBestCand((*probeMuIt), nullptr, m_propagator, m_useAtVtxCoord);

    // Fill the control histograms with the probe muon kinematic variables used
    m_ControlHistos[kCtrlProbeEta]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffEta));
    m_ControlHistos[kCtrlProbePhi]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffPhi));
    m_ControlHistos[kCtrlProbePt]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffPt));

    gmtIt = gmtCands->begin(0);  // use only on L1T muons from BX 0

    for (; gmtIt != gmtEnd; ++gmtIt) {
      MuonGmtPair pairTmpCand((*probeMuIt), &(*gmtIt), m_propagator, m_useAtVtxCoord);

      if ((pairTmpCand.dR() < m_maxGmtMuonDR) && (pairTmpCand.dR() < pairBestCand.dR())) {
        pairBestCand = pairTmpCand;
      }
    }
    m_MuonGmtPairs.push_back(pairBestCand);
    m_ControlHistos[kCtrlMuonGmtDeltaR]->Fill(pairBestCand.dR());
  }
}

//_____________________________________________________________________
double L1TMuonDQMOffline::matchHlt(edm::Handle<TriggerEvent>& triggerEvent, const reco::Muon* mu) {
  double matchDeltaR = 9999;

  TriggerObjectCollection trigObjs = triggerEvent->getObjects();

  vector<int>::const_iterator trigIndexIt = m_trigIndices.begin();
  vector<int>::const_iterator trigIndexEnd = m_trigIndices.end();

  for (; trigIndexIt != trigIndexEnd; ++trigIndexIt) {
    const vector<string> moduleLabels(m_hltConfig.moduleLabels(*trigIndexIt));
    const unsigned moduleIndex = m_hltConfig.size((*trigIndexIt)) - 2;
    const unsigned hltFilterIndex = triggerEvent->filterIndex(InputTag(moduleLabels[moduleIndex], "", m_trigProcess));

    if (hltFilterIndex < triggerEvent->sizeFilters()) {
      const Keys triggerKeys(triggerEvent->filterKeys(hltFilterIndex));
      const Vids triggerVids(triggerEvent->filterIds(hltFilterIndex));
      const unsigned nTriggers = triggerVids.size();
      for (size_t iTrig = 0; iTrig < nTriggers; ++iTrig) {
        const TriggerObject trigObject = trigObjs[triggerKeys[iTrig]];
        double dRtmp = deltaR((*mu), trigObject);
        if (dRtmp < matchDeltaR)
          matchDeltaR = dRtmp;
      }
    }
  }
  return matchDeltaR;
}

std::vector<float> L1TMuonDQMOffline::getHistBinsEff(EffType eff) {
  if (eff == kEffPt) {
    std::vector<float> effVsPtBins(m_effVsPtBins.begin(), m_effVsPtBins.end());
    return effVsPtBins;
  }
  if (eff == kEffPhi) {
    std::vector<float> effVsPhiBins(m_effVsPhiBins.begin(), m_effVsPhiBins.end());
    return effVsPhiBins;
  }
  if (eff == kEffEta) {
    std::vector<float> effVsEtaBins(m_effVsEtaBins.begin(), m_effVsEtaBins.end());
    return effVsEtaBins;
  }
  if (eff == kEffVtx) {
    std::vector<float> effVsVtxBins(m_effVsVtxBins.begin(), m_effVsVtxBins.end());
    return effVsVtxBins;
  }
  return {0., 1.};
}

std::tuple<int, double, double> L1TMuonDQMOffline::getHistBinsRes(ResType res) {
  if (res == kResPt)
    return {50, -2., 2.};
  if (res == kRes1OverPt)
    return {50, -2., 2.};
  if (res == kResQOverPt)
    return {50, -2., 2.};
  if (res == kResPhi)
    return {96, -0.2, 0.2};
  if (res == kResEta)
    return {100, -0.1, 0.1};
  if (res == kResCh)
    return {5, -2, 3};
  return {1, 0, 1};
}

//define this as a plug-in
DEFINE_FWK_MODULE(L1TMuonDQMOffline);