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   std::vector<Plot1DInBinVariable> genJetPlots_matched;
0118   std::vector<Plot1DInBinVariable> genJetPlots_unmatched;
0119   std::vector<Plot1DInBinVariable> recoJetPlots;
0120   std::vector<Plot1DInBinVariable> recoJetPlots_matched;
0121   std::vector<Plot1DInBinVariable> recoJetPlots_unmatched;
0122   // Is this data or MC?
0123   bool isMC;
0124
0125   float jetDeltaR;
0126
0127   bool genJetsOn, recoJetsOn;
0128
0129   std::string jetCollectionName;
0130
0131   edm::InputTag recoJetsLabel;
0132   edm::InputTag genJetsLabel;
0133   edm::EDGetTokenT<edm::View<pat::Jet>> recoJetsToken;
0134   edm::EDGetTokenT<edm::View<reco::Jet>> genJetsToken;
0135   edm::EDGetTokenT<reco::CandViewMatchMap> srcRefToJetMap;
0136
0137   void fillJetResponse(edm::View<pat::Jet>& recoJetCollection, edm::View<reco::Jet>& genJetCollection);
0138   void prepareJetResponsePlots(const std::vector<edm::ParameterSet>& genjet_plots_pset);
0139   void prepareGenJetPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset);
0140   void prepareGenJetMatchedPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset);
0141   void prepareGenJetUnmatchedPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset);
0142   void prepareRecoJetPlots(const std::vector<edm::ParameterSet>& recojet_plots_pset);
0143   void prepareRecoJetMatchedPlots(const std::vector<edm::ParameterSet>& recojet_plots_pset);
0144   void prepareRecoJetUnmatchedPlots(const std::vector<edm::ParameterSet>& recojet_plots_pset);
0145 };
0146
0147 void PFJetAnalyzerDQM::prepareJetResponsePlots(const std::vector<edm::ParameterSet>& response_plots) {
0148   for (auto& pset : response_plots) {
0149     //Low and high edges of the pt and eta bins for jets to pass to be filled into this histogram
0150     const auto ptbin_low = pset.getParameter<double>("ptBinLow");
0151     const auto ptbin_high = pset.getParameter<double>("ptBinHigh");
0152     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0153     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0154
0155     const auto response_nbins = pset.getParameter<uint32_t>("responseNbins");
0156     const auto response_low = pset.getParameter<double>("responseLow");
0157     const auto response_high = pset.getParameter<double>("responseHigh");
0158
0159     const auto name = pset.getParameter<std::string>("name");
0160     const auto title = pset.getParameter<std::string>("title");
0161
0162     // for title of raw jet response histograms
0163     auto rawTitle = title;
0164     rawTitle = rawTitle.replace(rawTitle.begin(), rawTitle.begin(), "Raw ");
0165
0166     jetResponsePlots.push_back(Plot1DInBin(
0167         name, title, response_nbins, response_low, response_high, ptbin_low, ptbin_high, etabin_low, etabin_high));
0168
0169     jetResponsePlots_noJEC.push_back(Plot1DInBin(
0170         name, rawTitle, response_nbins, response_low, response_high, ptbin_low, ptbin_high, etabin_low, etabin_high));
0171   }
0172   if (jetResponsePlots.size() > 200) {
0173     throw std::runtime_error("Requested too many jet response plots, aborting as this seems unusual.");
0174   }
0175 }
0176
0177 void PFJetAnalyzerDQM::prepareGenJetPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset) {
0178   for (auto& pset : genjet_plots_pset) {
0179     const auto name = pset.getParameter<std::string>("name");
0180     const auto title = pset.getParameter<std::string>("title");
0181
0182     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0183     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0184     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0185
0186     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0187     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0188
0189     genJetPlots.push_back(Plot1DInBinVariable(
0190         name,
0191         title,
0192         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0193         0.0,
0194         0.0,
0195         etabin_low,
0196         etabin_high));
0197   }
0198 }
0199
0200 void PFJetAnalyzerDQM::prepareGenJetMatchedPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset) {
0201   for (auto& pset : genjet_plots_pset) {
0202     const auto name = pset.getParameter<std::string>("name") + "_matched";
0203     const auto title = "Matched " + pset.getParameter<std::string>("title");
0204
0205     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0206     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0207     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0208
0209     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0210     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0211
0212     genJetPlots_matched.push_back(Plot1DInBinVariable(
0213         name,
0214         title,
0215         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0216         0.0,
0217         0.0,
0218         etabin_low,
0219         etabin_high));
0220   }
0221 }
0222
0223 void PFJetAnalyzerDQM::prepareGenJetUnmatchedPlots(const std::vector<edm::ParameterSet>& genjet_plots_pset) {
0224   for (auto& pset : genjet_plots_pset) {
0225     const auto name = pset.getParameter<std::string>("name") + "_unmatched";
0226     const auto title = "Unmatched " + pset.getParameter<std::string>("title");
0227
0228     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0229     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0230     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0231
0232     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0233     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0234
0235     genJetPlots_unmatched.push_back(Plot1DInBinVariable(
0236         name,
0237         title,
0238         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0239         0.0,
0240         0.0,
0241         etabin_low,
0242         etabin_high));
0243   }
0244 }
0245
0246 void PFJetAnalyzerDQM::prepareRecoJetPlots(const std::vector<edm::ParameterSet>& recojet_plots_pset) {
0247   for (auto& pset : recojet_plots_pset) {
0248     const auto name = pset.getParameter<std::string>("name");
0249     const auto title = pset.getParameter<std::string>("title");
0250
0251     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0252     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0253     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0254
0255     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0256     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0257
0258     recoJetPlots.push_back(Plot1DInBinVariable(
0259         name,
0260         title,
0261         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0262         0.0,
0263         0.0,
0264         etabin_low,
0265         etabin_high));
0266   }
0267 }
0268
0269 void PFJetAnalyzerDQM::prepareRecoJetMatchedPlots(const std::vector<edm::ParameterSet>& recojet_plots_pset) {
0270   for (auto& pset : recojet_plots_pset) {
0271     const auto name = pset.getParameter<std::string>("name") + "_matched";
0272     const auto title = "Matched " + pset.getParameter<std::string>("title");
0273
0274     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0275     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0276     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0277
0278     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0279     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0280
0281     recoJetPlots_matched.push_back(Plot1DInBinVariable(
0282         name,
0283         title,
0284         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0285         0.0,
0286         0.0,
0287         etabin_low,
0288         etabin_high));
0289   }
0290 }
0291
0292 void PFJetAnalyzerDQM::prepareRecoJetUnmatchedPlots(const std::vector<edm::ParameterSet>& recojet_plots_pset) {
0293   for (auto& pset : recojet_plots_pset) {
0294     const auto name = pset.getParameter<std::string>("name") + "_unmatched";
0295     const auto title = "Unmatched " + pset.getParameter<std::string>("title");
0296
0297     //Low and high edges of the eta bins for jets to pass to be filled into this histogram
0298     const auto ptbins_d = pset.getParameter<std::vector<double>>("ptBins");
0299     std::vector<float> ptbins(ptbins_d.begin(), ptbins_d.end());
0300
0301     const auto etabin_low = pset.getParameter<double>("etaBinLow");
0302     const auto etabin_high = pset.getParameter<double>("etaBinHigh");
0303
0304     recoJetPlots_unmatched.push_back(Plot1DInBinVariable(
0305         name,
0306         title,
0307         std::make_unique<TH1F>(name.c_str(), title.c_str(), static_cast<int>(ptbins.size()) - 1, ptbins.data()),
0308         0.0,
0309         0.0,
0310         etabin_low,
0311         etabin_high));
0312   }
0313 }
0314
0315 PFJetAnalyzerDQM::PFJetAnalyzerDQM(const edm::ParameterSet& iConfig) {
0316   recoJetsLabel = iConfig.getParameter<edm::InputTag>("recoJetCollection");
0317   genJetsLabel = iConfig.getParameter<edm::InputTag>("genJetCollection");
0318
0319   //label for making new folder
0320   jetCollectionName = recoJetsLabel.label();
0321
0322   //DeltaR for reco to gen jet matching
0323   jetDeltaR = iConfig.getParameter<double>("jetDeltaR");
0324
0325   //for turn genJet on/off
0326   genJetsOn = iConfig.getParameter<bool>("genJetsOn");
0327   recoJetsOn = iConfig.getParameter<bool>("recoJetsOn");
0328
0329   //Create all jet response plots in bins of genjet pt and eta
0330   const auto& response_plots = iConfig.getParameter<std::vector<edm::ParameterSet>>("responsePlots");
0331   prepareJetResponsePlots(response_plots);
0332
0333   const auto& genjet_plots = iConfig.getParameter<std::vector<edm::ParameterSet>>("genJetPlots");
0334   prepareGenJetPlots(genjet_plots);
0335   prepareGenJetMatchedPlots(genjet_plots);
0336   prepareGenJetUnmatchedPlots(genjet_plots);
0337
0338   const auto& recojet_plots = iConfig.getParameter<std::vector<edm::ParameterSet>>("recoJetPlots");
0339   prepareRecoJetPlots(recojet_plots);
0340   prepareRecoJetMatchedPlots(recojet_plots);
0341   prepareRecoJetUnmatchedPlots(recojet_plots);
0342
0343   recoJetsToken = consumes<edm::View<pat::Jet>>(recoJetsLabel);
0344   genJetsToken = consumes<edm::View<reco::Jet>>(genJetsLabel);
0345 }
0346
0347 void PFJetAnalyzerDQM::fillJetResponse(edm::View<pat::Jet>& recoJetCollection, edm::View<reco::Jet>& genJetCollection) {
0348   //match gen jets to reco jets, require minimum jetDeltaR, choose closest, do not try to match charge
0349   std::vector<int> matchIndices;
0350   std::vector<int> matchIndicesReco;
0351   PFB::match(genJetCollection, recoJetCollection, matchIndices, false, jetDeltaR);
0352   PFB::match(recoJetCollection, genJetCollection, matchIndicesReco, false, jetDeltaR);
0353
0354   //Fill recojet pt if recoJetOn
0355   for (unsigned int i = 0; i < recoJetCollection.size(); i++) {
0356     const auto& recoJet = recoJetCollection.at(i);
0357     const auto pt_reco = recoJet.pt();
0358     const auto eta_reco = abs(recoJet.eta());
0359     const int iMatch_reco = matchIndicesReco[i];
0360     if (recoJetsOn) {
0361       for (auto& plot : recoJetPlots) {
0362         if (plot.isInEtaBin(eta_reco)) {
0363           plot.fill(pt_reco);
0364         }
0365       }
0366       if (iMatch_reco != -1) {
0367         for (auto& plot : recoJetPlots_matched) {
0368           if (plot.isInEtaBin(eta_reco)) {
0369             plot.fill(pt_reco);
0370           }
0371         }
0372       } else {
0373         for (auto& plot : recoJetPlots_unmatched) {
0374           if (plot.isInEtaBin(eta_reco)) {
0375             plot.fill(pt_reco);
0376           }
0377         }
0378       }
0379     }
0380   }
0381
0382   for (unsigned int i = 0; i < genJetCollection.size(); i++) {
0383     const auto& genJet = genJetCollection.at(i);
0384     const auto pt_gen = genJet.pt();
0385     const auto eta_gen = abs(genJet.eta());
0386     const int iMatch = matchIndices[i];
0387
0388     //Fill genjet pt if genJetOn
0389     if (genJetsOn) {
0390       for (auto& plot : genJetPlots) {
0391         if (plot.isInEtaBin(eta_gen)) {
0392           plot.fill(pt_gen);
0393         }
0394       }
0395     }
0396     if (recoJetsOn) {
0397       if (iMatch != -1) {
0398         for (auto& plot : genJetPlots_matched) {
0399           if (plot.isInEtaBin(eta_gen)) {
0400             plot.fill(pt_gen);
0401           }
0402         }
0403       } else {
0404         for (auto& plot : genJetPlots_unmatched) {
0405           if (plot.isInEtaBin(eta_gen)) {
0406             plot.fill(pt_gen);
0407           }
0408         }
0409       }
0410     }
0411
0412     //If gen jet had a matched reco jet
0413     if (iMatch != -1) {
0414       const auto& recoJet = recoJetCollection[iMatch];
0415       auto pt_reco = recoJet.pt();
0416
0417       const auto response = pt_reco / pt_gen;
0418       const auto response_raw = pt_reco * recoJet.jecFactor("Uncorrected") / pt_gen;
0419
0420       //Loop linearly through all plots and check if they match the pt and eta bin
0421       //this is not algorithmically optimal but we don't expect to more than a few hundred plots
0422       //If this turns out to be a problem, can easily make a 2D-map for indices
0423       for (auto& plot : jetResponsePlots) {
0424         if (plot.isInPtEtaBin(pt_gen, eta_gen)) {
0425           plot.fill(response);
0426         }
0427       }
0428       // this loop should be for NoJEC plots
0429       for (auto& plot : jetResponsePlots_noJEC) {
0430         if (plot.isInPtEtaBin(pt_gen, eta_gen)) {
0431           plot.fill(response_raw);
0432         }
0433       }
0434     }
0435   }
0436 }
0437
0438 void PFJetAnalyzerDQM::bookHistograms(DQMStore::IBooker& booker, edm::Run const&, edm::EventSetup const&) {
0439   booker.setCurrentFolder("ParticleFlow/JetResponse/" + jetCollectionName + "/JEC/");
0440   for (auto& plot : jetResponsePlots) {
0441     plot.book(booker);
0442   }
0443   //Book plots for noJEC
0444   booker.setCurrentFolder("ParticleFlow/JetResponse/" + jetCollectionName + "/noJEC/");
0445   for (auto& plot : jetResponsePlots_noJEC) {
0446     plot.book(booker);
0447   }
0448
0449   if (recoJetsOn) {
0450     booker.setCurrentFolder("ParticleFlow/JetResponse/" + jetCollectionName + "/noJEC/");
0451     for (auto& plot : genJetPlots_matched) {
0452       plot.book(booker);
0453     }
0454     for (auto& plot : genJetPlots_unmatched) {
0455       plot.book(booker);
0456     }
0457     booker.setCurrentFolder("ParticleFlow/JetResponse/" + jetCollectionName + "/JEC/");
0458     for (auto& plot : recoJetPlots) {
0459       plot.book(booker);
0460     }
0461     for (auto& plot : recoJetPlots_matched) {
0462       plot.book(booker);
0463     }
0464     for (auto& plot : recoJetPlots_unmatched) {
0465       plot.book(booker);
0466     }
0467   }
0468
0469   //Book plots for gen-jet pt spectra
0470   if (genJetsOn) {
0471     booker.setCurrentFolder("ParticleFlow/GenJets/");
0472     for (auto& plot : genJetPlots) {
0473       plot.book(booker);
0474     }
0475   }
0476 }
0477
0478 void PFJetAnalyzerDQM::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0479   edm::Handle<edm::View<pat::Jet>> recoJetCollectionHandle;
0480   iEvent.getByToken(recoJetsToken, recoJetCollectionHandle);
0481   auto recoJetCollection = *recoJetCollectionHandle;
0482
0483   isMC = !iEvent.isRealData();
0484
0485   if (isMC) {
0486     edm::Handle<edm::View<reco::Jet>> genJetCollectionHandle;
0487     iEvent.getByToken(genJetsToken, genJetCollectionHandle);
0488     auto genJetCollection = *genJetCollectionHandle;
0489
0490     fillJetResponse(recoJetCollection, genJetCollection);
0491   }
0492 }
0493
0494 #include "FWCore/Framework/interface/MakerMacros.h"
0495 DEFINE_FWK_MODULE(PFJetAnalyzerDQM);