Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:35 

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     edm::LogInfo("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     edm::LogInfo("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         if (m_verbose)
0205           edm::LogInfo("L1TMuonDQMOffline") << "Found trigger " << tmpName << " with index " << tIndex << endl;
0206       }
0207     }
0208     if (tIndex < 0 && m_verbose)
0209       edm::LogInfo("L1TMuonDQMOffline") << "Warning: Could not find trigger " << (*trigNamesIt) << endl;
0210   }
0211 }
0212 
0213 //_____________________________________________________________________
0214 void L1TMuonDQMOffline::analyze(const Event& iEvent, const EventSetup& eventSetup) {
0215   if (m_verbose)
0216     edm::LogInfo("L1TMuonDQMOffline") << "\nProcessing new event in analyze" << endl;
0217 
0218   auto const propagator = m_propagatorSetup.init(eventSetup);
0219 
0220   Handle<reco::MuonCollection> muons;
0221   iEvent.getByToken(m_MuonInputTag, muons);
0222   Handle<BeamSpot> beamSpot;
0223   iEvent.getByToken(m_BsInputTag, beamSpot);
0224   Handle<VertexCollection> vertex;
0225   iEvent.getByToken(m_VtxInputTag, vertex);
0226   Handle<l1t::MuonBxCollection> gmtCands;
0227   iEvent.getByToken(m_GmtInputTag, gmtCands);
0228   Handle<edm::TriggerResults> trigResults;
0229   iEvent.getByToken(m_trigProcess_token, trigResults);
0230   edm::Handle<trigger::TriggerEvent> trigEvent;
0231   iEvent.getByToken(m_trigInputTag, trigEvent);
0232 
0233   // Only process event if at least one trigger fired
0234   bool pass_HLT = false;
0235   for (const auto& tIndex : m_trigIndices) {
0236     if (trigResults->accept(tIndex)) {
0237       if (m_verbose)
0238         edm::LogInfo("L1TMuonDQMOffline") << "Fired trigger " << m_hltConfig.triggerName(tIndex) << endl;
0239       pass_HLT = true;
0240       break;
0241     }
0242   }
0243   if (not pass_HLT) {
0244     if (m_verbose)
0245       edm::LogInfo("L1TMuonDQMOffline") << "Did not fire any triggers - quitting" << endl;
0246     return;
0247   }
0248 
0249   const auto nVtx = getNVertices(vertex);
0250   const Vertex primaryVertex = getPrimaryVertex(vertex, beamSpot);
0251 
0252   getTightMuons(muons, primaryVertex);
0253   getProbeMuons(trigResults, trigEvent);  // CB add flag to run on orthogonal datasets (no T&P)
0254 
0255   getMuonGmtPairs(gmtCands, propagator);
0256 
0257   if (m_verbose)
0258     edm::LogInfo("L1TMuonDQMOffline") << "Computing efficiencies with " << m_MuonGmtPairs.size() << " muonGmtPairs, "
0259                                       << m_TightMuons.size() << " tight muons, " << m_ProbeMuons.size()
0260                                       << " probe muons" << endl;
0261 
0262   vector<MuonGmtPair>::const_iterator muonGmtPairsIt = m_MuonGmtPairs.begin();
0263   vector<MuonGmtPair>::const_iterator muonGmtPairsEnd = m_MuonGmtPairs.end();
0264 
0265   // To fill once for global eta and once for TF eta region of the L1T muon.
0266   // The second entry is a placeholder and will be replaced by the TF eta region of the L1T muon.
0267   std::array<EtaRegion, 2> regsToFill{{kEtaRegionAll, kEtaRegionAll}};
0268 
0269   for (; muonGmtPairsIt != muonGmtPairsEnd; ++muonGmtPairsIt) {
0270     // Fill the resolution histograms
0271     if ((muonGmtPairsIt->etaRegion() != kEtaRegionOut) && (muonGmtPairsIt->gmtPt() > 0)) {
0272       regsToFill[1] = muonGmtPairsIt->etaRegion();
0273       m_histoKeyResType histoKeyRes = {kResPt, kEtaRegionAll, kQualAll};
0274       for (const auto var : m_resTypes) {
0275         const auto varToFill = muonGmtPairsIt->getDeltaVar(var);
0276         std::get<0>(histoKeyRes) = var;
0277         // Fill for the global eta and for TF eta region that the probe muon is in
0278         for (const auto regToFill : regsToFill) {
0279           std::get<1>(histoKeyRes) = regToFill;
0280           for (const auto qualLevel : m_qualLevelsRes) {
0281             // This assumes that the qualLevel enum has increasing qualities
0282             // HW quality levels can be 0, 4, 8, or 12
0283             int qualCut = qualLevel * 4;
0284             if (muonGmtPairsIt->gmtQual() >= qualCut) {
0285               std::get<2>(histoKeyRes) = qualLevel;
0286               m_ResolutionHistos[histoKeyRes]->Fill(varToFill);
0287               if (m_verbose)
0288                 edm::LogInfo("L1TMuonDQMOffline")
0289                     << "Filled resolution histo[" << std::get<0>(histoKeyRes) << ", " << std::get<1>(histoKeyRes)
0290                     << ", " << std::get<2>(histoKeyRes) << "] with " << varToFill << endl;
0291             }
0292           }
0293         }
0294       }
0295     }
0296 
0297     // Fill the efficiency numerator and denominator histograms
0298     if (muonGmtPairsIt->etaRegion() != kEtaRegionOut) {
0299       unsigned int cutsCounter = 0;
0300       for (const auto& cut : m_cuts) {
0301         const auto gmtPtCut = cut.first;
0302         const auto qualLevel = cut.second;
0303         const bool gmtAboveCut = (muonGmtPairsIt->gmtPt() > gmtPtCut);
0304 
0305         // default keys
0306         m_histoKeyEffDenVarType histoKeyEffDenVar = {kEffPt, gmtPtCut, kEtaRegionAll};
0307         m_histoKeyEffNumVarType histoKeyEffNumVar = {kEffPt, gmtPtCut, kEtaRegionAll, qualLevel};
0308 
0309         regsToFill[1] = muonGmtPairsIt->etaRegion();
0310         for (const auto var : m_effTypes) {
0311           if (var != kEffPt) {
0312             if (muonGmtPairsIt->pt() < m_recoToL1PtCutFactor * gmtPtCut)
0313               break;  // efficiency at plateau
0314           }
0315           double varToFill;
0316           if (var == kEffVtx) {
0317             varToFill = static_cast<double>(nVtx);
0318           } else {
0319             varToFill = muonGmtPairsIt->getVar(var);
0320           }
0321           // Fill denominators
0322           if (var == kEffEta) {
0323             m_EfficiencyDenEtaHistos[gmtPtCut]->Fill(varToFill);
0324             if (m_verbose)
0325               edm::LogInfo("L1TMuonDQMOffline")
0326                   << "Filled eff denom eta histo[" << gmtPtCut << "] with " << varToFill << endl;
0327           } else {
0328             std::get<0>(histoKeyEffDenVar) = var;
0329             // Fill for the global eta and for TF eta region that the probe muon is in
0330             for (const auto regToFill : regsToFill) {
0331               if (var == kEffPt) {
0332                 if (cutsCounter == 0) {
0333                   m_EfficiencyDenPtHistos[regToFill]->Fill(varToFill);
0334                   if (m_verbose)
0335                     edm::LogInfo("L1TMuonDQMOffline")
0336                         << "Filled eff denom pT histo[" << regToFill << "] with " << varToFill << endl;
0337                 }
0338               } else {
0339                 std::get<2>(histoKeyEffDenVar) = regToFill;
0340                 m_EfficiencyDenVarHistos[histoKeyEffDenVar]->Fill(varToFill);
0341                 if (m_verbose)
0342                   edm::LogInfo("L1TMuonDQMOffline") << "Filled eff denom histo[" << std::get<0>(histoKeyEffDenVar)
0343                                                     << ", " << std::get<1>(histoKeyEffDenVar) << ", "
0344                                                     << std::get<2>(histoKeyEffDenVar) << "] with " << varToFill << endl;
0345               }
0346             }
0347           }
0348           // Fill numerators
0349           std::get<0>(histoKeyEffNumVar) = var;
0350           // This assumes that the qualLevel enum has increasing qualities
0351           if (gmtAboveCut && muonGmtPairsIt->gmtQual() >= qualLevel * 4) {
0352             if (var == kEffEta) {
0353               m_histoKeyEffNumEtaType histoKeyEffNumEta = {gmtPtCut, qualLevel};
0354               m_EfficiencyNumEtaHistos[histoKeyEffNumEta]->Fill(varToFill);
0355               if (m_verbose)
0356                 edm::LogInfo("L1TMuonDQMOffline")
0357                     << "Filled eff num eta histo[" << std::get<0>(histoKeyEffNumEta) << ", "
0358                     << std::get<1>(histoKeyEffNumEta) << "] with " << varToFill << endl;
0359             } else {
0360               std::get<3>(histoKeyEffNumVar) = qualLevel;
0361               // Fill for the global eta and for TF eta region that the probe muon is in
0362               for (const auto regToFill : regsToFill) {
0363                 std::get<2>(histoKeyEffNumVar) = regToFill;
0364                 m_EfficiencyNumVarHistos[histoKeyEffNumVar]->Fill(varToFill);
0365                 if (m_verbose)
0366                   edm::LogInfo("L1TMuonDQMOffline")
0367                       << "Filled eff num histo[" << std::get<0>(histoKeyEffNumVar) << ", "
0368                       << std::get<1>(histoKeyEffNumVar) << ", " << std::get<2>(histoKeyEffNumVar) << ", "
0369                       << std::get<3>(histoKeyEffNumVar) << "] with " << varToFill << endl;
0370               }
0371             }
0372           }
0373         }
0374         ++cutsCounter;
0375       }
0376     }
0377   }
0378 
0379   if (m_verbose)
0380     edm::LogInfo("L1TMuonDQMOffline") << "Computation finished" << endl;
0381 }
0382 
0383 //_____________________________________________________________________
0384 void L1TMuonDQMOffline::bookControlHistos(DQMStore::IBooker& ibooker) {
0385   if (m_verbose)
0386     edm::LogInfo("L1TMuonDQMOffline") << "Booking Control Plot Histos" << endl;
0387 
0388   ibooker.setCurrentFolder(m_HistFolder + "/control_variables");
0389 
0390   m_ControlHistos[kCtrlMuonGmtDeltaR] = ibooker.book1D("MuonGmtDeltaR", "MuonGmtDeltaR; #DeltaR", 50, 0., 0.5);
0391   m_ControlHistos[kCtrlNTightVsAll] =
0392       ibooker.book2D("NTightVsAll", "NTightVsAll; # muons; # tight muons", 20, -0.5, 19.5, 16, -0.5, 15.5);
0393   m_ControlHistos[kCtrlNProbesVsTight] =
0394       ibooker.book2D("NProbesVsTight", "NProbesVsTight; # tight muons; # probe muons", 8, -0.5, 7.5, 8, -0.5, 7.5);
0395 
0396   m_ControlHistos[kCtrlTagPt] = ibooker.book1D("TagMuonPt", "TagMuonPt; p_{T}", 50, 0., 100.);
0397   m_ControlHistos[kCtrlTagPhi] = ibooker.book1D("TagMuonPhi", "TagMuonPhi; #phi", 66, -3.3, 3.3);
0398   m_ControlHistos[kCtrlTagEta] = ibooker.book1D("TagMuonEta", "TagMuonEta; #eta", 50, -2.5, 2.5);
0399 
0400   m_ControlHistos[kCtrlProbePt] = ibooker.book1D("ProbeMuonPt", "ProbeMuonPt; p_{T}", 50, 0., 100.);
0401   m_ControlHistos[kCtrlProbePhi] = ibooker.book1D("ProbeMuonPhi", "ProbeMuonPhi; #phi", 66, -3.3, 3.3);
0402   m_ControlHistos[kCtrlProbeEta] = ibooker.book1D("ProbeMuonEta", "ProbeMuonEta; #eta", 50, -2.5, 2.5);
0403 
0404   m_ControlHistos[kCtrlTagProbeDr] =
0405       ibooker.book1D("TagMuonProbeMuonDeltaR", "TagMuonProbeMuonDeltaR; #DeltaR", 50, 0., 5.0);
0406   m_ControlHistos[kCtrlTagHltDr] = ibooker.book1D("TagMuonHltDeltaR", "TagMuonHltDeltaR;#DeltaR", 55, 0., 0.11);
0407 }
0408 
0409 //_____________________________________________________________________
0410 void L1TMuonDQMOffline::bookEfficiencyHistos(DQMStore::IBooker& ibooker) {
0411   ibooker.setCurrentFolder(m_HistFolder + "/numerators_and_denominators");
0412 
0413   for (const auto var : m_effTypes) {
0414     auto histBins = getHistBinsEff(var);
0415     // histograms for eta variable get a special treatment
0416     if (var == kEffEta) {
0417       for (const auto& cut : m_cuts) {
0418         const auto gmtPtCut = cut.first;
0419         const auto qualLevel = cut.second;
0420         std::string name = "effDen_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut);
0421         m_EfficiencyDenEtaHistos[gmtPtCut] =
0422             ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0423         name = "effNum_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" + m_qualStrings[qualLevel];
0424         m_histoKeyEffNumEtaType histoKeyEffNumEta = {gmtPtCut, qualLevel};
0425         m_EfficiencyNumEtaHistos[histoKeyEffNumEta] =
0426             ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0427       }
0428     } else {
0429       for (const auto etaReg : m_etaRegions) {
0430         // denominator histograms for pt variable get a special treatment
0431         if (var == kEffPt) {
0432           std::string name = "effDen_" + m_effStrings[var] + "_" + m_etaStrings[etaReg];
0433           m_EfficiencyDenPtHistos[etaReg] =
0434               ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0435         } else {
0436           for (const auto& cut : m_cuts) {
0437             const int gmtPtCut = cut.first;
0438             std::string name =
0439                 "effDen_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" + m_etaStrings[etaReg];
0440             m_histoKeyEffDenVarType histoKeyEffDenVar = {var, gmtPtCut, etaReg};
0441             m_EfficiencyDenVarHistos[histoKeyEffDenVar] =
0442                 ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0443           }
0444         }
0445         for (const auto& cut : m_cuts) {
0446           const auto gmtPtCut = cut.first;
0447           const auto qualLevel = cut.second;
0448           std::string name = "effNum_" + m_effStrings[var] + "_" + std::to_string(gmtPtCut) + "_" +
0449                              m_etaStrings[etaReg] + "_" + m_qualStrings[qualLevel];
0450           m_histoKeyEffNumVarType histoKeyEffNum = {var, gmtPtCut, etaReg, qualLevel};
0451           m_EfficiencyNumVarHistos[histoKeyEffNum] =
0452               ibooker.book1D(name, name + ";" + m_effLabelStrings[var], histBins.size() - 1, &histBins[0]);
0453         }
0454       }
0455     }
0456   }
0457 }
0458 
0459 void L1TMuonDQMOffline::bookResolutionHistos(DQMStore::IBooker& ibooker) {
0460   if (m_verbose)
0461     edm::LogInfo("L1TMuonDQMOffline") << "Booking Resolution Plot Histos" << endl;
0462   ibooker.setCurrentFolder(m_HistFolder + "/resolution");
0463 
0464   for (const auto var : m_resTypes) {
0465     auto nbins = std::get<0>(getHistBinsRes(var));
0466     auto xmin = std::get<1>(getHistBinsRes(var));
0467     auto xmax = std::get<2>(getHistBinsRes(var));
0468     for (const auto etaReg : m_etaRegions) {
0469       for (const auto qualLevel : m_qualLevelsRes) {
0470         m_histoKeyResType histoKeyRes = {var, etaReg, qualLevel};
0471         std::string name =
0472             "resolution_" + m_resStrings[var] + "_" + m_etaStrings[etaReg] + "_" + m_qualStrings[qualLevel];
0473         m_ResolutionHistos[histoKeyRes] = ibooker.book1D(name, name + ";" + m_resLabelStrings[var], nbins, xmin, xmax);
0474       }
0475     }
0476   }
0477 }
0478 
0479 //_____________________________________________________________________
0480 const unsigned int L1TMuonDQMOffline::getNVertices(Handle<VertexCollection>& vertex) {
0481   unsigned int nVtx = 0;
0482 
0483   if (vertex.isValid()) {
0484     for (const auto& vertexIt : *vertex) {
0485       if (vertexIt.isValid() && !vertexIt.isFake()) {
0486         ++nVtx;
0487       }
0488     }
0489   }
0490   return nVtx;
0491 }
0492 
0493 //_____________________________________________________________________
0494 const reco::Vertex L1TMuonDQMOffline::getPrimaryVertex(Handle<VertexCollection>& vertex, Handle<BeamSpot>& beamSpot) {
0495   Vertex::Point posVtx;
0496   Vertex::Error errVtx;
0497 
0498   bool hasPrimaryVertex = false;
0499 
0500   if (vertex.isValid()) {
0501     vector<Vertex>::const_iterator vertexIt = vertex->begin();
0502     vector<Vertex>::const_iterator vertexEnd = vertex->end();
0503 
0504     for (; vertexIt != vertexEnd; ++vertexIt) {
0505       if (vertexIt->isValid() && !vertexIt->isFake()) {
0506         posVtx = vertexIt->position();
0507         errVtx = vertexIt->error();
0508         hasPrimaryVertex = true;
0509         break;
0510       }
0511     }
0512   }
0513 
0514   if (!hasPrimaryVertex) {
0515     posVtx = beamSpot->position();
0516     errVtx(0, 0) = beamSpot->BeamWidthX();
0517     errVtx(1, 1) = beamSpot->BeamWidthY();
0518     errVtx(2, 2) = beamSpot->sigmaZ();
0519   }
0520   const Vertex primaryVertex(posVtx, errVtx);
0521   return primaryVertex;
0522 }
0523 
0524 //_____________________________________________________________________
0525 void L1TMuonDQMOffline::getTightMuons(edm::Handle<reco::MuonCollection>& muons, const Vertex& vertex) {
0526   if (m_verbose)
0527     edm::LogInfo("L1TMuonDQMOffline") << "Getting tight muons" << endl;
0528   m_TightMuons.clear();
0529   MuonCollection::const_iterator muonIt = muons->begin();
0530   MuonCollection::const_iterator muonEnd = muons->end();
0531 
0532   for (; muonIt != muonEnd; ++muonIt) {
0533     if (muon::isTightMuon((*muonIt), vertex)) {
0534       m_TightMuons.push_back(&(*muonIt));
0535     }
0536   }
0537   m_ControlHistos[kCtrlNTightVsAll]->Fill(muons->size(), m_TightMuons.size());
0538 }
0539 
0540 //_____________________________________________________________________
0541 void L1TMuonDQMOffline::getProbeMuons(Handle<edm::TriggerResults>& trigResults,
0542                                       edm::Handle<trigger::TriggerEvent>& trigEvent) {
0543   if (m_verbose)
0544     edm::LogInfo("L1TMuonDQMOffline") << "getting probe muons" << endl;
0545   m_ProbeMuons.clear();
0546   std::vector<const reco::Muon*> tagMuonsInHist;
0547 
0548   tagMuonsInHist.clear();
0549 
0550   vector<const reco::Muon*>::const_iterator probeCandIt = m_TightMuons.begin();
0551   vector<const reco::Muon*>::const_iterator tightMuonsEnd = m_TightMuons.end();
0552 
0553   for (; probeCandIt != tightMuonsEnd; ++probeCandIt) {
0554     bool isProbe = false;
0555     vector<const reco::Muon*>::const_iterator tagCandIt = m_TightMuons.begin();
0556     float deltar = 0.;
0557 
0558     for (; tagCandIt != tightMuonsEnd; ++tagCandIt) {
0559       bool tagMuonAlreadyInHist = false;
0560       bool tagHasTrig = false;
0561       float eta = (*tagCandIt)->eta();
0562       float phi = (*tagCandIt)->phi();
0563       float pt = (*tagCandIt)->pt();
0564       float dEta = eta - (*probeCandIt)->eta();
0565       float dPhi = phi - (*probeCandIt)->phi();
0566       deltar = sqrt(dEta * dEta + dPhi * dPhi);
0567 
0568       if ((*tagCandIt) == (*probeCandIt) || deltar < m_minTagProbeDR)
0569         continue;  // CB has a little bias for closed-by muons
0570       auto matchHltDeltaR = matchHlt(trigEvent, (*tagCandIt));
0571       tagHasTrig = (matchHltDeltaR < m_maxHltMuonDR) && (pt > m_TagPtCut);
0572       isProbe |= tagHasTrig;
0573       if (tagHasTrig) {
0574         if (std::distance(m_TightMuons.begin(), m_TightMuons.end()) > 2) {
0575           for (vector<const reco::Muon*>::const_iterator tagMuonsInHistIt = tagMuonsInHist.begin();
0576                tagMuonsInHistIt != tagMuonsInHist.end();
0577                ++tagMuonsInHistIt) {
0578             if ((*tagCandIt) == (*tagMuonsInHistIt)) {
0579               tagMuonAlreadyInHist = true;
0580               break;
0581             }
0582           }
0583           if (tagMuonAlreadyInHist == false)
0584             tagMuonsInHist.push_back((*tagCandIt));
0585         }
0586         if (tagMuonAlreadyInHist == false) {
0587           m_ControlHistos[kCtrlTagEta]->Fill(eta);
0588           m_ControlHistos[kCtrlTagPhi]->Fill(phi);
0589           m_ControlHistos[kCtrlTagPt]->Fill(pt);
0590           m_ControlHistos[kCtrlTagProbeDr]->Fill(deltar);
0591           m_ControlHistos[kCtrlTagHltDr]->Fill(matchHltDeltaR);
0592         }
0593       }
0594     }
0595     if (isProbe)
0596       m_ProbeMuons.push_back((*probeCandIt));
0597   }
0598   m_ControlHistos[kCtrlNProbesVsTight]->Fill(m_TightMuons.size(), m_ProbeMuons.size());
0599 }
0600 
0601 //_____________________________________________________________________
0602 void L1TMuonDQMOffline::getMuonGmtPairs(edm::Handle<l1t::MuonBxCollection>& gmtCands,
0603                                         PropagateToMuon const& propagator) {
0604   m_MuonGmtPairs.clear();
0605   if (m_verbose)
0606     edm::LogInfo("L1TMuonDQMOffline") << "Getting muon GMT pairs" << endl;
0607 
0608   vector<const reco::Muon*>::const_iterator probeMuIt = m_ProbeMuons.begin();
0609   vector<const reco::Muon*>::const_iterator probeMuEnd = m_ProbeMuons.end();
0610 
0611   l1t::MuonBxCollection::const_iterator gmtIt;
0612   l1t::MuonBxCollection::const_iterator gmtEnd = gmtCands->end(0);
0613 
0614   for (; probeMuIt != probeMuEnd; ++probeMuIt) {
0615     MuonGmtPair pairBestCand((*probeMuIt), nullptr, propagator, m_useAtVtxCoord);
0616 
0617     // Fill the control histograms with the probe muon kinematic variables used
0618     m_ControlHistos[kCtrlProbeEta]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffEta));
0619     m_ControlHistos[kCtrlProbePhi]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffPhi));
0620     m_ControlHistos[kCtrlProbePt]->Fill(pairBestCand.getVar(L1TMuonDQMOffline::kEffPt));
0621 
0622     gmtIt = gmtCands->begin(0);  // use only on L1T muons from BX 0
0623 
0624     for (; gmtIt != gmtEnd; ++gmtIt) {
0625       MuonGmtPair pairTmpCand((*probeMuIt), &(*gmtIt), propagator, m_useAtVtxCoord);
0626 
0627       if ((pairTmpCand.dR() < m_maxGmtMuonDR) && (pairTmpCand.dR() < pairBestCand.dR())) {
0628         pairBestCand = pairTmpCand;
0629       }
0630     }
0631     m_MuonGmtPairs.push_back(pairBestCand);
0632     m_ControlHistos[kCtrlMuonGmtDeltaR]->Fill(pairBestCand.dR());
0633   }
0634 }
0635 
0636 //_____________________________________________________________________
0637 double L1TMuonDQMOffline::matchHlt(edm::Handle<TriggerEvent>& triggerEvent, const reco::Muon* mu) {
0638   double matchDeltaR = 9999;
0639 
0640   TriggerObjectCollection trigObjs = triggerEvent->getObjects();
0641 
0642   vector<int>::const_iterator trigIndexIt = m_trigIndices.begin();
0643   vector<int>::const_iterator trigIndexEnd = m_trigIndices.end();
0644 
0645   for (; trigIndexIt != trigIndexEnd; ++trigIndexIt) {
0646     const vector<string> moduleLabels(m_hltConfig.moduleLabels(*trigIndexIt));
0647     // For some reason various modules now come *after* "hltBoolEnd"
0648     // Really just want module one index before "hltBoolEnd" - AWB 2022.09.28
0649     // const unsigned moduleIndex = m_hltConfig.size((*trigIndexIt)) - 2;
0650     unsigned int moduleIndex = 999999;
0651     for (int ii = 0; ii < int(moduleLabels.size()); ii++) {
0652       if (moduleLabels[ii] == "hltBoolEnd") {
0653         moduleIndex = ii - 1;
0654         break;
0655       }
0656     }
0657     if (moduleIndex == 999999) {
0658       edm::LogError("L1TMuonDQMOffline") << "Found no module label in trigger " << (*trigIndexIt) << endl;
0659       continue;
0660     }
0661 
0662     const unsigned hltFilterIndex = triggerEvent->filterIndex(InputTag(moduleLabels[moduleIndex], "", m_trigProcess));
0663 
0664     if (hltFilterIndex < triggerEvent->sizeFilters()) {
0665       const Keys triggerKeys(triggerEvent->filterKeys(hltFilterIndex));
0666       const Vids triggerVids(triggerEvent->filterIds(hltFilterIndex));
0667       const unsigned nTriggers = triggerVids.size();
0668 
0669       for (size_t iTrig = 0; iTrig < nTriggers; ++iTrig) {
0670         const TriggerObject trigObject = trigObjs[triggerKeys[iTrig]];
0671         if (m_verbose)
0672           edm::LogInfo("L1TMuonDQMOffline") << "Found trigObject with pt = " << trigObject.pt()
0673                                             << ", eta = " << trigObject.eta() << ", phi = " << trigObject.phi() << endl;
0674         if (m_verbose)
0675           edm::LogInfo("L1TMuonDQMOffline") << "Compare to muon with pt = " << (*mu).pt() << ", eta = " << (*mu).eta()
0676                                             << ", phi = " << (*mu).phi() << endl;
0677         double dRtmp = deltaR((*mu), trigObject);
0678         if (dRtmp < matchDeltaR)
0679           matchDeltaR = dRtmp;
0680       }
0681     }
0682   }
0683   return matchDeltaR;
0684 }
0685 
0686 std::vector<float> L1TMuonDQMOffline::getHistBinsEff(EffType eff) {
0687   if (eff == kEffPt) {
0688     std::vector<float> effVsPtBins(m_effVsPtBins.begin(), m_effVsPtBins.end());
0689     return effVsPtBins;
0690   }
0691   if (eff == kEffPhi) {
0692     std::vector<float> effVsPhiBins(m_effVsPhiBins.begin(), m_effVsPhiBins.end());
0693     return effVsPhiBins;
0694   }
0695   if (eff == kEffEta) {
0696     std::vector<float> effVsEtaBins(m_effVsEtaBins.begin(), m_effVsEtaBins.end());
0697     return effVsEtaBins;
0698   }
0699   if (eff == kEffVtx) {
0700     std::vector<float> effVsVtxBins(m_effVsVtxBins.begin(), m_effVsVtxBins.end());
0701     return effVsVtxBins;
0702   }
0703   return {0., 1.};
0704 }
0705 
0706 std::tuple<int, double, double> L1TMuonDQMOffline::getHistBinsRes(ResType res) {
0707   if (res == kResPt)
0708     return {50, -2., 2.};
0709   if (res == kRes1OverPt)
0710     return {50, -2., 2.};
0711   if (res == kResQOverPt)
0712     return {50, -2., 2.};
0713   if (res == kResPhi)
0714     return {96, -0.2, 0.2};
0715   if (res == kResEta)
0716     return {100, -0.1, 0.1};
0717   if (res == kResCh)
0718     return {5, -2, 3};
0719   return {1, 0, 1};
0720 }
0721 
0722 //define this as a plug-in
0723 DEFINE_FWK_MODULE(L1TMuonDQMOffline);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!