RecoParticleFlow/plugins/PFJetAnalyzerDQM.cc

0001 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0002 #include "DQMOffline/PFTau/interface/Matchers.h"
0003 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0004 #include "DQMServices/Core/interface/DQMStore.h"
0005 #include "DataFormats/JetReco/interface/Jet.h"
0006 #include "DataFormats/JetReco/interface/PFJet.h"
0007 #include "DataFormats/PatCandidates/interface/Jet.h"
0008 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013
0014 #include <algorithm>
0015 #include <numeric>
0016 #include <regex>
0017 #include <sstream>
0018 #include <vector>
0019 #include <memory>
0020
0021 class PFJetAnalyzerDQM : public DQMEDAnalyzer {
0022 public:
0023   PFJetAnalyzerDQM(const edm::ParameterSet&);
0024   void analyze(const edm::Event&, const edm::EventSetup&) override;
0025
0026 protected:
0027   //Book histograms
0028   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0029
0030 private:
0031   class Plot1DInBin {
0032   public:
0033     const std::string name, title;
0034     const uint32_t nbins;
0035     const float min, max;
0036     const float ptbin_low, ptbin_high, etabin_low, etabin_high;
0037     MonitorElement* plot_;
0038
0039     Plot1DInBin(const std::string _name,
0040                 const std::string _title,
0041                 const uint32_t _nbins,
0042                 const float _min,
0043                 const float _max,
0044                 float _ptbin_low,
0045                 float _ptbin_high,
0046                 float _etabin_low,
0047                 float _etabin_high)
0048         : name(_name),
0049           title(_title),
0050           nbins(_nbins),
0051           min(_min),
0052           max(_max),
0053           ptbin_low(_ptbin_low),
0054           ptbin_high(_ptbin_high),
0055           etabin_low(_etabin_low),
0056           etabin_high(_etabin_high) {}
0057
0058     void book(DQMStore::IBooker& booker) { plot_ = booker.book1D(name, title, nbins, min, max); }
0059
0060     void fill(float value) {
0061       assert(plot_ != nullptr);
0062       plot_->Fill(value);
0063     }
0064
0065     //Check if a jet with a value v would be in the bin that applies to this plot
0066     bool isInBin(float v, float low, float high) { return v >= low && v < high; }
0067
0068     bool isInPtBin(float pt) { return isInBin(pt, ptbin_low, ptbin_high); }
0069
0070     bool isInEtaBin(float eta) { return isInBin(eta, etabin_low, etabin_high); }
0071
0072     bool isInPtEtaBin(float pt, float eta) { return isInPtBin(pt) && isInEtaBin(eta); }
0073   };
0074
0075   class Plot1DInBinVariable {
0076   public:
0077     const std::string name, title;
0078     std::unique_ptr<TH1F> base_hist;
0079     const float ptbin_low, ptbin_high, etabin_low, etabin_high;
0080     MonitorElement* plot_;
0081
0082     Plot1DInBinVariable(const std::string _name,
0083                         const std::string _title,
0084                         std::unique_ptr<TH1F> _base_hist,
0085                         float _ptbin_low,
0086                         float _ptbin_high,
0087                         float _etabin_low,
0088                         float _etabin_high)
0089         : name(_name),
0090           title(_title),
0091           base_hist(std::move(_base_hist)),
0092           ptbin_low(_ptbin_low),
0093           ptbin_high(_ptbin_high),
0094           etabin_low(_etabin_low),
0095           etabin_high(_etabin_high) {}
0096
0097     void book(DQMStore::IBooker& booker) { plot_ = booker.book1D(name.c_str(), base_hist.get()); }
0098
0099     void fill(float value) {
0100       assert(plot_ != nullptr);
0101       plot_->Fill(value);
0102     }
0103
0104     //Check if a jet with a value v would be in the bin that applies to this plot
0105     bool isInBin(float v, float low, float high) { return v >= low && v < high; }
0106
0107     bool isInPtBin(float pt) { return isInBin(pt, ptbin_low, ptbin_high); }
0108
0109     bool isInEtaBin(float eta) { return isInBin(eta, etabin_low, etabin_high); }
0110
0111     bool isInPtEtaBin(float pt, float eta) { return isInPtBin(pt) && isInEtaBin(eta); }
0112   };
0113
0114   std::vector<Plot1DInBin> jetResponsePlots;
0115   std::vector<Plot1DInBin> jetResponsePlots_noJEC;
0116   std::vector<Plot1DInBinVariable> genJetPlots;
0117
0118   // Is this data or MC?
0119   bool isMC;
0120
0121   float jetDeltaR;
0122
0123   bool genJetsOn;
0124
0125   std::string jetCollectionName;
0126
0127   edm::InputTag recoJetsLabel;
0128   edm::InputTag genJetsLabel;
0129   edm::EDGetTokenT<edm::View<pat::Jet>> recoJetsToken;
0130   edm::EDGetTokenT<edm::View<reco::Jet>> genJetsToken;
0131   edm::EDGetTokenT<reco::CandViewMatchMap> srcRefToJetMap;
0132
0133   void fillJetResponse(edm::View<pat::Jet>& recoJetCollection, edm::View<reco::Jet>& genJetCollection);
0134   void prepareJetResponsePlots(const std::vector<edm::ParameterSet>& genjet_plots_pset);
0135   void prepareGenJetPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset);
0136 };
0137
0138 void PFJetAnalyzerDQM::prepareJetResponsePlots(const std::vector<edm::ParameterSet>& response_plots) {
0139   for (auto& pset : response_plots) {
0140     //Low and high edges of the pt and eta bins for jets to pass to be filled into this histogram
0141     const auto ptbin_low = pset.getParameter<double>("ptBinLow");
0142     const auto ptbin_high = pset.getParameter<double>("ptBinHigh");
0143     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0144     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0145
0146     const auto response_nbins = pset.getParameter<uint32_t>("responseNbins");
0147     const auto response_low = pset.getParameter<double>("responseLow");
0148     const auto response_high = pset.getParameter<double>("responseHigh");
0149
0150     const auto name = pset.getParameter<std::string>("name");
0151     const auto title = pset.getParameter<std::string>("title");
0152
0153     // for title of raw jet response histograms
0154     auto rawTitle = title;
0155     rawTitle = rawTitle.replace(rawTitle.begin(), rawTitle.begin(), "Raw ");
0156
0157     jetResponsePlots.push_back(Plot1DInBin(
0158         name, title, response_nbins, response_low, response_high, ptbin_low, ptbin_high, etabin_low, etabin_high));
0159
0160     jetResponsePlots_noJEC.push_back(Plot1DInBin(
0161         name, rawTitle, response_nbins, response_low, response_high, ptbin_low, ptbin_high, etabin_low, etabin_high));
0162   }
0163   if (jetResponsePlots.size() > 200) {
0164     throw std::runtime_error("Requested too many jet response plots, aborting as this seems unusual.");
0165   }
0166 }
0167
0168 void PFJetAnalyzerDQM::prepareGenJetPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset) {
0169   for (auto& pset : genjet_plots_pset) {
0170     const auto name = pset.getParameter<std::string>("name");
0171     const auto title = pset.getParameter<std::string>("title");
0172
0173     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0174     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0175     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0176
0177     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0178     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0179
0180     genJetPlots.push_back(Plot1DInBinVariable(
0181         name,
0182         title,
0183         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0184         0.0,
0185         0.0,
0186         etabin_low,
0187         etabin_high));
0188   }
0189 }
0190
0191 PFJetAnalyzerDQM::PFJetAnalyzerDQM(const edm::ParameterSet& iConfig) {
0192   recoJetsLabel = iConfig.getParameter<edm::InputTag>("recoJetCollection");
0193   genJetsLabel = iConfig.getParameter<edm::InputTag>("genJetCollection");
0194
0195   //label for making new folder
0196   jetCollectionName = recoJetsLabel.label();
0197
0198   //DeltaR for reco to gen jet matching
0199   jetDeltaR = iConfig.getParameter<double>("jetDeltaR");
0200
0201   //for turn genJet on/off
0202   genJetsOn = iConfig.getParameter<bool>("genJetsOn");
0203
0204   //Create all jet response plots in bins of genjet pt and eta
0205   const auto& response_plots = iConfig.getParameter<std::vector<edm::ParameterSet>>("responsePlots");
0206   prepareJetResponsePlots(response_plots);
0207
0208   const auto& genjet_plots = iConfig.getParameter<std::vector<edm::ParameterSet>>("genJetPlots");
0209   prepareGenJetPlots(genjet_plots);
0210
0211   recoJetsToken = consumes<edm::View<pat::Jet>>(recoJetsLabel);
0212   genJetsToken = consumes<edm::View<reco::Jet>>(genJetsLabel);
0213 }
0214
0215 void PFJetAnalyzerDQM::fillJetResponse(edm::View<pat::Jet>& recoJetCollection, edm::View<reco::Jet>& genJetCollection) {
0216   //match gen jets to reco jets, require minimum jetDeltaR, choose closest, do not try to match charge
0217   std::vector<int> matchIndices;
0218   PFB::match(genJetCollection, recoJetCollection, matchIndices, false, jetDeltaR);
0219
0220   for (unsigned int i = 0; i < genJetCollection.size(); i++) {
0221     const auto& genJet = genJetCollection.at(i);
0222     const auto pt_gen = genJet.pt();
0223     const auto eta_gen = abs(genJet.eta());
0224     const int iMatch = matchIndices[i];
0225
0226     //Fill genjet pt if genJetOn
0227     if (genJetsOn) {
0228       for (auto& plot : genJetPlots) {
0229         if (plot.isInEtaBin(eta_gen)) {
0230           plot.fill(pt_gen);
0231         }
0232       }
0233     }
0234
0235     //If gen jet had a matched reco jet
0236     if (iMatch != -1) {
0237       const auto& recoJet = recoJetCollection[iMatch];
0238       auto pt_reco = recoJet.pt();
0239
0240       const auto response = pt_reco / pt_gen;
0241       const auto response_raw = pt_reco * recoJet.jecFactor("Uncorrected") / pt_gen;
0242
0243       //Loop linearly through all plots and check if they match the pt and eta bin
0244       //this is not algorithmically optimal but we don't expect to more than a few hundred plots
0245       //If this turns out to be a problem, can easily make a 2D-map for indices
0246       for (auto& plot : jetResponsePlots) {
0247         if (plot.isInPtEtaBin(pt_gen, eta_gen)) {
0248           plot.fill(response);
0249         }
0250       }
0251       // this loop should be for NoJEC plots
0252       for (auto& plot : jetResponsePlots_noJEC) {
0253         if (plot.isInPtEtaBin(pt_gen, eta_gen)) {
0254           plot.fill(response_raw);
0255         }
0256       }
0257     }
0258   }
0259 }
0260
0261 void PFJetAnalyzerDQM::bookHistograms(DQMStore::IBooker& booker, edm::Run const&, edm::EventSetup const&) {
0262   booker.setCurrentFolder("ParticleFlow/JetResponse/" + jetCollectionName + "/JEC/");
0263   for (auto& plot : jetResponsePlots) {
0264     plot.book(booker);
0265   }
0266   //Book plots for noJEC
0267   booker.setCurrentFolder("ParticleFlow/JetResponse/" + jetCollectionName + "/noJEC/");
0268   for (auto& plot : jetResponsePlots_noJEC) {
0269     plot.book(booker);
0270   }
0271
0272   //Book plots for gen-jet pt spectra
0273   if (genJetsOn) {
0274     booker.setCurrentFolder("ParticleFlow/GenJets/");
0275     for (auto& plot : genJetPlots) {
0276       plot.book(booker);
0277     }
0278   }
0279 }
0280
0281 void PFJetAnalyzerDQM::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0282   edm::Handle<edm::View<pat::Jet>> recoJetCollectionHandle;
0283   iEvent.getByToken(recoJetsToken, recoJetCollectionHandle);
0284   auto recoJetCollection = *recoJetCollectionHandle;
0285
0286   isMC = !iEvent.isRealData();
0287
0288   if (isMC) {
0289     edm::Handle<edm::View<reco::Jet>> genJetCollectionHandle;
0290     iEvent.getByToken(genJetsToken, genJetCollectionHandle);
0291     auto genJetCollection = *genJetCollectionHandle;
0292
0293     fillJetResponse(recoJetCollection, genJetCollection);
0294   }
0295 }
0296
0297 #include "FWCore/Framework/interface/MakerMacros.h"
0298 DEFINE_FWK_MODULE(PFJetAnalyzerDQM);