File indexing completed on 2022-03-09 00:51:02
0001
0002
0003
0004
0005 #include "DQMOffline/Trigger/interface/HLTMuonMatchAndPlot.h"
0006
0007 #include "FWCore/Framework/interface/Run.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/Common/interface/TriggerNames.h"
0014 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0015 #include <boost/algorithm/string/replace.hpp>
0016
0017 #include <iostream>
0018 #include <string>
0019 #include <utility>
0020 #include <utility>
0021
0022
0023
0024
0025 using namespace std;
0026 using namespace edm;
0027 using namespace reco;
0028 using namespace trigger;
0029 using namespace l1extra;
0030
0031 using vstring = std::vector<std::string>;
0032
0033
0034
0035
0036
0037 HLTMuonMatchAndPlot::HLTMuonMatchAndPlot(const ParameterSet& pset, string hltPath, string moduleLabel, bool islastfilter)
0038 : hltProcessName_(pset.getParameter<string>("hltProcessName")),
0039 destination_(pset.getUntrackedParameter<string>("destination")),
0040 requiredTriggers_(pset.getUntrackedParameter<vstring>("requiredTriggers")),
0041 targetParams_(pset.getParameterSet("targetParams")),
0042 probeParams_(pset.getParameterSet("probeParams")),
0043 hltPath_(std::move(hltPath)),
0044 moduleLabel_(std::move(moduleLabel)),
0045 isLastFilter_(islastfilter),
0046 hasTargetRecoCuts(targetParams_.exists("recoCuts")),
0047 hasProbeRecoCuts(probeParams_.exists("recoCuts")),
0048 targetMuonSelector_(targetParams_.getUntrackedParameter<string>("recoCuts", "")),
0049 targetZ0Cut_(targetParams_.getUntrackedParameter<double>("z0Cut", 0.)),
0050 targetD0Cut_(targetParams_.getUntrackedParameter<double>("d0Cut", 0.)),
0051 targetptCutZ_(targetParams_.getUntrackedParameter<double>("ptCut_Z", 20.)),
0052 targetptCutJpsi_(targetParams_.getUntrackedParameter<double>("ptCut_Jpsi", 20.)),
0053 probeMuonSelector_(probeParams_.getUntrackedParameter<string>("recoCuts", "")),
0054 probeZ0Cut_(probeParams_.getUntrackedParameter<double>("z0Cut", 0.)),
0055 probeD0Cut_(probeParams_.getUntrackedParameter<double>("d0Cut", 0.)),
0056 triggerSelector_(targetParams_.getUntrackedParameter<string>("hltCuts", "")),
0057 hasTriggerCuts_(targetParams_.exists("hltCuts")) {
0058
0059 fillMapFromPSet(binParams_, pset, "binParams");
0060 fillMapFromPSet(plotCuts_, pset, "plotCuts");
0061
0062
0063 triggerLevel_ = "L3";
0064 TPRegexp levelRegexp("L[1-3]");
0065
0066
0067 TObjArray* levelArray = levelRegexp.MatchS(moduleLabel_);
0068 if (levelArray->GetEntriesFast() > 0) {
0069 triggerLevel_ = static_cast<const char*>(((TObjString*)levelArray->At(0))->GetString());
0070 }
0071 delete levelArray;
0072
0073
0074 cutMinPt_ = 3;
0075 TPRegexp ptRegexp("Mu([0-9]*)");
0076 TObjArray* objArray = ptRegexp.MatchS(hltPath_);
0077 if (objArray->GetEntriesFast() >= 2) {
0078 auto* ptCutString = (TObjString*)objArray->At(1);
0079 cutMinPt_ = atoi(ptCutString->GetString());
0080 cutMinPt_ = ceil(cutMinPt_ * plotCuts_["minPtFactor"]);
0081 }
0082 delete objArray;
0083 }
0084
0085 void HLTMuonMatchAndPlot::beginRun(DQMStore::IBooker& iBooker, const edm::Run& iRun, const edm::EventSetup& iSetup) {
0086 TPRegexp suffixPtCut("Mu[0-9]+$");
0087
0088 string baseDir = destination_;
0089 if (baseDir[baseDir.size() - 1] != '/')
0090 baseDir += '/';
0091 string pathSansSuffix = hltPath_;
0092 if (hltPath_.rfind("_v") < hltPath_.length())
0093 pathSansSuffix = hltPath_.substr(0, hltPath_.rfind("_v"));
0094
0095 if (isLastFilter_)
0096 iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0097 else
0098 iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0099
0100
0101
0102
0103 for (const auto& suffix : EFFICIENCY_SUFFIXES) {
0104 if (isLastFilter_)
0105 iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0106 else
0107 iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0108
0109 book1D(iBooker, "efficiencyEta_" + suffix, "eta", ";#eta;");
0110 book1D(iBooker, "efficiencyPhi_" + suffix, "phi", ";#phi;");
0111 book1D(iBooker, "efficiencyTurnOn_" + suffix, "pt", ";p_{T};");
0112
0113 if (isLastFilter_)
0114 iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0115 else
0116 iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0117
0118 if (!isLastFilter_)
0119 continue;
0120
0121 book1D(iBooker, "efficiencyCharge_" + suffix, "charge", ";charge;");
0122 }
0123 }
0124
0125 void HLTMuonMatchAndPlot::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0126
0127 void HLTMuonMatchAndPlot::analyze(Handle<MuonCollection>& allMuons,
0128 Handle<BeamSpot>& beamSpot,
0129 Handle<VertexCollection>& vertices,
0130 Handle<TriggerEvent>& triggerSummary,
0131 Handle<TriggerResults>& triggerResults,
0132 const edm::TriggerNames& trigNames) {
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173 MuonCollection targetMuons =
0174 selectedMuons(*allMuons, *beamSpot, hasTargetRecoCuts, targetMuonSelector_, targetD0Cut_, targetZ0Cut_);
0175 MuonCollection probeMuons =
0176 selectedMuons(*allMuons, *beamSpot, hasProbeRecoCuts, probeMuonSelector_, probeD0Cut_, probeZ0Cut_);
0177 TriggerObjectCollection allTriggerObjects = triggerSummary->getObjects();
0178 TriggerObjectCollection hltMuons =
0179 selectedTriggerObjects(allTriggerObjects, *triggerSummary, hasTriggerCuts_, triggerSelector_);
0180
0181
0182
0183 vector<size_t> matches = matchByDeltaR(targetMuons, hltMuons, plotCuts_[triggerLevel_ + "DeltaR"]);
0184
0185
0186 for (size_t i = 0; i < targetMuons.size(); i++) {
0187 Muon& muon = targetMuons[i];
0188
0189
0190 for (const auto& suffix : EFFICIENCY_SUFFIXES) {
0191
0192 if (suffix == "numer" && matches[i] >= targetMuons.size())
0193 continue;
0194
0195 if (muon.pt() > cutMinPt_) {
0196 hists_["efficiencyEta_" + suffix]->Fill(muon.eta());
0197 }
0198
0199 if (fabs(muon.eta()) < plotCuts_["maxEta"]) {
0200 hists_["efficiencyTurnOn_" + suffix]->Fill(muon.pt());
0201 }
0202
0203 if (muon.pt() > cutMinPt_ && fabs(muon.eta()) < plotCuts_["maxEta"]) {
0204 const Track* track = nullptr;
0205 if (muon.isTrackerMuon())
0206 track = &*muon.innerTrack();
0207 else if (muon.isStandAloneMuon())
0208 track = &*muon.outerTrack();
0209 if (track) {
0210 hists_["efficiencyPhi_" + suffix]->Fill(muon.phi());
0211
0212 if (isLastFilter_) {
0213 hists_["efficiencyCharge_" + suffix]->Fill(muon.charge());
0214 }
0215 }
0216 }
0217 }
0218
0219 if (!isLastFilter_)
0220 continue;
0221
0222
0223 if (matches[i] >= targetMuons.size())
0224 continue;
0225 }
0226
0227
0228
0229
0230
0231 if (!isLastFilter_)
0232 return;
0233 unsigned int numTriggers = trigNames.size();
0234 bool passTrigger = false;
0235 if (requiredTriggers_.empty())
0236 passTrigger = true;
0237 for (auto const& requiredTrigger : requiredTriggers_) {
0238 for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0239 passTrigger = (trigNames.triggerName(hltIndex).find(requiredTrigger) != std::string::npos &&
0240 triggerResults->wasrun(hltIndex) && triggerResults->accept(hltIndex));
0241 if (passTrigger)
0242 break;
0243 }
0244 }
0245
0246 int nMatched = 0;
0247 for (unsigned long matche : matches) {
0248 if (matche < targetMuons.size())
0249 nMatched++;
0250 }
0251
0252 string nonSameSignPath = hltPath_;
0253 bool ssPath = false;
0254 if (nonSameSignPath.rfind("_SameSign") < nonSameSignPath.length()) {
0255 ssPath = true;
0256 nonSameSignPath = boost::replace_all_copy<string>(nonSameSignPath, "_SameSign", "");
0257 nonSameSignPath = nonSameSignPath.substr(0, nonSameSignPath.rfind("_v") + 2);
0258 }
0259 bool passTriggerSS = false;
0260 if (ssPath) {
0261 for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0262 passTriggerSS =
0263 passTriggerSS || (trigNames.triggerName(hltIndex).substr(0, nonSameSignPath.size()) == nonSameSignPath &&
0264 triggerResults->wasrun(hltIndex) && triggerResults->accept(hltIndex));
0265 }
0266 }
0267
0268 string nonDZPath = hltPath_;
0269 bool dzPath = false;
0270 if (nonDZPath.rfind("_DZ") < nonDZPath.length()) {
0271 dzPath = true;
0272 nonDZPath = boost::replace_all_copy<string>(nonDZPath, "_DZ", "");
0273 nonDZPath = nonDZPath.substr(0, nonDZPath.rfind("_v") + 2);
0274 }
0275 bool passTriggerDZ = false;
0276 if (dzPath) {
0277 for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0278 passTriggerDZ = passTriggerDZ || (trigNames.triggerName(hltIndex).find(nonDZPath) != std::string::npos &&
0279 triggerResults->wasrun(hltIndex) && triggerResults->accept(hltIndex));
0280 }
0281 }
0282
0283 }
0284
0285
0286 void HLTMuonMatchAndPlot::fillEdges(size_t& nBins, float*& edges, const vector<double>& binning) {
0287 if (binning.size() < 3) {
0288 LogWarning("HLTMuonVal") << "Invalid binning parameters!";
0289 return;
0290 }
0291
0292
0293 if (binning.size() == 3) {
0294 nBins = binning[0];
0295 edges = new float[nBins + 1];
0296 const double min = binning[1];
0297 const double binwidth = (binning[2] - binning[1]) / nBins;
0298 for (size_t i = 0; i <= nBins; i++)
0299 edges[i] = min + (binwidth * i);
0300 }
0301
0302
0303 else {
0304 nBins = binning.size() - 1;
0305 edges = new float[nBins + 1];
0306 for (size_t i = 0; i <= nBins; i++)
0307 edges[i] = binning[i];
0308 }
0309 }
0310
0311
0312
0313
0314 template <class T>
0315 void HLTMuonMatchAndPlot::fillMapFromPSet(map<string, T>& m, const ParameterSet& pset, const string& target) {
0316
0317 ParameterSet targetPset;
0318 if (pset.existsAs<ParameterSet>(target, true))
0319 targetPset = pset.getParameterSet(target);
0320 else if (pset.existsAs<ParameterSet>(target, false))
0321 targetPset = pset.getUntrackedParameterSet(target);
0322
0323
0324 vector<string> names = targetPset.getParameterNames();
0325 vector<string>::const_iterator iter;
0326
0327 for (iter = names.begin(); iter != names.end(); ++iter) {
0328 if (targetPset.existsAs<T>(*iter, true))
0329 m[*iter] = targetPset.getParameter<T>(*iter);
0330 else if (targetPset.existsAs<T>(*iter, false))
0331 m[*iter] = targetPset.getUntrackedParameter<T>(*iter);
0332 }
0333 }
0334
0335
0336 template <class T1, class T2>
0337 vector<size_t> HLTMuonMatchAndPlot::matchByDeltaR(const vector<T1>& collection1,
0338 const vector<T2>& collection2,
0339 const double maxDeltaR) {
0340 const size_t n1 = collection1.size();
0341 const size_t n2 = collection2.size();
0342
0343 vector<size_t> result(n1, -1);
0344 vector<vector<double> > deltaRMatrix(n1, vector<double>(n2, NOMATCH));
0345
0346 for (size_t i = 0; i < n1; i++)
0347 for (size_t j = 0; j < n2; j++) {
0348 deltaRMatrix[i][j] = deltaR(collection1[i], collection2[j]);
0349 }
0350
0351
0352 for (size_t k = 0; k < n1; k++) {
0353 size_t i_min = -1;
0354 size_t j_min = -1;
0355 double minDeltaR = maxDeltaR;
0356
0357 for (size_t i = 0; i < n1; i++)
0358 for (size_t j = 0; j < n2; j++)
0359 if (deltaRMatrix[i][j] < minDeltaR) {
0360 i_min = i;
0361 j_min = j;
0362 minDeltaR = deltaRMatrix[i][j];
0363 }
0364
0365 if (minDeltaR < maxDeltaR) {
0366 result[i_min] = j_min;
0367 deltaRMatrix[i_min] = vector<double>(n2, NOMATCH);
0368 for (size_t i = 0; i < n1; i++)
0369 deltaRMatrix[i][j_min] = NOMATCH;
0370 }
0371 }
0372
0373 return result;
0374 }
0375
0376 MuonCollection HLTMuonMatchAndPlot::selectedMuons(const MuonCollection& allMuons,
0377 const BeamSpot& beamSpot,
0378 bool hasRecoCuts,
0379 const StringCutObjectSelector<reco::Muon>& selector,
0380 double d0Cut,
0381 double z0Cut) {
0382
0383 if (!hasRecoCuts)
0384 return MuonCollection();
0385
0386 MuonCollection reducedMuons;
0387 for (auto const& mu : allMuons) {
0388 const Track* track = nullptr;
0389 if (mu.isTrackerMuon())
0390 track = &*mu.innerTrack();
0391 else if (mu.isStandAloneMuon())
0392 track = &*mu.outerTrack();
0393 if (track && selector(mu) && fabs(track->dxy(beamSpot.position())) < d0Cut &&
0394 fabs(track->dz(beamSpot.position())) < z0Cut)
0395 reducedMuons.push_back(mu);
0396 }
0397
0398 return reducedMuons;
0399 }
0400
0401 TriggerObjectCollection HLTMuonMatchAndPlot::selectedTriggerObjects(
0402 const TriggerObjectCollection& triggerObjects,
0403 const TriggerEvent& triggerSummary,
0404 bool hasTriggerCuts,
0405 const StringCutObjectSelector<TriggerObject>& triggerSelector) {
0406 if (!hasTriggerCuts)
0407 return TriggerObjectCollection();
0408
0409 InputTag filterTag(moduleLabel_, "", hltProcessName_);
0410 size_t filterIndex = triggerSummary.filterIndex(filterTag);
0411
0412 TriggerObjectCollection selectedObjects;
0413
0414 if (filterIndex < triggerSummary.sizeFilters()) {
0415 const Keys& keys = triggerSummary.filterKeys(filterIndex);
0416 for (unsigned short key : keys) {
0417 TriggerObject foundObject = triggerObjects[key];
0418 if (triggerSelector(foundObject))
0419 selectedObjects.push_back(foundObject);
0420 }
0421 }
0422
0423 return selectedObjects;
0424 }
0425
0426 void HLTMuonMatchAndPlot::book1D(DQMStore::IBooker& iBooker, string name, const string& binningType, string title) {
0427
0428
0429
0430
0431
0432
0433 size_t nBins;
0434 float* edges = nullptr;
0435 fillEdges(nBins, edges, binParams_[binningType]);
0436 hists_[name] = iBooker.book1D(name, title, nBins, edges);
0437 if (hists_[name])
0438 if (hists_[name]->getTH1F()->GetSumw2N())
0439 hists_[name]->enableSumw2();
0440
0441 if (edges)
0442 delete[] edges;
0443 }
0444
0445 void HLTMuonMatchAndPlot::book2D(DQMStore::IBooker& iBooker,
0446 const string& name,
0447 const string& binningTypeX,
0448 const string& binningTypeY,
0449 const string& title) {
0450
0451
0452
0453
0454
0455
0456 size_t nBinsX;
0457 float* edgesX = nullptr;
0458 fillEdges(nBinsX, edgesX, binParams_[binningTypeX]);
0459
0460 size_t nBinsY;
0461 float* edgesY = nullptr;
0462 fillEdges(nBinsY, edgesY, binParams_[binningTypeY]);
0463
0464 hists_[name] = iBooker.book2D(name.c_str(), title.c_str(), nBinsX, edgesX, nBinsY, edgesY);
0465 if (hists_[name])
0466 if (hists_[name]->getTH2F()->GetSumw2N())
0467 hists_[name]->enableSumw2();
0468
0469 if (edgesX)
0470 delete[] edgesX;
0471 if (edgesY)
0472 delete[] edgesY;
0473 }