Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-23 12:04:15

0001 /**
0002  * \file L1TMuonDQMOffline.cc
0003  *
0004  * \author J. Pela, C. Battilana
0005  *
0006  * Stage2 Muons implementation: Anna Stakia
0007  *
0008  */
0009 
0010 #include "DQMOffline/L1Trigger/interface/L1TMuonDQMOffline.h"
0011 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0012 
0013 #include <array>
0014 
0015 using namespace reco;
0016 using namespace trigger;
0017 using namespace edm;
0018 using namespace std;
0019 using namespace l1t;
0020 
0021 //__________RECO-GMT Muon Pair Helper Class____________________________
0022 
0023 MuonGmtPair::MuonGmtPair(const reco::Muon* muon,
0024                          const l1t::Muon* regMu,
0025                          const PropagateToMuon& propagator,
0026                          bool useAtVtxCoord)
0027     : m_muon(muon), m_regMu(regMu) {
0028   if (m_regMu) {
0029     if (useAtVtxCoord) {
0030       m_gmtEta = m_regMu->etaAtVtx();
0031       m_gmtPhi = m_regMu->phiAtVtx();
0032     } else {
0033       m_gmtEta = m_regMu->eta();
0034       m_gmtPhi = m_regMu->phi();
0035     }
0036   } else {
0037     m_gmtEta = -5.;
0038     m_gmtPhi = -5.;
0039   }
0040   if (m_muon) {
0041     TrajectoryStateOnSurface trajectory = propagator.extrapolate(*m_muon);
0042     if (trajectory.isValid()) {
0043       m_eta = trajectory.globalPosition().eta();
0044       m_phi = trajectory.globalPosition().phi();
0045     } else {
0046       m_eta = 999.;
0047       m_phi = 999.;
0048     }
0049   } else {
0050     m_eta = 999.;
0051     m_phi = 999.;
0052   }
0053 };
0054 
0055 MuonGmtPair::MuonGmtPair(const MuonGmtPair& muonGmtPair) {
0056   m_muon = muonGmtPair.m_muon;
0057   m_regMu = muonGmtPair.m_regMu;
0058 
0059   m_gmtEta = muonGmtPair.m_gmtEta;
0060   m_gmtPhi = muonGmtPair.m_gmtPhi;
0061 
0062   m_eta = muonGmtPair.m_eta;
0063   m_phi = muonGmtPair.m_phi;
0064 }
0065 
0066 double MuonGmtPair::dR() {
0067   float dEta = m_regMu ? (m_gmtEta - m_eta) : 999.;
0068   float dPhi = m_regMu ? reco::deltaPhi(m_gmtPhi, m_phi) : 999.;
0069   return sqrt(dEta * dEta + dPhi * dPhi);
0070 }
0071 
0072 L1TMuonDQMOffline::EtaRegion MuonGmtPair::etaRegion() const {
0073   if (std::abs(m_eta) < 0.83)
0074     return L1TMuonDQMOffline::kEtaRegionBmtf;
0075   if (std::abs(m_eta) < 1.24)
0076     return L1TMuonDQMOffline::kEtaRegionOmtf;
0077   if (std::abs(m_eta) < 2.4)
0078     return L1TMuonDQMOffline::kEtaRegionEmtf;
0079   return L1TMuonDQMOffline::kEtaRegionOut;
0080 }
0081 
0082 double MuonGmtPair::getDeltaVar(const L1TMuonDQMOffline::ResType type) const {
0083   if (type == L1TMuonDQMOffline::kResPt)
0084     return (gmtPt() - pt()) / pt();
0085   if (type == L1TMuonDQMOffline::kRes1OverPt)
0086     return (pt() - gmtPt()) / gmtPt();  // (1/gmtPt - 1/pt) / (1/pt)
0087   if (type == L1TMuonDQMOffline::kResQOverPt)
0088     return (gmtCharge() * charge() * pt() - gmtPt()) /
0089            gmtPt();  // (gmtCharge/gmtPt - charge/pt) / (charge/pt) with gmtCharge/charge = gmtCharge*charge
0090   if (type == L1TMuonDQMOffline::kResPhi)
0091     return reco::deltaPhi(gmtPhi(), m_phi);
0092   if (type == L1TMuonDQMOffline::kResEta)
0093     return gmtEta() - m_eta;
0094   if (type == L1TMuonDQMOffline::kResCh)
0095     return gmtCharge() - charge();
0096   return -999.;
0097 }
0098 
0099 double MuonGmtPair::getVar(const L1TMuonDQMOffline::EffType type) const {
0100   if (type == L1TMuonDQMOffline::kEffPt)
0101     return pt();
0102   if (type == L1TMuonDQMOffline::kEffPhi)
0103     return m_phi;
0104   if (type == L1TMuonDQMOffline::kEffEta)
0105     return m_eta;
0106   return -999.;
0107 }
0108 
0109 //__________DQM_base_class_______________________________________________
0110 L1TMuonDQMOffline::L1TMuonDQMOffline(const ParameterSet& ps)
0111     : m_propagatorSetup(ps.getParameter<edm::ParameterSet>("muProp"), consumesCollector()),
0112       m_effTypes({kEffPt, kEffPhi, kEffEta, kEffVtx}),
0113       m_resTypes({kResPt, kResQOverPt, kResPhi, kResEta}),
0114       m_etaRegions({kEtaRegionAll, kEtaRegionBmtf, kEtaRegionOmtf, kEtaRegionEmtf}),
0115       m_qualLevelsRes({kQualAll}),
0116       m_effStrings({{kEffPt, "pt"}, {kEffPhi, "phi"}, {kEffEta, "eta"}, {kEffVtx, "vtx"}}),
0117       m_effLabelStrings({{kEffPt, "p_{T} (GeV)"}, {kEffPhi, "#phi"}, {kEffEta, "#eta"}, {kEffVtx, "# vertices"}}),
0118       m_resStrings({{kResPt, "pt"},
0119                     {kRes1OverPt, "1overpt"},
0120                     {kResQOverPt, "qoverpt"},
0121                     {kResPhi, "phi"},
0122                     {kResEta, "eta"},
0123                     {kResCh, "charge"}}),
0124       m_resLabelStrings({{kResPt, "(p_{T}^{L1} - p_{T}^{reco})/p_{T}^{reco}"},
0125                          {kRes1OverPt, "(p_{T}^{reco} - p_{T}^{L1})/p_{T}^{L1}"},
0126                          {kResQOverPt, "(q^{L1}*q^{reco}*p_{T}^{reco} - p_{T}^{L1})/p_{T}^{L1}"},
0127                          {kResPhi, "#phi_{L1} - #phi_{reco}"},
0128                          {kResEta, "#eta_{L1} - #eta_{reco}"},
0129                          {kResCh, "charge^{L1} - charge^{reco}"}}),
0130       m_etaStrings({{kEtaRegionAll, "etaMin0_etaMax2p4"},
0131                     {kEtaRegionBmtf, "etaMin0_etaMax0p83"},
0132                     {kEtaRegionOmtf, "etaMin0p83_etaMax1p24"},
0133                     {kEtaRegionEmtf, "etaMin1p24_etaMax2p4"}}),
0134       m_qualStrings(
0135           {{kQualAll, "qualAll"}, {kQualOpen, "qualOpen"}, {kQualDouble, "qualDouble"}, {kQualSingle, "qualSingle"}}),
0136       m_verbose(ps.getUntrackedParameter<bool>("verbose")),
0137       m_HistFolder(ps.getUntrackedParameter<string>("histFolder")),
0138       m_TagPtCut(ps.getUntrackedParameter<double>("tagPtCut")),
0139       m_recoToL1PtCutFactor(ps.getUntrackedParameter<double>("recoToL1PtCutFactor")),
0140       m_cutsVPSet(ps.getUntrackedParameter<std::vector<edm::ParameterSet>>("cuts")),
0141       m_MuonInputTag(consumes<reco::MuonCollection>(ps.getUntrackedParameter<InputTag>("muonInputTag"))),
0142       m_GmtInputTag(consumes<l1t::MuonBxCollection>(ps.getUntrackedParameter<InputTag>("gmtInputTag"))),
0143       m_VtxInputTag(consumes<VertexCollection>(ps.getUntrackedParameter<InputTag>("vtxInputTag"))),
0144       m_BsInputTag(consumes<BeamSpot>(ps.getUntrackedParameter<InputTag>("bsInputTag"))),
0145       m_trigInputTag(consumes<trigger::TriggerEvent>(ps.getUntrackedParameter<InputTag>("trigInputTag"))),
0146       m_trigProcess(ps.getUntrackedParameter<string>("trigProcess")),
0147       m_trigProcess_token(consumes<edm::TriggerResults>(ps.getUntrackedParameter<InputTag>("trigProcess_token"))),
0148       m_trigNames(ps.getUntrackedParameter<vector<string>>("triggerNames")),
0149       m_effVsPtBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsPtBins")),
0150       m_effVsPhiBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsPhiBins")),
0151       m_effVsEtaBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsEtaBins")),
0152       m_effVsVtxBins(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsVtxBins")),
0153       m_useAtVtxCoord(ps.getUntrackedParameter<bool>("useL1AtVtxCoord")),
0154       m_maxGmtMuonDR(0.3),
0155       m_minTagProbeDR(0.5),
0156       m_maxHltMuonDR(0.1) {
0157   if (m_verbose)
0158     cout << "[L1TMuonDQMOffline:] ____________ Storage initialization ____________ " << endl;
0159 
0160   for (const auto& cutsPSet : m_cutsVPSet) {
0161     const auto qCut = cutsPSet.getUntrackedParameter<int>("qualCut");
0162     QualLevel qLevel = kQualAll;
0163     if (qCut > 11) {
0164       qLevel = kQualSingle;
0165     } else if (qCut > 7) {
0166       qLevel = kQualDouble;
0167     } else if (qCut > 3) {
0168       qLevel = kQualOpen;
0169     }
0170     m_cuts.emplace_back(std::make_pair(cutsPSet.getUntrackedParameter<int>("ptCut"), qLevel));
0171   }
0172 }
0173 
0174 //_____________________________________________________________________
0175 L1TMuonDQMOffline::~L1TMuonDQMOffline() {}
0176 //----------------------------------------------------------------------
0177 void L1TMuonDQMOffline::dqmBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0178   if (m_verbose)
0179     cout << "[L1TMuonDQMOffline:] Called beginRun." << endl;
0180   bool changed = true;
0181   m_hltConfig.init(run, iSetup, m_trigProcess, changed);
0182 }
0183 
0184 //_____________________________________________________________________
0185 void L1TMuonDQMOffline::bookHistograms(DQMStore::IBooker& ibooker, const edm::Run& run, const edm::EventSetup& iSetup) {
0186   //book histos
0187   bookControlHistos(ibooker);
0188   bookEfficiencyHistos(ibooker);
0189   bookResolutionHistos(ibooker);
0190 
0191   vector<string>::const_iterator trigNamesIt = m_trigNames.begin();
0192   vector<string>::const_iterator trigNamesEnd = m_trigNames.end();
0193 
0194   for (; trigNamesIt != trigNamesEnd; ++trigNamesIt) {
0195     TString tNameTmp = TString(*trigNamesIt);  // use TString as it handles regex
0196     TRegexp tNamePattern = TRegexp(tNameTmp, true);
0197     int tIndex = -1;
0198 
0199     for (unsigned ipath = 0; ipath < m_hltConfig.size(); ++ipath) {
0200       TString tmpName = TString(m_hltConfig.triggerName(ipath));
0201       if (tmpName.Contains(tNamePattern)) {
0202         tIndex = int(ipath);
0203         m_trigIndices.push_back(tIndex);
0204       }
0205     }
0206     if (tIndex < 0 && m_verbose)
0207       cout << "[L1TMuonDQMOffline:] Warning: Could not find trigger " << (*trigNamesIt) << endl;
0208   }
0209 }
0210 
0211 //_____________________________________________________________________
0212 void L1TMuonDQMOffline::analyze(const Event& iEvent, const EventSetup& eventSetup) {
0213   auto const propagator = m_propagatorSetup.init(eventSetup);
0214 
0215   Handle<reco::MuonCollection> muons;
0216   iEvent.getByToken(m_MuonInputTag, muons);
0217   Handle<BeamSpot> beamSpot;
0218   iEvent.getByToken(m_BsInputTag, beamSpot);
0219   Handle<VertexCollection> vertex;
0220   iEvent.getByToken(m_VtxInputTag, vertex);
0221   Handle<l1t::MuonBxCollection> gmtCands;
0222   iEvent.getByToken(m_GmtInputTag, gmtCands);
0223   Handle<edm::TriggerResults> trigResults;
0224   iEvent.getByToken(m_trigProcess_token, trigResults);
0225   edm::Handle<trigger::TriggerEvent> trigEvent;
0226   iEvent.getByToken(m_trigInputTag, trigEvent);
0227 
0228   const auto nVtx = getNVertices(vertex);
0229   const Vertex primaryVertex = getPrimaryVertex(vertex, beamSpot);
0230 
0231   getTightMuons(muons, primaryVertex);
0232   getProbeMuons(trigResults, trigEvent);  // CB add flag to run on orthogonal datasets (no T&P)
0233 
0234   getMuonGmtPairs(gmtCands, propagator);
0235 
0236   if (m_verbose)
0237     cout << "[L1TMuonDQMOffline:] Computing efficiencies" << endl;
0238 
0239   vector<MuonGmtPair>::const_iterator muonGmtPairsIt = m_MuonGmtPairs.begin();
0240   vector<MuonGmtPair>::const_iterator muonGmtPairsEnd = m_MuonGmtPairs.end();
0241 
0242   // To fill once for global eta and once for TF eta region of the L1T muon.
0243   // The second entry is a placeholder and will be replaced by the TF eta region of the L1T muon.
0244   std::array<EtaRegion, 2> regsToFill{{kEtaRegionAll, kEtaRegionAll}};
0245 
0246   for (; muonGmtPairsIt != muonGmtPairsEnd; ++muonGmtPairsIt) {
0247     // Fill the resolution histograms
0248     if ((muonGmtPairsIt->etaRegion() != kEtaRegionOut) && (muonGmtPairsIt->gmtPt() > 0)) {
0249       regsToFill[1] = muonGmtPairsIt->etaRegion();
0250       m_histoKeyResType histoKeyRes = {kResPt, kEtaRegionAll, kQualAll};
0251       for (const auto var : m_resTypes) {
0252         const auto varToFill = muonGmtPairsIt->getDeltaVar(var);
0253         std::get<0>(histoKeyRes) = var;
0254         // Fill for the global eta and for TF eta region that the probe muon is in
0255         for (const auto regToFill : regsToFill) {
0256           std::get<1>(histoKeyRes) = regToFill;
0257           for (const auto qualLevel : m_qualLevelsRes) {
0258             // This assumes that the qualLevel enum has increasing qualities
0259             // HW quality levels can be 0, 4, 8, or 12
0260             int qualCut = qualLevel * 4;
0261             if (muonGmtPairsIt->gmtQual() >= qualCut) {
0262               std::get<2>(histoKeyRes) = qualLevel;
0263               m_ResolutionHistos[histoKeyRes]->Fill(varToFill);
0264             }
0265           }
0266         }
0267       }
0268     }
0269 
0270     // Fill the efficiency numerator and denominator histograms
0271     if (muonGmtPairsIt->etaRegion() != kEtaRegionOut) {
0272       unsigned int cutsCounter = 0;
0273       for (const auto& cut : m_cuts) {
0274         const auto gmtPtCut = cut.first;
0275         const auto qualLevel = cut.second;
0276         const bool gmtAboveCut = (muonGmtPairsIt->gmtPt() > gmtPtCut);
0277 
0278         // default keys
0279         m_histoKeyEffDenVarType histoKeyEffDenVar = {kEffPt, gmtPtCut, kEtaRegionAll};
0280         m_histoKeyEffNumVarType histoKeyEffNumVar = {kEffPt, gmtPtCut, kEtaRegionAll, qualLevel};
0281 
0282         regsToFill[1] = muonGmtPairsIt->etaRegion();
0283         for (const auto var : m_effTypes) {
0284           if (var != kEffPt) {
0285             if (muonGmtPairsIt->pt() < m_recoToL1PtCutFactor * gmtPtCut)
0286               break;  // efficiency at plateau
0287           }
0288           double varToFill;
0289           if (var == kEffVtx) {
0290             varToFill = static_cast<double>(nVtx);
0291           } else {
0292             varToFill = muonGmtPairsIt->getVar(var);
0293           }
0294           // Fill denominators
0295           if (var == kEffEta) {
0296             m_EfficiencyDenEtaHistos[gmtPtCut]->Fill(varToFill);
0297           } else {
0298             std::get<0>(histoKeyEffDenVar) = var;
0299             // Fill for the global eta and for TF eta region that the probe muon is in
0300             for (const auto regToFill : regsToFill) {
0301               if (var == kEffPt) {
0302                 if (cutsCounter == 0) {
0303                   m_EfficiencyDenPtHistos[regToFill]->Fill(varToFill);
0304                 }
0305               } else {
0306                 std::get<2>(histoKeyEffDenVar) = regToFill;
0307                 m_EfficiencyDenVarHistos[histoKeyEffDenVar]->Fill(varToFill);
0308               }
0309             }
0310           }
0311           // Fill numerators
0312           std::get<0>(histoKeyEffNumVar) = var;
0313           // This assumes that the qualLevel enum has increasing qualities
0314           if (gmtAboveCut && muonGmtPairsIt->gmtQual() >= qualLevel * 4) {
0315             if (var == kEffEta) {
0316               m_histoKeyEffNumEtaType histoKeyEffNumEta = {gmtPtCut, qualLevel};
0317               m_EfficiencyNumEtaHistos[histoKeyEffNumEta]->Fill(varToFill);
0318             } else {
0319               std::get<3>(histoKeyEffNumVar) = qualLevel;
0320               // Fill for the global eta and for TF eta region that the probe muon is in
0321               for (const auto regToFill : regsToFill) {
0322                 std::get<2>(histoKeyEffNumVar) = regToFill;
0323                 m_EfficiencyNumVarHistos[histoKeyEffNumVar]->Fill(varToFill);
0324               }
0325             }
0326           }
0327         }
0328         ++cutsCounter;
0329       }
0330     }
0331   }
0332 
0333   if (m_verbose)
0334     cout << "[L1TMuonDQMOffline:] Computation finished" << endl;
0335 }
0336 
0337 //_____________________________________________________________________
0338 void L1TMuonDQMOffline::bookControlHistos(DQMStore::IBooker& ibooker) {
0339   if (m_verbose)
0340     cout << "[L1TMuonDQMOffline:] Booking Control Plot Histos" << endl;
0341 
0342   ibooker.setCurrentFolder(m_HistFolder + "/control_variables");
0343 
0344   m_ControlHistos[kCtrlMuonGmtDeltaR] = ibooker.book1D("MuonGmtDeltaR", "MuonGmtDeltaR; #DeltaR", 50, 0., 0.5);
0345   m_ControlHistos[kCtrlNTightVsAll] =
0346       ibooker.book2D("NTightVsAll", "NTightVsAll; # muons; # tight muons", 20, -0.5, 19.5, 16, -0.5, 15.5);
0347   m_ControlHistos[kCtrlNProbesVsTight] =
0348       ibooker.book2D("NProbesVsTight", "NProbesVsTight; # tight muons; # probe muons", 8, -0.5, 7.5, 8, -0.5, 7.5);
0349 
0350   m_ControlHistos[kCtrlTagPt] = ibooker.book1D("TagMuonPt", "TagMuonPt; p_{T}", 50, 0., 100.);
0351   m_ControlHistos[kCtrlTagPhi] = ibooker.book1D("TagMuonPhi", "TagMuonPhi; #phi", 66, -3.3, 3.3);
0352   m_ControlHistos[kCtrlTagEta] = ibooker.book1D("TagMuonEta", "TagMuonEta; #eta", 50, -2.5, 2.5);
0353 
0354   m_ControlHistos[kCtrlProbePt] = ibooker.book1D("ProbeMuonPt", "ProbeMuonPt; p_{T}", 50, 0., 100.);
0355   m_ControlHistos[kCtrlProbePhi] = ibooker.book1D("ProbeMuonPhi", "ProbeMuonPhi; #phi", 66, -3.3, 3.3);
0356   m_ControlHistos[kCtrlProbeEta] = ibooker.book1D("ProbeMuonEta", "ProbeMuonEta; #eta", 50, -2.5, 2.5);
0357 
0358   m_ControlHistos[kCtrlTagProbeDr] =
0359       ibooker.book1D("TagMuonProbeMuonDeltaR", "TagMuonProbeMuonDeltaR; #DeltaR", 50, 0., 5.0);
0360   m_ControlHistos[kCtrlTagHltDr] = ibooker.book1D("TagMuonHltDeltaR", "TagMuonHltDeltaR;#DeltaR", 55, 0., 0.11);
0361 }
0362 
0363 //_____________________________________________________________________
0364 void L1TMuonDQMOffline::bookEfficiencyHistos(DQMStore::IBooker& ibooker) {
0365   ibooker.setCurrentFolder(m_HistFolder + "/numerators_and_denominators");
0366 
0367   for (const auto var : m_effTypes) {
0368     auto histBins = getHistBinsEff(var);
0369     // histograms for eta variable get a special treatment
0370     if (var == kEffEta) {
0371       for (const auto& cut : m_cuts) {
0372         const auto gmtPtCut = cut.first;
0373         const auto qualLevel = cut.second;
0374         std::string name = "effDen_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut);
0375         m_EfficiencyDenEtaHistos[gmtPtCut] =
0376             ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0377         name = "effNum_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" + m_qualStrings[qualLevel];
0378         m_histoKeyEffNumEtaType histoKeyEffNumEta = {gmtPtCut, qualLevel};
0379         m_EfficiencyNumEtaHistos[histoKeyEffNumEta] =
0380             ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0381       }
0382     } else {
0383       for (const auto etaReg : m_etaRegions) {
0384         // denominator histograms for pt variable get a special treatment
0385         if (var == kEffPt) {
0386           std::string name = "effDen_" + m_effStrings[var] + "_" + m_etaStrings[etaReg];
0387           m_EfficiencyDenPtHistos[etaReg] =
0388               ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0389         } else {
0390           for (const auto& cut : m_cuts) {
0391             const int gmtPtCut = cut.first;
0392             std::string name =
0393                 "effDen_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" + m_etaStrings[etaReg];
0394             m_histoKeyEffDenVarType histoKeyEffDenVar = {var, gmtPtCut, etaReg};
0395             m_EfficiencyDenVarHistos[histoKeyEffDenVar] =
0396                 ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0397           }
0398         }
0399         for (const auto& cut : m_cuts) {
0400           const auto gmtPtCut = cut.first;
0401           const auto qualLevel = cut.second;
0402           std::string name = "effNum_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" +
0403                              m_etaStrings[etaReg] + "_" + m_qualStrings[qualLevel];
0404           m_histoKeyEffNumVarType histoKeyEffNum = {var, gmtPtCut, etaReg, qualLevel};
0405           m_EfficiencyNumVarHistos[histoKeyEffNum] =
0406               ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0407         }
0408       }
0409     }
0410   }
0411 }
0412 
0413 void L1TMuonDQMOffline::bookResolutionHistos(DQMStore::IBooker& ibooker) {
0414   if (m_verbose)
0415     cout << "[L1TMuonOffline:] Booking Resolution Plot Histos" << endl;
0416   ibooker.setCurrentFolder(m_HistFolder + "/resolution");
0417 
0418   for (const auto var : m_resTypes) {
0419     auto nbins = std::get<0>(getHistBinsRes(var));
0420     auto xmin = std::get<1>(getHistBinsRes(var));
0421     auto xmax = std::get<2>(getHistBinsRes(var));
0422     for (const auto etaReg : m_etaRegions) {
0423       for (const auto qualLevel : m_qualLevelsRes) {
0424         m_histoKeyResType histoKeyRes = {var, etaReg, qualLevel};
0425         std::string name =
0426             "resolution_" + m_resStrings[var] + "_" + m_etaStrings[etaReg] + "_" + m_qualStrings[qualLevel];
0427         m_ResolutionHistos[histoKeyRes] = ibooker.book1D(name, name + ";" + m_resLabelStrings[var], nbins, xmin, xmax);
0428       }
0429     }
0430   }
0431 }
0432 
0433 //_____________________________________________________________________
0434 const unsigned int L1TMuonDQMOffline::getNVertices(Handle<VertexCollection>& vertex) {
0435   unsigned int nVtx = 0;
0436 
0437   if (vertex.isValid()) {
0438     for (const auto& vertexIt : *vertex) {
0439       if (vertexIt.isValid() && !vertexIt.isFake()) {
0440         ++nVtx;
0441       }
0442     }
0443   }
0444   return nVtx;
0445 }
0446 
0447 //_____________________________________________________________________
0448 const reco::Vertex L1TMuonDQMOffline::getPrimaryVertex(Handle<VertexCollection>& vertex, Handle<BeamSpot>& beamSpot) {
0449   Vertex::Point posVtx;
0450   Vertex::Error errVtx;
0451 
0452   bool hasPrimaryVertex = false;
0453 
0454   if (vertex.isValid()) {
0455     vector<Vertex>::const_iterator vertexIt = vertex->begin();
0456     vector<Vertex>::const_iterator vertexEnd = vertex->end();
0457 
0458     for (; vertexIt != vertexEnd; ++vertexIt) {
0459       if (vertexIt->isValid() && !vertexIt->isFake()) {
0460         posVtx = vertexIt->position();
0461         errVtx = vertexIt->error();
0462         hasPrimaryVertex = true;
0463         break;
0464       }
0465     }
0466   }
0467 
0468   if (!hasPrimaryVertex) {
0469     posVtx = beamSpot->position();
0470     errVtx(0, 0) = beamSpot->BeamWidthX();
0471     errVtx(1, 1) = beamSpot->BeamWidthY();
0472     errVtx(2, 2) = beamSpot->sigmaZ();
0473   }
0474   const Vertex primaryVertex(posVtx, errVtx);
0475   return primaryVertex;
0476 }
0477 
0478 //_____________________________________________________________________
0479 void L1TMuonDQMOffline::getTightMuons(edm::Handle<reco::MuonCollection>& muons, const Vertex& vertex) {
0480   if (m_verbose)
0481     cout << "[L1TMuonDQMOffline:] Getting tight muons" << endl;
0482   m_TightMuons.clear();
0483   MuonCollection::const_iterator muonIt = muons->begin();
0484   MuonCollection::const_iterator muonEnd = muons->end();
0485 
0486   for (; muonIt != muonEnd; ++muonIt) {
0487     if (muon::isTightMuon((*muonIt), vertex)) {
0488       m_TightMuons.push_back(&(*muonIt));
0489     }
0490   }
0491   m_ControlHistos[kCtrlNTightVsAll]->Fill(muons->size(), m_TightMuons.size());
0492 }
0493 
0494 //_____________________________________________________________________
0495 void L1TMuonDQMOffline::getProbeMuons(Handle<edm::TriggerResults>& trigResults,
0496                                       edm::Handle<trigger::TriggerEvent>& trigEvent) {
0497   if (m_verbose)
0498     cout << "[L1TMuonDQMOffline:] getting probe muons" << endl;
0499   m_ProbeMuons.clear();
0500   std::vector<const reco::Muon*> tagMuonsInHist;
0501 
0502   tagMuonsInHist.clear();
0503 
0504   vector<const reco::Muon*>::const_iterator probeCandIt = m_TightMuons.begin();
0505   vector<const reco::Muon*>::const_iterator tightMuonsEnd = m_TightMuons.end();
0506 
0507   for (; probeCandIt != tightMuonsEnd; ++probeCandIt) {
0508     bool isProbe = false;
0509     vector<const reco::Muon*>::const_iterator tagCandIt = m_TightMuons.begin();
0510     float deltar = 0.;
0511 
0512     for (; tagCandIt != tightMuonsEnd; ++tagCandIt) {
0513       bool tagMuonAlreadyInHist = false;
0514       bool tagHasTrig = false;
0515       float eta = (*tagCandIt)->eta();
0516       float phi = (*tagCandIt)->phi();
0517       float pt = (*tagCandIt)->pt();
0518       float dEta = eta - (*probeCandIt)->eta();
0519       float dPhi = phi - (*probeCandIt)->phi();
0520       deltar = sqrt(dEta * dEta + dPhi * dPhi);
0521 
0522       if ((*tagCandIt) == (*probeCandIt) || deltar < m_minTagProbeDR)
0523         continue;  // CB has a little bias for closed-by muons
0524       auto matchHltDeltaR = matchHlt(trigEvent, (*tagCandIt));
0525       tagHasTrig = (matchHltDeltaR < m_maxHltMuonDR) && (pt > m_TagPtCut);
0526       isProbe |= tagHasTrig;
0527       if (tagHasTrig) {
0528         if (std::distance(m_TightMuons.begin(), m_TightMuons.end()) > 2) {
0529           for (vector<const reco::Muon*>::const_iterator tagMuonsInHistIt = tagMuonsInHist.begin();
0530                tagMuonsInHistIt != tagMuonsInHist.end();
0531                ++tagMuonsInHistIt) {
0532             if ((*tagCandIt) == (*tagMuonsInHistIt)) {
0533               tagMuonAlreadyInHist = true;
0534               break;
0535             }
0536           }
0537           if (tagMuonAlreadyInHist == false)
0538             tagMuonsInHist.push_back((*tagCandIt));
0539         }
0540         if (tagMuonAlreadyInHist == false) {
0541           m_ControlHistos[kCtrlTagEta]->Fill(eta);
0542           m_ControlHistos[kCtrlTagPhi]->Fill(phi);
0543           m_ControlHistos[kCtrlTagPt]->Fill(pt);
0544           m_ControlHistos[kCtrlTagProbeDr]->Fill(deltar);
0545           m_ControlHistos[kCtrlTagHltDr]->Fill(matchHltDeltaR);
0546         }
0547       }
0548     }
0549     if (isProbe)
0550       m_ProbeMuons.push_back((*probeCandIt));
0551   }
0552   m_ControlHistos[kCtrlNProbesVsTight]->Fill(m_TightMuons.size(), m_ProbeMuons.size());
0553 }
0554 
0555 //_____________________________________________________________________
0556 void L1TMuonDQMOffline::getMuonGmtPairs(edm::Handle<l1t::MuonBxCollection>& gmtCands,
0557                                         PropagateToMuon const& propagator) {
0558   m_MuonGmtPairs.clear();
0559   if (m_verbose)
0560     cout << "[L1TMuonDQMOffline:] Getting muon GMT pairs" << endl;
0561 
0562   vector<const reco::Muon*>::const_iterator probeMuIt = m_ProbeMuons.begin();
0563   vector<const reco::Muon*>::const_iterator probeMuEnd = m_ProbeMuons.end();
0564 
0565   l1t::MuonBxCollection::const_iterator gmtIt;
0566   l1t::MuonBxCollection::const_iterator gmtEnd = gmtCands->end(0);
0567 
0568   for (; probeMuIt != probeMuEnd; ++probeMuIt) {
0569     MuonGmtPair pairBestCand((*probeMuIt), nullptr, propagator, m_useAtVtxCoord);
0570 
0571     // Fill the control histograms with the probe muon kinematic variables used
0572     m_ControlHistos[kCtrlProbeEta]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffEta));
0573     m_ControlHistos[kCtrlProbePhi]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffPhi));
0574     m_ControlHistos[kCtrlProbePt]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffPt));
0575 
0576     gmtIt = gmtCands->begin(0);  // use only on L1T muons from BX 0
0577 
0578     for (; gmtIt != gmtEnd; ++gmtIt) {
0579       MuonGmtPair pairTmpCand((*probeMuIt), &(*gmtIt), propagator, m_useAtVtxCoord);
0580 
0581       if ((pairTmpCand.dR() < m_maxGmtMuonDR) && (pairTmpCand.dR() < pairBestCand.dR())) {
0582         pairBestCand = pairTmpCand;
0583       }
0584     }
0585     m_MuonGmtPairs.push_back(pairBestCand);
0586     m_ControlHistos[kCtrlMuonGmtDeltaR]->Fill(pairBestCand.dR());
0587   }
0588 }
0589 
0590 //_____________________________________________________________________
0591 double L1TMuonDQMOffline::matchHlt(edm::Handle<TriggerEvent>& triggerEvent, const reco::Muon* mu) {
0592   double matchDeltaR = 9999;
0593 
0594   TriggerObjectCollection trigObjs = triggerEvent->getObjects();
0595 
0596   vector<int>::const_iterator trigIndexIt = m_trigIndices.begin();
0597   vector<int>::const_iterator trigIndexEnd = m_trigIndices.end();
0598 
0599   for (; trigIndexIt != trigIndexEnd; ++trigIndexIt) {
0600     const vector<string> moduleLabels(m_hltConfig.moduleLabels(*trigIndexIt));
0601     const unsigned moduleIndex = m_hltConfig.size((*trigIndexIt)) - 2;
0602     const unsigned hltFilterIndex = triggerEvent->filterIndex(InputTag(moduleLabels[moduleIndex], "", m_trigProcess));
0603 
0604     if (hltFilterIndex < triggerEvent->sizeFilters()) {
0605       const Keys triggerKeys(triggerEvent->filterKeys(hltFilterIndex));
0606       const Vids triggerVids(triggerEvent->filterIds(hltFilterIndex));
0607       const unsigned nTriggers = triggerVids.size();
0608       for (size_t iTrig = 0; iTrig < nTriggers; ++iTrig) {
0609         const TriggerObject trigObject = trigObjs[triggerKeys[iTrig]];
0610         double dRtmp = deltaR((*mu), trigObject);
0611         if (dRtmp < matchDeltaR)
0612           matchDeltaR = dRtmp;
0613       }
0614     }
0615   }
0616   return matchDeltaR;
0617 }
0618 
0619 std::vector<float> L1TMuonDQMOffline::getHistBinsEff(EffType eff) {
0620   if (eff == kEffPt) {
0621     std::vector<float> effVsPtBins(m_effVsPtBins.begin(), m_effVsPtBins.end());
0622     return effVsPtBins;
0623   }
0624   if (eff == kEffPhi) {
0625     std::vector<float> effVsPhiBins(m_effVsPhiBins.begin(), m_effVsPhiBins.end());
0626     return effVsPhiBins;
0627   }
0628   if (eff == kEffEta) {
0629     std::vector<float> effVsEtaBins(m_effVsEtaBins.begin(), m_effVsEtaBins.end());
0630     return effVsEtaBins;
0631   }
0632   if (eff == kEffVtx) {
0633     std::vector<float> effVsVtxBins(m_effVsVtxBins.begin(), m_effVsVtxBins.end());
0634     return effVsVtxBins;
0635   }
0636   return {0., 1.};
0637 }
0638 
0639 std::tuple<int, double, double> L1TMuonDQMOffline::getHistBinsRes(ResType res) {
0640   if (res == kResPt)
0641     return {50, -2., 2.};
0642   if (res == kRes1OverPt)
0643     return {50, -2., 2.};
0644   if (res == kResQOverPt)
0645     return {50, -2., 2.};
0646   if (res == kResPhi)
0647     return {96, -0.2, 0.2};
0648   if (res == kResEta)
0649     return {100, -0.1, 0.1};
0650   if (res == kResCh)
0651     return {5, -2, 3};
0652   return {1, 0, 1};
0653 }
0654 
0655 //define this as a plug-in
0656 DEFINE_FWK_MODULE(L1TMuonDQMOffline);