[DijetAnalysis] Compute weights for PU reweighting as a function of H…

…LT path.

DijetAnalysis/data/pileup/2016/byHLT/pileupCalc_ratio.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+
+for t in 60 80 140 200 260 320 400 500; do
+    ../../../../scripts/pu_reweighting.py datamc_pu_weights \
+      -d "nPUMean_data_HLT_DiPFJetAve${t}.root" \
+      -m "../nPUMean_mc.root" \
+      -o "nPUMean_ratio_HLT_DiPFJetAve${t}.root"
+done

DijetAnalysis/interface/NtupleProducer.h

@@ -256,6 +256,8 @@ namespace dijet {
         std::unique_ptr<karma::PileupWeightProvider> m_puWeightProvider;
         std::unique_ptr<karma::PileupWeightProvider> m_puWeightProviderAlt;
 
+        std::vector<karma::PileupWeightProvider*> m_puWeightProvidersByHLT;
+
         // -- handles and tokens
         typename edm::Handle<karma::Event> karmaEventHandle;
         edm::EDGetTokenT<karma::Event> karmaEventToken;

DijetAnalysis/python/Configuration/dijetAnalysis_94X_Run2016_17Jul2018_cff.py

@@ -491,6 +491,10 @@ def setup_pipeline(process, options, pipeline_name, jet_algo_name, jec_shift=Non
                 os.getenv('CMSSW_BASE'),
                 YEAR="2016",
             ),
+            pileupWeightByHLTFileBasename = "{}/src/Karma/DijetAnalysis/data/pileup/{YEAR}/byHLT/nPUMean_ratio".format(
+                os.getenv('CMSSW_BASE'),
+                YEAR="2016",
+            ),
             pileupWeightHistogramName = "pileup",
 
             # jet ID (for event-based jet ID in PostProcessing using branches 'jet1id', 'jet2id')

DijetAnalysis/python/NtupleProducer_cfi.py

@@ -29,6 +29,8 @@
         pileupWeightFileAlt = cms.string("nPUMean_ratio_alternative.root"),
         pileupWeightHistogramName = cms.string("pileup"),
 
+        pileupWeightByHLTFileBasename = cms.string("nPUMean_ratio"),
+
         weightForStitching = cms.double(1.0),
 
         # YAML files specifying analysis binning

DijetAnalysis/src/NtupleFlatOutput.cc

@@ -123,6 +123,7 @@ void dijet::NtupleFlatOutput::setUpTTree(TTree* tree, dijet::NtupleEntry* produc
         ADD_BRANCH(tree, productForFill, weightForStitching, D);
         ADD_BRANCH(tree, productForFill, pileupWeight, D);
         ADD_BRANCH(tree, productForFill, pileupWeightAlt, D);
+        ADD_BRANCH(tree, productForFill, pileupWeightByActiveHLT, D);
         // binning values
         ADD_BRANCH(tree, productForFill, binningValue, D);
     }

DijetAnalysis/src/NtupleProducer.cc

@@ -5,6 +5,8 @@
 
 #include "Karma/DijetAnalysis/interface/NtupleProducer.h"
 
+#include "boost/filesystem/path.hpp"
+#include "boost/filesystem/operations.hpp"
 
 // -- constructor
 dijet::NtupleProducer::NtupleProducer(const edm::ParameterSet& config, const dijet::NtupleProducerGlobalCache* globalCache) : m_configPSet(config) {
@@ -56,6 +58,25 @@ dijet::NtupleProducer::NtupleProducer(const edm::ParameterSet& config, const dij
                 )
             );
         }
+        // can provide pileup weight files for each HLT path
+        auto pileupWeightByHLTFileBasename = m_configPSet.getParameter<std::string>("pileupWeightByHLTFileBasename");
+        if (!pileupWeightByHLTFileBasename.empty()) {
+            m_puWeightProvidersByHLT.resize(globalCache->hltPaths_.size());
+            for (size_t iHLTPath = 0; iHLTPath < globalCache->hltPaths_.size(); ++iHLTPath) {
+                std::string pileupWeightFileName = pileupWeightByHLTFileBasename + "_" + globalCache->hltPaths_.at(iHLTPath) + ".root";
+
+                if (!boost::filesystem::exists(pileupWeightFileName)) {
+                    std::cout << "No HLT-dependent pileup weight information found for trigger path: " << globalCache->hltPaths_.at(iHLTPath) << std::endl;
+                    continue;
+                }
+
+                std::cout << "Reading HLT-dependent pileup weight information from file: " << pileupWeightFileName << std::endl;
+                m_puWeightProvidersByHLT[iHLTPath] = new karma::PileupWeightProvider(
+                    pileupWeightByHLTFileBasename + "_" + globalCache->hltPaths_[iHLTPath] + ".root",
+                    m_configPSet.getParameter<std::string>("pileupWeightHistogramName")
+                );
+            }
+        }
     }
 
 
@@ -267,9 +288,11 @@ void dijet::NtupleProducer::produce(edm::Event& event, const edm::EventSetup& se
     outputNtupleEntry->weightForStitching = m_weightForStitching;
 
     // prefiring weights
-    outputNtupleEntry->prefiringWeight = *(this->karmaPrefiringWeightHandle);
-    outputNtupleEntry->prefiringWeightUp = *(this->karmaPrefiringWeightUpHandle);
-    outputNtupleEntry->prefiringWeightDown = *(this->karmaPrefiringWeightDownHandle);
+    if (m_isData) {
+        outputNtupleEntry->prefiringWeight = *(this->karmaPrefiringWeightHandle);
+        outputNtupleEntry->prefiringWeightUp = *(this->karmaPrefiringWeightUpHandle);
+        outputNtupleEntry->prefiringWeightDown = *(this->karmaPrefiringWeightDownHandle);
+    }
 
     // -- trigger results
     dijet::TriggerBits bitsetHLTBits;
@@ -509,6 +532,19 @@ void dijet::NtupleProducer::produce(edm::Event& event, const edm::EventSetup& se
     outputNtupleEntry->metRaw = this->karmaMETCollectionHandle->at(0).uncorP4.Pt();
     outputNtupleEntry->sumEtRaw = this->karmaMETCollectionHandle->at(0).uncorSumEt;
 
+    // determine PU weight as a function of the active HLT path
+    if (!m_isData) {
+        if (outputNtupleEntry->indexActiveTriggerPathJet12PtAve >= 0) {
+            auto* puWeightByHLTProvider = m_puWeightProvidersByHLT.at(outputNtupleEntry->indexActiveTriggerPathJet12PtAve);
+            if (puWeightByHLTProvider) {
+                outputNtupleEntry->pileupWeightByActiveHLT = puWeightByHLTProvider->getPileupWeight(outputNtupleEntry->nPUMean);
+            }
+        }
+        else {
+            outputNtupleEntry->pileupWeightByActiveHLT = -1.0;
+        }
+    }
+
     // move outputs to event tree
     event.put(std::move(outputNtupleEntry));
 }

DijetAnalysisFormats/interface/Ntuple.h

@@ -120,6 +120,7 @@ namespace dijet {
         double weightForStitching = 1.0;
         double pileupWeight =  1.0;
         double pileupWeightAlt =  1.0;
+        double pileupWeightByActiveHLT = 1.0;
 
         // binning values
         double binningValue = UNDEFINED_DOUBLE;