Add new weight for PU reweighting

cp3-llbb · Nov 17, 2015 · b7aec0e · b7aec0e
1 parent 0bd3ae4
commit b7aec0e
Show file tree

Hide file tree

Showing 7 changed files with 162 additions and 1 deletion.
diff --git a/data/PUProfiles/pu_histogram_all_2015.root b/data/PUProfiles/pu_histogram_all_2015.root
diff --git a/interface/EventProducer.h b/interface/EventProducer.h
@@ -2,6 +2,7 @@
 #define EVENT_PRODUCER
 
 #include <cp3_llbb/Framework/interface/Producer.h>
+#include <cp3_llbb/Framework/interface/PUReweighter.h>
 
 #include <SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h>
 #include <SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h>
@@ -13,6 +14,9 @@ class EventProducer: public Framework::Producer {
             Producer(name, tree, config)
         {
 
+            if (config.getUntrackedParameter<bool>("compute_pu_weights", true))
+                m_pu_reweighter = std::make_shared<Framework::PUReweighter>(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<Framework::PUReweighter> 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
@@ -0,0 +1,86 @@
+#pragma once
+
+#include <cstdint>
+#include <map>
+#include <memory>
+
+#include <TH1.h>
+
+#include <FWCore/ParameterSet/interface/ParameterSet.h>
+
+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<TH1> m_pu_weights;
+
+            std::map<PUProfile, std::vector<float>> 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
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <algorithm>
+
 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
@@ -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
@@ -11,6 +11,8 @@ void EventProducer::produce(edm::Event& event_, const edm::EventSetup& eventSetu
 
     rho = *rho_handle;
 
+    pu_weight = 1;
+
     edm::Handle<std::vector<PileupSummaryInfo>> 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
@@ -0,0 +1,53 @@
+#include <iostream>
+#include <memory>
+
+#include <TFile.h>
+#include <TH1F.h>
+
+#include <cp3_llbb/Framework/interface/PUReweighter.h>
+
+namespace Framework {
+
+    PUReweighter::PUReweighter(const edm::ParameterSet& config, PUProfile mc_pu_profile) {
+
+        std::string file = config.getUntrackedParameter<edm::FileInPath>("data_pu_profile").fullPath();
+
+        // Retrieve PU histogram from file
+        std::unique_ptr<TFile> f(TFile::Open(file.c_str()));
+        TH1* data_pu_histogram = static_cast<TH1*>(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<TH1F> mc_pu_histogram = std::make_shared<TH1F>("pu_mc", "", n_bins, 0, x_max);
+        mc_pu_histogram->SetDirectory(nullptr);
+
+        const std::vector<float>& 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<TH1> pu_weights(static_cast<TH1*>(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));
+    }
+};