#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;
 }