User Tools

Site Tools


users:ngadiuba:tools:puweightrecipe:puweighth

Differences

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

Link to this comparison view

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