diff --git a/data/PUProfiles/pu_histogram_all_2015.root b/data/PUProfiles/pu_histogram_all_2015.root new file mode 100644 index 0000000..b6fdda0 Binary files /dev/null and b/data/PUProfiles/pu_histogram_all_2015.root differ diff --git a/interface/EventProducer.h b/interface/EventProducer.h index af3aabb..d56acc7 100644 --- a/interface/EventProducer.h +++ b/interface/EventProducer.h @@ -2,6 +2,7 @@ #define EVENT_PRODUCER #include +#include #include #include @@ -13,6 +14,9 @@ class EventProducer: public Framework::Producer { Producer(name, tree, config) { + if (config.getUntrackedParameter("compute_pu_weights", true)) + m_pu_reweighter = std::make_shared(config.getParameterSet("pu_reweighter"), Framework::PUProfile::Run2015_25ns); + } virtual ~EventProducer() {} @@ -38,6 +42,8 @@ class EventProducer: public Framework::Producer { float m_event_weight_sum = 0; + std::shared_ptr m_pu_reweighter; + public: // Tree members @@ -49,6 +55,7 @@ class EventProducer: public Framework::Producer { BRANCH(npu, int); BRANCH(true_interactions, float); + BRANCH(pu_weight, float); BRANCH(pt_hat, float); BRANCH(weight, float); diff --git a/interface/PUReweighter.h b/interface/PUReweighter.h new file mode 100644 index 0000000..cb51fca --- /dev/null +++ b/interface/PUReweighter.h @@ -0,0 +1,86 @@ +#pragma once + +#include +#include +#include + +#include + +#include + +namespace Framework { + enum PUProfile : uint8_t { + Run2015_50ns, + Run2015_25ns + }; + + class PUReweighter { + public: + PUReweighter(const edm::ParameterSet&, PUProfile mc_pu_profile); + + float getWeight(float n_interactions); + + private: + std::shared_ptr m_pu_weights; + + std::map> m_pu_profiles { + { + Run2015_25ns, + { + 4.8551E-07, + 1.74806E-06, + 3.30868E-06, + 1.62972E-05, + 4.95667E-05, + 0.000606966, + 0.003307249, + 0.010340741, + 0.022852296, + 0.041948781, + 0.058609363, + 0.067475755, + 0.072817826, + 0.075931405, + 0.076782504, + 0.076202319, + 0.074502547, + 0.072355135, + 0.069642102, + 0.064920999, + 0.05725576, + 0.047289348, + 0.036528446, + 0.026376131, + 0.017806872, + 0.011249422, + 0.006643385, + 0.003662904, + 0.001899681, + 0.00095614, + 0.00050028, + 0.000297353, + 0.000208717, + 0.000165856, + 0.000139974, + 0.000120481, + 0.000103826, + 8.88868E-05, + 7.53323E-05, + 6.30863E-05, + 5.21356E-05, + 4.24754E-05, + 3.40876E-05, + 2.69282E-05, + 2.09267E-05, + 1.5989E-05, + 4.8551E-06, + 2.42755E-06, + 4.8551E-07, + 2.42755E-07, + 1.21378E-07, + 4.8551E-08 + } + } + }; + }; +}; diff --git a/interface/Tools.h b/interface/Tools.h index 45c19ea..91dfce0 100644 --- a/interface/Tools.h +++ b/interface/Tools.h @@ -1,5 +1,7 @@ #pragma once +#include + namespace pat { class Jet; } @@ -23,4 +25,8 @@ namespace Tools { bool passTightId(const pat::Jet& jet); bool passTightLeptonVetoId(const pat::Jet& jet); }; + + inline bool caseInsensitiveEquals(const std::string& a, const std::string& b) { + return std::equal(a.begin(), a.end(), b.begin(), [](char a, char b) { return std::tolower(a) == std::tolower(b); }); + } }; diff --git a/python/EventProducer.py b/python/EventProducer.py index 506e72d..c7f71f7 100644 --- a/python/EventProducer.py +++ b/python/EventProducer.py @@ -5,6 +5,10 @@ prefix = cms.string('event_'), enable = cms.bool(True), parameters = cms.PSet( - pu_summary = cms.untracked.InputTag('slimmedAddPileupInfo') + pu_summary = cms.untracked.InputTag('slimmedAddPileupInfo'), + compute_pu_weights = cms.untracked.bool(True), + pu_reweighter = cms.PSet( + data_pu_profile = cms.untracked.FileInPath('cp3_llbb/Framework/data/PUProfiles/pu_histogram_all_2015.root') + ) ) ) diff --git a/src/EventProducer.cc b/src/EventProducer.cc index 6219686..5e6eea5 100644 --- a/src/EventProducer.cc +++ b/src/EventProducer.cc @@ -11,6 +11,8 @@ void EventProducer::produce(edm::Event& event_, const edm::EventSetup& eventSetu rho = *rho_handle; + pu_weight = 1; + edm::Handle> pu_infos; if (event_.getByToken(m_pu_info_token, pu_infos)) { for (const auto& pu_info: *pu_infos) { @@ -19,6 +21,9 @@ void EventProducer::produce(edm::Event& event_, const edm::EventSetup& eventSetu npu = pu_info.getPU_NumInteractions(); true_interactions = pu_info.getTrueNumInteractions(); + + if (m_pu_reweighter.get()) + pu_weight = m_pu_reweighter->getWeight(true_interactions); } } diff --git a/src/PUReweighter.cc b/src/PUReweighter.cc new file mode 100644 index 0000000..49014dc --- /dev/null +++ b/src/PUReweighter.cc @@ -0,0 +1,53 @@ +#include +#include + +#include +#include + +#include + +namespace Framework { + + PUReweighter::PUReweighter(const edm::ParameterSet& config, PUProfile mc_pu_profile) { + + std::string file = config.getUntrackedParameter("data_pu_profile").fullPath(); + + // Retrieve PU histogram from file + std::unique_ptr f(TFile::Open(file.c_str())); + TH1* data_pu_histogram = static_cast(f->Get("pileup")); + + size_t x_max = (size_t) data_pu_histogram->GetXaxis()->GetXmax(); + size_t n_bins = data_pu_histogram->GetNbinsX(); + + std::shared_ptr mc_pu_histogram = std::make_shared("pu_mc", "", n_bins, 0, x_max); + mc_pu_histogram->SetDirectory(nullptr); + + const std::vector& coeffs = m_pu_profiles[mc_pu_profile]; + for (size_t i = 1; i <= (size_t) data_pu_histogram->GetNbinsX(); i++) { + size_t coef_index = (size_t) data_pu_histogram->GetBinCenter(i); + float coef = (coef_index) < coeffs.size() ? coeffs[coef_index] : 0.; + + mc_pu_histogram->SetBinContent(i, coef); + } + + data_pu_histogram->Scale(1. / data_pu_histogram->Integral()); + mc_pu_histogram->Scale(1. / mc_pu_histogram->Integral()); + + std::shared_ptr pu_weights(static_cast(data_pu_histogram->Clone())); + pu_weights->SetDirectory(nullptr); + + pu_weights->Divide(mc_pu_histogram.get()); + + m_pu_weights = pu_weights; + } + + float PUReweighter::getWeight(float n_interactions) { + + TH1* pu_weights = m_pu_weights.get(); + + if (! pu_weights) + return 1; + + return pu_weights->GetBinContent(pu_weights->GetXaxis()->FindBin(n_interactions)); + } +};