LumiReweightingStandAlone.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_interface_LumiReWeighting_h
#define PhysicsTools_Utilities_interface_LumiReWeighting_h

#include "TH1F.h"
#include "TH3.h"
#include "TFile.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include <string>
#include <vector>
#include <cmath>
#include <algorithm>
#include <iostream>

namespace reweight {

  // add a class to shift the mean of a poisson-like luminosity distribution by an arbitrary amount.
  // Only valid for small (<1.5) shifts about the 2011 lumi distribution for now, because shifts are non-linear
  // Each PoissonMeanShifter does one shift, so defining multiples can give you an arbitrary collection

  class PoissonMeanShifter {
  public:
    PoissonMeanShifter(){};

    PoissonMeanShifter(float Shift) {
      // these are the polynomial or exponential coefficients for each bin of a 25-bin sequence that
      // convert the Distribution of the 2011 luminosity to something with a lower or higher peak luminosity.
      // The distributions aren't quite poisson because they model luminosity decreasing during a fill. This implies that
      // they do get wider as the mean increases, so the weights are not linear with increasing mean.

      static const double p0_minus[20] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
                                          1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
      static const double p1_minus[20] = {-0.677786, -0.619614, -0.49465, -0.357963, -0.238359, -0.110002, 0.0348629,
                                          0.191263,  0.347648,  0.516615, 0.679646,  0.836673,  0.97764,   1.135,
                                          1.29922,   1.42467,   1.55901,  1.61762,   1.67275,   1.96008};
      static const double p2_minus[20] = {
          0.526164, 0.251816, 0.11049,  0.026917, -0.0464692, -0.087022, -0.0931581, -0.0714295, -0.0331772, 0.0347473,
          0.108658, 0.193048, 0.272314, 0.376357, 0.4964,     0.58854,   0.684959,   0.731063,   0.760044,   1.02386};

      static const double p1_expoM[5] = {1.63363e-03, 6.79290e-04, 3.69900e-04, 2.24349e-04, 9.87156e-06};

      static const double p2_expoM[5] = {2.64692, 3.26585, 3.53229, 4.18035, 5.64027};

      static const double p0_plus[20] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
                                         1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
      static const double p1_plus[20] = {-0.739059, -0.594445, -0.477276, -0.359707, -0.233573, -0.103458, 0.0373401,
                                         0.176571,  0.337617,  0.499074,  0.675126,  0.840522,  1.00917,   1.15847,
                                         1.23816,   1.44271,   1.52982,   1.46385,   1.5802,    0.988689};
      static const double p2_plus[20] = {0.208068,    0.130033,   0.0850356,  0.0448344,  0.000749832,
                                         -0.0331347,  -0.0653281, -0.0746009, -0.0800667, -0.0527636,
                                         -0.00402649, 0.103338,   0.261261,   0.491084,   0.857966,
                                         1.19495,     1.75071,    2.65559,    3.35433,    5.48835};

      static const double p1_expoP[5] = {1.42463e-01, 4.18966e-02, 1.12697e-01, 1.66197e-01, 1.50768e-01};

      static const double p2_expoP[5] = {1.98758, 2.27217, 2.26799, 2.38455, 2.52428};

      // initialize weights based on desired Shift

      for (int ibin = 0; ibin < 20; ibin++) {
        if (Shift < .0) {
          Pweight_[ibin] = p0_minus[ibin] + p1_minus[ibin] * Shift + p2_minus[ibin] * Shift * Shift;
        } else {
          Pweight_[ibin] = p0_plus[ibin] + p1_plus[ibin] * Shift + p2_plus[ibin] * Shift * Shift;
        }
      }

      // last few bins fit better to an exponential...

      for (int ibin = 20; ibin < 25; ibin++) {
        if (Shift < 0.) {
          Pweight_[ibin] = p1_expoM[ibin - 20] * exp(p2_expoM[ibin - 20] * Shift);
        } else {
          Pweight_[ibin] = p1_expoP[ibin - 20] * exp(p2_expoP[ibin - 20] * Shift);
        }
      }
    };

    double ShiftWeight(int ibin) {
      if (ibin < 25 && ibin >= 0) {
        return Pweight_[ibin];
      } else {
        return 0;
      }
    };

