Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-07-12 22:34:35 

0001 
0002 /** \file HLTMuonPlotter.cc
0003  */
0004 
0005 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/Math/interface/deltaPhi.h"
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "DataFormats/MuonReco/interface/Muon.h"
0010 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "HLTriggerOffline/Muon/interface/HLTMuonPlotter.h"
0014 
0015 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
0016 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0017 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
0018 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeedCollection.h"
0019 
0020 using namespace std;
0021 using namespace edm;
0022 using namespace reco;
0023 using namespace trigger;
0024 
0025 namespace {
0026   const unsigned int kNull = (unsigned int)-1;
0027 }
0028 
0029 typedef vector<ParameterSet> Parameters;
0030 
0031 HLTMuonPlotter::HLTMuonPlotter(const ParameterSet &pset,
0032                                string hltPath,
0033                                const std::vector<string> &moduleLabels,
0034                                const std::vector<string> &stepLabels,
0035                                const edm::EDGetTokenT<trigger::TriggerEventWithRefs> &triggerEventToken,
0036                                const edm::EDGetTokenT<reco::GenParticleCollection> &genParticlesToken,
0037                                const edm::EDGetTokenT<reco::MuonCollection> &recoMuonsToken,
0038                                const L1MuonMatcherAlgoForDQM &l1Matcher)
0039     : l1Matcher_(l1Matcher) {
0040   hltPath_ = hltPath;
0041   moduleLabels_ = moduleLabels;
0042   stepLabels_ = stepLabels;
0043   hltProcessName_ = pset.getParameter<string>("hltProcessName");
0044 
0045   cutsDr_ = pset.getParameter<vector<double>>("cutsDr");
0046 
0047   parametersEta_ = pset.getParameter<vector<double>>("parametersEta");
0048   parametersPhi_ = pset.getParameter<vector<double>>("parametersPhi");
0049   parametersTurnOn_ = pset.getParameter<vector<double>>("parametersTurnOn");
0050 
0051   genMuonCut_ = pset.getParameter<string>("genMuonCut");
0052   recMuonCut_ = pset.getParameter<string>("recMuonCut");
0053 
0054   genMuonSelector_ = nullptr;
0055   recMuonSelector_ = nullptr;
0056 
0057   hltTriggerSummaryRAW_ = triggerEventToken;
0058   genParticleLabel_ = genParticlesToken;
0059   recMuonLabel_ = recoMuonsToken;
0060 }
0061 
0062 void HLTMuonPlotter::beginJob() {}
0063 
0064 void HLTMuonPlotter::beginRun(DQMStore::IBooker &iBooker, const Run &iRun, const EventSetup &iSetup) {
0065   l1Matcher_.init(iSetup);
0066 
0067   cutMaxEta_ = 2.4;
0068   if (hltPath_.find("eta2p1") != string::npos)
0069     cutMaxEta_ = 2.1;
0070 
0071   // Choose a pT cut for gen/rec muons based on the pT cut in the hltPath_
0072   unsigned int threshold = 0;
0073   TPRegexp ptRegexp("Mu([0-9]+)");
0074   TObjArray *regexArray = ptRegexp.MatchS(hltPath_);
0075   if (regexArray->GetEntriesFast() == 2) {
0076     threshold = atoi(((TObjString *)regexArray->At(1))->GetString());
0077   }
0078   delete regexArray;
0079   // We select a whole number min pT cut slightly above the hltPath_'s final
0080   // pt threshold, then subtract a bit to let through particle gun muons with
0081   // exact integer pT:
0082   cutMinPt_ = ceil(threshold * 1.1) - 0.01;
0083   if (cutMinPt_ < 0.)
0084     cutMinPt_ = 0.;
0085 
0086   string baseDir = "HLT/Muon/Distributions/";
0087   iBooker.setCurrentFolder(baseDir + hltPath_);
0088 
0089   vector<string> sources(2);
0090   sources[0] = "gen";
0091   sources[1] = "rec";
0092 
0093   elements_["CutMinPt"] = iBooker.bookFloat("CutMinPt");
0094   elements_["CutMaxEta"] = iBooker.bookFloat("CutMaxEta");
0095   elements_["CutMinPt"]->Fill(cutMinPt_);
0096   elements_["CutMaxEta"]->Fill(cutMaxEta_);
0097 
0098   for (size_t i = 0; i < sources.size(); i++) {
0099     string source = sources[i];
0100     for (size_t j = 0; j < stepLabels_.size(); j++) {
0101       bookHist(iBooker, hltPath_, stepLabels_[j], source, "Eta");
0102       bookHist(iBooker, hltPath_, stepLabels_[j], source, "Phi");
0103       bookHist(iBooker, hltPath_, stepLabels_[j], source, "MaxPt1");
0104       bookHist(iBooker, hltPath_, stepLabels_[j], source, "MaxPt2");
0105     }
0106   }
0107 }
0108 
0109 void HLTMuonPlotter::analyze(const Event &iEvent, const EventSetup &iSetup) {
0110   LogTrace("HLTMuonVal") << "In HLTMuonPlotter::analyze,  "
0111                          << "Event: " << iEvent.id();
0112 
0113   // cout << hltPath_ << endl;
0114   // for (size_t i = 0; i < moduleLabels_.size(); i++)
0115   //   cout << "    " << moduleLabels_[i] << endl;
0116 
0117   Handle<TriggerEventWithRefs> rawTriggerEvent;
0118   Handle<MuonCollection> recMuons;
0119   Handle<GenParticleCollection> genParticles;
0120 
0121   iEvent.getByToken(hltTriggerSummaryRAW_, rawTriggerEvent);
0122   if (rawTriggerEvent.failedToGet()) {
0123     LogError("HLTMuonVal") << "No trigger summary found";
0124     return;
0125   }
0126   iEvent.getByToken(recMuonLabel_, recMuons);
0127   iEvent.getByToken(genParticleLabel_, genParticles);
0128 
0129   vector<string> sources;
0130   if (genParticles.isValid())
0131     sources.push_back("gen");
0132   if (recMuons.isValid())
0133     sources.push_back("rec");
0134 
0135   for (size_t sourceNo = 0; sourceNo < sources.size(); sourceNo++) {
0136     string source = sources[sourceNo];
0137 
0138     // If this is the first event, initialize selectors
0139     if (!genMuonSelector_)
0140       genMuonSelector_ = new StringCutObjectSelector<reco::GenParticle>(genMuonCut_);
0141     if (!recMuonSelector_)
0142       recMuonSelector_ = new StringCutObjectSelector<reco::Muon>(recMuonCut_);
0143 
0144     // Make each good gen/rec muon into the base cand for a MatchStruct
0145     vector<MatchStruct> matches;
0146     if (source == "gen" && genParticles.isValid())
0147       for (size_t i = 0; i < genParticles->size(); i++)
0148         if ((*genMuonSelector_)(genParticles->at(i)))
0149           matches.push_back(MatchStruct(&genParticles->at(i)));
0150     if (source == "rec" && recMuons.isValid())
0151       for (size_t i = 0; i < recMuons->size(); i++)
0152         if ((*recMuonSelector_)(recMuons->at(i)))
0153           matches.push_back(MatchStruct(&recMuons->at(i)));
0154 
0155     // Sort the MatchStructs by pT for later filling of turn-on curve
0156     sort(matches.begin(), matches.end(), matchesByDescendingPt());
0157 
0158     const bool isDoubleMuonPath = (hltPath_.find("Double") != string::npos);
0159     const size_t nFilters = moduleLabels_.size();
0160     const size_t nSteps = stepLabels_.size();
0161     const size_t nStepsHlt = nSteps - 2;
0162     const int nObjectsToPassPath = (isDoubleMuonPath) ? 2 : 1;
0163     l1t::MuonVectorRef candsL1;
0164     vector<vector<RecoChargedCandidateRef>> refsHlt(nStepsHlt);
0165     vector<vector<const RecoChargedCandidate *>> candsHlt(nStepsHlt);
0166 
0167     for (size_t i = 0; i < nFilters; i++) {
0168       const int hltStep = i - 1;
0169       InputTag tag = InputTag(moduleLabels_[i], "", hltProcessName_);
0170       size_t iFilter = rawTriggerEvent->filterIndex(tag);
0171       if (iFilter < rawTriggerEvent->size()) {
0172         if (i == 0)
0173           rawTriggerEvent->getObjects(iFilter, TriggerL1Mu, candsL1);
0174         else
0175           rawTriggerEvent->getObjects(iFilter, TriggerMuon, refsHlt[hltStep]);
0176       } else
0177         LogTrace("HLTMuonVal") << "No collection with label " << tag;
0178     }
0179     for (size_t i = 0; i < nStepsHlt; i++)
0180       for (size_t j = 0; j < refsHlt[i].size(); j++)
0181         if (refsHlt[i][j].isAvailable()) {
0182           candsHlt[i].push_back(&*refsHlt[i][j]);
0183         } else {
0184           LogWarning("HLTMuonPlotter") << "Ref refsHlt[i][j]: product not available " << i << " " << j;
0185         }
0186 
0187     // Add trigger objects to the MatchStructs
0188     findMatches(matches, candsL1, candsHlt);
0189 
0190     vector<size_t> matchesInEtaRange;
0191     vector<bool> hasMatch(matches.size(), true);
0192 
0193     for (size_t step = 0; step < nSteps; step++) {
0194       size_t hltStep = (step >= 2) ? step - 2 : 0;
0195       if (nSteps == 6)
0196         hltStep = hltStep - 1;  // case of the tracker muon (it has no L2)
0197       size_t level = 0;
0198       if ((stepLabels_[step].find("L3TkIso") != string::npos) || (stepLabels_[step].find("TkTkIso") != string::npos))
0199         level = 6;
0200       else if ((stepLabels_[step].find("L3HcalIso") != string::npos) ||
0201                (stepLabels_[step].find("TkEcalIso") != string::npos))
0202         level = 5;
0203       else if ((stepLabels_[step].find("L3EcalIso") != string::npos) ||
0204                (stepLabels_[step].find("TkEcalIso") != string::npos))
0205         level = 4;
0206       else if ((stepLabels_[step].find("L3") != string::npos) || (stepLabels_[step].find("Tk") != string::npos))
0207         level = 3;
0208       else if (stepLabels_[step].find("L2") != string::npos)
0209         level = 2;
0210       else if (stepLabels_[step].find("L1") != string::npos)
0211         level = 1;
0212 
0213       for (size_t j = 0; j < matches.size(); j++) {
0214         if (level == 0) {
0215           if (fabs(matches[j].candBase->eta()) < cutMaxEta_)
0216             matchesInEtaRange.push_back(j);
0217         } else if (level == 1) {
0218           if (matches[j].candL1 == nullptr)
0219             hasMatch[j] = false;
0220         } else if (level >= 2) {
0221           if (matches[j].candHlt[hltStep] == nullptr)
0222             hasMatch[j] = false;
0223           else if (!hasMatch[j]) {
0224             LogTrace("HLTMuonVal") << "Match found for HLT step " << hltStep << " of " << nStepsHlt
0225                                    << " without previous match!";
0226             break;
0227           }
0228         }
0229       }
0230 
0231       if (std::count(hasMatch.begin(), hasMatch.end(), true) < nObjectsToPassPath)
0232         break;
0233 
0234       string pre = source + "Pass";
0235       string post = "_" + stepLabels_[step];
0236 
0237       for (size_t j = 0; j < matches.size(); j++) {
0238         float pt = matches[j].candBase->pt();
0239         float eta = matches[j].candBase->eta();
0240         float phi = matches[j].candBase->phi();
0241         if (hasMatch[j]) {
0242           if (!matchesInEtaRange.empty() && j == matchesInEtaRange[0])
0243             elements_[pre + "MaxPt1" + post]->Fill(pt);
0244           if (matchesInEtaRange.size() >= 2 && j == matchesInEtaRange[1])
0245             elements_[pre + "MaxPt2" + post]->Fill(pt);
0246           if (pt > cutMinPt_) {
0247             elements_[pre + "Eta" + post]->Fill(eta);
0248             if (fabs(eta) < cutMaxEta_)
0249               elements_[pre + "Phi" + post]->Fill(phi);
0250           }
0251         }
0252       }
0253     }
0254 
0255   }  // End loop over sources
0256 }
0257 
0258 void HLTMuonPlotter::findMatches(vector<MatchStruct> &matches,
0259                                  const l1t::MuonVectorRef &candsL1,
0260                                  const std::vector<vector<const RecoChargedCandidate *>> &candsHlt) {
0261   set<size_t>::iterator it;
0262 
0263   set<size_t> indicesL1;
0264   for (size_t i = 0; i < candsL1.size(); i++)
0265     indicesL1.insert(i);
0266 
0267   vector<set<size_t>> indicesHlt(candsHlt.size());
0268   for (size_t i = 0; i < candsHlt.size(); i++)
0269     for (size_t j = 0; j < candsHlt[i].size(); j++)
0270       indicesHlt[i].insert(j);
0271 
0272   for (size_t i = 0; i < matches.size(); i++) {
0273     const Candidate *cand = matches[i].candBase;
0274 
0275     double bestDeltaR = cutsDr_[0];
0276     size_t bestMatch = kNull;
0277     for (it = indicesL1.begin(); it != indicesL1.end(); it++) {
0278       if (candsL1[*it].isAvailable()) {
0279         double dR = deltaR(cand->eta(), cand->phi(), candsL1[*it]->eta(), candsL1[*it]->phi());
0280         if (dR < bestDeltaR) {
0281           bestMatch = *it;
0282           bestDeltaR = dR;
0283         }
0284         // TrajectoryStateOnSurface propagated;
0285         // float dR = 999., dPhi = 999.;
0286         // bool isValid = l1Matcher_.match(* cand, * candsL1[*it],
0287         //                                 dR, dPhi, propagated);
0288         // if (isValid && dR < bestDeltaR) {
0289         //   bestMatch = *it;
0290         //   bestDeltaR = dR;
0291         // }
0292       } else {
0293         LogWarning("HLTMuonPlotter") << "Ref candsL1[*it]: product not available " << *it;
0294       }
0295     }
0296 
0297     if (bestMatch != kNull)
0298       matches[i].candL1 = &*candsL1[bestMatch];
0299     indicesL1.erase(bestMatch);
0300 
0301     matches[i].candHlt.assign(candsHlt.size(), nullptr);
0302     for (size_t j = 0; j < candsHlt.size(); j++) {
0303       size_t level = (candsHlt.size() == 4) ? (j < 2) ? 2 : 3 : (candsHlt.size() == 2) ? (j < 1) ? 2 : 3 : 2;
0304       bestDeltaR = cutsDr_[level - 2];
0305       bestMatch = kNull;
0306       for (it = indicesHlt[j].begin(); it != indicesHlt[j].end(); it++) {
0307         double dR = deltaR(cand->eta(), cand->phi(), candsHlt[j][*it]->eta(), candsHlt[j][*it]->phi());
0308         if (dR < bestDeltaR) {
0309           bestMatch = *it;
0310           bestDeltaR = dR;
0311         }
0312       }
0313       if (bestMatch != kNull)
0314         matches[i].candHlt[j] = candsHlt[j][bestMatch];
0315       indicesHlt[j].erase(bestMatch);
0316     }
0317 
0318     //     cout << "    Muon: " << cand->eta() << ", ";
0319     //     if (matches[i].candL1) cout << matches[i].candL1->eta() << ", ";
0320     //     else cout << "none, ";
0321     //     for (size_t j = 0; j < candsHlt.size(); j++)
0322     //       if (matches[i].candHlt[j]) cout << matches[i].candHlt[j]->eta() <<
0323     //       ", "; else cout << "none, ";
0324     //     cout << endl;
0325   }
0326 }
0327 
0328 void HLTMuonPlotter::bookHist(DQMStore::IBooker &iBooker, string path, string label, string source, string type) {
0329   string sourceUpper = source;
0330   sourceUpper[0] = toupper(sourceUpper[0]);
0331   string name = source + "Pass" + type + "_" + label;
0332   TH1F *h;
0333 
0334   if (type.find("MaxPt") != string::npos) {
0335     string desc = (type == "MaxPt1") ? "Leading" : "Next-to-Leading";
0336     string title = "pT of " + desc + " " + sourceUpper + " Muon " + "matched to " + label;
0337     const size_t nBins = parametersTurnOn_.size() - 1;
0338     float *edges = new float[nBins + 1];
0339     for (size_t i = 0; i < nBins + 1; i++)
0340       edges[i] = parametersTurnOn_[i];
0341     h = new TH1F(name.c_str(), title.c_str(), nBins, edges);
0342   }
0343 
0344   else {
0345     string symbol = (type == "Eta") ? "#eta" : "#phi";
0346     string title = symbol + " of " + sourceUpper + " Muons " + "matched to " + label;
0347     vector<double> params = (type == "Eta") ? parametersEta_ : parametersPhi_;
0348     int nBins = (int)params[0];
0349     double min = params[1];
0350     double max = params[2];
0351     h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
0352   }
0353 
0354   h->Sumw2();
0355   elements_[name] = iBooker.book1D(name, h);
0356   delete h;
0357 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!