File indexing completed on 2022-07-12 22:34:35
0001
0002
0003
0004
0005 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/Math/interface/deltaPhi.h"
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "DataFormats/MuonReco/interface/Muon.h"
0010 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "HLTriggerOffline/Muon/interface/HLTMuonPlotter.h"
0014
0015 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
0016 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0017 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
0018 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeedCollection.h"
0019
0020 using namespace std;
0021 using namespace edm;
0022 using namespace reco;
0023 using namespace trigger;
0024
0025 namespace {
0026 const unsigned int kNull = (unsigned int)-1;
0027 }
0028
0029 typedef vector<ParameterSet> Parameters;
0030
0031 HLTMuonPlotter::HLTMuonPlotter(const ParameterSet &pset,
0032 string hltPath,
0033 const std::vector<string> &moduleLabels,
0034 const std::vector<string> &stepLabels,
0035 const edm::EDGetTokenT<trigger::TriggerEventWithRefs> &triggerEventToken,
0036 const edm::EDGetTokenT<reco::GenParticleCollection> &genParticlesToken,
0037 const edm::EDGetTokenT<reco::MuonCollection> &recoMuonsToken,
0038 const L1MuonMatcherAlgoForDQM &l1Matcher)
0039 : l1Matcher_(l1Matcher) {
0040 hltPath_ = hltPath;
0041 moduleLabels_ = moduleLabels;
0042 stepLabels_ = stepLabels;
0043 hltProcessName_ = pset.getParameter<string>("hltProcessName");
0044
0045 cutsDr_ = pset.getParameter<vector<double>>("cutsDr");
0046
0047 parametersEta_ = pset.getParameter<vector<double>>("parametersEta");
0048 parametersPhi_ = pset.getParameter<vector<double>>("parametersPhi");
0049 parametersTurnOn_ = pset.getParameter<vector<double>>("parametersTurnOn");
0050
0051 genMuonCut_ = pset.getParameter<string>("genMuonCut");
0052 recMuonCut_ = pset.getParameter<string>("recMuonCut");
0053
0054 genMuonSelector_ = nullptr;
0055 recMuonSelector_ = nullptr;
0056
0057 hltTriggerSummaryRAW_ = triggerEventToken;
0058 genParticleLabel_ = genParticlesToken;
0059 recMuonLabel_ = recoMuonsToken;
0060 }
0061
0062 void HLTMuonPlotter::beginJob() {}
0063
0064 void HLTMuonPlotter::beginRun(DQMStore::IBooker &iBooker, const Run &iRun, const EventSetup &iSetup) {
0065 l1Matcher_.init(iSetup);
0066
0067 cutMaxEta_ = 2.4;
0068 if (hltPath_.find("eta2p1") != string::npos)
0069 cutMaxEta_ = 2.1;
0070
0071
0072 unsigned int threshold = 0;
0073 TPRegexp ptRegexp("Mu([0-9]+)");
0074 TObjArray *regexArray = ptRegexp.MatchS(hltPath_);
0075 if (regexArray->GetEntriesFast() == 2) {
0076 threshold = atoi(((TObjString *)regexArray->At(1))->GetString());
0077 }
0078 delete regexArray;
0079
0080
0081
0082 cutMinPt_ = ceil(threshold * 1.1) - 0.01;
0083 if (cutMinPt_ < 0.)
0084 cutMinPt_ = 0.;
0085
0086 string baseDir = "HLT/Muon/Distributions/";
0087 iBooker.setCurrentFolder(baseDir + hltPath_);
0088
0089 vector<string> sources(2);
0090 sources[0] = "gen";
0091 sources[1] = "rec";
0092
0093 elements_["CutMinPt"] = iBooker.bookFloat("CutMinPt");
0094 elements_["CutMaxEta"] = iBooker.bookFloat("CutMaxEta");
0095 elements_["CutMinPt"]->Fill(cutMinPt_);
0096 elements_["CutMaxEta"]->Fill(cutMaxEta_);
0097
0098 for (size_t i = 0; i < sources.size(); i++) {
0099 string source = sources[i];
0100 for (size_t j = 0; j < stepLabels_.size(); j++) {
0101 bookHist(iBooker, hltPath_, stepLabels_[j], source, "Eta");
0102 bookHist(iBooker, hltPath_, stepLabels_[j], source, "Phi");
0103 bookHist(iBooker, hltPath_, stepLabels_[j], source, "MaxPt1");
0104 bookHist(iBooker, hltPath_, stepLabels_[j], source, "MaxPt2");
0105 }
0106 }
0107 }
0108
0109 void HLTMuonPlotter::analyze(const Event &iEvent, const EventSetup &iSetup) {
0110 LogTrace("HLTMuonVal") << "In HLTMuonPlotter::analyze, "
0111 << "Event: " << iEvent.id();
0112
0113
0114
0115
0116
0117 Handle<TriggerEventWithRefs> rawTriggerEvent;
0118 Handle<MuonCollection> recMuons;
0119 Handle<GenParticleCollection> genParticles;
0120
0121 iEvent.getByToken(hltTriggerSummaryRAW_, rawTriggerEvent);
0122 if (rawTriggerEvent.failedToGet()) {
0123 LogError("HLTMuonVal") << "No trigger summary found";
0124 return;
0125 }
0126 iEvent.getByToken(recMuonLabel_, recMuons);
0127 iEvent.getByToken(genParticleLabel_, genParticles);
0128
0129 vector<string> sources;
0130 if (genParticles.isValid())
0131 sources.push_back("gen");
0132 if (recMuons.isValid())
0133 sources.push_back("rec");
0134
0135 for (size_t sourceNo = 0; sourceNo < sources.size(); sourceNo++) {
0136 string source = sources[sourceNo];
0137
0138
0139 if (!genMuonSelector_)
0140 genMuonSelector_ = new StringCutObjectSelector<reco::GenParticle>(genMuonCut_);
0141 if (!recMuonSelector_)
0142 recMuonSelector_ = new StringCutObjectSelector<reco::Muon>(recMuonCut_);
0143
0144
0145 vector<MatchStruct> matches;
0146 if (source == "gen" && genParticles.isValid())
0147 for (size_t i = 0; i < genParticles->size(); i++)
0148 if ((*genMuonSelector_)(genParticles->at(i)))
0149 matches.push_back(MatchStruct(&genParticles->at(i)));
0150 if (source == "rec" && recMuons.isValid())
0151 for (size_t i = 0; i < recMuons->size(); i++)
0152 if ((*recMuonSelector_)(recMuons->at(i)))
0153 matches.push_back(MatchStruct(&recMuons->at(i)));
0154
0155
0156 sort(matches.begin(), matches.end(), matchesByDescendingPt());
0157
0158 const bool isDoubleMuonPath = (hltPath_.find("Double") != string::npos);
0159 const size_t nFilters = moduleLabels_.size();
0160 const size_t nSteps = stepLabels_.size();
0161 const size_t nStepsHlt = nSteps - 2;
0162 const int nObjectsToPassPath = (isDoubleMuonPath) ? 2 : 1;
0163 l1t::MuonVectorRef candsL1;
0164 vector<vector<RecoChargedCandidateRef>> refsHlt(nStepsHlt);
0165 vector<vector<const RecoChargedCandidate *>> candsHlt(nStepsHlt);
0166
0167 for (size_t i = 0; i < nFilters; i++) {
0168 const int hltStep = i - 1;
0169 InputTag tag = InputTag(moduleLabels_[i], "", hltProcessName_);
0170 size_t iFilter = rawTriggerEvent->filterIndex(tag);
0171 if (iFilter < rawTriggerEvent->size()) {
0172 if (i == 0)
0173 rawTriggerEvent->getObjects(iFilter, TriggerL1Mu, candsL1);
0174 else
0175 rawTriggerEvent->getObjects(iFilter, TriggerMuon, refsHlt[hltStep]);
0176 } else
0177 LogTrace("HLTMuonVal") << "No collection with label " << tag;
0178 }
0179 for (size_t i = 0; i < nStepsHlt; i++)
0180 for (size_t j = 0; j < refsHlt[i].size(); j++)
0181 if (refsHlt[i][j].isAvailable()) {
0182 candsHlt[i].push_back(&*refsHlt[i][j]);
0183 } else {
0184 LogWarning("HLTMuonPlotter") << "Ref refsHlt[i][j]: product not available " << i << " " << j;
0185 }
0186
0187
0188 findMatches(matches, candsL1, candsHlt);
0189
0190 vector<size_t> matchesInEtaRange;
0191 vector<bool> hasMatch(matches.size(), true);
0192
0193 for (size_t step = 0; step < nSteps; step++) {
0194 size_t hltStep = (step >= 2) ? step - 2 : 0;
0195 if (nSteps == 6)
0196 hltStep = hltStep - 1;
0197 size_t level = 0;
0198 if ((stepLabels_[step].find("L3TkIso") != string::npos) || (stepLabels_[step].find("TkTkIso") != string::npos))
0199 level = 6;
0200 else if ((stepLabels_[step].find("L3HcalIso") != string::npos) ||
0201 (stepLabels_[step].find("TkEcalIso") != string::npos))
0202 level = 5;
0203 else if ((stepLabels_[step].find("L3EcalIso") != string::npos) ||
0204 (stepLabels_[step].find("TkEcalIso") != string::npos))
0205 level = 4;
0206 else if ((stepLabels_[step].find("L3") != string::npos) || (stepLabels_[step].find("Tk") != string::npos))
0207 level = 3;
0208 else if (stepLabels_[step].find("L2") != string::npos)
0209 level = 2;
0210 else if (stepLabels_[step].find("L1") != string::npos)
0211 level = 1;
0212
0213 for (size_t j = 0; j < matches.size(); j++) {
0214 if (level == 0) {
0215 if (fabs(matches[j].candBase->eta()) < cutMaxEta_)
0216 matchesInEtaRange.push_back(j);
0217 } else if (level == 1) {
0218 if (matches[j].candL1 == nullptr)
0219 hasMatch[j] = false;
0220 } else if (level >= 2) {
0221 if (matches[j].candHlt[hltStep] == nullptr)
0222 hasMatch[j] = false;
0223 else if (!hasMatch[j]) {
0224 LogTrace("HLTMuonVal") << "Match found for HLT step " << hltStep << " of " << nStepsHlt
0225 << " without previous match!";
0226 break;
0227 }
0228 }
0229 }
0230
0231 if (std::count(hasMatch.begin(), hasMatch.end(), true) < nObjectsToPassPath)
0232 break;
0233
0234 string pre = source + "Pass";
0235 string post = "_" + stepLabels_[step];
0236
0237 for (size_t j = 0; j < matches.size(); j++) {
0238 float pt = matches[j].candBase->pt();
0239 float eta = matches[j].candBase->eta();
0240 float phi = matches[j].candBase->phi();
0241 if (hasMatch[j]) {
0242 if (!matchesInEtaRange.empty() && j == matchesInEtaRange[0])
0243 elements_[pre + "MaxPt1" + post]->Fill(pt);
0244 if (matchesInEtaRange.size() >= 2 && j == matchesInEtaRange[1])
0245 elements_[pre + "MaxPt2" + post]->Fill(pt);
0246 if (pt > cutMinPt_) {
0247 elements_[pre + "Eta" + post]->Fill(eta);
0248 if (fabs(eta) < cutMaxEta_)
0249 elements_[pre + "Phi" + post]->Fill(phi);
0250 }
0251 }
0252 }
0253 }
0254
0255 }
0256 }
0257
0258 void HLTMuonPlotter::findMatches(vector<MatchStruct> &matches,
0259 const l1t::MuonVectorRef &candsL1,
0260 const std::vector<vector<const RecoChargedCandidate *>> &candsHlt) {
0261 set<size_t>::iterator it;
0262
0263 set<size_t> indicesL1;
0264 for (size_t i = 0; i < candsL1.size(); i++)
0265 indicesL1.insert(i);
0266
0267 vector<set<size_t>> indicesHlt(candsHlt.size());
0268 for (size_t i = 0; i < candsHlt.size(); i++)
0269 for (size_t j = 0; j < candsHlt[i].size(); j++)
0270 indicesHlt[i].insert(j);
0271
0272 for (size_t i = 0; i < matches.size(); i++) {
0273 const Candidate *cand = matches[i].candBase;
0274
0275 double bestDeltaR = cutsDr_[0];
0276 size_t bestMatch = kNull;
0277 for (it = indicesL1.begin(); it != indicesL1.end(); it++) {
0278 if (candsL1[*it].isAvailable()) {
0279 double dR = deltaR(cand->eta(), cand->phi(), candsL1[*it]->eta(), candsL1[*it]->phi());
0280 if (dR < bestDeltaR) {
0281 bestMatch = *it;
0282 bestDeltaR = dR;
0283 }
0284
0285
0286
0287
0288
0289
0290
0291
0292 } else {
0293 LogWarning("HLTMuonPlotter") << "Ref candsL1[*it]: product not available " << *it;
0294 }
0295 }
0296
0297 if (bestMatch != kNull)
0298 matches[i].candL1 = &*candsL1[bestMatch];
0299 indicesL1.erase(bestMatch);
0300
0301 matches[i].candHlt.assign(candsHlt.size(), nullptr);
0302 for (size_t j = 0; j < candsHlt.size(); j++) {
0303 size_t level = (candsHlt.size() == 4) ? (j < 2) ? 2 : 3 : (candsHlt.size() == 2) ? (j < 1) ? 2 : 3 : 2;
0304 bestDeltaR = cutsDr_[level - 2];
0305 bestMatch = kNull;
0306 for (it = indicesHlt[j].begin(); it != indicesHlt[j].end(); it++) {
0307 double dR = deltaR(cand->eta(), cand->phi(), candsHlt[j][*it]->eta(), candsHlt[j][*it]->phi());
0308 if (dR < bestDeltaR) {
0309 bestMatch = *it;
0310 bestDeltaR = dR;
0311 }
0312 }
0313 if (bestMatch != kNull)
0314 matches[i].candHlt[j] = candsHlt[j][bestMatch];
0315 indicesHlt[j].erase(bestMatch);
0316 }
0317
0318
0319
0320
0321
0322
0323
0324
0325 }
0326 }
0327
0328 void HLTMuonPlotter::bookHist(DQMStore::IBooker &iBooker, string path, string label, string source, string type) {
0329 string sourceUpper = source;
0330 sourceUpper[0] = toupper(sourceUpper[0]);
0331 string name = source + "Pass" + type + "_" + label;
0332 TH1F *h;
0333
0334 if (type.find("MaxPt") != string::npos) {
0335 string desc = (type == "MaxPt1") ? "Leading" : "Next-to-Leading";
0336 string title = "pT of " + desc + " " + sourceUpper + " Muon " + "matched to " + label;
0337 const size_t nBins = parametersTurnOn_.size() - 1;
0338 float *edges = new float[nBins + 1];
0339 for (size_t i = 0; i < nBins + 1; i++)
0340 edges[i] = parametersTurnOn_[i];
0341 h = new TH1F(name.c_str(), title.c_str(), nBins, edges);
0342 }
0343
0344 else {
0345 string symbol = (type == "Eta") ? "#eta" : "#phi";
0346 string title = symbol + " of " + sourceUpper + " Muons " + "matched to " + label;
0347 vector<double> params = (type == "Eta") ? parametersEta_ : parametersPhi_;
0348 int nBins = (int)params[0];
0349 double min = params[1];
0350 double max = params[2];
0351 h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
0352 }
0353
0354 h->Sumw2();
0355 elements_[name] = iBooker.book1D(name, h);
0356 delete h;
0357 }