    double ShiftWeight(float pvnum) {
      int ibin = int(pvnum);

      if (ibin < 25 && ibin >= 0) {
        return Pweight_[ibin];
      } else {
        return 0;
      }
    };

  private:
    double Pweight_[25];
  };

  class LumiReWeighting {
  public:
    LumiReWeighting(){};

    LumiReWeighting(std::string generatedFile, std::string dataFile, std::string GenHistName, std::string DataHistName)
        : generatedFileName_(generatedFile),
          dataFileName_(dataFile),
          GenHistName_(GenHistName),
          DataHistName_(DataHistName) {
      generatedFile_ = new TFile(generatedFileName_.c_str());  //MC distribution
      dataFile_ = new TFile(dataFileName_.c_str());            //Data distribution

      Data_distr_ = new TH1F(*(static_cast<TH1F*>(dataFile_->Get(DataHistName_.c_str())->Clone())));
      MC_distr_ = new TH1F(*(static_cast<TH1F*>(generatedFile_->Get(GenHistName_.c_str())->Clone())));

      // normalize both histograms first

      Data_distr_->Scale(1.0 / Data_distr_->Integral());
      MC_distr_->Scale(1.0 / MC_distr_->Integral());

      weights_ = new TH1F(*(Data_distr_));

      // MC * data/MC = data, so the weights are data/MC:

      weights_->SetName("lumiWeights");

      TH1F* den = new TH1F(*(MC_distr_));

      weights_->Divide(den);  // so now the average weight should be 1.0

      std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;

      int NBins = weights_->GetNbinsX();

      for (int ibin = 1; ibin < NBins + 1; ++ibin) {
        std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
      }

      weightOOT_init();

      FirstWarning_ = true;
    }

    LumiReWeighting(const std::vector<float>& MC_distr, const std::vector<float>& Lumi_distr) {
      // no histograms for input: use vectors

      // now, make histograms out of them:

      // first, check they are the same size...

      if (MC_distr.size() != Lumi_distr.size()) {
        std::cerr << "ERROR: LumiReWeighting: input vectors have different sizes. Quitting... \n";
        return;
      }

      Int_t NBins = MC_distr.size();

      MC_distr_ = new TH1F("MC_distr", "MC dist", NBins, -0.5, float(NBins) - 0.5);
      Data_distr_ = new TH1F("Data_distr", "Data dist", NBins, -0.5, float(NBins) - 0.5);

      weights_ = new TH1F("luminumer", "luminumer", NBins, -0.5, float(NBins) - 0.5);
      TH1F* den = new TH1F("lumidenom", "lumidenom", NBins, -0.5, float(NBins) - 0.5);

      for (int ibin = 1; ibin < NBins + 1; ++ibin) {
        weights_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
        Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
        den->SetBinContent(ibin, MC_distr[ibin - 1]);
        MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
      }

      // check integrals, make sure things are normalized

      float deltaH = weights_->Integral();
      if (fabs(1.0 - deltaH) > 0.02) {  //*OOPS*...
        weights_->Scale(1.0 / deltaH);
        Data_distr_->Scale(1.0 / deltaH);
      }
      float deltaMC = den->Integral();
      if (fabs(1.0 - deltaMC) > 0.02) {
        den->Scale(1.0 / deltaMC);
        MC_distr_->Scale(1.0 / deltaMC);
      }

      weights_->Divide(den);  // so now the average weight should be 1.0

      std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;

      for (int ibin = 1; ibin < NBins + 1; ++ibin) {
        std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
      }

      weightOOT_init();

      FirstWarning_ = true;
    }

