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       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   // Create std::map<string, T> from ParameterSets.
0059   fillMapFromPSet(binParams_, pset, "binParams");
0060   fillMapFromPSet(plotCuts_, pset, "plotCuts");
0061
0062   // Get the trigger level.
0063   triggerLevel_ = "L3";
0064   TPRegexp levelRegexp("L[1-3]");
0065   //  size_t nModules = moduleLabels_.size();
0066   //  cout << moduleLabel_ << " " << hltPath_ << endl;
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   // Get the pT cut by parsing the name of the HLT path.
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   // Form is book1D(name, binningType, title) where 'binningType' is used
0101   // to fetch the bin settings from binParams_.
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;  //this will be plotted only for the last filter
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   if(gen != 0) {
0138     for(g_part = gen->begin(); g_part != gen->end(); g_part++){
0139       if( abs(g_part->pdgId()) ==  13) {
0140         if(!( g_part->status() ==1 || (g_part->status() ==2 && abs(g_part->pdgId())==5))) continue;
0141         bool GenMomExists  (true);
0142         bool GenGrMomExists(true);
0143         if( g_part->pt() < 10.0 )  continue;
0144         if( fabs(g_part->eta()) > 2.4 ) continue;
0145         int gen_id= g_part->pdgId();
0146         const GenParticle* gen_lept = &(*g_part);
0147         // get mother of gen_lept
0148         const GenParticle* gen_mom = static_cast<const GenParticle*> (gen_lept->mother());
0149         if(gen_mom==NULL) GenMomExists=false;
0150         int m_id=-999;
0151         if(GenMomExists) m_id = gen_mom->pdgId();
0152         if(m_id != gen_id || !GenMomExists);
0153         else{
0154           int id= m_id;
0155           while(id == gen_id && GenMomExists){
0156             gen_mom = static_cast<const GenParticle*> (gen_mom->mother());
0157             if(gen_mom==NULL){
0158               GenMomExists=false;
0159             }
0160             if(GenMomExists) id=gen_mom->pdgId();
0161           }
0162         }
0163         if(GenMomExists) m_id = gen_mom->pdgId();
0164         gen_leptsp.push_back(gen_lept);
0165         gen_momsp.push_back(gen_mom);
0166       }
0167     }
0168   }
0169   std::vector<GenParticle> gen_lepts;
0170   for(size_t i = 0; i < gen_leptsp.size(); i++) {
0171     gen_lepts.push_back(*gen_leptsp[i]);
0172   }
0173   */
0174
0175   // Select objects based on the configuration.
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   // Fill plots for HLT muons.
0184
0185   // Find the best trigger object matches for the targetMuons.
0186   vector<size_t> matches = matchByDeltaR(targetMuons, hltMuons, plotCuts_[triggerLevel_ + "DeltaR"]);
0187
0188   // Fill plots for matched muons.
0189   for (size_t i = 0; i < targetMuons.size(); i++) {
0190     Muon& muon = targetMuons[i];
0191
0192     // Fill numerators and denominators for efficiency plots.
0193     for (const auto& suffix : EFFICIENCY_SUFFIXES) {
0194       // If no match was found, then the numerator plots don't get filled.
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     }  // finish loop numerator / denominator...
0223
0224     if (!isLastFilter_)
0225       continue;
0226     // Fill plots for tag and probe
0227     // Muon cannot be a tag because doesn't match an hlt muon
0228     if (matches[i] >= targetMuons.size())
0229       continue;
0230   }  // End loop over targetMuons.
0231
0232   // fill eff histograms for reference trigger method
0233   // Denominator: events passing reference trigger and two target muons
0234   // Numerator:   events in the denominator with two target muons
0235   // matched to hlt muons
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 }  // End analyze() method.
0282
0283 // Method to fill binning parameters from a vector of doubles.
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   // Fixed-width binning.
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   // Variable-width binning.
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 // This is an unorthodox method of getting parameters, but cleaner in my mind
0311 // It looks for parameter 'target' in the pset, and assumes that all entries
0312 // have the same class (T), filling the values into a std::map.
0313 template <class T>
0314 void HLTMuonMatchAndPlot::fillMapFromPSet(map<string, T>& m, const ParameterSet& pset, const string& target) {
0315   // Get the ParameterSet with name 'target' from 'pset'
0316   ParameterSet targetPset;
0317   if (pset.existsAs<ParameterSet>(target, true))  // target is tracked
0318     targetPset = pset.getParameterSet(target);
0319   else if (pset.existsAs<ParameterSet>(target, false))  // target is untracked
0320     targetPset = pset.getUntrackedParameterSet(target);
0321
0322   // Get the parameter names from targetPset
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))  // target is tracked
0328       m[*iter] = targetPset.getParameter<T>(*iter);
0329     else if (targetPset.existsAs<T>(*iter, false))  // target is untracked
0330       m[*iter] = targetPset.getUntrackedParameter<T>(*iter);
0331   }
0332 }
0333
0334 // A generic method to find the best deltaR matches from 2 collections.
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   // Run through the matrix n1 times to make sure we've found all matches.
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     // find the smallest deltaR
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     // If a match has been made, save it and make those candidates unavailable.
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   // If there is no selector (recoCuts does not exists), return an empty collection.
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   /* Properly delete the array of floats that has been allocated on
0427    * the heap by fillEdges.  Avoid multiple copies and internal ROOT
0428    * clones by simply creating the histograms directly in the DQMStore
0429    * using the appropriate book1D method to handle the variable bins
0430    * case. */
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   /* Properly delete the arrays of floats that have been allocated on
0450    * the heap by fillEdges.  Avoid multiple copies and internal ROOT
0451    * clones by simply creating the histograms directly in the DQMStore
0452    * using the appropriate book2D method to handle the variable bins
0453    * case. */
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 }