This is an old revision of the document!
#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;
}