    void weight3D_init(float ScaleFactor, std::string WeightOutputFile = "") {
      //create histogram to write output weights, save pain of generating them again...

      TH3D* WHist = new TH3D("WHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
      TH3D* DHist = new TH3D("DHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
      TH3D* MHist = new TH3D("MHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);

      using std::min;

      if (MC_distr_->GetEntries() == 0) {
        std::cout << " MC and Data distributions are not initialized! You must call the LumiReWeighting constructor. "
                  << std::endl;
      }

      // arrays for storing number of interactions

      double MC_ints[50][50][50];
      double Data_ints[50][50][50];

      for (int i = 0; i < 50; i++) {
        for (int j = 0; j < 50; j++) {
          for (int k = 0; k < 50; k++) {
            MC_ints[i][j][k] = 0.;
            Data_ints[i][j][k] = 0.;
          }
        }
      }

      double factorial[50];
      double PowerSer[50];
      double base = 1.;

      factorial[0] = 1.;
      PowerSer[0] = 1.;

      for (int i = 1; i < 50; ++i) {
        base = base * float(i);
        factorial[i] = base;
      }

      double x;
      double xweight;
      double probi, probj, probk;
      double Expval, mean;
      int xi;

      // Get entries for Data, MC, fill arrays:
      int NMCbin = MC_distr_->GetNbinsX();

      for (int jbin = 1; jbin < NMCbin + 1; jbin++) {
        x = MC_distr_->GetBinCenter(jbin);
        xweight = MC_distr_->GetBinContent(jbin);  //use as weight for matrix

        //for Summer 11, we have this int feature:
        xi = int(x);

        // Generate Poisson distribution for each value of the mean
        mean = double(xi);

        if (mean < 0.) {
          std::cout << "LumiReweighting:BadInputValue"
                    << " Your histogram generates MC luminosity values less than zero!"
                    << " Please Check.  Terminating." << std::endl;
        }

        if (mean == 0.) {
          Expval = 1.;
        } else {
          Expval = exp(-1. * mean);
        }

        base = 1.;

        for (int i = 1; i < 50; ++i) {
          base = base * mean;
          PowerSer[i] = base;  // PowerSer is mean^i
        }

        // compute poisson probability for each Nvtx in weight matrix
        for (int i = 0; i < 50; i++) {
          probi = PowerSer[i] / factorial[i] * Expval;
          for (int j = 0; j < 50; j++) {
            probj = PowerSer[j] / factorial[j] * Expval;
            for (int k = 0; k < 50; k++) {
              probk = PowerSer[k] / factorial[k] * Expval;
              // joint probability is product of event weights multiplied by weight of input distribution bin
              MC_ints[i][j][k] = MC_ints[i][j][k] + probi * probj * probk * xweight;
            }
          }
        }
      }

      int NDatabin = Data_distr_->GetNbinsX();

      for (int jbin = 1; jbin < NDatabin + 1; jbin++) {
        mean = (Data_distr_->GetBinCenter(jbin)) * ScaleFactor;
        xweight = Data_distr_->GetBinContent(jbin);

        // Generate poisson distribution for each value of the mean
        if (mean < 0.) {
          std::cout << "LumiReweighting:BadInputValue"
                    << " Your histogram generates MC luminosity values less than zero!"
                    << " Please Check.  Terminating." << std::endl;
        }

        if (mean == 0.) {
          Expval = 1.;
        } else {
          Expval = exp(-1. * mean);
        }

        base = 1.;

        for (int i = 1; i < 50; ++i) {
          base = base * mean;
          PowerSer[i] = base;
        }

        // compute poisson probability for each Nvtx in weight matrix

        for (int i = 0; i < 50; i++) {
          probi = PowerSer[i] / factorial[i] * Expval;
          for (int j = 0; j < 50; j++) {
            probj = PowerSer[j] / factorial[j] * Expval;
            for (int k = 0; k < 50; k++) {
              probk = PowerSer[k] / factorial[k] * Expval;
              // joint probability is product of event weights multiplied by weight of input distribution bin
              Data_ints[i][j][k] = Data_ints[i][j][k] + probi * probj * probk * xweight;
            }
          }
        }
      }

      for (int i = 0; i < 50; i++) {
        //if(i<5) std::cout << "i = " << i << std::endl;
        for (int j = 0; j < 50; j++) {
          for (int k = 0; k < 50; k++) {
            if ((MC_ints[i][j][k]) > 0.) {
              Weight3D_[i][j][k] = Data_ints[i][j][k] / MC_ints[i][j][k];
            } else {
              Weight3D_[i][j][k] = 0.;
            }
            WHist->SetBinContent(i + 1, j + 1, k + 1, Weight3D_[i][j][k]);
            DHist->SetBinContent(i + 1, j + 1, k + 1, Data_ints[i][j][k]);
            MHist->SetBinContent(i + 1, j + 1, k + 1, MC_ints[i][j][k]);
            //      if(i<5 && j<5 && k<5) std::cout << Weight3D_[i][j][k] << " " ;
          }
          //      if(i<5 && j<5) std::cout << std::endl;
        }
      }

      if (!WeightOutputFile.empty()) {
        std::cout << " 3D Weight Matrix initialized! " << std::endl;
        std::cout << " Writing weights to file " << WeightOutputFile << " for re-use...  " << std::endl;

        TFile* outfile = new TFile(WeightOutputFile.c_str(), "RECREATE");
        WHist->Write();
        MHist->Write();
        DHist->Write();
        outfile->Write();
        outfile->Close();
        outfile->Delete();
      }

      return;
    }

    void weight3D_set(std::string WeightFileName) {
      TFile* infile = new TFile(WeightFileName.c_str());
      TH1F* WHist = (TH1F*)infile->Get("WHist");

      // Check if the histogram exists
      if (!WHist) {
        std::cout << " Could not find the histogram WHist in the file "
                  << "in the file " << WeightFileName << "." << std::endl;
        return;
      }

      for (int i = 0; i < 50; i++) {
        for (int j = 0; j < 50; j++) {
          for (int k = 0; k < 50; k++) {
            Weight3D_[i][j][k] = WHist->GetBinContent(i, j, k);
          }
        }
      }

      std::cout << " 3D Weight Matrix initialized! " << std::endl;

      return;
    }

    void weightOOT_init() {
      // The following are poisson distributions with different means, where the maximum
      // of the function has been normalized to weight 1.0
      // These are used to reweight the out-of-time pileup to match the in-time distribution.
      // The total event weight is the product of the in-time weight, the out-of-time weight,
      // and a residual correction to fix the distortions caused by the fact that the out-of-time
      // distribution is not flat.

      static const double weight_24[25] = {0,           0,           0,           0,          2.46277e-06,
                                           2.95532e-05, 0.000104668, 0.000401431, 0.00130034, 0.00342202,
                                           0.00818132,  0.0175534,   0.035784,    0.0650836,  0.112232,
                                           0.178699,    0.268934,    0.380868,    0.507505,   0.640922,
                                           0.768551,    0.877829,    0.958624,    0.99939,    1};

      static const double weight_23[25] = {0,           1.20628e-06, 1.20628e-06, 2.41255e-06, 1.20628e-05,
                                           6.39326e-05, 0.000252112, 0.000862487, 0.00244995,  0.00616527,
                                           0.0140821,   0.0293342,   0.0564501,   0.100602,    0.164479,
                                           0.252659,    0.36268,     0.491427,    0.627979,    0.75918,
                                           0.873185,    0.957934,    0.999381,    1,           0.957738};

      static const double weight_22[25] = {
          0,         0,         0,         5.88636e-06, 3.0609e-05, 0.000143627, 0.000561558, 0.00173059, 0.00460078,
          0.0110616, 0.0238974, 0.0475406, 0.0875077,   0.148682,   0.235752,    0.343591,    0.473146,   0.611897,
          0.748345,  0.865978,  0.953199,  0.997848,    1,          0.954245,    0.873688};

      static const double weight_21[25] = {
          0,         0,         1.15381e-06, 8.07665e-06, 7.1536e-05, 0.000280375, 0.00107189, 0.00327104, 0.00809396,
          0.0190978, 0.0401894, 0.0761028,   0.13472,     0.216315,   0.324649,    0.455125,   0.598241,   0.739215,
          0.861866,  0.953911,  0.998918,    1,           0.956683,   0.872272,    0.76399};

      static const double weight_20[25] = {
          0,         0,         1.12532e-06, 2.58822e-05, 0.000145166, 0.000633552, 0.00215048, 0.00592816, 0.0145605,
          0.0328367, 0.0652649, 0.11893,     0.19803,     0.305525,    0.436588,    0.581566,   0.727048,   0.8534,
          0.949419,  0.999785,  1,           0.953008,    0.865689,    0.753288,    0.62765};
      static const double weight_19[25] = {
          0,         0,        1.20714e-05, 5.92596e-05, 0.000364337, 0.00124994, 0.00403953, 0.0108149, 0.025824,
          0.0544969, 0.103567, 0.17936,     0.283532,    0.416091,    0.562078,   0.714714,   0.846523,  0.947875,
          1,         0.999448, 0.951404,    0.859717,    0.742319,    0.613601,   0.48552};

      static const double weight_18[25] = {
          0,         3.20101e-06, 2.88091e-05, 0.000164319, 0.000719161, 0.00250106, 0.00773685, 0.0197513, 0.0443693,
          0.0885998, 0.159891,    0.262607,    0.392327,    0.543125,    0.69924,    0.837474,   0.943486,  0.998029,
          1,         0.945937,    0.851807,    0.729309,    0.596332,    0.467818,   0.350434};

      static const double weight_17[25] = {
          1.03634e-06, 7.25437e-06, 4.97443e-05, 0.000340956, 0.00148715, 0.00501485, 0.0143067, 0.034679, 0.0742009,
          0.140287,    0.238288,    0.369416,    0.521637,    0.682368,   0.828634,   0.939655,  1,        0.996829,
          0.94062,     0.841575,    0.716664,    0.582053,    0.449595,   0.331336,   0.234332};

      static const double weight_16[25] = {
          4.03159e-06, 2.41895e-05, 0.000141106, 0.00081942, 0.00314565, 0.00990662, 0.026293, 0.0603881, 0.120973,
          0.214532,    0.343708,    0.501141,    0.665978,   0.820107,   0.938149,   1,        0.99941,   0.940768,
          0.837813,    0.703086,    0.564023,    0.42928,    0.312515,   0.216251,   0.14561};

      static const double weight_15[25] = {
          9.76084e-07, 5.07564e-05, 0.000303562, 0.00174036, 0.00617959, 0.0188579, 0.047465, 0.101656, 0.189492,
          0.315673,    0.474383,    0.646828,    0.809462,   0.934107,   0.998874,  1,        0.936163, 0.827473,
          0.689675,    0.544384,    0.40907,     0.290648,   0.198861,   0.12951,   0.0808051};

      static const double weight_14[25] = {
          1.13288e-05, 0.000124617, 0.000753365, 0.00345056, 0.0123909, 0.0352712, 0.0825463, 0.16413,  0.287213,
          0.44615,     0.625826,    0.796365,    0.930624,   0.999958,  1,         0.934414,  0.816456, 0.672939,
          0.523033,    0.386068,    0.269824,    0.180342,   0.114669,  0.0698288, 0.0406496};

      static const double weight_13[25] = {
          2.54296e-05, 0.000261561, 0.00167018, 0.00748083, 0.0241308, 0.0636801, 0.138222, 0.255814, 0.414275,
          0.600244,    0.779958,    0.92256,    0.999155,   1,         0.927126,  0.804504, 0.651803, 0.497534,
          0.35976,     0.245834,    0.160904,   0.0991589,  0.0585434, 0.0332437, 0.0180159};

      static const double weight_12[25] = {
          5.85742e-05, 0.000627706, 0.00386677, 0.0154068, 0.0465892, 0.111683,  0.222487,  0.381677, 0.5719,
          0.765001,    0.915916,    1,          0.999717,  0.921443,  0.791958,  0.632344,  0.475195, 0.334982,
          0.223666,    0.141781,    0.0851538,  0.048433,  0.0263287, 0.0133969, 0.00696683};

      static const double weight_11[25] = {
          0.00015238, 0.00156064, 0.00846044, 0.0310939, 0.0856225, 0.187589,   0.343579,  0.541892, 0.74224,
          0.909269,   0.998711,   1,          0.916889,  0.77485,   0.608819,   0.447016,  0.307375, 0.198444,
          0.121208,   0.070222,   0.0386492,  0.0201108, 0.0100922, 0.00484937, 0.00222458};

      static const double weight_10[25] = {
          0.000393044, 0.00367001, 0.0179474, 0.060389,   0.151477,   0.302077,   0.503113,   0.720373, 0.899568,
          1,           0.997739,   0.909409,  0.75728,    0.582031,   0.415322,   0.277663,   0.174147, 0.102154,
          0.0566719,   0.0298642,  0.0147751, 0.00710995, 0.00319628, 0.00140601, 0.000568796};

      static const double weight_9[25] = {
          0.00093396, 0.00854448, 0.0380306,  0.113181,  0.256614,    0.460894,    0.690242,   0.888781,  1,
          0.998756,   0.899872,   0.735642,   0.552532,  0.382726,    0.246114,    0.147497,   0.0825541, 0.0441199,
          0.0218157,  0.0103578,  0.00462959, 0.0019142, 0.000771598, 0.000295893, 0.000111529};

      static const double weight_8[25] = {
          0.00240233, 0.0192688, 0.0768653,  0.205008,    0.410958,    0.65758,     0.875657,   0.999886,  1,
          0.889476,   0.711446,  0.517781,   0.345774,    0.212028,    0.121208,    0.0644629,  0.0324928, 0.0152492,
          0.00673527, 0.0028547, 0.00117213, 0.000440177, 0.000168471, 5.80689e-05, 1.93563e-05};

      static const double weight_7[25] = {0.00617233,  0.0428714,   0.150018,    0.350317,    0.612535,
                                          0.856525,    0.999923,    1,           0.87544,     0.679383,
                                          0.478345,    0.303378,    0.176923,    0.0950103,   0.0476253,
                                          0.0222211,   0.00972738,  0.00392962,  0.0015258,   0.000559168,
                                          0.000183928, 6.77983e-05, 1.67818e-05, 7.38398e-06, 6.71271e-07};

      static const double weight_6[25] = {0.0154465,   0.0923472,   0.277322,    0.55552,     0.833099,
                                          0.999035,    1,           0.855183,    0.641976,    0.428277,
                                          0.256804,    0.139798,    0.0700072,   0.0321586,   0.0137971,
                                          0.00544756,  0.00202316,  0.000766228, 0.000259348, 8.45836e-05,
                                          1.80362e-05, 8.70713e-06, 3.73163e-06, 6.21938e-07, 0};

      static const double weight_5[25] = {0.0382845,   0.191122,    0.478782,    0.797314,    1,
                                          0.997148,    0.831144,    0.59461,     0.371293,    0.205903,
                                          0.103102,    0.0471424,   0.0194997,   0.00749415,  0.00273709,
                                          0.000879189, 0.000286049, 0.000102364, 1.70606e-05, 3.98081e-06,
                                          2.27475e-06, 0,           0,           0,           0};

      static const double weight_4[25] = {0.0941305,   0.373824,    0.750094,    1,           0.997698,
                                          0.800956,    0.532306,    0.304597,    0.152207,    0.0676275,
                                          0.0270646,   0.00975365,  0.00326077,  0.00101071,  0.000301781,
                                          7.41664e-05, 1.58563e-05, 3.58045e-06, 1.02299e-06, 0,
                                          5.11493e-07, 0,           0,           0,           0};

      static const double weight_3[25] = {0.222714,    0.667015,    1,           0.999208,    0.750609,
                                          0.449854,    0.224968,    0.0965185,   0.0361225,   0.012084,
                                          0.00359618,  0.000977166, 0.000239269, 6.29422e-05, 1.16064e-05,
                                          1.78559e-06, 0,           4.46398e-07, 0,           0,
                                          0,           0,           0,           0,           0};

      static const double weight_2[25] = {
          0.499541,   0.999607,    1,           0.666607,    0.333301,    0.13279, 0.0441871, 0.0127455, 0.00318434,
          0.00071752, 0.000132204, 2.69578e-05, 5.16999e-06, 2.21571e-06, 0,       0,         0,         0,
          0,          0,           0,           0,           0,           0,       0};

      static const double weight_1[25] = {
          0.999165,    1,           0.499996,    0.166868, 0.0414266, 0.00831053, 0.00137472, 0.000198911, 2.66302e-05,
          2.44563e-06, 2.71737e-07, 2.71737e-07, 0,        0,         0,          0,          0,           0,
          0,           0,           0,           0,        0,         0,          0};

      static const double weight_0[25] = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

      //WeightOOTPU_ = {0};

      const double* WeightPtr = 0;

      for (int iint = 0; iint < 25; ++iint) {
        if (iint == 0)
          WeightPtr = weight_0;
        if (iint == 1)
          WeightPtr = weight_1;
        if (iint == 2)
          WeightPtr = weight_2;
        if (iint == 3)
          WeightPtr = weight_3;
        if (iint == 4)
          WeightPtr = weight_4;
        if (iint == 5)
          WeightPtr = weight_5;
        if (iint == 6)
          WeightPtr = weight_6;
        if (iint == 7)
          WeightPtr = weight_7;
        if (iint == 8)
          WeightPtr = weight_8;
        if (iint == 9)
          WeightPtr = weight_9;
        if (iint == 10)
          WeightPtr = weight_10;
        if (iint == 11)
          WeightPtr = weight_11;
        if (iint == 12)
          WeightPtr = weight_12;
        if (iint == 13)
          WeightPtr = weight_13;
        if (iint == 14)
          WeightPtr = weight_14;
        if (iint == 15)
          WeightPtr = weight_15;
        if (iint == 16)
          WeightPtr = weight_16;
        if (iint == 17)
          WeightPtr = weight_17;
        if (iint == 18)
          WeightPtr = weight_18;
        if (iint == 19)
          WeightPtr = weight_19;
        if (iint == 20)
          WeightPtr = weight_20;
        if (iint == 21)
          WeightPtr = weight_21;
        if (iint == 22)
          WeightPtr = weight_22;
        if (iint == 23)
          WeightPtr = weight_23;
        if (iint == 24)
          WeightPtr = weight_24;

        for (int ibin = 0; ibin < 25; ++ibin) {
          WeightOOTPU_[iint][ibin] = *(WeightPtr + ibin);
        }
      }
    }

    double ITweight(int npv) {
      int bin = weights_->GetXaxis()->FindBin(npv);
      return weights_->GetBinContent(bin);
    }

    double ITweight3BX(float ave_int) {
      int bin = weights_->GetXaxis()->FindBin(ave_int);
      return weights_->GetBinContent(bin);
    }

    double weight(float n_int) {
      int bin = weights_->GetXaxis()->FindBin(n_int);
      return weights_->GetBinContent(bin);
    }

    double weight3D(int pv1, int pv2, int pv3) {
      using std::min;

      int npm1 = min(pv1, 34);
      int np0 = min(pv2, 34);
      int npp1 = min(pv3, 34);

      return Weight3D_[npm1][np0][npp1];
    }

    double weightOOT(int npv_in_time, int npv_m50nsBX) {
      static const double Correct_Weights2011[25] = {
          // residual correction to match lumi spectrum
          5.30031,  2.07903,  1.40729,  1.27687, 1.0702, 0.902094, 0.902345, 0.931449, 0.78202,
          0.824686, 0.837735, 0.910261, 1.01394, 1.1599, 1.12778,  1.58423,  1.78868,  1.58296,
          2.3291,   3.86641,  0,        0,       0,      0,        0};

      if (FirstWarning_) {
        std::cout << " **** Warning: Out-of-time pileup reweighting appropriate only for PU_S3 **** " << std::endl;
        std::cout << " ****                          will be applied                           **** " << std::endl;

        FirstWarning_ = false;
      }

      // Note: for the "uncorrelated" out-of-time pileup, reweighting is only done on the 50ns
      // "late" bunch (BX=+1), since that is basically the only one that matters in terms of
      // energy deposition.

      if (npv_in_time < 0) {
        std::cerr << " no in-time beam crossing found\n! ";
        std::cerr << " Returning event weight=0\n! ";
        return 0.;
      }
      if (npv_m50nsBX < 0) {
        std::cerr << " no out-of-time beam crossing found\n! ";
        std::cerr << " Returning event weight=0\n! ";
        return 0.;
      }

      int bin = weights_->GetXaxis()->FindBin(npv_in_time);

      double inTimeWeight = weights_->GetBinContent(bin);

      double TotalWeight = 1.0;

      TotalWeight = inTimeWeight * WeightOOTPU_[bin - 1][npv_m50nsBX] * Correct_Weights2011[bin - 1];

      return TotalWeight;
    }

  protected:
    std::string generatedFileName_;
    std::string dataFileName_;
    std::string GenHistName_;
    std::string DataHistName_;
    TFile* generatedFile_;
    TFile* dataFile_;
    TH1F* weights_;

    //keep copies of normalized distributions:
    TH1F* MC_distr_;
    TH1F* Data_distr_;

    double WeightOOTPU_[25][25];
    double Weight3D_[50][50][50];

    bool FirstWarning_;
  };
}  // namespace reweight

#endif