New PU weight retrival; New weights with uncertainties

AnalysisStep/interface/PileUpWeight.h

@@ -0,0 +1,29 @@
+#ifndef PileUpWeight_h
+#define PileUpWeight_h
+
+/** \class PileUpWeight
+ *
+ *  Implement PU reweighting.
+ *
+ */
+
+#include <TH1F.h>
+#include <TFile.h>
+#include <string>
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/FileInPath.h"
+
+class PileUpWeight {
+public:
+  enum class PUvar : int {NOMINAL=0, VARUP=1, VARDOWN=2};
+
+  PileUpWeight(int MC, int target); 
+
+  float weight(float input, PUvar var = PUvar::NOMINAL);
+
+  std::unique_ptr<TH1> h_nominal;
+  std::unique_ptr<TH1> h_up;
+  std::unique_ptr<TH1> h_down;
+};
+#endif

AnalysisStep/src/PileUpWeight.cc

@@ -0,0 +1,35 @@
+#include "../interface/PileUpWeight.h"
+
+float PileUpWeight::weight(float input, PileUpWeight::PUvar var) {
+
+  if (h_nominal != nullptr && var == PUvar::NOMINAL) {
+    return h_nominal->GetBinContent(h_nominal->FindBin(input));
+  } else if (h_down != nullptr && var == PUvar::VARDOWN) {
+    return h_down->GetBinContent(h_down->FindBin(input));
+  } else if (h_up != nullptr && var == PUvar::VARUP) {
+    return h_up->GetBinContent(h_up->FindBin(input));
+  } else {
+    edm::LogError("PU reweight") << "Did not find reweighting histogram to reweight variation=";// << var;
+    return -99.;
+  }
+}
+
+
+PileUpWeight::PileUpWeight(int MC, int target) { 
+
+ if (MC==2017&&target==2017) {
+    edm::FileInPath fip("ZZAnalysis/AnalysisStep/data/PileUpWeights/puWeightsMoriond17_v2.root");
+
+    TFile *fPUWeight = TFile::Open(fip.fullPath().data(),"READ");
+
+    h_nominal.reset((TH1*)fPUWeight->Get("weights")->Clone());
+    h_up.reset((TH1*)fPUWeight->Get("weights_varUp")->Clone());
+    h_down.reset((TH1*)fPUWeight->Get("weights_varDn")->Clone());
+
+    fPUWeight->Close();
+  }
+  if(h_nominal == nullptr) {
+     edm::LogError("PU reweight") << "Did not find reweighting histogram to reweight MC=" << MC << " to data=" << target;
+     edm::LogError("PU reweight") << "Legacy configurations are in PUReweight.cc";
+  }
+}

AnalysisStep/test/Ntuplizers/HZZ4lNtupleMaker.cc

@@ -46,7 +46,9 @@
 #include <ZZAnalysis/AnalysisStep/interface/DaughterDataHelpers.h>
 #include <ZZAnalysis/AnalysisStep/interface/FinalStates.h>
 #include <ZZAnalysis/AnalysisStep/interface/MCHistoryTools.h>
-#include <ZZAnalysis/AnalysisStep/interface/PUReweight.h>
+//#include <ZZAnalysis/AnalysisStep/interface/PileUpWeight.h>
+#include <ZZAnalysis/AnalysisStep/interface/PileUpWeight.h>
+
 #include "ZZAnalysis/AnalysisStep/interface/EwkCorrections.h"
 #include "ZZAnalysis/AnalysisStep/interface/LHEHandler.h"
 #include "ZZAnalysis/AnalysisStep/src/kFactors.C"
@@ -90,6 +92,9 @@ namespace {
   Short_t NObsInt  = 0;
   Float_t NTrueInt  = 0;
   Float_t PUWeight  = 0;
+  Float_t PUWeight_Up  = 0;
+  Float_t PUWeight_Dn  = 0;
+
   Float_t KFactor_QCD_ggZZ_Nominal = 0;
   Float_t KFactor_QCD_ggZZ_PDFScaleDn = 0;
   Float_t KFactor_QCD_ggZZ_PDFScaleUp = 0;
@@ -162,6 +167,12 @@ namespace {
   //std::vector<float> LepPhotonIso;
   std::vector<float> LepCombRelIsoPF;
   std::vector<short> LepisLoose;
+  std::vector<float> LepRecoSF;
+  std::vector<float> LepRecoSF_Unc;
+  std::vector<float> LepSelSF;
+  std::vector<float> LepSelSF_Unc;
+
+
   std::vector<float> fsrPt;
   std::vector<float> fsrEta;
   std::vector<float> fsrPhi;
@@ -328,6 +339,7 @@ class HZZ4lNtupleMaker : public edm::EDAnalyzer {
     const math::XYZTLorentzVector Lep1, const math::XYZTLorentzVector Lep2, const math::XYZTLorentzVector Lep3, const math::XYZTLorentzVector Lep4);
   void FillAssocLepGenInfo(std::vector<const reco::Candidate *>& AssocLeps);
 
+
   Float_t getAllWeight(const reco::Candidate* Lep) const;
   Float_t getHqTWeight(double mH, double genPt) const;
   Float_t getFakeWeight(Float_t LepPt, Float_t LepEta, Int_t LepID, Int_t LepZ1ID);
@@ -405,7 +417,7 @@ class HZZ4lNtupleMaker : public edm::EDAnalyzer {
 
   edm::EDGetTokenT<edm::MergeableCounter> preSkimToken;
 
-  PUReweight reweight;
+  PileUpWeight pileUpReweight;
 
   //counters
   Float_t Nevt_Gen;
@@ -474,7 +486,7 @@ HZZ4lNtupleMaker::HZZ4lNtupleMaker(const edm::ParameterSet& pset) :
   apply_K_NNLOQCD_ZZQQB(pset.getParameter<bool>("Apply_K_NNLOQCD_ZZQQB")),
   apply_K_NLOEW_ZZQQB(pset.getParameter<bool>("Apply_K_NLOEW_ZZQQB")),
 
-  reweight(),
+  pileUpReweight(myHelper.sampleType(), myHelper.setup()),
   sampleName(pset.getParameter<string>("sampleName")),
   hTH2D_Mu_All(0),
   hTH2F_El_Reco(0),
@@ -692,7 +704,9 @@ void HZZ4lNtupleMaker::analyze(const edm::Event& event, const edm::EventSetup& e
     }
 
     // get PU weight
-    PUWeight = reweight.weight(myHelper.sampleType(), myHelper.setup(), NTrueInt);
+    PUWeight = pileUpReweight.weight(NTrueInt);
+    PUWeight_Up = pileUpReweight.weight(NTrueInt, PileUpWeight::PUvar::VARUP);
+    PUWeight_Dn = pileUpReweight.weight(NTrueInt, PileUpWeight::PUvar::VARDOWN);
 
     event.getByToken(genParticleToken, genParticles);
     event.getByToken(genInfoToken, genInfo);
@@ -1666,7 +1680,11 @@ Float_t HZZ4lNtupleMaker::getAllWeight(const reco::Candidate* Lep) const
 
     Float_t lookup_pT = 50.;  // FIXME: the histogram contains 1 bin only, and overflows/underflows are intended to be included (?)
 
-    weight *= hTH2F_El_Reco->GetBinContent(hTH2F_El_Reco->GetXaxis()->FindBin(SCeta),hTH2F_El_Reco->GetYaxis()->FindBin(lookup_pT));
+    Float_t RecoSF = hTH2F_El_Reco->GetBinContent(hTH2F_El_Reco->GetXaxis()->FindBin(SCeta),hTH2F_El_Reco->GetYaxis()->FindBin(lookup_pT));
+    //RecoSF_Unc = = hTH2F_El_Reco->GetBinError(hTH2F_El_Reco->GetXaxis()->FindBin(SCeta),hTH2F_El_Reco->GetYaxis()->FindBin(lookup_pT));
+
+    //LepRecoSF.push_back(RecoSF)
+    weight *= RecoSF;
 
     if(mySIP >= 4.0 ) { // FIXME: use a better way to find RSE electrons!
         // No SF for RSE yet
@@ -2004,6 +2022,8 @@ void HZZ4lNtupleMaker::BookAllBranches(){
     myTree->Book("genProcessId", genProcessId, failedTreeLevel >= minimalFailedTree);
     myTree->Book("genHEPMCweight", genHEPMCweight, failedTreeLevel >= minimalFailedTree);
     myTree->Book("PUWeight", PUWeight, failedTreeLevel >= minimalFailedTree);
+    myTree->Book("PUWeight_Dn", PUWeight_Dn, failedTreeLevel >= minimalFailedTree);
+    myTree->Book("PUWeight_Up", PUWeight_Up, failedTreeLevel >= minimalFailedTree);
     myTree->Book("dataMCWeight", dataMCWeight, false);
     myTree->Book("trigEffWeight", trigEffWeight, false);
     myTree->Book("overallEventWeight", overallEventWeight, false);