Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:18:53 

0001 /**\class HeavyFlavorValidation HeavyFlavorValidation.cc HLTriggerOfflineHeavyFlavor/src/HeavyFlavorValidation.cc
0002 
0003  Description: Analyzer to fill Monitoring Elements for muon, dimuon and trigger path efficiency studies (HLT/RECO, RECO/GEN)
0004 
0005  Implementation:
0006      matching is based on closest in delta R, no duplicates allowed. Generated to Global based on momentum at IP; L1, L2, L2v to Global based on position in muon system, L3 to Global based on momentum at IP.
0007 */
0008 // Original Author:  Zoltan Gecse
0009 
0010 #include <memory>
0011 #include <initializer_list>
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "FWCore/Common/interface/TriggerNames.h"
0019 
0020 #include "DataFormats/Common/interface/RefToBase.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0023 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0024 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0025 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0026 #include "DataFormats/Candidate/interface/Candidate.h"
0027 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0028 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0029 #include "DataFormats/MuonReco/interface/Muon.h"
0030 #include "DataFormats/Math/interface/deltaR.h"
0031 #include "DataFormats/Common/interface/TriggerResults.h"
0032 #include "DataFormats/Candidate/interface/Particle.h"
0033 #include "DataFormats/HLTReco/interface/TriggerEventWithRefs.h"
0034 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0035 
0036 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0037 
0038 #include "DQMServices/Core/interface/DQMStore.h"
0039 #include <DQMServices/Core/interface/DQMEDAnalyzer.h>
0040 
0041 #include "CommonTools/Utils/interface/PtComparator.h"
0042 
0043 #include "TLorentzVector.h"
0044 
0045 using namespace std;
0046 using namespace edm;
0047 using namespace reco;
0048 using namespace l1extra;
0049 using namespace trigger;
0050 
0051 class HeavyFlavorValidation : public DQMEDAnalyzer {
0052 public:
0053   explicit HeavyFlavorValidation(const edm::ParameterSet &);
0054   ~HeavyFlavorValidation() override;
0055 
0056 protected:
0057   void dqmBeginRun(const edm::Run &, const edm::EventSetup &) override;
0058   void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
0059   void analyze(const edm::Event &, const edm::EventSetup &) override;
0060 
0061 private:
0062   int getMotherId(const Candidate *p);
0063   void match(MonitorElement *me,
0064              vector<LeafCandidate> &from,
0065              vector<LeafCandidate> &to,
0066              double deltaRMatchingCut,
0067              vector<int> &map);
0068   void myBook2D(DQMStore::IBooker &ibooker,
0069                 TString name,
0070                 vector<double> &xBins,
0071                 TString xLabel,
0072                 vector<double> &yBins,
0073                 TString yLabel,
0074                 TString title);
0075   void myBook2D(DQMStore::IBooker &ibooker,
0076                 TString name,
0077                 vector<double> &xBins,
0078                 TString xLabel,
0079                 vector<double> &yBins,
0080                 TString yLabel) {
0081     myBook2D(ibooker, name, xBins, xLabel, yBins, yLabel, name);
0082   }
0083   void myBookProfile2D(DQMStore::IBooker &ibooker,
0084                        TString name,
0085                        vector<double> &xBins,
0086                        TString xLabel,
0087                        vector<double> &yBins,
0088                        TString yLabel,
0089                        TString title);
0090   void myBookProfile2D(DQMStore::IBooker &ibooker,
0091                        TString name,
0092                        vector<double> &xBins,
0093                        TString xLabel,
0094                        vector<double> &yBins,
0095                        TString yLabel) {
0096     myBookProfile2D(ibooker, name, xBins, xLabel, yBins, yLabel, name);
0097   }
0098   void myBook1D(DQMStore::IBooker &ibooker, TString name, vector<double> &xBins, TString label, TString title);
0099   void myBook1D(DQMStore::IBooker &ibooker, TString name, vector<double> &xBins, TString label) {
0100     myBook1D(ibooker, name, xBins, label, name);
0101   }
0102 
0103   /**
0104      * Get the filter "level" (as it is defined for the use of this module and its corresponding
0105      * harvesting module).
0106      *
0107      * level 1 - 3 -> more or less synonymously to the the trigger levels
0108      * level 4 and 5 -> vertex, dz, track, etc.. filters
0109      *
0110      * See the comments in the definition for some more details.
0111      */
0112   int getFilterLevel(const std::string &moduleName, const HLTConfigProvider &hltConfig);
0113 
0114   string dqmFolder;
0115   string triggerProcessName;
0116   string triggerPathName;
0117 
0118   EDGetTokenT<TriggerEventWithRefs> triggerSummaryRAWTag;
0119   EDGetTokenT<TriggerEvent> triggerSummaryAODTag;
0120   InputTag triggerResultsTag;
0121   EDGetTokenT<TriggerResults> triggerResultsToken;
0122   InputTag recoMuonsTag;
0123   EDGetTokenT<MuonCollection> recoMuonsToken;
0124   InputTag genParticlesTag;
0125   EDGetTokenT<GenParticleCollection> genParticlesToken;
0126 
0127   vector<int> motherIDs;
0128   double genGlobDeltaRMatchingCut;
0129   double globL1DeltaRMatchingCut;
0130   double globL2DeltaRMatchingCut;
0131   double globL3DeltaRMatchingCut;
0132   vector<double> deltaEtaBins;
0133   vector<double> deltaPhiBins;
0134   vector<double> muonPtBins;
0135   vector<double> muonEtaBins;
0136   vector<double> muonPhiBins;
0137   vector<double> dimuonPtBins;
0138   vector<double> dimuonEtaBins;
0139   vector<double> dimuonDRBins;
0140   map<TString, MonitorElement *> ME;
0141   vector<pair<string, int> > filterNamesLevels;
0142   const double muonMass;
0143 };
0144 
0145 HeavyFlavorValidation::HeavyFlavorValidation(const ParameterSet &pset)
0146     :  //get parameters
0147       dqmFolder(pset.getUntrackedParameter<string>("DQMFolder")),
0148       triggerProcessName(pset.getUntrackedParameter<string>("TriggerProcessName")),
0149       triggerPathName(pset.getUntrackedParameter<string>("TriggerPathName")),
0150       motherIDs(pset.getUntrackedParameter<vector<int> >("MotherIDs")),
0151       genGlobDeltaRMatchingCut(pset.getUntrackedParameter<double>("GenGlobDeltaRMatchingCut")),
0152       globL1DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL1DeltaRMatchingCut")),
0153       globL2DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL2DeltaRMatchingCut")),
0154       globL3DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL3DeltaRMatchingCut")),
0155       deltaEtaBins(pset.getUntrackedParameter<vector<double> >("DeltaEtaBins")),
0156       deltaPhiBins(pset.getUntrackedParameter<vector<double> >("DeltaPhiBins")),
0157       muonPtBins(pset.getUntrackedParameter<vector<double> >("MuonPtBins")),
0158       muonEtaBins(pset.getUntrackedParameter<vector<double> >("MuonEtaBins")),
0159       muonPhiBins(pset.getUntrackedParameter<vector<double> >("MuonPhiBins")),
0160       dimuonPtBins(pset.getUntrackedParameter<vector<double> >("DimuonPtBins")),
0161       dimuonEtaBins(pset.getUntrackedParameter<vector<double> >("DimuonEtaBins")),
0162       dimuonDRBins(pset.getUntrackedParameter<vector<double> >("DimuonDRBins")),
0163       muonMass(0.106) {
0164   triggerSummaryRAWTag = consumes<TriggerEventWithRefs>(
0165       InputTag(pset.getUntrackedParameter<string>("TriggerSummaryRAW"), "", triggerProcessName));
0166   triggerSummaryAODTag =
0167       consumes<TriggerEvent>(InputTag(pset.getUntrackedParameter<string>("TriggerSummaryAOD"), "", triggerProcessName));
0168   triggerResultsTag = InputTag(pset.getUntrackedParameter<string>("TriggerResults"), "", triggerProcessName);
0169   triggerResultsToken = consumes<TriggerResults>(triggerResultsTag);
0170   recoMuonsTag = pset.getParameter<InputTag>("RecoMuons");
0171   recoMuonsToken = consumes<MuonCollection>(recoMuonsTag);
0172   genParticlesTag = pset.getParameter<InputTag>("GenParticles");
0173   genParticlesToken = consumes<GenParticleCollection>(genParticlesTag);
0174 }
0175 
0176 void HeavyFlavorValidation::dqmBeginRun(const edm::Run &iRun, const edm::EventSetup &iSetup) {
0177   //discover HLT configuration
0178   HLTConfigProvider hltConfig;
0179   bool isChanged;
0180   if (hltConfig.init(iRun, iSetup, triggerProcessName, isChanged)) {
0181     LogDebug("HLTriggerOfflineHeavyFlavor")
0182         << "Successfully initialized HLTConfigProvider with process name: " << triggerProcessName << endl;
0183   } else {
0184     LogWarning("HLTriggerOfflineHeavyFlavor")
0185         << "Could not initialize HLTConfigProvider with process name: " << triggerProcessName << endl;
0186     return;
0187   }
0188   vector<string> triggerNames = hltConfig.triggerNames();
0189   for (const auto &trigName : triggerNames) {
0190     // TString triggerName = trigName;
0191     if (trigName.find(triggerPathName) != std::string::npos) {
0192       vector<string> moduleNames = hltConfig.moduleLabels(trigName);
0193       for (const auto &moduleName : moduleNames) {
0194         const int level = getFilterLevel(moduleName, hltConfig);
0195         if (level > 0) {
0196           filterNamesLevels.push_back({moduleName, level});
0197         }
0198       }
0199       break;
0200     }
0201   }
0202 
0203   if (filterNamesLevels.empty()) {
0204     LogDebug("HLTriggerOfflineHeavyFlavor") << "Bad Trigger Path: " << triggerPathName << endl;
0205     return;
0206   } else {
0207     std::string str;
0208     str.reserve(
0209         512);  // avoid too many realloctions in the following loop (allows for filter names with roughly 100 chars each)
0210     for (const auto &filters : filterNamesLevels)
0211       str = str + " " + filters.first;
0212     LogDebug("HLTriggerOfflineHeavyFlavor") << "Trigger Path: " << triggerPathName << " has filters:" << str;
0213   }
0214 }
0215 
0216 void HeavyFlavorValidation::bookHistograms(DQMStore::IBooker &ibooker,
0217                                            edm::Run const &iRun,
0218                                            edm::EventSetup const &iSetup) {
0219   ibooker.cd();
0220   ibooker.setCurrentFolder((dqmFolder + "/") + triggerProcessName + "/" + triggerPathName);
0221 
0222   // create Monitor Elements
0223   // Eta Pt Single
0224   myBook2D(ibooker, "genMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
0225   myBook2D(ibooker, "globMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
0226   myBook2D(ibooker, "globMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
0227   myBook2D(ibooker, "pathMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
0228   myBook2D(ibooker, "resultMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
0229   // Eta Pt Single Resolution
0230   myBookProfile2D(ibooker, "resGlobGen_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
0231   myBookProfile2D(
0232       ibooker, "resPathGlob_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
0233   // Eta Pt Double
0234   myBook2D(ibooker, "genDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
0235   myBook2D(ibooker, "globDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
0236   myBook2D(ibooker, "globDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
0237   myBook2D(
0238       ibooker, "pathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
0239   myBook2D(ibooker, "resultDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
0240   myBook2D(
0241       ibooker, "diPathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
0242   // Eta Phi Single
0243   myBook2D(ibooker, "genMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
0244   myBook2D(ibooker, "globMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
0245   myBook2D(ibooker, "globMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
0246   myBook2D(ibooker, "pathMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", triggerPathName);
0247   myBook2D(ibooker, "resultMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
0248   // Rap Pt Double
0249   myBook2D(ibooker, "genDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
0250   myBook2D(ibooker, "globDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
0251   myBook2D(ibooker, "globDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
0252   myBook2D(ibooker,
0253            "pathDimuon_recoRapPt",
0254            dimuonEtaBins,
0255            "#mu#mu rapidity",
0256            dimuonPtBins,
0257            " #mu#mu pT (GeV)",
0258            triggerPathName);
0259   myBook2D(ibooker, "resultDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
0260   myBook2D(ibooker,
0261            "diPathDimuon_recoRapPt",
0262            dimuonEtaBins,
0263            "#mu#mu rapidity",
0264            dimuonPtBins,
0265            " #mu#mu pT (GeV)",
0266            triggerPathName);
0267   // Pt DR Double
0268   myBook2D(ibooker, "genDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
0269   myBook2D(ibooker, "globDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
0270   myBook2D(ibooker, "globDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
0271   myBook2D(ibooker,
0272            "pathDimuon_recoPtDR",
0273            dimuonPtBins,
0274            " #mu#mu pT (GeV)",
0275            dimuonDRBins,
0276            "#mu#mu #Delta R at IP",
0277            triggerPathName);
0278   myBook2D(ibooker,
0279            "diPathDimuon_recoPtDR",
0280            dimuonPtBins,
0281            " #mu#mu pT (GeV)",
0282            dimuonDRBins,
0283            "#mu#mu #Delta R at IP",
0284            triggerPathName);
0285   myBook2D(ibooker, "resultDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
0286   // Pt DRpos Double
0287   myBook2D(ibooker, "globDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
0288   myBook2D(ibooker,
0289            "pathDimuon_recoPtDRpos",
0290            dimuonPtBins,
0291            " #mu#mu pT (GeV)",
0292            dimuonDRBins,
0293            "#mu#mu #Delta R in MS",
0294            triggerPathName);
0295   myBook2D(ibooker,
0296            "diPathDimuon_recoPtDRpos",
0297            dimuonPtBins,
0298            " #mu#mu pT (GeV)",
0299            dimuonDRBins,
0300            "#mu#mu #Delta R in MS",
0301            triggerPathName);
0302   myBook2D(
0303       ibooker, "resultDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
0304 
0305   // Matching
0306   myBook2D(ibooker, "globGen_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi");
0307   myBook2D(
0308       ibooker, "pathGlob_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", triggerPathName);
0309   // Size of containers
0310   vector<double> sizeBins;
0311   sizeBins.push_back(10);
0312   sizeBins.push_back(0);
0313   sizeBins.push_back(10);
0314   myBook1D(ibooker, "genMuon_size", sizeBins, "container size");
0315   myBook1D(ibooker, "globMuon_size", sizeBins, "container size");
0316   myBook1D(ibooker, "pathMuon_size", sizeBins, "container size", triggerPathName);
0317 }
0318 
0319 void HeavyFlavorValidation::analyze(const Event &iEvent, const EventSetup &iSetup) {
0320   if (filterNamesLevels.empty()) {
0321     return;
0322   }
0323   //access the containers and create LeafCandidate copies
0324   vector<LeafCandidate> genMuons;
0325   Handle<GenParticleCollection> genParticles;
0326   iEvent.getByToken(genParticlesToken, genParticles);
0327   if (genParticles.isValid()) {
0328     for (GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++p) {
0329       if (p->status() == 1 && std::abs(p->pdgId()) == 13 &&
0330           (find(motherIDs.begin(), motherIDs.end(), -1) != motherIDs.end() ||
0331            find(motherIDs.begin(), motherIDs.end(), getMotherId(&(*p))) != motherIDs.end())) {
0332         genMuons.push_back(*p);
0333       }
0334     }
0335   } else {
0336     LogDebug("HLTriggerOfflineHeavyFlavor") << "Could not access GenParticleCollection" << endl;
0337   }
0338   sort(genMuons.begin(), genMuons.end(), GreaterByPt<LeafCandidate>());
0339   ME["genMuon_size"]->Fill(genMuons.size());
0340   LogDebug("HLTriggerOfflineHeavyFlavor")
0341       << "GenParticleCollection from " << genParticlesTag << " has size: " << genMuons.size() << endl;
0342 
0343   vector<LeafCandidate> globMuons;
0344   vector<LeafCandidate> globMuons_position;
0345   Handle<MuonCollection> recoMuonsHandle;
0346   iEvent.getByToken(recoMuonsToken, recoMuonsHandle);
0347   if (recoMuonsHandle.isValid()) {
0348     for (MuonCollection::const_iterator p = recoMuonsHandle->begin(); p != recoMuonsHandle->end(); ++p) {
0349       if (p->isGlobalMuon()) {
0350         globMuons.push_back(*p);
0351         globMuons_position.push_back(LeafCandidate(p->charge(),
0352                                                    math::XYZTLorentzVector(p->outerTrack()->innerPosition().x(),
0353                                                                            p->outerTrack()->innerPosition().y(),
0354                                                                            p->outerTrack()->innerPosition().z(),
0355                                                                            0.)));
0356       }
0357     }
0358   } else {
0359     LogDebug("HLTriggerOfflineHeavyFlavor") << "Could not access reco Muons" << endl;
0360   }
0361   ME["globMuon_size"]->Fill(globMuons.size());
0362   LogDebug("HLTriggerOfflineHeavyFlavor")
0363       << "Global Muons from " << recoMuonsTag << " has size: " << globMuons.size() << endl;
0364 
0365   // access RAW trigger event
0366   vector<vector<LeafCandidate> > muonsAtFilter;
0367   vector<vector<LeafCandidate> > muonPositionsAtFilter;
0368   for (size_t i = 0; i < filterNamesLevels.size(); i++) {
0369     muonsAtFilter.push_back(vector<LeafCandidate>());
0370     muonPositionsAtFilter.push_back(vector<LeafCandidate>());
0371   }
0372   Handle<TriggerEventWithRefs> rawTriggerEvent;
0373   iEvent.getByToken(triggerSummaryRAWTag, rawTriggerEvent);
0374   if (rawTriggerEvent.isValid()) {
0375     for (size_t i = 0; i < filterNamesLevels.size(); i++) {
0376       size_t index = rawTriggerEvent->filterIndex(InputTag(filterNamesLevels[i].first, "", triggerProcessName));
0377       if (index < rawTriggerEvent->size()) {
0378         if (filterNamesLevels[i].second == 1) {
0379           vector<L1MuonParticleRef> l1Cands;
0380           rawTriggerEvent->getObjects(index, TriggerL1Mu, l1Cands);
0381           for (size_t j = 0; j < l1Cands.size(); j++) {
0382             muonsAtFilter[i].push_back(*l1Cands[j]);
0383           }
0384         } else {
0385           vector<RecoChargedCandidateRef> hltCands;
0386           rawTriggerEvent->getObjects(index, TriggerMuon, hltCands);
0387           for (size_t j = 0; j < hltCands.size(); j++) {
0388             muonsAtFilter[i].push_back(*hltCands[j]);
0389             if (filterNamesLevels[i].second == 2) {
0390               muonPositionsAtFilter[i].push_back(
0391                   LeafCandidate(hltCands[j]->charge(),
0392                                 math::XYZTLorentzVector(hltCands[j]->track()->innerPosition().x(),
0393                                                         hltCands[j]->track()->innerPosition().y(),
0394                                                         hltCands[j]->track()->innerPosition().z(),
0395                                                         0.)));
0396             }
0397           }
0398         }
0399       }
0400       //ME[TString::Format("filt%dMuon_size", int(i + 1))]->Fill(muonsAtFilter[i].size());
0401       LogDebug("HLTriggerOfflineHeavyFlavor")
0402           << "Filter \"" << filterNamesLevels[i].first << "\" has " << muonsAtFilter[i].size() << " muons" << endl;
0403     }
0404   } else {
0405     LogDebug("HLTriggerOfflineHeavyFlavor") << "Could not access RAWTriggerEvent" << endl;
0406   }
0407 
0408   // access AOD trigger event
0409   vector<LeafCandidate> pathMuons;
0410   Handle<TriggerEvent> aodTriggerEvent;
0411   iEvent.getByToken(triggerSummaryAODTag, aodTriggerEvent);
0412   if (aodTriggerEvent.isValid()) {
0413     TriggerObjectCollection allObjects = aodTriggerEvent->getObjects();
0414     for (int i = 0; i < aodTriggerEvent->sizeFilters(); i++) {
0415       if (aodTriggerEvent->filterTag(i) == InputTag((filterNamesLevels.end() - 1)->first, "", triggerProcessName)) {
0416         Keys keys = aodTriggerEvent->filterKeys(i);
0417         for (size_t j = 0; j < keys.size(); j++) {
0418           pathMuons.push_back(LeafCandidate(
0419               allObjects[keys[j]].id() > 0 ? 1 : -1,
0420               math::PtEtaPhiMLorentzVector(
0421                   allObjects[keys[j]].pt(), allObjects[keys[j]].eta(), allObjects[keys[j]].phi(), muonMass)));
0422         }
0423       }
0424     }
0425     ME["pathMuon_size"]->Fill(pathMuons.size());
0426     LogDebug("HLTriggerOfflineHeavyFlavor")
0427         << "Path \"" << triggerPathName << "\" has " << pathMuons.size() << " muons at last filter \""
0428         << (filterNamesLevels.end() - 1)->first << "\"" << endl;
0429   } else {
0430     LogDebug("HLTriggerOfflineHeavyFlavor") << "Could not access AODTriggerEvent" << endl;
0431   }
0432 
0433   // access Trigger Results
0434   bool triggerFired = false;
0435   Handle<TriggerResults> triggerResults;
0436   iEvent.getByToken(triggerResultsToken, triggerResults);
0437   if (triggerResults.isValid()) {
0438     LogDebug("HLTriggerOfflineHeavyFlavor") << "Successfully initialized " << triggerResultsTag << endl;
0439     const edm::TriggerNames &triggerNames = iEvent.triggerNames(*triggerResults);
0440     bool hlt_exists = false;
0441     for (unsigned int i = 0; i != triggerNames.size(); i++) {
0442       TString hlt_name = triggerNames.triggerName(i);
0443       if (hlt_name.Contains(triggerPathName)) {
0444         triggerFired = triggerResults->accept(i);
0445         hlt_exists = true;
0446         break;
0447       }
0448     }
0449     if (!hlt_exists) {
0450       LogDebug("HLTriggerOfflineHeavyFlavor") << triggerResultsTag << " has no trigger: " << triggerPathName << endl;
0451     }
0452   } else {
0453     LogDebug("HLTriggerOfflineHeavyFlavor") << "Could not initialize " << triggerResultsTag << endl;
0454   }
0455 
0456   //create matching maps
0457   vector<int> glob_gen(genMuons.size(), -1);
0458   match(ME["globGen_deltaEtaDeltaPhi"], genMuons, globMuons, genGlobDeltaRMatchingCut, glob_gen);
0459   vector<vector<int> > filt_glob;
0460   vector<int> path_glob(globMuons.size(), -1);
0461   if ((filterNamesLevels.end() - 1)->second == 1) {
0462     match(ME["pathGlob_deltaEtaDeltaPhi"], globMuons_position, pathMuons, globL1DeltaRMatchingCut, path_glob);
0463   } else if ((filterNamesLevels.end() - 1)->second == 2) {
0464     match(ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons, globL2DeltaRMatchingCut, path_glob);
0465   } else if ((filterNamesLevels.end() - 1)->second > 2) {
0466     match(ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons, globL3DeltaRMatchingCut, path_glob);
0467   }
0468 
0469   //fill histos
0470   bool first = true;
0471   for (size_t i = 0; i < genMuons.size(); i++) {
0472     ME["genMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
0473     ME["genMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
0474     if (glob_gen[i] != -1) {
0475       ME["resGlobGen_genEtaPt"]->Fill(
0476           genMuons[i].eta(), genMuons[i].pt(), (globMuons[glob_gen[i]].pt() - genMuons[i].pt()) / genMuons[i].pt());
0477       ME["globMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
0478       ME["globMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
0479       ME["globMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
0480       ME["globMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
0481       if (path_glob[glob_gen[i]] != -1) {
0482         ME["resPathGlob_recoEtaPt"]->Fill(
0483             globMuons[glob_gen[i]].eta(),
0484             globMuons[glob_gen[i]].pt(),
0485             (pathMuons[path_glob[glob_gen[i]]].pt() - globMuons[glob_gen[i]].pt()) / globMuons[glob_gen[i]].pt());
0486         ME["pathMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
0487         ME["pathMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
0488       }
0489       //highest pt muon
0490       if (first) {
0491         first = false;
0492         if (triggerFired) {
0493           ME["resultMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
0494           ME["resultMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
0495         }
0496       }
0497     }
0498   }
0499 
0500   //fill dimuon histograms (highest pT, opposite charge)
0501   int secondMuon = 0;
0502   for (size_t j = 1; j < genMuons.size(); j++) {
0503     if (genMuons[0].charge() * genMuons[j].charge() == -1) {
0504       secondMuon = j;
0505       break;
0506     }
0507   }
0508   if (secondMuon > 0) {
0509     //two generated
0510     double genDimuonPt = (genMuons[0].p4() + genMuons[secondMuon].p4()).pt();
0511     double genDimuonEta = (genMuons[0].p4() + genMuons[secondMuon].p4()).eta();
0512     double genDimuonRap = (genMuons[0].p4() + genMuons[secondMuon].p4()).Rapidity();
0513     double genDimuonDR = deltaR<LeafCandidate, LeafCandidate>(genMuons[0], genMuons[secondMuon]);
0514     bool highPt = genMuons[0].pt() > 7. && genMuons[secondMuon].pt() > 7;
0515     ME["genDimuon_genEtaPt"]->Fill(genDimuonEta, genDimuonPt);
0516     ME["genDimuon_genRapPt"]->Fill(genDimuonRap, genDimuonPt);
0517     if (highPt)
0518       ME["genDimuon_genPtDR"]->Fill(genDimuonPt, genDimuonDR);
0519     //two global
0520     if (glob_gen[0] != -1 && glob_gen[secondMuon] != -1) {
0521       ME["globDimuon_genEtaPt"]->Fill(genDimuonEta, genDimuonPt);
0522       ME["globDimuon_genRapPt"]->Fill(genDimuonRap, genDimuonPt);
0523       if (highPt)
0524         ME["globDimuon_genPtDR"]->Fill(genDimuonPt, genDimuonDR);
0525       double globDimuonPt = (globMuons[glob_gen[0]].p4() + globMuons[glob_gen[secondMuon]].p4()).pt();
0526       double globDimuonEta = (globMuons[glob_gen[0]].p4() + globMuons[glob_gen[secondMuon]].p4()).eta();
0527       double globDimuonRap = (globMuons[glob_gen[0]].p4() + globMuons[glob_gen[secondMuon]].p4()).Rapidity();
0528       double globDimuonDR =
0529           deltaR<LeafCandidate, LeafCandidate>(globMuons[glob_gen[0]], globMuons[glob_gen[secondMuon]]);
0530       double globDimuonDRpos = deltaR<LeafCandidate, LeafCandidate>(globMuons_position[glob_gen[0]],
0531                                                                     globMuons_position[glob_gen[secondMuon]]);
0532       ME["globDimuon_recoEtaPt"]->Fill(globDimuonEta, globDimuonPt);
0533       ME["globDimuon_recoRapPt"]->Fill(globDimuonRap, globDimuonPt);
0534       if (highPt)
0535         ME["globDimuon_recoPtDR"]->Fill(globDimuonPt, globDimuonDR);
0536       if (highPt)
0537         ME["globDimuon_recoPtDRpos"]->Fill(globDimuonPt, globDimuonDRpos);
0538       if (path_glob[glob_gen[0]] != -1 && path_glob[glob_gen[secondMuon]] != -1) {
0539         ME["diPathDimuon_recoEtaPt"]->Fill(globDimuonEta, globDimuonPt);
0540         ME["diPathDimuon_recoRapPt"]->Fill(globDimuonRap, globDimuonPt);
0541         if (highPt)
0542           ME["diPathDimuon_recoPtDR"]->Fill(globDimuonPt, globDimuonDR);
0543         if (highPt)
0544           ME["diPathDimuon_recoPtDRpos"]->Fill(globDimuonPt, globDimuonDRpos);
0545       }
0546       //one path object
0547       if (path_glob[glob_gen[0]] != -1 || path_glob[glob_gen[secondMuon]] != -1) {
0548         ME["pathDimuon_recoEtaPt"]->Fill(globDimuonEta, globDimuonPt);
0549         ME["pathDimuon_recoRapPt"]->Fill(globDimuonRap, globDimuonPt);
0550         if (highPt)
0551           ME["pathDimuon_recoPtDR"]->Fill(globDimuonPt, globDimuonDR);
0552         if (highPt)
0553           ME["pathDimuon_recoPtDRpos"]->Fill(globDimuonPt, globDimuonDRpos);
0554       }
0555       //trigger result
0556       if (triggerFired) {
0557         ME["resultDimuon_recoEtaPt"]->Fill(globDimuonEta, globDimuonPt);
0558         ME["resultDimuon_recoRapPt"]->Fill(globDimuonRap, globDimuonPt);
0559         if (highPt)
0560           ME["resultDimuon_recoPtDR"]->Fill(globDimuonPt, globDimuonDR);
0561         if (highPt)
0562           ME["resultDimuon_recoPtDRpos"]->Fill(globDimuonPt, globDimuonDRpos);
0563       }
0564     }
0565   }
0566 }
0567 
0568 int HeavyFlavorValidation::getMotherId(const Candidate *p) {
0569   const Candidate *mother = p->mother();
0570   if (mother) {
0571     if (mother->pdgId() == p->pdgId()) {
0572       return getMotherId(mother);
0573     } else {
0574       return mother->pdgId();
0575     }
0576   } else {
0577     return 0;
0578   }
0579 }
0580 
0581 void HeavyFlavorValidation::match(MonitorElement *me,
0582                                   vector<LeafCandidate> &from,
0583                                   vector<LeafCandidate> &to,
0584                                   double dRMatchingCut,
0585                                   vector<int> &map) {
0586   vector<double> dR(from.size());
0587   for (size_t i = 0; i < from.size(); i++) {
0588     map[i] = -1;
0589     dR[i] = 10.;
0590     //find closest
0591     for (size_t j = 0; j < to.size(); j++) {
0592       double dRtmp = deltaR<double>(from[i].eta(), from[i].phi(), to[j].eta(), to[j].phi());
0593       if (dRtmp < dR[i]) {
0594         dR[i] = dRtmp;
0595         map[i] = j;
0596       }
0597     }
0598     //fill matching histo
0599     if (map[i] != -1) {
0600       me->Fill(to[map[i]].eta() - from[i].eta(), deltaPhi<double>(to[map[i]].phi(), from[i].phi()));
0601     }
0602     //apply matching cut
0603     if (dR[i] > dRMatchingCut) {
0604       map[i] = -1;
0605     }
0606     //remove duplication
0607     if (map[i] != -1) {
0608       for (size_t k = 0; k < i; k++) {
0609         if (map[k] != -1 && map[i] == map[k]) {
0610           if (dR[i] < dR[k]) {
0611             map[k] = -1;
0612           } else {
0613             map[i] = -1;
0614           }
0615           break;
0616         }
0617       }
0618     }
0619   }
0620 }
0621 
0622 void HeavyFlavorValidation::myBook2D(DQMStore::IBooker &ibooker,
0623                                      TString name,
0624                                      vector<double> &ptBins,
0625                                      TString ptLabel,
0626                                      vector<double> &etaBins,
0627                                      TString etaLabel,
0628                                      TString title) {
0629   //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
0630   int ptN = ptBins.size() == 3 ? (int)ptBins[0] + 1 : ptBins.size();
0631   Double_t *pt = new Double_t[ptN];
0632   for (int i = 0; i < ptN; i++) {
0633     pt[i] = ptBins.size() == 3 ? ptBins[1] + i * (ptBins[2] - ptBins[1]) / ptBins[0] : ptBins[i];
0634   }
0635   int etaN = etaBins.size() == 3 ? (int)etaBins[0] + 1 : etaBins.size();
0636   Double_t *eta = new Double_t[etaN];
0637   for (int i = 0; i < etaN; i++) {
0638     eta[i] = etaBins.size() == 3 ? etaBins[1] + i * (etaBins[2] - etaBins[1]) / etaBins[0] : etaBins[i];
0639   }
0640   TH2F *h = new TH2F(name, name, ptN - 1, pt, etaN - 1, eta);
0641   h->SetXTitle(ptLabel);
0642   h->SetYTitle(etaLabel);
0643   h->SetTitle(title);
0644   ME[name] = ibooker.book2D(name.Data(), h);
0645   delete h;
0646 }
0647 
0648 void HeavyFlavorValidation::myBookProfile2D(DQMStore::IBooker &ibooker,
0649                                             TString name,
0650                                             vector<double> &ptBins,
0651                                             TString ptLabel,
0652                                             vector<double> &etaBins,
0653                                             TString etaLabel,
0654                                             TString title) {
0655   //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
0656   int ptN = ptBins.size() == 3 ? (int)ptBins[0] + 1 : ptBins.size();
0657   Double_t *pt = new Double_t[ptN];
0658   for (int i = 0; i < ptN; i++) {
0659     pt[i] = ptBins.size() == 3 ? ptBins[1] + i * (ptBins[2] - ptBins[1]) / ptBins[0] : ptBins[i];
0660   }
0661   int etaN = etaBins.size() == 3 ? (int)etaBins[0] + 1 : etaBins.size();
0662   Double_t *eta = new Double_t[etaN];
0663   for (int i = 0; i < etaN; i++) {
0664     eta[i] = etaBins.size() == 3 ? etaBins[1] + i * (etaBins[2] - etaBins[1]) / etaBins[0] : etaBins[i];
0665   }
0666   TProfile2D *h = new TProfile2D(name, name, ptN - 1, pt, etaN - 1, eta);
0667   h->SetXTitle(ptLabel);
0668   h->SetYTitle(etaLabel);
0669   h->SetTitle(title);
0670   ME[name] = ibooker.bookProfile2D(name.Data(), h);
0671   delete h;
0672 }
0673 
0674 void HeavyFlavorValidation::myBook1D(
0675     DQMStore::IBooker &ibooker, TString name, vector<double> &bins, TString label, TString title) {
0676   //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
0677   int binsN = bins.size() == 3 ? (int)bins[0] + 1 : bins.size();
0678   Double_t *myBins = new Double_t[binsN];
0679   for (int i = 0; i < binsN; i++) {
0680     myBins[i] = bins.size() == 3 ? bins[1] + i * (bins[2] - bins[1]) / bins[0] : bins[i];
0681   }
0682   TH1F *h = new TH1F(name, name, binsN - 1, myBins);
0683   h->SetXTitle(label);
0684   h->SetTitle(title);
0685   ME[name] = ibooker.book1D(name.Data(), h);
0686   delete h;
0687 }
0688 
0689 int HeavyFlavorValidation::getFilterLevel(const std::string &moduleName, const HLTConfigProvider &hltConfig) {
0690   // helper lambda to check if a string contains a substring
0691   const auto contains = [](const std::string &s, const std::string &sub) -> bool {
0692     return s.find(sub) != std::string::npos;
0693   };
0694 
0695   // helper lambda to check if a string contains any of a list of substrings
0696   const auto containsAny = [](const std::string &s, const std::vector<std::string> &subs) -> bool {
0697     for (const auto &sub : subs) {
0698       if (s.find(sub) != std::string::npos)
0699         return true;
0700     }
0701     return false;
0702   };
0703 
0704   // helper lambda to check if string s is any of the strings in vector ms
0705   const auto isAnyOf = [](const std::string &s, const std::vector<std::string> &ms) -> bool {
0706     for (const auto &m : ms) {
0707       if (s == m)
0708         return true;
0709     }
0710     return false;
0711   };
0712 
0713   // tmadlener, 20.08.2017:
0714   // define the valid module names for the different "levels", to add a little bit more stability
0715   // to the checking compared to just doing some name matching.
0716   // Note, that the name matching is not completely remved, since at level 4 and 5 some of the
0717   // valid modules are the same, so that the name matching is still needed.
0718   // With the current definition this yields the exact same levels as before, but weeds out some
0719   // of the "false" positives at level 3 (naming matches also to some HLTMuonL1TFilter modules due to
0720   // the 'forIterL3' in the name)
0721   const std::string l1Filter = "HLTMuonL1TFilter";
0722   const std::string l2Filter = "HLTMuonL2FromL1TPreFilter";
0723   const std::vector<std::string> l3Filters = {"HLTMuonDimuonL3Filter", "HLTMuonL3PreFilter"};
0724   const std::vector<std::string> l4Filters = {
0725       "HLTDisplacedmumuFilter", "HLTDiMuonGlbTrkFilter", "HLTMuonTrackMassFilter"};
0726   const std::vector<std::string> l5Filters = {"HLTmumutkFilter", "HLT2MuonMuonDZ", "HLTDisplacedmumuFilter"};
0727 
0728   if (contains(moduleName, "Filter") && hltConfig.moduleEDMType(moduleName) == "EDFilter") {
0729     if (contains(moduleName, "L1") && !contains(moduleName, "ForIterL3") &&
0730         hltConfig.moduleType(moduleName) == l1Filter) {
0731       return 1;
0732     }
0733     if (contains(moduleName, "L2") && hltConfig.moduleType(moduleName) == l2Filter) {
0734       return 2;
0735     }
0736     if (contains(moduleName, "L3") && isAnyOf(hltConfig.moduleType(moduleName), l3Filters)) {
0737       return 3;
0738     }
0739     if (containsAny(moduleName, {"DisplacedmumuFilter", "DiMuon", "MuonL3Filtered", "TrackMassFiltered"}) &&
0740         isAnyOf(hltConfig.moduleType(moduleName), l4Filters)) {
0741       return 4;
0742     }
0743     if (containsAny(moduleName, {"Vertex", "Dz"}) && isAnyOf(hltConfig.moduleType(moduleName), l5Filters)) {
0744       return 5;
0745     }
0746   }
0747 
0748   return -1;
0749 }
0750 
0751 HeavyFlavorValidation::~HeavyFlavorValidation() {}
0752 
0753 //define this as a plug-in
0754 DEFINE_FWK_MODULE(HeavyFlavorValidation);
[ 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!