Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:53:31 

0001 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "DataFormats/Math/interface/deltaPhi.h"
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 #include "DataFormats/MuonReco/interface/Muon.h"
0006 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0007 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
0008 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0009 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
0010 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeedCollection.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/path_configuration.h"
0014 #include "HLTriggerOffline/Muon/interface/HLTMuonPlotter.h"
0015 
0016 #include <algorithm>
0017 
0018 #include "TPRegexp.h"
0019 #include "TObjArray.h"
0020 #include "TObjString.h"
0021 
0022 namespace {
0023   const unsigned int kNull = (unsigned int)-1;
0024 }
0025 
0026 HLTMuonPlotter::HLTMuonPlotter(const edm::ParameterSet &pset,
0027                                const std::string &hltPath,
0028                                const std::vector<std::string> &moduleLabels,
0029                                const std::vector<std::string> &stepLabels,
0030                                const edm::EDGetTokenT<trigger::TriggerEventWithRefs> &triggerEventWithRefsToken,
0031                                const edm::EDGetTokenT<reco::GenParticleCollection> &genParticlesToken,
0032                                const edm::EDGetTokenT<reco::MuonCollection> &recoMuonsToken,
0033                                const L1MuonMatcherAlgoForDQM &l1Matcher)
0034     : hltPath_(hltPath),
0035       hltProcessName_(pset.getParameter<std::string>("hltProcessName")),
0036       moduleLabels_(moduleLabels),
0037       stepLabels_(stepLabels),
0038       triggerEventWithRefsToken_(triggerEventWithRefsToken),
0039       genParticleToken_(genParticlesToken),
0040       recMuonToken_(recoMuonsToken),
0041       genMuonSelector_(pset.getParameter<std::string>("genMuonCut")),
0042       recMuonSelector_(pset.getParameter<std::string>("recMuonCut")),
0043       cutsDr_(pset.getParameter<std::vector<double>>("cutsDr")),
0044       parametersEta_(pset.getParameter<std::vector<double>>("parametersEta")),
0045       parametersPhi_(pset.getParameter<std::vector<double>>("parametersPhi")),
0046       parametersTurnOn_(pset.getParameter<std::vector<double>>("parametersTurnOn")),
0047       l1Matcher_(l1Matcher),
0048       isInvalid_(false) {
0049   if (moduleLabels_.empty()) {
0050     edm::LogError("HLTMuonPlotter") << "Invalid inputs: 'moduleLabels_' is empty."
0051                                     << "\nMonitorElements for HLT path '" << hltPath_ << "' will not be produced.";
0052     isInvalid_ = true;
0053   } else if (stepLabels_.size() != moduleLabels_.size() + 1) {
0054     edm::LogError err("HLTMuonPlotter");
0055     err << "Invalid inputs: 'stepLabels_.size()' must equal 'moduleLabels_.size() + 1'.";
0056     err << "\nMonitorElements for HLT path '" << hltPath_ << "' will not be produced.";
0057     err << "\n  stepLabels_ = (";
0058     for (auto const &foo : stepLabels_)
0059       err << " " << foo;
0060     err << " )";
0061     err << "\n  moduleLabels_ = (";
0062     for (auto const &foo : moduleLabels_)
0063       err << " " << foo;
0064     err << " )";
0065     isInvalid_ = true;
0066   }
0067 }
0068 
0069 void HLTMuonPlotter::beginRun(DQMStore::IBooker &iBooker, const edm::Run &iRun, const edm::EventSetup &iSetup) {
0070   if (isInvalid_)
0071     return;
0072 
0073   l1Matcher_.init(iSetup);
0074 
0075   cutMaxEta_ = 2.4;
0076   if (hltPath_.find("eta2p1") != std::string::npos)
0077     cutMaxEta_ = 2.1;
0078 
0079   // Choose a pT cut for gen/rec muons based on the pT cut in the hltPath_
0080   unsigned int threshold = 0;
0081   TPRegexp ptRegexp("Mu([0-9]+)");
0082   TObjArray *regexArray = ptRegexp.MatchS(hltPath_);
0083   if (regexArray->GetEntriesFast() == 2) {
0084     threshold = atoi(((TObjString *)regexArray->At(1))->GetString());
0085   }
0086   delete regexArray;
0087   // We select a whole number min pT cut slightly above the hltPath_'s final
0088   // pt threshold, then subtract a bit to let through particle gun muons with
0089   // exact integer pT:
0090   cutMinPt_ = ceil(threshold * 1.1) - 0.01;
0091   if (cutMinPt_ < 0.)
0092     cutMinPt_ = 0.;
0093 
0094   std::string baseDir = "HLT/Muon/Distributions/";
0095   iBooker.setCurrentFolder(baseDir + hltPath_);
0096 
0097   std::vector<std::string> sources(2);
0098   sources[0] = "gen";
0099   sources[1] = "rec";
0100 
0101   elements_["CutMinPt"] = iBooker.bookFloat("CutMinPt");
0102   elements_["CutMaxEta"] = iBooker.bookFloat("CutMaxEta");
0103   elements_["CutMinPt"]->Fill(cutMinPt_);
0104   elements_["CutMaxEta"]->Fill(cutMaxEta_);
0105 
0106   for (size_t i = 0; i < sources.size(); i++) {
0107     std::string source = sources[i];
0108     for (size_t j = 0; j < stepLabels_.size(); j++) {
0109       bookHist(iBooker, hltPath_, stepLabels_[j], source, "Eta");
0110       bookHist(iBooker, hltPath_, stepLabels_[j], source, "Phi");
0111       bookHist(iBooker, hltPath_, stepLabels_[j], source, "MaxPt1");
0112       bookHist(iBooker, hltPath_, stepLabels_[j], source, "MaxPt2");
0113     }
0114   }
0115 }
0116 
0117 void HLTMuonPlotter::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0118   if (isInvalid_)
0119     return;
0120 
0121   auto const triggerEventWithRefsHandle = iEvent.getHandle(triggerEventWithRefsToken_);
0122   if (triggerEventWithRefsHandle.failedToGet()) {
0123     edm::LogError("HLTMuonPlotter") << "No trigger summary found";
0124     return;
0125   }
0126 
0127   auto const recoMuonsHandle = iEvent.getHandle(recMuonToken_);
0128   auto const genParticlesHandle = iEvent.getHandle(genParticleToken_);
0129 
0130   LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter] --------------";
0131   LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter] Event: " << iEvent.id();
0132 
0133   const int nFilters = moduleLabels_.size();
0134   const int nSteps = stepLabels_.size();
0135   const bool hasStepL1 = (stepLabels_.size() > 1 and stepLabels_[1] == "L1");
0136   const int nStepsHlt = hasStepL1 ? nSteps - 2 : nSteps - 1;
0137   const bool isDoubleMuonPath = (hltPath_.find("Double") != std::string::npos);
0138   const int nObjectsToPassPath = (isDoubleMuonPath) ? 2 : 1;
0139 
0140   LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter] HLTPath=" << hltPath_ << " nFilters=" << nFilters
0141                              << " nSteps=" << nSteps << " hasStepL1=" << hasStepL1 << " nStepsHlt=" << nStepsHlt
0142                              << " isDoubleMuonPath=" << isDoubleMuonPath
0143                              << " nObjectsToPassPath=" << nObjectsToPassPath;
0144 
0145   if (nFilters + 1 == nSteps) {
0146     LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]   stepLabels | moduleLabels";
0147     for (int istep = 0; istep < nSteps; ++istep) {
0148       if (istep == 0)
0149         LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]   " << stepLabels_[istep] << " | [N/A]";
0150       else
0151         LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]   " << stepLabels_[istep] << " | " << moduleLabels_[istep - 1];
0152     }
0153   }
0154 
0155   std::vector<std::string> sources;
0156   if (genParticlesHandle.isValid())
0157     sources.push_back("gen");
0158   if (recoMuonsHandle.isValid())
0159     sources.push_back("rec");
0160 
0161   for (size_t sourceNo = 0; sourceNo < sources.size(); sourceNo++) {
0162     std::string const &source = sources[sourceNo];
0163 
0164     LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter] source=" << source;
0165 
0166     // Make each good gen/rec muon into the base cand for a MatchStruct
0167     std::vector<MatchStruct> matches;
0168 
0169     if (source == "gen") {
0170       if (genParticlesHandle.isValid()) {
0171         matches.reserve(matches.size() + genParticlesHandle->size());
0172         for (auto const &genp : *genParticlesHandle)
0173           if (genMuonSelector_(genp))
0174             matches.emplace_back(MatchStruct(&genp));
0175       }
0176     } else if (source == "rec") {
0177       if (recoMuonsHandle.isValid()) {
0178         matches.reserve(matches.size() + recoMuonsHandle->size());
0179         for (auto const &recomu : *recoMuonsHandle)
0180           if (recMuonSelector_(recomu))
0181             matches.emplace_back(MatchStruct(&recomu));
0182       }
0183     }
0184 
0185     // Sort the MatchStructs by pT for later filling of turn-on curve
0186     std::sort(matches.begin(), matches.end(), matchesByDescendingPt());
0187 
0188     l1t::MuonVectorRef candsL1;
0189     std::vector<std::vector<reco::RecoChargedCandidateRef>> refsHlt(nStepsHlt);
0190     std::vector<std::vector<const reco::RecoChargedCandidate *>> candsHlt(nStepsHlt);
0191 
0192     for (int idx = 0; idx < nFilters; ++idx) {
0193       auto const moduleLabelStripped =
0194           edm::path_configuration::removeSchedulingTokensFromModuleLabel(moduleLabels_[idx]);
0195       auto const iTag = edm::InputTag(moduleLabelStripped, "", hltProcessName_);
0196       auto const iFilter = triggerEventWithRefsHandle->filterIndex(iTag);
0197       auto const iFilterValid = (iFilter < triggerEventWithRefsHandle->size());
0198 
0199       LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]   InputTag[" << idx << "]: " << moduleLabels_[idx]
0200                                  << "::" << hltProcessName_ << " (filterIndex = " << iFilter
0201                                  << ", valid = " << iFilterValid << ")";
0202 
0203       if (iFilterValid) {
0204         if (idx == 0 and hasStepL1)
0205           triggerEventWithRefsHandle->getObjects(iFilter, trigger::TriggerL1Mu, candsL1);
0206         else {
0207           auto const hltStep = hasStepL1 ? idx - 1 : idx;
0208           triggerEventWithRefsHandle->getObjects(iFilter, trigger::TriggerMuon, refsHlt[hltStep]);
0209         }
0210       } else
0211         LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]    No collection with " << iTag;
0212     }
0213 
0214     for (int i = 0; i < nStepsHlt; i++) {
0215       for (size_t j = 0; j < refsHlt[i].size(); j++) {
0216         if (refsHlt[i][j].isAvailable())
0217           candsHlt[i].push_back(&*refsHlt[i][j]);
0218         else
0219           edm::LogWarning("HLTMuonPlotter") << "Ref refsHlt[i][j]: product not available " << i << " " << j;
0220       }
0221     }
0222 
0223     // Add trigger objects to the MatchStructs
0224     findMatches(matches, candsL1, candsHlt);
0225 
0226     LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]   Number of Candidates = " << matches.size();
0227 
0228     for (auto const &match_i : matches) {
0229       if (!match_i.candBase)
0230         continue;
0231       LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]    CandBase: pt=" << match_i.candBase->pt()
0232                                  << " eta=" << match_i.candBase->eta() << " phi=" << match_i.candBase->phi();
0233       if (match_i.candL1)
0234         LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]      CandL1: pt=" << match_i.candL1->pt()
0235                                    << " eta=" << match_i.candL1->eta() << " phi=" << match_i.candL1->phi();
0236       else
0237         LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]      CandL1: NULL";
0238 
0239       int ihlt = -1;
0240       for (auto const *chlt : match_i.candHlt) {
0241         ++ihlt;
0242         if (chlt)
0243           LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]      CandHLT[" << ihlt << "]: pt=" << chlt->pt()
0244                                      << " eta=" << chlt->eta() << " phi=" << chlt->phi();
0245         else
0246           LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]      CandHLT[" << ihlt << "]: NULL";
0247       }
0248     }
0249 
0250     std::vector<size_t> matchesInEtaRange;
0251     std::vector<bool> hasMatch(matches.size(), true);
0252 
0253     for (int step = 0; step < nSteps; step++) {
0254       int const hltStep = hasStepL1 ? step - 2 : step - 1;
0255       size_t level = 0;
0256       if ((stepLabels_[step].find("L3TkIso") != std::string::npos) ||
0257           (stepLabels_[step].find("TkTkIso") != std::string::npos))
0258         level = 6;
0259       else if ((stepLabels_[step].find("L3HcalIso") != std::string::npos) ||
0260                (stepLabels_[step].find("TkEcalIso") != std::string::npos))
0261         level = 5;
0262       else if ((stepLabels_[step].find("L3EcalIso") != std::string::npos) ||
0263                (stepLabels_[step].find("TkEcalIso") != std::string::npos))
0264         level = 4;
0265       else if ((stepLabels_[step].find("L3") != std::string::npos) ||
0266                (stepLabels_[step].find("Tk") != std::string::npos))
0267         level = 3;
0268       else if (stepLabels_[step].find("L2") != std::string::npos)
0269         level = 2;
0270       else if (stepLabels_[step].find("L1") != std::string::npos)
0271         level = 1;
0272 
0273       for (size_t j = 0; j < matches.size(); j++) {
0274         if (level == 0) {
0275           if (std::abs(matches[j].candBase->eta()) < cutMaxEta_)
0276             matchesInEtaRange.push_back(j);
0277         } else if (level == 1) {
0278           if (matches[j].candL1 == nullptr)
0279             hasMatch[j] = false;
0280         } else if (level >= 2) {
0281           if (matches[j].candHlt.at(hltStep) == nullptr)
0282             hasMatch[j] = false;
0283           else if (!hasMatch[j]) {
0284             LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]     match found for " << source << " candidate " << j
0285                                        << " in HLT step " << hltStep << " of " << nStepsHlt
0286                                        << " without previous match!";
0287             break;
0288           }
0289         }
0290       }
0291 
0292       LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]    (step=" << step << ", level=" << level
0293                                  << ", hltStep=" << hltStep << ") matchesInEtaRange: [ "
0294                                  << this->vector_to_string(matchesInEtaRange) << " ]";
0295 
0296       LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]    (step=" << step << ", level=" << level
0297                                  << ", hltStep=" << hltStep << ") hasMatch: [ " << this->vector_to_string(hasMatch)
0298                                  << " ]";
0299 
0300       if (std::count(hasMatch.begin(), hasMatch.end(), true) < nObjectsToPassPath)
0301         break;
0302 
0303       std::string const pre = source + "Pass";
0304       std::string const post = "_" + stepLabels_[step];
0305 
0306       for (size_t j = 0; j < matches.size(); j++) {
0307         float const pt = matches[j].candBase->pt();
0308         float const eta = matches[j].candBase->eta();
0309         float const phi = matches[j].candBase->phi();
0310         if (hasMatch[j]) {
0311           if (!matchesInEtaRange.empty() && j == matchesInEtaRange[0]) {
0312             elements_[pre + "MaxPt1" + post]->Fill(pt);
0313             LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]     FILL(" << pre + "MaxPt1" + post << ") value = " << pt;
0314           }
0315           if (matchesInEtaRange.size() >= 2 && j == matchesInEtaRange[1]) {
0316             elements_[pre + "MaxPt2" + post]->Fill(pt);
0317             LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]     FILL(" << pre + "MaxPt2" + post << ") value = " << pt;
0318           }
0319           if (pt > cutMinPt_) {
0320             elements_[pre + "Eta" + post]->Fill(eta);
0321             LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]     FILL(" << pre + "Eta" + post << ") value = " << eta;
0322             if (std::abs(eta) < cutMaxEta_) {
0323               elements_[pre + "Phi" + post]->Fill(phi);
0324               LogTrace("HLTMuonPlotter") << "[HLTMuonPlotter]     FILL(" << pre + "Phi" + post << ") value = " << phi;
0325             }
0326           }
0327         }
0328       }
0329     }
0330   }  // End loop over sources
0331 }
0332 
0333 void HLTMuonPlotter::findMatches(std::vector<MatchStruct> &matches,
0334                                  const l1t::MuonVectorRef &candsL1,
0335                                  const std::vector<std::vector<const reco::RecoChargedCandidate *>> &candsHlt) {
0336   std::set<size_t>::iterator it;
0337 
0338   std::set<size_t> indicesL1;
0339   for (size_t i = 0; i < candsL1.size(); i++)
0340     indicesL1.insert(i);
0341 
0342   std::vector<set<size_t>> indicesHlt(candsHlt.size());
0343   for (size_t i = 0; i < candsHlt.size(); i++)
0344     for (size_t j = 0; j < candsHlt[i].size(); j++)
0345       indicesHlt[i].insert(j);
0346 
0347   for (size_t i = 0; i < matches.size(); i++) {
0348     const reco::Candidate *cand = matches[i].candBase;
0349 
0350     double bestDeltaR = cutsDr_[0];
0351     size_t bestMatch = kNull;
0352     for (it = indicesL1.begin(); it != indicesL1.end(); it++) {
0353       if (candsL1[*it].isAvailable()) {
0354         double dR = deltaR(cand->eta(), cand->phi(), candsL1[*it]->eta(), candsL1[*it]->phi());
0355         if (dR < bestDeltaR) {
0356           bestMatch = *it;
0357           bestDeltaR = dR;
0358         }
0359         // TrajectoryStateOnSurface propagated;
0360         // float dR = 999., dPhi = 999.;
0361         // bool isValid = l1Matcher_.match(* cand, * candsL1[*it],
0362         //                                 dR, dPhi, propagated);
0363         // if (isValid && dR < bestDeltaR) {
0364         //   bestMatch = *it;
0365         //   bestDeltaR = dR;
0366         // }
0367       } else {
0368         edm::LogWarning("HLTMuonPlotter") << "Ref candsL1[*it]: product not available " << *it;
0369       }
0370     }
0371 
0372     if (bestMatch != kNull)
0373       matches[i].candL1 = &*candsL1[bestMatch];
0374     indicesL1.erase(bestMatch);
0375 
0376     matches[i].candHlt.assign(candsHlt.size(), nullptr);
0377     for (size_t j = 0; j < candsHlt.size(); j++) {
0378       size_t level = (candsHlt.size() == 4) ? (j < 2) ? 2 : 3 : (candsHlt.size() == 2) ? (j < 1) ? 2 : 3 : 2;
0379       bestDeltaR = cutsDr_[level - 2];
0380       bestMatch = kNull;
0381       for (it = indicesHlt[j].begin(); it != indicesHlt[j].end(); it++) {
0382         double dR = deltaR(cand->eta(), cand->phi(), candsHlt[j][*it]->eta(), candsHlt[j][*it]->phi());
0383         if (dR < bestDeltaR) {
0384           bestMatch = *it;
0385           bestDeltaR = dR;
0386         }
0387       }
0388       if (bestMatch != kNull)
0389         matches[i].candHlt[j] = candsHlt[j][bestMatch];
0390       indicesHlt[j].erase(bestMatch);
0391     }
0392   }
0393 }
0394 
0395 void HLTMuonPlotter::bookHist(DQMStore::IBooker &iBooker,
0396                               std::string const &path,
0397                               std::string const &label,
0398                               std::string const &source,
0399                               std::string const &type) {
0400   std::string sourceUpper = source;
0401   sourceUpper[0] = toupper(sourceUpper[0]);
0402   std::string name = source + "Pass" + type + "_" + label;
0403   TH1F *h;
0404 
0405   if (type.find("MaxPt") != std::string::npos) {
0406     std::string desc = (type == "MaxPt1") ? "Leading" : "Next-to-Leading";
0407     std::string title = "pT of " + desc + " " + sourceUpper + " Muon " + "matched to " + label;
0408     const size_t nBins = parametersTurnOn_.size() - 1;
0409     float *edges = new float[nBins + 1];
0410     for (size_t i = 0; i < nBins + 1; i++)
0411       edges[i] = parametersTurnOn_[i];
0412     h = new TH1F(name.c_str(), title.c_str(), nBins, edges);
0413   } else {
0414     std::string symbol = (type == "Eta") ? "#eta" : "#phi";
0415     std::string title = symbol + " of " + sourceUpper + " Muons " + "matched to " + label;
0416     std::vector<double> params = (type == "Eta") ? parametersEta_ : parametersPhi_;
0417     int nBins = (int)params[0];
0418     double min = params[1];
0419     double max = params[2];
0420     h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
0421   }
0422 
0423   h->Sumw2();
0424   elements_[name] = iBooker.book1D(name, h);
0425   delete h;
0426 }
[ 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!