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