/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/PhysicsTools/Utilities/interface/LumiReweightingStandAlone.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_interface_LumiReWeighting_h
#define PhysicsTools_Utilities_interface_LumiReWeighting_h


#include "TH1.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>

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 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 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 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 double p1_expoM[5] = {
        1.63363e-03,
        6.79290e-04,
        3.69900e-04,
        2.24349e-04,
        9.87156e-06
      };

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


      static 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 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 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 double p1_expoP[5] = {
        1.42463e-01,
        4.18966e-02,
        1.12697e-01,
        1.66197e-01,
        1.50768e-01
      };

      static 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 TH1(  *(static_cast<TH1*>(dataFile_->Get( DataHistName_.c_str() )->Clone() )) );
          MC_distr_ = new TH1(  *(static_cast<TH1*>(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 TH1( *(Data_distr_)) ;

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

          weights_->SetName("lumiWeights");

          TH1* den = new TH1(*(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( std::vector< float > MC_distr, 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);
        TH1* 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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};

        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 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_;
      TH1  *weights_;

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

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

      bool FirstWarning_;


  };
}



#endif