Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:28

0001 #ifndef PhysicsTools_Utilities_interface_LumiReWeighting_h
0002 #define PhysicsTools_Utilities_interface_LumiReWeighting_h
0003 
0004 /**
0005   \class    LumiReWeighting LumiReWeighting.h "PhysicsTools/Utilities/interface/LumiReWeighting.h"
0006   \brief    Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data
0007 
0008   This class will trivially take two histograms:
0009   1. The generated "flat-to-N" distributions from a given processing
0010   2. A histogram generated from the "estimatePileup" macro here:
0011 
0012    -- This is the Stand-Alone version that doesn't use any CMS classes --
0013 
0014   https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup
0015 
0016   \authors Salvatore Rappoccio, Mike Hildreth
0017 */
0018 
0019 #include "TH1F.h"
0020 #include "TH3.h"
0021 #include "TFile.h"
0022 #include "TRandom1.h"
0023 #include "TRandom2.h"
0024 #include "TRandom3.h"
0025 #include "TStopwatch.h"
0026 #include <string>
0027 #include <vector>
0028 #include <cmath>
0029 #include <algorithm>
0030 #include <iostream>
0031 
0032 namespace reweight {
0033 
0034   // add a class to shift the mean of a poisson-like luminosity distribution by an arbitrary amount.
0035   // Only valid for small (<1.5) shifts about the 2011 lumi distribution for now, because shifts are non-linear
0036   // Each PoissonMeanShifter does one shift, so defining multiples can give you an arbitrary collection
0037 
0038   class PoissonMeanShifter {
0039   public:
0040     PoissonMeanShifter(){};
0041 
0042     PoissonMeanShifter(float Shift) {
0043       // these are the polynomial or exponential coefficients for each bin of a 25-bin sequence that
0044       // convert the Distribution of the 2011 luminosity to something with a lower or higher peak luminosity.
0045       // The distributions aren't quite poisson because they model luminosity decreasing during a fill. This implies that
0046       // they do get wider as the mean increases, so the weights are not linear with increasing mean.
0047 
0048       static const double p0_minus[20] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
0049                                           1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
0050       static const double p1_minus[20] = {-0.677786, -0.619614, -0.49465, -0.357963, -0.238359, -0.110002, 0.0348629,
0051                                           0.191263,  0.347648,  0.516615, 0.679646,  0.836673,  0.97764,   1.135,
0052                                           1.29922,   1.42467,   1.55901,  1.61762,   1.67275,   1.96008};
0053       static const double p2_minus[20] = {
0054           0.526164, 0.251816, 0.11049,  0.026917, -0.0464692, -0.087022, -0.0931581, -0.0714295, -0.0331772, 0.0347473,
0055           0.108658, 0.193048, 0.272314, 0.376357, 0.4964,     0.58854,   0.684959,   0.731063,   0.760044,   1.02386};
0056 
0057       static const double p1_expoM[5] = {1.63363e-03, 6.79290e-04, 3.69900e-04, 2.24349e-04, 9.87156e-06};
0058 
0059       static const double p2_expoM[5] = {2.64692, 3.26585, 3.53229, 4.18035, 5.64027};
0060 
0061       static const double p0_plus[20] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
0062                                          1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
0063       static const double p1_plus[20] = {-0.739059, -0.594445, -0.477276, -0.359707, -0.233573, -0.103458, 0.0373401,
0064                                          0.176571,  0.337617,  0.499074,  0.675126,  0.840522,  1.00917,   1.15847,
0065                                          1.23816,   1.44271,   1.52982,   1.46385,   1.5802,    0.988689};
0066       static const double p2_plus[20] = {0.208068,    0.130033,   0.0850356,  0.0448344,  0.000749832,
0067                                          -0.0331347,  -0.0653281, -0.0746009, -0.0800667, -0.0527636,
0068                                          -0.00402649, 0.103338,   0.261261,   0.491084,   0.857966,
0069                                          1.19495,     1.75071,    2.65559,    3.35433,    5.48835};
0070 
0071       static const double p1_expoP[5] = {1.42463e-01, 4.18966e-02, 1.12697e-01, 1.66197e-01, 1.50768e-01};
0072 
0073       static const double p2_expoP[5] = {1.98758, 2.27217, 2.26799, 2.38455, 2.52428};
0074 
0075       // initialize weights based on desired Shift
0076 
0077       for (int ibin = 0; ibin < 20; ibin++) {
0078         if (Shift < .0) {
0079           Pweight_[ibin] = p0_minus[ibin] + p1_minus[ibin] * Shift + p2_minus[ibin] * Shift * Shift;
0080         } else {
0081           Pweight_[ibin] = p0_plus[ibin] + p1_plus[ibin] * Shift + p2_plus[ibin] * Shift * Shift;
0082         }
0083       }
0084 
0085       // last few bins fit better to an exponential...
0086 
0087       for (int ibin = 20; ibin < 25; ibin++) {
0088         if (Shift < 0.) {
0089           Pweight_[ibin] = p1_expoM[ibin - 20] * exp(p2_expoM[ibin - 20] * Shift);
0090         } else {
0091           Pweight_[ibin] = p1_expoP[ibin - 20] * exp(p2_expoP[ibin - 20] * Shift);
0092         }
0093       }
0094     };
0095 
0096     double ShiftWeight(int ibin) {
0097       if (ibin < 25 && ibin >= 0) {
0098         return Pweight_[ibin];
0099       } else {
0100         return 0;
0101       }
0102     };
0103 
0104     double ShiftWeight(float pvnum) {
0105       int ibin = int(pvnum);
0106 
0107       if (ibin < 25 && ibin >= 0) {
0108         return Pweight_[ibin];
0109       } else {
0110         return 0;
0111       }
0112     };
0113 
0114   private:
0115     double Pweight_[25];
0116   };
0117 
0118   class LumiReWeighting {
0119   public:
0120     LumiReWeighting(){};
0121 
0122     LumiReWeighting(std::string generatedFile, std::string dataFile, std::string GenHistName, std::string DataHistName)
0123         : generatedFileName_(generatedFile),
0124           dataFileName_(dataFile),
0125           GenHistName_(GenHistName),
0126           DataHistName_(DataHistName) {
0127       generatedFile_ = new TFile(generatedFileName_.c_str());  //MC distribution
0128       dataFile_ = new TFile(dataFileName_.c_str());            //Data distribution
0129 
0130       Data_distr_ = new TH1F(*(static_cast<TH1F*>(dataFile_->Get(DataHistName_.c_str())->Clone())));
0131       MC_distr_ = new TH1F(*(static_cast<TH1F*>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
0132 
0133       // normalize both histograms first
0134 
0135       Data_distr_->Scale(1.0 / Data_distr_->Integral());
0136       MC_distr_->Scale(1.0 / MC_distr_->Integral());
0137 
0138       weights_ = new TH1F(*(Data_distr_));
0139 
0140       // MC * data/MC = data, so the weights are data/MC:
0141 
0142       weights_->SetName("lumiWeights");
0143 
0144       TH1F* den = new TH1F(*(MC_distr_));
0145 
0146       weights_->Divide(den);  // so now the average weight should be 1.0
0147 
0148       std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;
0149 
0150       int NBins = weights_->GetNbinsX();
0151 
0152       for (int ibin = 1; ibin < NBins + 1; ++ibin) {
0153         std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
0154       }
0155 
0156       weightOOT_init();
0157 
0158       FirstWarning_ = true;
0159     }
0160 
0161     LumiReWeighting(const std::vector<float>& MC_distr, const std::vector<float>& Lumi_distr) {
0162       // no histograms for input: use vectors
0163 
0164       // now, make histograms out of them:
0165 
0166       // first, check they are the same size...
0167 
0168       if (MC_distr.size() != Lumi_distr.size()) {
0169         std::cerr << "ERROR: LumiReWeighting: input vectors have different sizes. Quitting... \n";
0170         return;
0171       }
0172 
0173       Int_t NBins = MC_distr.size();
0174 
0175       MC_distr_ = new TH1F("MC_distr", "MC dist", NBins, -0.5, float(NBins) - 0.5);
0176       Data_distr_ = new TH1F("Data_distr", "Data dist", NBins, -0.5, float(NBins) - 0.5);
0177 
0178       weights_ = new TH1F("luminumer", "luminumer", NBins, -0.5, float(NBins) - 0.5);
0179       TH1F* den = new TH1F("lumidenom", "lumidenom", NBins, -0.5, float(NBins) - 0.5);
0180 
0181       for (int ibin = 1; ibin < NBins + 1; ++ibin) {
0182         weights_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
0183         Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
0184         den->SetBinContent(ibin, MC_distr[ibin - 1]);
0185         MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
0186       }
0187 
0188       // check integrals, make sure things are normalized
0189 
0190       float deltaH = weights_->Integral();
0191       if (fabs(1.0 - deltaH) > 0.02) {  //*OOPS*...
0192         weights_->Scale(1.0 / deltaH);
0193         Data_distr_->Scale(1.0 / deltaH);
0194       }
0195       float deltaMC = den->Integral();
0196       if (fabs(1.0 - deltaMC) > 0.02) {
0197         den->Scale(1.0 / deltaMC);
0198         MC_distr_->Scale(1.0 / deltaMC);
0199       }
0200 
0201       weights_->Divide(den);  // so now the average weight should be 1.0
0202 
0203       std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;
0204 
0205       for (int ibin = 1; ibin < NBins + 1; ++ibin) {
0206         std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
0207       }
0208 
0209       weightOOT_init();
0210 
0211       FirstWarning_ = true;
0212     }
0213 
0214     void weight3D_init(float ScaleFactor, std::string WeightOutputFile = "") {
0215       //create histogram to write output weights, save pain of generating them again...
0216 
0217       TH3D* WHist = new TH3D("WHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
0218       TH3D* DHist = new TH3D("DHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
0219       TH3D* MHist = new TH3D("MHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
0220 
0221       using std::min;
0222 
0223       if (MC_distr_->GetEntries() == 0) {
0224         std::cout << " MC and Data distributions are not initialized! You must call the LumiReWeighting constructor. "
0225                   << std::endl;
0226       }
0227 
0228       // arrays for storing number of interactions
0229 
0230       double MC_ints[50][50][50];
0231       double Data_ints[50][50][50];
0232 
0233       for (int i = 0; i < 50; i++) {
0234         for (int j = 0; j < 50; j++) {
0235           for (int k = 0; k < 50; k++) {
0236             MC_ints[i][j][k] = 0.;
0237             Data_ints[i][j][k] = 0.;
0238           }
0239         }
0240       }
0241 
0242       double factorial[50];
0243       double PowerSer[50];
0244       double base = 1.;
0245 
0246       factorial[0] = 1.;
0247       PowerSer[0] = 1.;
0248 
0249       for (int i = 1; i < 50; ++i) {
0250         base = base * float(i);
0251         factorial[i] = base;
0252       }
0253 
0254       double x;
0255       double xweight;
0256       double probi, probj, probk;
0257       double Expval, mean;
0258       int xi;
0259 
0260       // Get entries for Data, MC, fill arrays:
0261       int NMCbin = MC_distr_->GetNbinsX();
0262 
0263       for (int jbin = 1; jbin < NMCbin + 1; jbin++) {
0264         x = MC_distr_->GetBinCenter(jbin);
0265         xweight = MC_distr_->GetBinContent(jbin);  //use as weight for matrix
0266 
0267         //for Summer 11, we have this int feature:
0268         xi = int(x);
0269 
0270         // Generate Poisson distribution for each value of the mean
0271         mean = double(xi);
0272 
0273         if (mean < 0.) {
0274           std::cout << "LumiReweighting:BadInputValue"
0275                     << " Your histogram generates MC luminosity values less than zero!"
0276                     << " Please Check.  Terminating." << std::endl;
0277         }
0278 
0279         if (mean == 0.) {
0280           Expval = 1.;
0281         } else {
0282           Expval = exp(-1. * mean);
0283         }
0284 
0285         base = 1.;
0286 
0287         for (int i = 1; i < 50; ++i) {
0288           base = base * mean;
0289           PowerSer[i] = base;  // PowerSer is mean^i
0290         }
0291 
0292         // compute poisson probability for each Nvtx in weight matrix
0293         for (int i = 0; i < 50; i++) {
0294           probi = PowerSer[i] / factorial[i] * Expval;
0295           for (int j = 0; j < 50; j++) {
0296             probj = PowerSer[j] / factorial[j] * Expval;
0297             for (int k = 0; k < 50; k++) {
0298               probk = PowerSer[k] / factorial[k] * Expval;
0299               // joint probability is product of event weights multiplied by weight of input distribution bin
0300               MC_ints[i][j][k] = MC_ints[i][j][k] + probi * probj * probk * xweight;
0301             }
0302           }
0303         }
0304       }
0305 
0306       int NDatabin = Data_distr_->GetNbinsX();
0307 
0308       for (int jbin = 1; jbin < NDatabin + 1; jbin++) {
0309         mean = (Data_distr_->GetBinCenter(jbin)) * ScaleFactor;
0310         xweight = Data_distr_->GetBinContent(jbin);
0311 
0312         // Generate poisson distribution for each value of the mean
0313         if (mean < 0.) {
0314           std::cout << "LumiReweighting:BadInputValue"
0315                     << " Your histogram generates MC luminosity values less than zero!"
0316                     << " Please Check.  Terminating." << std::endl;
0317         }
0318 
0319         if (mean == 0.) {
0320           Expval = 1.;
0321         } else {
0322           Expval = exp(-1. * mean);
0323         }
0324 
0325         base = 1.;
0326 
0327         for (int i = 1; i < 50; ++i) {
0328           base = base * mean;
0329           PowerSer[i] = base;
0330         }
0331 
0332         // compute poisson probability for each Nvtx in weight matrix
0333 
0334         for (int i = 0; i < 50; i++) {
0335           probi = PowerSer[i] / factorial[i] * Expval;
0336           for (int j = 0; j < 50; j++) {
0337             probj = PowerSer[j] / factorial[j] * Expval;
0338             for (int k = 0; k < 50; k++) {
0339               probk = PowerSer[k] / factorial[k] * Expval;
0340               // joint probability is product of event weights multiplied by weight of input distribution bin
0341               Data_ints[i][j][k] = Data_ints[i][j][k] + probi * probj * probk * xweight;
0342             }
0343           }
0344         }
0345       }
0346 
0347       for (int i = 0; i < 50; i++) {
0348         //if(i<5) std::cout << "i = " << i << std::endl;
0349         for (int j = 0; j < 50; j++) {
0350           for (int k = 0; k < 50; k++) {
0351             if ((MC_ints[i][j][k]) > 0.) {
0352               Weight3D_[i][j][k] = Data_ints[i][j][k] / MC_ints[i][j][k];
0353             } else {
0354               Weight3D_[i][j][k] = 0.;
0355             }
0356             WHist->SetBinContent(i + 1, j + 1, k + 1, Weight3D_[i][j][k]);
0357             DHist->SetBinContent(i + 1, j + 1, k + 1, Data_ints[i][j][k]);
0358             MHist->SetBinContent(i + 1, j + 1, k + 1, MC_ints[i][j][k]);
0359             //      if(i<5 && j<5 && k<5) std::cout << Weight3D_[i][j][k] << " " ;
0360           }
0361           //      if(i<5 && j<5) std::cout << std::endl;
0362         }
0363       }
0364 
0365       if (!WeightOutputFile.empty()) {
0366         std::cout << " 3D Weight Matrix initialized! " << std::endl;
0367         std::cout << " Writing weights to file " << WeightOutputFile << " for re-use...  " << std::endl;
0368 
0369         TFile* outfile = new TFile(WeightOutputFile.c_str(), "RECREATE");
0370         WHist->Write();
0371         MHist->Write();
0372         DHist->Write();
0373         outfile->Write();
0374         outfile->Close();
0375         outfile->Delete();
0376       }
0377 
0378       return;
0379     }
0380 
0381     void weight3D_set(std::string WeightFileName) {
0382       TFile* infile = new TFile(WeightFileName.c_str());
0383       TH1F* WHist = (TH1F*)infile->Get("WHist");
0384 
0385       // Check if the histogram exists
0386       if (!WHist) {
0387         std::cout << " Could not find the histogram WHist in the file "
0388                   << "in the file " << WeightFileName << "." << std::endl;
0389         return;
0390       }
0391 
0392       for (int i = 0; i < 50; i++) {
0393         for (int j = 0; j < 50; j++) {
0394           for (int k = 0; k < 50; k++) {
0395             Weight3D_[i][j][k] = WHist->GetBinContent(i, j, k);
0396           }
0397         }
0398       }
0399 
0400       std::cout << " 3D Weight Matrix initialized! " << std::endl;
0401 
0402       return;
0403     }
0404 
0405     void weightOOT_init() {
0406       // The following are poisson distributions with different means, where the maximum
0407       // of the function has been normalized to weight 1.0
0408       // These are used to reweight the out-of-time pileup to match the in-time distribution.
0409       // The total event weight is the product of the in-time weight, the out-of-time weight,
0410       // and a residual correction to fix the distortions caused by the fact that the out-of-time
0411       // distribution is not flat.
0412 
0413       static const double weight_24[25] = {0,           0,           0,           0,          2.46277e-06,
0414                                            2.95532e-05, 0.000104668, 0.000401431, 0.00130034, 0.00342202,
0415                                            0.00818132,  0.0175534,   0.035784,    0.0650836,  0.112232,
0416                                            0.178699,    0.268934,    0.380868,    0.507505,   0.640922,
0417                                            0.768551,    0.877829,    0.958624,    0.99939,    1};
0418 
0419       static const double weight_23[25] = {0,           1.20628e-06, 1.20628e-06, 2.41255e-06, 1.20628e-05,
0420                                            6.39326e-05, 0.000252112, 0.000862487, 0.00244995,  0.00616527,
0421                                            0.0140821,   0.0293342,   0.0564501,   0.100602,    0.164479,
0422                                            0.252659,    0.36268,     0.491427,    0.627979,    0.75918,
0423                                            0.873185,    0.957934,    0.999381,    1,           0.957738};
0424 
0425       static const double weight_22[25] = {
0426           0,         0,         0,         5.88636e-06, 3.0609e-05, 0.000143627, 0.000561558, 0.00173059, 0.00460078,
0427           0.0110616, 0.0238974, 0.0475406, 0.0875077,   0.148682,   0.235752,    0.343591,    0.473146,   0.611897,
0428           0.748345,  0.865978,  0.953199,  0.997848,    1,          0.954245,    0.873688};
0429 
0430       static const double weight_21[25] = {
0431           0,         0,         1.15381e-06, 8.07665e-06, 7.1536e-05, 0.000280375, 0.00107189, 0.00327104, 0.00809396,
0432           0.0190978, 0.0401894, 0.0761028,   0.13472,     0.216315,   0.324649,    0.455125,   0.598241,   0.739215,
0433           0.861866,  0.953911,  0.998918,    1,           0.956683,   0.872272,    0.76399};
0434 
0435       static const double weight_20[25] = {
0436           0,         0,         1.12532e-06, 2.58822e-05, 0.000145166, 0.000633552, 0.00215048, 0.00592816, 0.0145605,
0437           0.0328367, 0.0652649, 0.11893,     0.19803,     0.305525,    0.436588,    0.581566,   0.727048,   0.8534,
0438           0.949419,  0.999785,  1,           0.953008,    0.865689,    0.753288,    0.62765};
0439       static const double weight_19[25] = {
0440           0,         0,        1.20714e-05, 5.92596e-05, 0.000364337, 0.00124994, 0.00403953, 0.0108149, 0.025824,
0441           0.0544969, 0.103567, 0.17936,     0.283532,    0.416091,    0.562078,   0.714714,   0.846523,  0.947875,
0442           1,         0.999448, 0.951404,    0.859717,    0.742319,    0.613601,   0.48552};
0443 
0444       static const double weight_18[25] = {
0445           0,         3.20101e-06, 2.88091e-05, 0.000164319, 0.000719161, 0.00250106, 0.00773685, 0.0197513, 0.0443693,
0446           0.0885998, 0.159891,    0.262607,    0.392327,    0.543125,    0.69924,    0.837474,   0.943486,  0.998029,
0447           1,         0.945937,    0.851807,    0.729309,    0.596332,    0.467818,   0.350434};
0448 
0449       static const double weight_17[25] = {
0450           1.03634e-06, 7.25437e-06, 4.97443e-05, 0.000340956, 0.00148715, 0.00501485, 0.0143067, 0.034679, 0.0742009,
0451           0.140287,    0.238288,    0.369416,    0.521637,    0.682368,   0.828634,   0.939655,  1,        0.996829,
0452           0.94062,     0.841575,    0.716664,    0.582053,    0.449595,   0.331336,   0.234332};
0453 
0454       static const double weight_16[25] = {
0455           4.03159e-06, 2.41895e-05, 0.000141106, 0.00081942, 0.00314565, 0.00990662, 0.026293, 0.0603881, 0.120973,
0456           0.214532,    0.343708,    0.501141,    0.665978,   0.820107,   0.938149,   1,        0.99941,   0.940768,
0457           0.837813,    0.703086,    0.564023,    0.42928,    0.312515,   0.216251,   0.14561};
0458 
0459       static const double weight_15[25] = {
0460           9.76084e-07, 5.07564e-05, 0.000303562, 0.00174036, 0.00617959, 0.0188579, 0.047465, 0.101656, 0.189492,
0461           0.315673,    0.474383,    0.646828,    0.809462,   0.934107,   0.998874,  1,        0.936163, 0.827473,
0462           0.689675,    0.544384,    0.40907,     0.290648,   0.198861,   0.12951,   0.0808051};
0463 
0464       static const double weight_14[25] = {
0465           1.13288e-05, 0.000124617, 0.000753365, 0.00345056, 0.0123909, 0.0352712, 0.0825463, 0.16413,  0.287213,
0466           0.44615,     0.625826,    0.796365,    0.930624,   0.999958,  1,         0.934414,  0.816456, 0.672939,
0467           0.523033,    0.386068,    0.269824,    0.180342,   0.114669,  0.0698288, 0.0406496};
0468 
0469       static const double weight_13[25] = {
0470           2.54296e-05, 0.000261561, 0.00167018, 0.00748083, 0.0241308, 0.0636801, 0.138222, 0.255814, 0.414275,
0471           0.600244,    0.779958,    0.92256,    0.999155,   1,         0.927126,  0.804504, 0.651803, 0.497534,
0472           0.35976,     0.245834,    0.160904,   0.0991589,  0.0585434, 0.0332437, 0.0180159};
0473 
0474       static const double weight_12[25] = {
0475           5.85742e-05, 0.000627706, 0.00386677, 0.0154068, 0.0465892, 0.111683,  0.222487,  0.381677, 0.5719,
0476           0.765001,    0.915916,    1,          0.999717,  0.921443,  0.791958,  0.632344,  0.475195, 0.334982,
0477           0.223666,    0.141781,    0.0851538,  0.048433,  0.0263287, 0.0133969, 0.00696683};
0478 
0479       static const double weight_11[25] = {
0480           0.00015238, 0.00156064, 0.00846044, 0.0310939, 0.0856225, 0.187589,   0.343579,  0.541892, 0.74224,
0481           0.909269,   0.998711,   1,          0.916889,  0.77485,   0.608819,   0.447016,  0.307375, 0.198444,
0482           0.121208,   0.070222,   0.0386492,  0.0201108, 0.0100922, 0.00484937, 0.00222458};
0483 
0484       static const double weight_10[25] = {
0485           0.000393044, 0.00367001, 0.0179474, 0.060389,   0.151477,   0.302077,   0.503113,   0.720373, 0.899568,
0486           1,           0.997739,   0.909409,  0.75728,    0.582031,   0.415322,   0.277663,   0.174147, 0.102154,
0487           0.0566719,   0.0298642,  0.0147751, 0.00710995, 0.00319628, 0.00140601, 0.000568796};
0488 
0489       static const double weight_9[25] = {
0490           0.00093396, 0.00854448, 0.0380306,  0.113181,  0.256614,    0.460894,    0.690242,   0.888781,  1,
0491           0.998756,   0.899872,   0.735642,   0.552532,  0.382726,    0.246114,    0.147497,   0.0825541, 0.0441199,
0492           0.0218157,  0.0103578,  0.00462959, 0.0019142, 0.000771598, 0.000295893, 0.000111529};
0493 
0494       static const double weight_8[25] = {
0495           0.00240233, 0.0192688, 0.0768653,  0.205008,    0.410958,    0.65758,     0.875657,   0.999886,  1,
0496           0.889476,   0.711446,  0.517781,   0.345774,    0.212028,    0.121208,    0.0644629,  0.0324928, 0.0152492,
0497           0.00673527, 0.0028547, 0.00117213, 0.000440177, 0.000168471, 5.80689e-05, 1.93563e-05};
0498 
0499       static const double weight_7[25] = {0.00617233,  0.0428714,   0.150018,    0.350317,    0.612535,
0500                                           0.856525,    0.999923,    1,           0.87544,     0.679383,
0501                                           0.478345,    0.303378,    0.176923,    0.0950103,   0.0476253,
0502                                           0.0222211,   0.00972738,  0.00392962,  0.0015258,   0.000559168,
0503                                           0.000183928, 6.77983e-05, 1.67818e-05, 7.38398e-06, 6.71271e-07};
0504 
0505       static const double weight_6[25] = {0.0154465,   0.0923472,   0.277322,    0.55552,     0.833099,
0506                                           0.999035,    1,           0.855183,    0.641976,    0.428277,
0507                                           0.256804,    0.139798,    0.0700072,   0.0321586,   0.0137971,
0508                                           0.00544756,  0.00202316,  0.000766228, 0.000259348, 8.45836e-05,
0509                                           1.80362e-05, 8.70713e-06, 3.73163e-06, 6.21938e-07, 0};
0510 
0511       static const double weight_5[25] = {0.0382845,   0.191122,    0.478782,    0.797314,    1,
0512                                           0.997148,    0.831144,    0.59461,     0.371293,    0.205903,
0513                                           0.103102,    0.0471424,   0.0194997,   0.00749415,  0.00273709,
0514                                           0.000879189, 0.000286049, 0.000102364, 1.70606e-05, 3.98081e-06,
0515                                           2.27475e-06, 0,           0,           0,           0};
0516 
0517       static const double weight_4[25] = {0.0941305,   0.373824,    0.750094,    1,           0.997698,
0518                                           0.800956,    0.532306,    0.304597,    0.152207,    0.0676275,
0519                                           0.0270646,   0.00975365,  0.00326077,  0.00101071,  0.000301781,
0520                                           7.41664e-05, 1.58563e-05, 3.58045e-06, 1.02299e-06, 0,
0521                                           5.11493e-07, 0,           0,           0,           0};
0522 
0523       static const double weight_3[25] = {0.222714,    0.667015,    1,           0.999208,    0.750609,
0524                                           0.449854,    0.224968,    0.0965185,   0.0361225,   0.012084,
0525                                           0.00359618,  0.000977166, 0.000239269, 6.29422e-05, 1.16064e-05,
0526                                           1.78559e-06, 0,           4.46398e-07, 0,           0,
0527                                           0,           0,           0,           0,           0};
0528 
0529       static const double weight_2[25] = {
0530           0.499541,   0.999607,    1,           0.666607,    0.333301,    0.13279, 0.0441871, 0.0127455, 0.00318434,
0531           0.00071752, 0.000132204, 2.69578e-05, 5.16999e-06, 2.21571e-06, 0,       0,         0,         0,
0532           0,          0,           0,           0,           0,           0,       0};
0533 
0534       static const double weight_1[25] = {
0535           0.999165,    1,           0.499996,    0.166868, 0.0414266, 0.00831053, 0.00137472, 0.000198911, 2.66302e-05,
0536           2.44563e-06, 2.71737e-07, 2.71737e-07, 0,        0,         0,          0,          0,           0,
0537           0,           0,           0,           0,        0,         0,          0};
0538 
0539       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};
0540 
0541       //WeightOOTPU_ = {0};
0542 
0543       const double* WeightPtr = nullptr;
0544 
0545       for (int iint = 0; iint < 25; ++iint) {
0546         if (iint == 0)
0547           WeightPtr = weight_0;
0548         if (iint == 1)
0549           WeightPtr = weight_1;
0550         if (iint == 2)
0551           WeightPtr = weight_2;
0552         if (iint == 3)
0553           WeightPtr = weight_3;
0554         if (iint == 4)
0555           WeightPtr = weight_4;
0556         if (iint == 5)
0557           WeightPtr = weight_5;
0558         if (iint == 6)
0559           WeightPtr = weight_6;
0560         if (iint == 7)
0561           WeightPtr = weight_7;
0562         if (iint == 8)
0563           WeightPtr = weight_8;
0564         if (iint == 9)
0565           WeightPtr = weight_9;
0566         if (iint == 10)
0567           WeightPtr = weight_10;
0568         if (iint == 11)
0569           WeightPtr = weight_11;
0570         if (iint == 12)
0571           WeightPtr = weight_12;
0572         if (iint == 13)
0573           WeightPtr = weight_13;
0574         if (iint == 14)
0575           WeightPtr = weight_14;
0576         if (iint == 15)
0577           WeightPtr = weight_15;
0578         if (iint == 16)
0579           WeightPtr = weight_16;
0580         if (iint == 17)
0581           WeightPtr = weight_17;
0582         if (iint == 18)
0583           WeightPtr = weight_18;
0584         if (iint == 19)
0585           WeightPtr = weight_19;
0586         if (iint == 20)
0587           WeightPtr = weight_20;
0588         if (iint == 21)
0589           WeightPtr = weight_21;
0590         if (iint == 22)
0591           WeightPtr = weight_22;
0592         if (iint == 23)
0593           WeightPtr = weight_23;
0594         if (iint == 24)
0595           WeightPtr = weight_24;
0596 
0597         for (int ibin = 0; ibin < 25; ++ibin) {
0598           WeightOOTPU_[iint][ibin] = *(WeightPtr + ibin);
0599         }
0600       }
0601     }
0602 
0603     double ITweight(int npv) {
0604       int bin = weights_->GetXaxis()->FindBin(npv);
0605       return weights_->GetBinContent(bin);
0606     }
0607 
0608     double ITweight3BX(float ave_int) {
0609       int bin = weights_->GetXaxis()->FindBin(ave_int);
0610       return weights_->GetBinContent(bin);
0611     }
0612 
0613     double weight(float n_int) {
0614       int bin = weights_->GetXaxis()->FindBin(n_int);
0615       return weights_->GetBinContent(bin);
0616     }
0617 
0618     double weight3D(int pv1, int pv2, int pv3) {
0619       using std::min;
0620 
0621       int npm1 = min(pv1, 34);
0622       int np0 = min(pv2, 34);
0623       int npp1 = min(pv3, 34);
0624 
0625       return Weight3D_[npm1][np0][npp1];
0626     }
0627 
0628     double weightOOT(int npv_in_time, int npv_m50nsBX) {
0629       static const double Correct_Weights2011[25] = {
0630           // residual correction to match lumi spectrum
0631           5.30031,  2.07903,  1.40729,  1.27687, 1.0702, 0.902094, 0.902345, 0.931449, 0.78202,
0632           0.824686, 0.837735, 0.910261, 1.01394, 1.1599, 1.12778,  1.58423,  1.78868,  1.58296,
0633           2.3291,   3.86641,  0,        0,       0,      0,        0};
0634 
0635       if (FirstWarning_) {
0636         std::cout << " **** Warning: Out-of-time pileup reweighting appropriate only for PU_S3 **** " << std::endl;
0637         std::cout << " ****                          will be applied                           **** " << std::endl;
0638 
0639         FirstWarning_ = false;
0640       }
0641 
0642       // Note: for the "uncorrelated" out-of-time pileup, reweighting is only done on the 50ns
0643       // "late" bunch (BX=+1), since that is basically the only one that matters in terms of
0644       // energy deposition.
0645 
0646       if (npv_in_time < 0) {
0647         std::cerr << " no in-time beam crossing found\n! ";
0648         std::cerr << " Returning event weight=0\n! ";
0649         return 0.;
0650       }
0651       if (npv_m50nsBX < 0) {
0652         std::cerr << " no out-of-time beam crossing found\n! ";
0653         std::cerr << " Returning event weight=0\n! ";
0654         return 0.;
0655       }
0656 
0657       int bin = weights_->GetXaxis()->FindBin(npv_in_time);
0658 
0659       double inTimeWeight = weights_->GetBinContent(bin);
0660 
0661       double TotalWeight = 1.0;
0662 
0663       TotalWeight = inTimeWeight * WeightOOTPU_[bin - 1][npv_m50nsBX] * Correct_Weights2011[bin - 1];
0664 
0665       return TotalWeight;
0666     }
0667 
0668   protected:
0669     std::string generatedFileName_;
0670     std::string dataFileName_;
0671     std::string GenHistName_;
0672     std::string DataHistName_;
0673     TFile* generatedFile_;
0674     TFile* dataFile_;
0675     TH1F* weights_;
0676 
0677     //keep copies of normalized distributions:
0678     TH1F* MC_distr_;
0679     TH1F* Data_distr_;
0680 
0681     double WeightOOTPU_[25][25];
0682     double Weight3D_[50][50][50];
0683 
0684     bool FirstWarning_;
0685   };
0686 }  // namespace reweight
0687 
0688 #endif