UZH

Physik-Institut

CMS Wiki Pages

User Tools

Site Tools

Trace: • puweighth
users:ngadiuba:tools:puweightrecipe:puweighth

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision		Previous revision
users:ngadiuba:tools:puweightrecipe:puweighth [2014/04/15 14:27] – created ngadiubausers:ngadiuba:tools:puweightrecipe:puweighth [2014/04/15 14:34] (current) – ngadiuba
Line 1: Line 1:
-#include "include/PUWeight.h" +    #ifndef PUWeight_H 
- +    #define PUWeight_H 
-#include "TFile.h" +    
- +    #include <string> 
-#include <iostream> +    #include <vector> 
- +    
-//============================================================================================== +    #include "TH1.h" 
-// Get MC pile-up scenario from string representation +    #include "TString.h" 
-PUWeight::Scenario PUWeight::toScenario(const std::string& str) { +    
-  PUWeight::Scenario sc = Summer12S10; +    // Compute pile-up weights to match data distribtution 
-  if( str == "PUS10" ) sc = Summer12S10; +    class PUWeight { 
-  else { +    public: 
-    std::cerr << "\n\nERROR unknown scenario '" << str << "'" << std::endl; +      enum Scenario { Summer12S10 }; 
-    throw std::exception(); +    
-  } +      static Scenario toScenario(const std::string& str); 
- +      static std::string toString(const Scenario sc); 
-  return sc; +    
-} +      // Constructor. Initializes default behaviour to return PU weight of 1 
- +      PUWeight(); 
-//============================================================================================== +    
-// MC pile-up scenario to string representation +      // Initialise weights for a given MC pile-up scenario. Can only be 
-std::string PUWeight::toString(const PUWeight::Scenario sc) { +      // called once. 
-  std::string str; +      void initPUWeights(const TString& nameOfDataDistribution, const PUWeight::Scenario sc) { 
-  if( sc == Summer12S10 ) str = "PUS10"; +     initPUWeights( std::string( nameOfDataDistribution.Data() ), sc ); 
-  else { +      } 
-    std::cerr << "\n\nERROR unknown scenario '" << sc << "'" << std::endl; +      void initPUWeights(const std::string& nameOfDataDistribution, const PUWeight::Scenario sc); 
-    throw std::exception(); +    
-  } +      // Get weight factor dependent on number of added PU interactions 
- +      // nPU is the 'true' number of interactions as explained in 
-  return str; +      // https://twiki.cern.ch/twiki/bin/viewauth/CMS/PileupJSONFileforData#2012_Pileup_JSON_Files 
-} +      double getPUWeight(const int nPU) const; 
- +    
-//============================================================================================== +    
-// Constructor. Initializes default behaviour to return PU weight of 1 +    private: 
-PUWeight::PUWeight() +      bool isInit_; 
-  : isInit_(false), nPUMax_(0) {} +      int nPUMax_; 
- +      std::vector<double> puWeigths_; 
-//============================================================================================== +    
-// Initialise weights for a given MC pile-up scenario. Can only be +      // Generate weights for given data PU distribution 
-// called once. +      std::vector<double> generateWeights(const Scenario sc, const TH1* data_npu_estimated) const; 
-void PUWeight::initPUWeights(const std::string& nameOfDataDistribution, const PUWeight::Scenario sc) { +    }; 
- +    
-  if( isInit_ ) { +    #endif //PUWeight_H
-    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; +
-} +
- +
users/ngadiuba/tools/puweightrecipe/puweighth.1397564868.txt.gz · Last modified: by ngadiuba

Page Tools

Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 4.0 International
CC Attribution-Share Alike 4.0 International Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki