File indexing completed on 2023-03-17 10:58:50
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 book1D(iBooker, "efficiencyNVertex_" + suffix, "NVertex", ";NVertex;");
0113
0114 if (isLastFilter_)
0115 iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0116 else
0117 iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0118
0119 if (!isLastFilter_)
0120 continue;
0121
0122 book1D(iBooker, "efficiencyCharge_" + suffix, "charge", ";charge;");
0123 book1D(iBooker, "efficiencyZ0_" + suffix, "z0", ";z0;");
0124 book1D(iBooker, "efficiency_DZ_Mu_" + suffix, "z0", ";z0;");
0125 }
0126 }
0127
0128 void HLTMuonMatchAndPlot::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0129
0130 void HLTMuonMatchAndPlot::analyze(Handle<MuonCollection>& allMuons,
0131 Handle<BeamSpot>& beamSpot,
0132 Handle<VertexCollection>& vertices,
0133 Handle<TriggerEvent>& triggerSummary,
0134 Handle<TriggerResults>& triggerResults,
0135 const edm::TriggerNames& trigNames) {
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
0174
0175
0176 MuonCollection targetMuons =
0177 selectedMuons(*allMuons, *beamSpot, hasTargetRecoCuts, targetMuonSelector_, targetD0Cut_, targetZ0Cut_);
0178 MuonCollection probeMuons =
0179 selectedMuons(*allMuons, *beamSpot, hasProbeRecoCuts, probeMuonSelector_, probeD0Cut_, probeZ0Cut_);
0180 TriggerObjectCollection allTriggerObjects = triggerSummary->getObjects();
0181 TriggerObjectCollection hltMuons =
0182 selectedTriggerObjects(allTriggerObjects, *triggerSummary, hasTriggerCuts_, triggerSelector_);
0183
0184
0185
0186 vector<size_t> matches = matchByDeltaR(targetMuons, hltMuons, plotCuts_[triggerLevel_ + "DeltaR"]);
0187
0188
0189 for (size_t i = 0; i < targetMuons.size(); i++) {
0190 Muon& muon = targetMuons[i];
0191
0192
0193 for (const auto& suffix : EFFICIENCY_SUFFIXES) {
0194
0195 if (suffix == "numer" && matches[i] >= targetMuons.size())
0196 continue;
0197
0198 if (muon.pt() > cutMinPt_) {
0199 hists_["efficiencyEta_" + suffix]->Fill(muon.eta());
0200 }
0201
0202 if (fabs(muon.eta()) < plotCuts_["maxEta"]) {
0203 hists_["efficiencyTurnOn_" + suffix]->Fill(muon.pt());
0204 }
0205
0206 if (muon.pt() > cutMinPt_ && fabs(muon.eta()) < plotCuts_["maxEta"]) {
0207 const Track* track = nullptr;
0208 if (muon.isTrackerMuon())
0209 track = &*muon.innerTrack();
0210 else if (muon.isStandAloneMuon())
0211 track = &*muon.outerTrack();
0212 if (track) {
0213 hists_["efficiencyNVertex_" + suffix]->Fill(vertices->size());
0214 hists_["efficiencyPhi_" + suffix]->Fill(muon.phi());
0215
0216 if (isLastFilter_) {
0217 hists_["efficiencyZ0_" + suffix]->Fill(track->dz(beamSpot->position()));
0218 hists_["efficiencyCharge_" + suffix]->Fill(muon.charge());
0219 }
0220 }
0221 }
0222 }
0223
0224 if (!isLastFilter_)
0225 continue;
0226
0227
0228 if (matches[i] >= targetMuons.size())
0229 continue;
0230 }
0231
0232
0233
0234
0235
0236 if (!isLastFilter_)
0237 return;
0238 unsigned int numTriggers = trigNames.size();
0239
0240 int nMatched = 0;
0241 for (unsigned long matche : matches) {
0242 if (matche < targetMuons.size())
0243 nMatched++;
0244 }
0245
0246 string nonDZPath = hltPath_;
0247 bool dzPath = false;
0248 if (nonDZPath.rfind("_DZ") < nonDZPath.length()) {
0249 dzPath = true;
0250 nonDZPath = boost::replace_all_copy<string>(nonDZPath, "_DZ", "");
0251 nonDZPath = nonDZPath.substr(0, nonDZPath.rfind("_v") + 2);
0252 }
0253 bool passTriggerDZ = false;
0254
0255 if (dzPath) {
0256 for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0257 passTriggerDZ = passTriggerDZ || (trigNames.triggerName(hltIndex).find(nonDZPath) != std::string::npos &&
0258 triggerResults->wasrun(hltIndex) && triggerResults->accept(hltIndex));
0259 }
0260 }
0261 if (dzPath && targetMuons.size() > 1 && passTriggerDZ) {
0262 const Track* track0 = nullptr;
0263 const Track* track1 = nullptr;
0264 if (targetMuons.at(0).isTrackerMuon())
0265 track0 = &*targetMuons.at(0).innerTrack();
0266 else if (targetMuons.at(0).isStandAloneMuon())
0267 track0 = &*targetMuons.at(0).outerTrack();
0268 if (targetMuons.at(1).isTrackerMuon())
0269 track1 = &*targetMuons.at(1).innerTrack();
0270 else if (targetMuons.at(1).isStandAloneMuon())
0271 track1 = &*targetMuons.at(1).outerTrack();
0272
0273 if (track0 && track1) {
0274 hists_["efficiency_DZ_Mu_denom"]->Fill(track0->dz(beamSpot->position()) - track1->dz(beamSpot->position()));
0275 if (nMatched > 1) {
0276 hists_["efficiency_DZ_Mu_numer"]->Fill(track0->dz(beamSpot->position()) - track1->dz(beamSpot->position()));
0277 }
0278 }
0279 }
0280
0281 }
0282
0283
0284 bool HLTMuonMatchAndPlot::fillEdges(size_t& nBins, float*& edges, const vector<double>& binning) {
0285 if (binning.size() < 3) {
0286 LogWarning("HLTMuonVal") << "Invalid binning parameters!";
0287 return false;
0288 }
0289
0290
0291 if (binning.size() == 3) {
0292 nBins = binning[0];
0293 edges = new float[nBins + 1];
0294 const double min = binning[1];
0295 const double binwidth = (binning[2] - binning[1]) / nBins;
0296 for (size_t i = 0; i <= nBins; i++)
0297 edges[i] = min + (binwidth * i);
0298 }
0299
0300
0301 else {
0302 nBins = binning.size() - 1;
0303 edges = new float[nBins + 1];
0304 for (size_t i = 0; i <= nBins; i++)
0305 edges[i] = binning[i];
0306 }
0307 return true;
0308 }
0309
0310
0311
0312
0313 template <class T>
0314 void HLTMuonMatchAndPlot::fillMapFromPSet(map<string, T>& m, const ParameterSet& pset, const string& target) {
0315
0316 ParameterSet targetPset;
0317 if (pset.existsAs<ParameterSet>(target, true))
0318 targetPset = pset.getParameterSet(target);
0319 else if (pset.existsAs<ParameterSet>(target, false))
0320 targetPset = pset.getUntrackedParameterSet(target);
0321
0322
0323 vector<string> names = targetPset.getParameterNames();
0324 vector<string>::const_iterator iter;
0325
0326 for (iter = names.begin(); iter != names.end(); ++iter) {
0327 if (targetPset.existsAs<T>(*iter, true))
0328 m[*iter] = targetPset.getParameter<T>(*iter);
0329 else if (targetPset.existsAs<T>(*iter, false))
0330 m[*iter] = targetPset.getUntrackedParameter<T>(*iter);
0331 }
0332 }
0333
0334
0335 template <class T1, class T2>
0336 vector<size_t> HLTMuonMatchAndPlot::matchByDeltaR(const vector<T1>& collection1,
0337 const vector<T2>& collection2,
0338 const double maxDeltaR) {
0339 const size_t n1 = collection1.size();
0340 const size_t n2 = collection2.size();
0341
0342 vector<size_t> result(n1, -1);
0343 vector<vector<double> > deltaRMatrix(n1, vector<double>(n2, NOMATCH));
0344
0345 for (size_t i = 0; i < n1; i++)
0346 for (size_t j = 0; j < n2; j++) {
0347 deltaRMatrix[i][j] = deltaR(collection1[i], collection2[j]);
0348 }
0349
0350
0351 for (size_t k = 0; k < n1; k++) {
0352 size_t i_min = -1;
0353 size_t j_min = -1;
0354 double minDeltaR = maxDeltaR;
0355
0356 for (size_t i = 0; i < n1; i++)
0357 for (size_t j = 0; j < n2; j++)
0358 if (deltaRMatrix[i][j] < minDeltaR) {
0359 i_min = i;
0360 j_min = j;
0361 minDeltaR = deltaRMatrix[i][j];
0362 }
0363
0364 if (minDeltaR < maxDeltaR) {
0365 result[i_min] = j_min;
0366 deltaRMatrix[i_min] = vector<double>(n2, NOMATCH);
0367 for (size_t i = 0; i < n1; i++)
0368 deltaRMatrix[i][j_min] = NOMATCH;
0369 }
0370 }
0371
0372 return result;
0373 }
0374
0375 MuonCollection HLTMuonMatchAndPlot::selectedMuons(const MuonCollection& allMuons,
0376 const BeamSpot& beamSpot,
0377 bool hasRecoCuts,
0378 const StringCutObjectSelector<reco::Muon>& selector,
0379 double d0Cut,
0380 double z0Cut) {
0381
0382 if (!hasRecoCuts)
0383 return MuonCollection();
0384
0385 MuonCollection reducedMuons;
0386 for (auto const& mu : allMuons) {
0387 const Track* track = nullptr;
0388 if (mu.isTrackerMuon())
0389 track = &*mu.innerTrack();
0390 else if (mu.isStandAloneMuon())
0391 track = &*mu.outerTrack();
0392 if (track && selector(mu) && fabs(track->dxy(beamSpot.position())) < d0Cut &&
0393 fabs(track->dz(beamSpot.position())) < z0Cut)
0394 reducedMuons.push_back(mu);
0395 }
0396
0397 return reducedMuons;
0398 }
0399
0400 TriggerObjectCollection HLTMuonMatchAndPlot::selectedTriggerObjects(
0401 const TriggerObjectCollection& triggerObjects,
0402 const TriggerEvent& triggerSummary,
0403 bool hasTriggerCuts,
0404 const StringCutObjectSelector<TriggerObject>& triggerSelector) {
0405 if (!hasTriggerCuts)
0406 return TriggerObjectCollection();
0407
0408 InputTag filterTag(moduleLabel_, "", hltProcessName_);
0409 size_t filterIndex = triggerSummary.filterIndex(filterTag);
0410
0411 TriggerObjectCollection selectedObjects;
0412
0413 if (filterIndex < triggerSummary.sizeFilters()) {
0414 const Keys& keys = triggerSummary.filterKeys(filterIndex);
0415 for (unsigned short key : keys) {
0416 TriggerObject foundObject = triggerObjects[key];
0417 if (triggerSelector(foundObject))
0418 selectedObjects.push_back(foundObject);
0419 }
0420 }
0421
0422 return selectedObjects;
0423 }
0424
0425 void HLTMuonMatchAndPlot::book1D(DQMStore::IBooker& iBooker, string name, const string& binningType, string title) {
0426
0427
0428
0429
0430
0431
0432 size_t nBins;
0433 float* edges = nullptr;
0434 bool bookhist = fillEdges(nBins, edges, binParams_[binningType]);
0435 if (bookhist) {
0436 hists_[name] = iBooker.book1D(name, title, nBins, edges);
0437 if (hists_[name]->getTH1F()->GetSumw2N())
0438 hists_[name]->enableSumw2();
0439
0440 delete[] edges;
0441 }
0442 }
0443
0444 void HLTMuonMatchAndPlot::book2D(DQMStore::IBooker& iBooker,
0445 const string& name,
0446 const string& binningTypeX,
0447 const string& binningTypeY,
0448 const string& title) {
0449
0450
0451
0452
0453
0454
0455 size_t nBinsX;
0456 float* edgesX = nullptr;
0457 bool bookhist = fillEdges(nBinsX, edgesX, binParams_[binningTypeX]);
0458
0459 size_t nBinsY;
0460 float* edgesY = nullptr;
0461 bookhist &= fillEdges(nBinsY, edgesY, binParams_[binningTypeY]);
0462 if (bookhist) {
0463 hists_[name] = iBooker.book2D(name.c_str(), title.c_str(), nBinsX, edgesX, nBinsY, edgesY);
0464 if (hists_[name]->getTH2F()->GetSumw2N())
0465 hists_[name]->enableSumw2();
0466
0467 delete[] edgesX;
0468 delete[] edgesY;
0469 }
0470 }