#include "include/PUWeight.h" #include "TFile.h" #include <iostream> //============================================================================================== // Get MC pile-up scenario from string representation PUWeight::Scenario PUWeight::toScenario(const std::string& str) { PUWeight::Scenario sc = Summer12S10; if( str == "PUS10" ) sc = Summer12S10; else { std::cerr << "\n\nERROR unknown scenario '" << str << "'" << std::endl; throw std::exception(); } return sc; } //============================================================================================== // MC pile-up scenario to string representation std::string PUWeight::toString(const PUWeight::Scenario sc) { std::string str; if( sc == Summer12S10 ) str = "PUS10"; else { std::cerr << "\n\nERROR unknown scenario '" << sc << "'" << std::endl; throw std::exception(); } return str; } //============================================================================================== // Constructor. Initializes default behaviour to return PU weight of 1 PUWeight::PUWeight() : isInit_(false), nPUMax_(0) {} //============================================================================================== // Initialise weights for a given MC pile-up scenario. Can only be // called once. void PUWeight::initPUWeights(const std::string& nameOfDataDistribution, const PUWeight::Scenario sc) { if( isInit_ ) { std::cerr << "\n\nERROR in PUWeight: weights already initialised" << std::endl; throw std::exception(); } // Get data distribution from file TFile file(nameOfDataDistribution.c_str(), "READ"); TH1* h = NULL; file.GetObject("pileup",h); if( h == NULL ) { std::cerr << "\n\nERROR in PUWeight: Histogram 'pileup' does not exist in file '" << nameOfDataDistribution << "'\n."; throw std::exception(); } h->SetDirectory(0); file.Close(); // Computing weights puWeigths_ = generateWeights(sc,h); nPUMax_ = puWeigths_.size(); // Clean up delete h; isInit_ = true; } //============================================================================================== // Get weight factor dependent on number of added PU interactions double PUWeight::getPUWeight(const int nPU) const { double w = 1.; if( isInit_ ) { if( nPU >= nPUMax_ ) { //std::cerr << "WARNING: Number of PU vertices = " << nPU << " out of histogram binning." << std::endl; // In case nPU is out-of data-profile binning, // use weight from last bin w = puWeigths_.back(); } else { w = puWeigths_.at(nPU); } } return w; } //============================================================================================== // Generate weights for given data PU distribution // Scenarios from: https://twiki.cern.ch/twiki/bin/view/CMS/Pileup_MC_Gen_Scenarios // Code adapted from: https://twiki.cern.ch/twiki/bin/viewauth/CMS/PileupReweighting std::vector<double> PUWeight::generateWeights(const PUWeight::Scenario sc, const TH1* data_npu_estimated) const { // Store probabilites for each pu bin unsigned int nPUMax = 0; double *npuProbs = 0; if( sc == Summer12S10 ) { nPUMax = 60; double npuSummer12_S10[60] = { 2.560E-06, 5.239E-06, 1.420E-05, 5.005E-05, 1.001E-04, 2.705E-04, 1.999E-03, 6.097E-03, 1.046E-02, 1.383E-02, 1.685E-02, 2.055E-02, 2.572E-02, 3.262E-02, 4.121E-02, 4.977E-02, 5.539E-02, 5.725E-02, 5.607E-02, 5.312E-02, 5.008E-02, 4.763E-02, 4.558E-02, 4.363E-02, 4.159E-02, 3.933E-02, 3.681E-02, 3.406E-02, 3.116E-02, 2.818E-02, 2.519E-02, 2.226E-02, 1.946E-02, 1.682E-02, 1.437E-02, 1.215E-02, 1.016E-02, 8.400E-03, 6.873E-03, 5.564E-03, 4.457E-03, 3.533E-03, 2.772E-03, 2.154E-03, 1.656E-03, 1.261E-03, 9.513E-04, 7.107E-04, 5.259E-04, 3.856E-04, 2.801E-04, 2.017E-04, 1.439E-04, 1.017E-04, 7.126E-05, 4.948E-05, 3.405E-05, 2.322E-05, 1.570E-05, 5.005E-06}; npuProbs = npuSummer12_S10; } // Check that binning of data-profile matches MC scenario if( nPUMax != static_cast<unsigned int>(data_npu_estimated->GetNbinsX()) ) { std::cerr << "\n\nERROR number of bins (" << data_npu_estimated->GetNbinsX() << ") in data PU-profile does not match number of bins (" << nPUMax << ") in MC scenario " << toString(sc) << std::endl; throw std::exception(); } std::vector<double> result(nPUMax,0.); double s = 0.; for(unsigned int npu = 0; npu < nPUMax; ++npu) { const double npu_estimated = data_npu_estimated->GetBinContent(data_npu_estimated->GetXaxis()->FindBin(npu)); result[npu] = npu_estimated / npuProbs[npu]; s += npu_estimated; } // normalize weights such that the total sum of weights over thw whole sample is 1.0, i.e., sum_i result[i] * npu_probs[i] should be 1.0 (!) for(unsigned int npu = 0; npu < nPUMax; ++npu) { result[npu] /= s; } return result; }