Trigger/src/HLTMuonMatchAndPlot.cc

0001 /** \file DQMOffline/Trigger/HLTMuonMatchAndPlot.cc
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 //////// Namespaces and Typedefs /////////////////////////////////////////////
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 //////// HLTMuonMatchAndPlot Class Members ///////////////////////////////////
0035
0036 /// Constructor
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       targetMuonEtaMax_(targetParams_.getUntrackedParameter<double>("recoMaxEtaCut", 0.)),
0047       targetMuonEtaMin_(targetParams_.getUntrackedParameter<double>("recoMinEtaCut", 0.)),
0048       targetIsMuonGlb_(targetParams_.getUntrackedParameter<bool>("recoGlbMuCut", false)),
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       probeMuonEtaMax_(probeParams_.getUntrackedParameter<double>("recoMaxEtaCut", 0.)),
0054       probeMuonEtaMin_(probeParams_.getUntrackedParameter<double>("recoMinEtaCut", 0.)),
0055       probeIsMuonGlb_(probeParams_.getUntrackedParameter<bool>("recoGlbMuCut", false)),
0056       probeZ0Cut_(probeParams_.getUntrackedParameter<double>("z0Cut", 0.)),
0057       probeD0Cut_(probeParams_.getUntrackedParameter<double>("d0Cut", 0.)),
0058       triggerEtaMaxCut_(targetParams_.getUntrackedParameter<double>("hltMaxEtaCut", 0.)),
0059       triggerEtaMinCut_(targetParams_.getUntrackedParameter<double>("hltMinEtaCut", 0.)) {
0060   // Create std::map<string, T> from ParameterSets.
0061   fillMapFromPSet(binParams_, pset, "binParams");
0062   fillMapFromPSet(plotCuts_, pset, "plotCuts");
0063
0064   // Get the trigger level.
0065   triggerLevel_ = "L3";
0066   TPRegexp levelRegexp("L[1-3]");
0067   //  size_t nModules = moduleLabels_.size();
0068   //  cout << moduleLabel_ << " " << hltPath_ << endl;
0069   TObjArray* levelArray = levelRegexp.MatchS(moduleLabel_);
0070   if (levelArray->GetEntriesFast() > 0) {
0071     triggerLevel_ = static_cast<const char*>(((TObjString*)levelArray->At(0))->GetString());
0072   }
0073   delete levelArray;
0074
0075   // Get the pT cut by parsing the name of the HLT path.
0076   cutMinPt_ = 3;
0077   TPRegexp ptRegexp("Mu([0-9]*)");
0078   TObjArray* objArray = ptRegexp.MatchS(hltPath_);
0079   if (objArray->GetEntriesFast() >= 2) {
0080     auto* ptCutString = (TObjString*)objArray->At(1);
0081     cutMinPt_ = atoi(ptCutString->GetString());
0082     cutMinPt_ = ceil(cutMinPt_ * plotCuts_["minPtFactor"]);
0083   }
0084   delete objArray;
0085 }
0086
0087 void HLTMuonMatchAndPlot::beginRun(DQMStore::IBooker& iBooker, const edm::Run& iRun, const edm::EventSetup& iSetup) {
0088   TPRegexp suffixPtCut("Mu[0-9]+$");
0089
0090   string baseDir = destination_;
0091   if (baseDir[baseDir.size() - 1] != '/')
0092     baseDir += '/';
0093   string pathSansSuffix = hltPath_;
0094   if (hltPath_.rfind("_v") < hltPath_.length())
0095     pathSansSuffix = hltPath_.substr(0, hltPath_.rfind("_v"));
0096
0097   if (isLastFilter_)
0098     iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0099   else
0100     iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0101
0102   // Form is book1D(name, binningType, title) where 'binningType' is used
0103   // to fetch the bin settings from binParams_.
0104
0105   for (const auto& suffix : EFFICIENCY_SUFFIXES) {
0106     if (isLastFilter_)
0107       iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0108     else
0109       iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0110
0111     book1D(iBooker, "efficiencyEta_" + suffix, "eta", ";#eta;");
0112     book1D(iBooker, "efficiencyPhi_" + suffix, "phi", ";#phi;");
0113     book1D(iBooker, "efficiencyTurnOn_" + suffix, "pt", ";p_{T};");
0114     book1D(iBooker, "efficiencyNVertex_" + suffix, "NVertex", ";NVertex;");
0115
0116     if (isLastFilter_)
0117       iBooker.setCurrentFolder(baseDir + pathSansSuffix);
0118     else
0119       iBooker.setCurrentFolder(baseDir + pathSansSuffix + "/" + moduleLabel_);
0120
0121     if (!isLastFilter_)
0122       continue;  //this will be plotted only for the last filter
0123
0124     book1D(iBooker, "efficiencyCharge_" + suffix, "charge", ";charge;");
0125     book1D(iBooker, "efficiencyZ0_" + suffix, "z0", ";z0;");
0126     book1D(iBooker, "efficiency_DZ_Mu_" + suffix, "z0", ";z0;");
0127   }
0128 }
0129
0130 void HLTMuonMatchAndPlot::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0131
0132 void HLTMuonMatchAndPlot::analyze(Handle<MuonCollection>& allMuons,
0133                                   Handle<BeamSpot>& beamSpot,
0134                                   Handle<VertexCollection>& vertices,
0135                                   Handle<TriggerEvent>& triggerSummary,
0136                                   Handle<TriggerResults>& triggerResults,
0137                                   const edm::TriggerNames& trigNames) {
0138   /*
0139   if(gen != 0) {
0140     for(g_part = gen->begin(); g_part != gen->end(); g_part++){
0141       if( abs(g_part->pdgId()) ==  13) {
0142         if(!( g_part->status() ==1 || (g_part->status() ==2 && abs(g_part->pdgId())==5))) continue;
0143         bool GenMomExists  (true);
0144         bool GenGrMomExists(true);
0145         if( g_part->pt() < 10.0 )  continue;
0146         if( fabs(g_part->eta()) > 2.4 ) continue;
0147         int gen_id= g_part->pdgId();
0148         const GenParticle* gen_lept = &(*g_part);
0149         // get mother of gen_lept
0150         const GenParticle* gen_mom = static_cast<const GenParticle*> (gen_lept->mother());
0151         if(gen_mom==NULL) GenMomExists=false;
0152         int m_id=-999;
0153         if(GenMomExists) m_id = gen_mom->pdgId();
0154         if(m_id != gen_id || !GenMomExists);
0155         else{
0156           int id= m_id;
0157           while(id == gen_id && GenMomExists){
0158             gen_mom = static_cast<const GenParticle*> (gen_mom->mother());
0159             if(gen_mom==NULL){
0160               GenMomExists=false;
0161             }
0162             if(GenMomExists) id=gen_mom->pdgId();
0163           }
0164         }
0165         if(GenMomExists) m_id = gen_mom->pdgId();
0166         gen_leptsp.push_back(gen_lept);
0167         gen_momsp.push_back(gen_mom);
0168       }
0169     }
0170   }
0171   std::vector<GenParticle> gen_lepts;
0172   for(size_t i = 0; i < gen_leptsp.size(); i++) {
0173     gen_lepts.push_back(*gen_leptsp[i]);
0174   }
0175   */
0176
0177   // Select objects based on the configuration.
0178   MuonCollection targetMuons = selectedMuons(*allMuons, *beamSpot, true);
0179   MuonCollection probeMuons = selectedMuons(*allMuons, *beamSpot, false);
0180   TriggerObjectCollection allTriggerObjects = triggerSummary->getObjects();
0181   TriggerObjectCollection hltMuons = selectedTriggerObjects(allTriggerObjects, *triggerSummary);
0182   // Fill plots for HLT muons.
0183
0184   // Find the best trigger object matches for the targetMuons.
0185   vector<size_t> matches = matchByDeltaR(targetMuons, hltMuons, plotCuts_[triggerLevel_ + "DeltaR"]);
0186
0187   // Fill plots for matched muons.
0188   for (size_t i = 0; i < targetMuons.size(); i++) {
0189     Muon& muon = targetMuons[i];
0190
0191     // Fill numerators and denominators for efficiency plots.
0192     for (const auto& suffix : EFFICIENCY_SUFFIXES) {
0193       // If no match was found, then the numerator plots don't get filled.
0194       if (suffix == "numer" && matches[i] >= targetMuons.size())
0195         continue;
0196
0197       if (muon.pt() > cutMinPt_) {
0198         hists_["efficiencyEta_" + suffix]->Fill(muon.eta());
0199       }
0200
0201       if (fabs(muon.eta()) < plotCuts_["maxEta"]) {
0202         hists_["efficiencyTurnOn_" + suffix]->Fill(muon.pt());
0203       }
0204
0205       if (muon.pt() > cutMinPt_ && fabs(muon.eta()) < plotCuts_["maxEta"]) {
0206         const Track* track = nullptr;
0207         if (muon.isTrackerMuon())
0208           track = &*muon.innerTrack();
0209         else if (muon.isStandAloneMuon())
0210           track = &*muon.outerTrack();
0211         if (track) {
0212           hists_["efficiencyNVertex_" + suffix]->Fill(vertices->size());
0213           hists_["efficiencyPhi_" + suffix]->Fill(muon.phi());
0214
0215           if (isLastFilter_) {
0216             hists_["efficiencyZ0_" + suffix]->Fill(track->dz(beamSpot->position()));
0217             hists_["efficiencyCharge_" + suffix]->Fill(muon.charge());
0218           }
0219         }
0220       }
0221     }  // finish loop numerator / denominator...
0222
0223     if (!isLastFilter_)
0224       continue;
0225     // Fill plots for tag and probe
0226     // Muon cannot be a tag because doesn't match an hlt muon
0227     if (matches[i] >= targetMuons.size())
0228       continue;
0229   }  // End loop over targetMuons.
0230
0231   // fill eff histograms for reference trigger method
0232   // Denominator: events passing reference trigger and two target muons
0233   // Numerator:   events in the denominator with two target muons
0234   // matched to hlt muons
0235   if (!isLastFilter_)
0236     return;
0237   unsigned int numTriggers = trigNames.size();
0238
0239   int nMatched = 0;
0240   for (unsigned long matche : matches) {
0241     if (matche < targetMuons.size())
0242       nMatched++;
0243   }
0244
0245   string nonDZPath = hltPath_;
0246   bool dzPath = false;
0247   if (nonDZPath.rfind("_DZ") < nonDZPath.length()) {
0248     dzPath = true;
0249     nonDZPath = boost::replace_all_copy<string>(nonDZPath, "_DZ", "");
0250     nonDZPath = nonDZPath.substr(0, nonDZPath.rfind("_v") + 2);
0251   }
0252   bool passTriggerDZ = false;
0253
0254   if (dzPath) {
0255     for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0256       passTriggerDZ = passTriggerDZ || (trigNames.triggerName(hltIndex).find(nonDZPath) != std::string::npos &&
0257                                         triggerResults->wasrun(hltIndex) && triggerResults->accept(hltIndex));
0258     }
0259   }
0260   if (dzPath && targetMuons.size() > 1 && passTriggerDZ) {
0261     const Track* track0 = nullptr;
0262     const Track* track1 = nullptr;
0263     if (targetMuons.at(0).isTrackerMuon())
0264       track0 = &*targetMuons.at(0).innerTrack();
0265     else if (targetMuons.at(0).isStandAloneMuon())
0266       track0 = &*targetMuons.at(0).outerTrack();
0267     if (targetMuons.at(1).isTrackerMuon())
0268       track1 = &*targetMuons.at(1).innerTrack();
0269     else if (targetMuons.at(1).isStandAloneMuon())
0270       track1 = &*targetMuons.at(1).outerTrack();
0271
0272     if (track0 && track1) {
0273       hists_["efficiency_DZ_Mu_denom"]->Fill(track0->dz(beamSpot->position()) - track1->dz(beamSpot->position()));
0274       if (nMatched > 1) {
0275         hists_["efficiency_DZ_Mu_numer"]->Fill(track0->dz(beamSpot->position()) - track1->dz(beamSpot->position()));
0276       }
0277     }
0278   }
0279
0280 }  // End analyze() method.
0281
0282 // Method to fill binning parameters from a vector of doubles.
0283 bool HLTMuonMatchAndPlot::fillEdges(size_t& nBins, float*& edges, const vector<double>& binning) {
0284   if (binning.size() < 3) {
0285     LogWarning("HLTMuonVal") << "Invalid binning parameters!";
0286     return false;
0287   }
0288
0289   // Fixed-width binning.
0290   if (binning.size() == 3) {
0291     nBins = binning[0];
0292     edges = new float[nBins + 1];
0293     const double min = binning[1];
0294     const double binwidth = (binning[2] - binning[1]) / nBins;
0295     for (size_t i = 0; i <= nBins; i++)
0296       edges[i] = min + (binwidth * i);
0297   }
0298
0299   // Variable-width binning.
0300   else {
0301     nBins = binning.size() - 1;
0302     edges = new float[nBins + 1];
0303     for (size_t i = 0; i <= nBins; i++)
0304       edges[i] = binning[i];
0305   }
0306   return true;
0307 }
0308
0309 // This is an unorthodox method of getting parameters, but cleaner in my mind
0310 // It looks for parameter 'target' in the pset, and assumes that all entries
0311 // have the same class (T), filling the values into a std::map.
0312 template <class T>
0313 void HLTMuonMatchAndPlot::fillMapFromPSet(map<string, T>& m, const ParameterSet& pset, const string& target) {
0314   // Get the ParameterSet with name 'target' from 'pset'
0315   ParameterSet targetPset;
0316   if (pset.existsAs<ParameterSet>(target, true))  // target is tracked
0317     targetPset = pset.getParameterSet(target);
0318   else if (pset.existsAs<ParameterSet>(target, false))  // target is untracked
0319     targetPset = pset.getUntrackedParameterSet(target);
0320
0321   // Get the parameter names from targetPset
0322   vector<string> names = targetPset.getParameterNames();
0323   vector<string>::const_iterator iter;
0324
0325   for (iter = names.begin(); iter != names.end(); ++iter) {
0326     if (targetPset.existsAs<T>(*iter, true))  // target is tracked
0327       m[*iter] = targetPset.getParameter<T>(*iter);
0328     else if (targetPset.existsAs<T>(*iter, false))  // target is untracked
0329       m[*iter] = targetPset.getUntrackedParameter<T>(*iter);
0330   }
0331 }
0332
0333 // A generic method to find the best deltaR matches from 2 collections.
0334 template <class T1, class T2>
0335 vector<size_t> HLTMuonMatchAndPlot::matchByDeltaR(const vector<T1>& collection1,
0336                                                   const vector<T2>& collection2,
0337                                                   const double maxDeltaR) {
0338   const size_t n1 = collection1.size();
0339   const size_t n2 = collection2.size();
0340
0341   vector<size_t> result(n1, -1);
0342   vector<vector<double> > deltaRMatrix(n1, vector<double>(n2, NOMATCH));
0343
0344   for (size_t i = 0; i < n1; i++)
0345     for (size_t j = 0; j < n2; j++) {
0346       deltaRMatrix[i][j] = deltaR(collection1[i], collection2[j]);
0347     }
0348
0349   // Run through the matrix n1 times to make sure we've found all matches.
0350   for (size_t k = 0; k < n1; k++) {
0351     size_t i_min = -1;
0352     size_t j_min = -1;
0353     double minDeltaR = maxDeltaR;
0354     // find the smallest deltaR
0355     for (size_t i = 0; i < n1; i++)
0356       for (size_t j = 0; j < n2; j++)
0357         if (deltaRMatrix[i][j] < minDeltaR) {
0358           i_min = i;
0359           j_min = j;
0360           minDeltaR = deltaRMatrix[i][j];
0361         }
0362     // If a match has been made, save it and make those candidates unavailable.
0363     if (minDeltaR < maxDeltaR) {
0364       result[i_min] = j_min;
0365       deltaRMatrix[i_min] = vector<double>(n2, NOMATCH);
0366       for (size_t i = 0; i < n1; i++)
0367         deltaRMatrix[i][j_min] = NOMATCH;
0368     }
0369   }
0370
0371   return result;
0372 }
0373
0374 MuonCollection HLTMuonMatchAndPlot::selectedMuons(const MuonCollection& allMuons,
0375                                                   const BeamSpot& beamSpot,
0376                                                   bool isTargetMuons) {
0377   MuonCollection reducedMuons;
0378   double RecoMuonEtaMax = isTargetMuons ? targetMuonEtaMax_ : probeMuonEtaMax_;
0379   double RecoMuonEtaMin = isTargetMuons ? targetMuonEtaMin_ : probeMuonEtaMin_;
0380   bool IsMuonGlb = isTargetMuons ? targetIsMuonGlb_ : probeIsMuonGlb_;
0381   double d0Cut = isTargetMuons ? targetD0Cut_ : probeD0Cut_;
0382   double z0Cut = isTargetMuons ? targetZ0Cut_ : probeZ0Cut_;
0383
0384   for (auto const& mu : allMuons) {
0385     const Track* track = nullptr;
0386     if (mu.isTrackerMuon())
0387       track = &*mu.innerTrack();
0388     else if (mu.isStandAloneMuon())
0389       track = &*mu.outerTrack();
0390     // minimun ID (requested for cosmics) is being a StandAlone muon
0391     bool muID = IsMuonGlb ? mu.isGlobalMuon() : mu.isStandAloneMuon();
0392     if (track && muID && abs(mu.eta()) < RecoMuonEtaMax && abs(mu.eta()) >= RecoMuonEtaMin &&
0393         fabs(track->dxy(beamSpot.position())) < d0Cut && fabs(track->dz(beamSpot.position())) < z0Cut)
0394       reducedMuons.push_back(mu);
0395   }
0396
0397   return reducedMuons;
0398 }
0399
0400 TriggerObjectCollection HLTMuonMatchAndPlot::selectedTriggerObjects(const TriggerObjectCollection& triggerObjects,
0401                                                                     const TriggerEvent& triggerSummary) {
0402   InputTag filterTag(moduleLabel_, "", hltProcessName_);
0403   size_t filterIndex = triggerSummary.filterIndex(filterTag);
0404
0405   TriggerObjectCollection selectedObjects;
0406   if (filterIndex < triggerSummary.sizeFilters()) {
0407     const Keys& keys = triggerSummary.filterKeys(filterIndex);
0408     for (unsigned short key : keys) {
0409       TriggerObject foundObject = triggerObjects[key];
0410       if (abs(foundObject.eta()) < triggerEtaMaxCut_ && abs(foundObject.eta()) >= triggerEtaMinCut_)
0411         selectedObjects.push_back(foundObject);
0412     }
0413   }
0414
0415   return selectedObjects;
0416 }
0417
0418 void HLTMuonMatchAndPlot::book1D(DQMStore::IBooker& iBooker, string name, const string& binningType, string title) {
0419   /* Properly delete the array of floats that has been allocated on
0420    * the heap by fillEdges.  Avoid multiple copies and internal ROOT
0421    * clones by simply creating the histograms directly in the DQMStore
0422    * using the appropriate book1D method to handle the variable bins
0423    * case. */
0424
0425   size_t nBins;
0426   float* edges = nullptr;
0427   bool bookhist = fillEdges(nBins, edges, binParams_[binningType]);
0428   if (bookhist) {
0429     hists_[name] = iBooker.book1D(name, title, nBins, edges);
0430     if (hists_[name]->getTH1F()->GetSumw2N())
0431       hists_[name]->enableSumw2();
0432
0433     delete[] edges;
0434   }
0435 }
0436
0437 void HLTMuonMatchAndPlot::book2D(DQMStore::IBooker& iBooker,
0438                                  const string& name,
0439                                  const string& binningTypeX,
0440                                  const string& binningTypeY,
0441                                  const string& title) {
0442   /* Properly delete the arrays of floats that have been allocated on
0443    * the heap by fillEdges.  Avoid multiple copies and internal ROOT
0444    * clones by simply creating the histograms directly in the DQMStore
0445    * using the appropriate book2D method to handle the variable bins
0446    * case. */
0447
0448   size_t nBinsX;
0449   float* edgesX = nullptr;
0450   bool bookhist = fillEdges(nBinsX, edgesX, binParams_[binningTypeX]);
0451
0452   size_t nBinsY;
0453   float* edgesY = nullptr;
0454   bookhist &= fillEdges(nBinsY, edgesY, binParams_[binningTypeY]);
0455   if (bookhist) {
0456     hists_[name] = iBooker.book2D(name.c_str(), title.c_str(), nBinsX, edgesX, nBinsY, edgesY);
0457     if (hists_[name]->getTH2F()->GetSumw2N())
0458       hists_[name]->enableSumw2();
0459
0460     delete[] edgesX;
0461     delete[] edgesY;
0462   }
0463 }