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;

}