yangyong/zSelection/LumiReWeighting.cc

//// from  PhysicsTools_Utilities_LumiReWeighting_cc


//int nWarning = 0; 

/**
  \class    LumiReWeighting LumiReWeighting.h "PhysicsTools/Utilities/interface/LumiReWeighting.h"
  \brief    Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data

  This class will trivially take two histograms:
  1. The generated "flat-to-N" distributions from a given processing (or any other generated input)
  2. A histogram generated from the "estimatePileup" macro here:

  https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup

  and produce weights to convert the input distribution (1) to the latter (2).

  \author Salvatore Rappoccio, modified by Mike Hildreth
  
*/

TH1D *weights_; 
double WeightOOTPU_[25][25];

void LumiReWeighting_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);
    }
  }

}


void LumiReWeighting_setpileupWeights(TH1D *hdata, TH1D *hmc){
  TH1D *num = (TH1D*)hdata->Clone("datawt");

  cout<<"dd" <<endl;
  weights_ = (TH1D*)hdata->Clone("normwt");
  double total = weights_->Integral();
  cout<<"dd1" <<endl;

  num ->Scale( 1.0/total);
  TH1D *den = (TH1D*)hmc->Clone("mcwt");
  double total1 = den->Integral();
  den ->Scale(1.0/total1);
  ////weights_ ->Divide( den);
  int nbin = weights_->GetNbinsX();
  
  double sum_wt = 0; 
  for(int b=1; b<= nbin; b++){
    if( den->GetBinContent(b) > 0){
      weights_ ->SetBinContent(b, num->GetBinContent(b)/ den->GetBinContent(b));
    }else{
      weights_ ->SetBinContent(b, 0);
    }
    sum_wt += weights_ ->GetBinContent(b) *  den->GetBinContent(b); 
    cout<<"intimepileweight "<< b-1 <<" "<<  weights_ ->GetBinContent(b)<<endl;
  }
  cout<<"sum wt "<< sum_wt <<" "<< weights_->GetMean()<<endl;  
  
  LumiReWeighting_weightOOT_init();
  
}


///Now use this version
vector<double>  LumiReWeighting_getWeights(TH1D *hdata, TH1D *hmc){
  TH1D *num = (TH1D*)hdata->Clone("datawt");
  double total = num->Integral();
  num ->Scale( 1.0/total);
  
  TH1D *den = (TH1D*)hmc->Clone("mcwt");
  double total1 = den->Integral();
  den ->Scale(1.0/total1);
  
  int nbindata = num->GetNbinsX();
  int nbinmc = den->GetNbinsX();
  
  if( nbinmc != nbindata ){
    cout<<" waring! nbinmc != nbindata " << nbinmc <<" "<< nbindata <<endl; 
  }
  
  
  vector<double> wts; 
  int nbin = nbinmc < nbindata ? nbinmc: nbindata; 
  
  for(int b=1; b<= nbin; b++){
    if( den->GetBinContent(b) > 0){
      wts.push_back( num->GetBinContent(b)/ den->GetBinContent(b));
    }else{
      wts.push_back(0);
    }
    cout<<" LumiReWeighting_getWeights " << wts[b-1] <<endl; 
    
  }
  if( nbinmc > nbindata){
    for(int b= nbindata; b < nbinmc; b++){
      wts.push_back(0);
    }
  }
  
  //final check
  if( int(wts.size()) != nbinmc){
    cout<<" problem of seeting pileup wts for MC " <<  wts.size() <<" "<< nbinmc <<endl; 
    exit(1);
  }

  return wts;
  
  
  ///LumiReWeighting_weightOOT_init();
  
}


double LumiReWeighting_weightINT_v1(vector<double> weight){
  
  if( weight.size() <=0){
    cout<<" LumiReWeighting_weightINT_v1 bad weight vector " << endl; 
    exit(1);
  }
  
  int npv = -1; 
  for(unsigned int j =0; j < pileupBunchX->size(); j++){
    int BX = pileupBunchX->at(j);
    if(BX==0){
      npv = pileupNInteraction->at(j);
      break; 
    }
  }
  if(npv < 0) {
    cout << " no in-time beam crossing found\n! " <<endl; 
    return 0;
  }
  if(npv > int(weight.size()) ) {
    cout<<"warning! LumiReWeighting_weightINT_v1 " << npv <<" "<< weight.size() <<endl; 
    return 0; 
  }
  return weight[npv];
  
}


double LumiReWeighting_weightINT_v2(int npv, vector<double> weight){
  
  if( weight.size() <=0){
    cout<<" LumiReWeighting_weightINT_v2 bad weight vector " << endl; 
    exit(1);
  }
  
//   int npv = -1; 
//   for(unsigned int j =0; j < pileupBunchX->size(); j++){
//     int BX = pileupBunchX->at(j);
//     if(BX==0){
//       npv = pileupNInteraction->at(j);
//       break; 
//     }
//   }
  if(npv < 0) {
    cout << " no in-time beam crossing found\n! " <<endl; 
    return 0;
  }
  if(npv > int(weight.size()) ) {
    cout<<"warning! LumiReWeighting_weightINT_v2 " << npv <<" "<< weight.size() <<endl; 
    return 0; 
  }
  return weight[npv];
  
}


double LumiReWeighting_weightINT(){
  
  int npv = -1; 
  for(unsigned int j =0; j < pileupBunchX->size(); j++){
    int BX = pileupBunchX->at(j);
    if(BX==0){
      npv = pileupNInteraction->at(j);
      break; 
    }
  }
  if(npv < 0) cout << " no in-time beam crossing found\n! " <<endl; 

  if( weights_ == NULL){
    cout<<" weights histogram not filled"<<endl;
    return 1; 
  }
  int bin = weights_->GetXaxis()->FindBin( npv );
  double inTimeWeight = weights_->GetBinContent( bin );
  return inTimeWeight; 
  
}


double LumiReWeighting_weightOOT(string dataset) {
    
  if( dataset.find("SIM") == string::npos){
    cout<<" weight not on MC ? " << dataset.c_str() <<endl; 
    exit(1);
  }
  
  
  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
  };                        
  
  
  int npv = -1;
  int npv50ns = -1;
  for(unsigned int j =0; j < pileupBunchX->size(); j++){
    int BX = pileupBunchX->at(j);
    if(BX==0){
      npv = pileupNInteraction->at(j);
    }
    if(BX == 1) { 
      npv50ns = pileupNInteraction->at(j);
    }
    
  }
  
  // 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 < 0) {
    std::cerr << " no in-time beam crossing found\n! " ;
    std::cerr << " Returning event weight=0\n! ";
    return 0.;
  }
 //  if(npv50ns < 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 );

  double inTimeWeight = weights_->GetBinContent( bin );
  
  double TotalWeight = 1.0;
  
  
  if( dataset.find("S3") != string::npos){    //if(Reweight_4_2_2p2_) {
    TotalWeight = inTimeWeight * WeightOOTPU_[bin-1][npv50ns] * Correct_Weights2011[bin-1];
  }
  else {
    TotalWeight = inTimeWeight;
  }
  
//   if( TotalWeight < 1E-5 && nWarning < 10 ){
//     cout<<"warnign 0 weight " << TotalWeight <<" "<< inTimeWeight <<" "<< bin <<" "<< npv<<" "<< npv50ns <<" "<< WeightOOTPU_[bin-1][npv50ns] <<" "<<Correct_Weights2011[bin-1] <<endl; 
//     nWarning ++; 
//   }
  
  
  return TotalWeight;
  
}


double LumiReWeighting_weightOOT_givenWeightVector(string dataset, vector<double> puwts) {
  
  if( dataset.find("SIM") == string::npos){
    cout<<" weight not on MC ? " << dataset.c_str() <<endl; 
    exit(1);
  }
  
  
  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
  };                        
  
  
  int npv = -1;
  int npv50ns = -1;
  for(unsigned int j =0; j < pileupBunchX->size(); j++){
    int BX = pileupBunchX->at(j);
    if(BX==0){
      npv = pileupNInteraction->at(j);
    }
    if(BX == 1) { 
      npv50ns = pileupNInteraction->at(j);
    }
    
  }
  
  // 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 < 0) {
    std::cerr << " no in-time beam crossing found\n! " ;
    std::cerr << " Returning event weight=0\n! ";
    return 0.;
  }
 //  if(npv50ns < 0) {
//     std::cerr << " no out-of-time beam crossing found\n! " ;
//     std::cerr << " Returning event weight=0\n! ";
//     return 0.;
//   }
  
  
  double inTimeWeight = puwts[npv];
  
  double TotalWeight = 1.0;
  
  
  if( dataset.find("S3") != string::npos){    //if(Reweight_4_2_2p2_) {
    TotalWeight = inTimeWeight * WeightOOTPU_[npv][npv50ns] * Correct_Weights2011[npv];
  }
  else {
    TotalWeight = inTimeWeight;
  }
  
//   if( TotalWeight < 1E-5 && nWarning < 10 ){
//     cout<<"warnign 0 weight " << TotalWeight <<" "<< inTimeWeight <<" "<< bin <<" "<< npv<<" "<< npv50ns <<" "<< WeightOOTPU_[bin-1][npv50ns] <<" "<<Correct_Weights2011[bin-1] <<endl; 
//     nWarning ++; 
//   }
  
  
  return TotalWeight;
  
}

Revision:	1.1
Committed:	Fri Sep 30 15:23:48 2011 UTC (13 years, 7 months ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-00-01, V01-00-00, HEAD
Log Message:	z selection code
#	User	Rev	Content
1	yangyong	1.1	//// from PhysicsTools_Utilities_LumiReWeighting_cc
2
3
4			//int nWarning = 0;
5
6			/**
7			\class LumiReWeighting LumiReWeighting.h "PhysicsTools/Utilities/interface/LumiReWeighting.h"
8			\brief Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data
9
10			This class will trivially take two histograms:
11			1. The generated "flat-to-N" distributions from a given processing (or any other generated input)
12			2. A histogram generated from the "estimatePileup" macro here:
13
14			https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup
15
16			and produce weights to convert the input distribution (1) to the latter (2).
17
18			\author Salvatore Rappoccio, modified by Mike Hildreth
19
20			*/
21
22			TH1D *weights_;
23			double WeightOOTPU_[25][25];
24
25			void LumiReWeighting_weightOOT_init() {
26
27			// The following are poisson distributions with different means, where the maximum
28			// of the function has been normalized to weight 1.0
29			// These are used to reweight the out-of-time pileup to match the in-time distribution.
30			// The total event weight is the product of the in-time weight, the out-of-time weight,
31			// and a residual correction to fix the distortions caused by the fact that the out-of-time
32			// distribution is not flat.
33
34			static double weight_24[25] = {
35			0,
36			0,
37			0,
38			0,
39			2.46277e-06,
40			2.95532e-05,
41			0.000104668,
42			0.000401431,
43			0.00130034,
44			0.00342202,
45			0.00818132,
46			0.0175534,
47			0.035784,
48			0.0650836,
49			0.112232,
50			0.178699,
51			0.268934,
52			0.380868,
53			0.507505,
54			0.640922,
55			0.768551,
56			0.877829,
57			0.958624,
58			0.99939,
59			1
60			};
61
62			static double weight_23[25] = {
63			0,
64			1.20628e-06,
65			1.20628e-06,
66			2.41255e-06,
67			1.20628e-05,
68			6.39326e-05,
69			0.000252112,
70			0.000862487,
71			0.00244995,
72			0.00616527,
73			0.0140821,
74			0.0293342,
75			0.0564501,
76			0.100602,
77			0.164479,
78			0.252659,
79			0.36268,
80			0.491427,
81			0.627979,
82			0.75918,
83			0.873185,
84			0.957934,
85			0.999381,
86			1,
87			0.957738
88			};
89
90			static double weight_22[25] = {
91			0,
92			0,
93			0,
94			5.88636e-06,
95			3.0609e-05,
96			0.000143627,
97			0.000561558,
98			0.00173059,
99			0.00460078,
100			0.0110616,
101			0.0238974,
102			0.0475406,
103			0.0875077,
104			0.148682,
105			0.235752,
106			0.343591,
107			0.473146,
108			0.611897,
109			0.748345,
110			0.865978,
111			0.953199,
112			0.997848,
113			1,
114			0.954245,
115			0.873688
116			};
117
118			static double weight_21[25] = {
119			0,
120			0,
121			1.15381e-06,
122			8.07665e-06,
123			7.1536e-05,
124			0.000280375,
125			0.00107189,
126			0.00327104,
127			0.00809396,
128			0.0190978,
129			0.0401894,
130			0.0761028,
131			0.13472,
132			0.216315,
133			0.324649,
134			0.455125,
135			0.598241,
136			0.739215,
137			0.861866,
138			0.953911,
139			0.998918,
140			1,
141			0.956683,
142			0.872272,
143			0.76399
144			};
145
146
147			static double weight_20[25] = {
148			0,
149			0,
150			1.12532e-06,
151			2.58822e-05,
152			0.000145166,
153			0.000633552,
154			0.00215048,
155			0.00592816,
156			0.0145605,
157			0.0328367,
158			0.0652649,
159			0.11893,
160			0.19803,
161			0.305525,
162			0.436588,
163			0.581566,
164			0.727048,
165			0.8534,
166			0.949419,
167			0.999785,
168			1,
169			0.953008,
170			0.865689,
171			0.753288,
172			0.62765
173			};
174			static double weight_19[25] = {
175			0,
176			0,
177			1.20714e-05,
178			5.92596e-05,
179			0.000364337,
180			0.00124994,
181			0.00403953,
182			0.0108149,
183			0.025824,
184			0.0544969,
185			0.103567,
186			0.17936,
187			0.283532,
188			0.416091,
189			0.562078,
190			0.714714,
191			0.846523,
192			0.947875,
193			1,
194			0.999448,
195			0.951404,
196			0.859717,
197			0.742319,
198			0.613601,
199			0.48552
200			};
201
202			static double weight_18[25] = {
203			0,
204			3.20101e-06,
205			2.88091e-05,
206			0.000164319,
207			0.000719161,
208			0.00250106,
209			0.00773685,
210			0.0197513,
211			0.0443693,
212			0.0885998,
213			0.159891,
214			0.262607,
215			0.392327,
216			0.543125,
217			0.69924,
218			0.837474,
219			0.943486,
220			0.998029,
221			1,
222			0.945937,
223			0.851807,
224			0.729309,
225			0.596332,
226			0.467818,
227			0.350434
228			};
229
230
231			static double weight_17[25] = {
232			1.03634e-06,
233			7.25437e-06,
234			4.97443e-05,
235			0.000340956,
236			0.00148715,
237			0.00501485,
238			0.0143067,
239			0.034679,
240			0.0742009,
241			0.140287,
242			0.238288,
243			0.369416,
244			0.521637,
245			0.682368,
246			0.828634,
247			0.939655,
248			1,
249			0.996829,
250			0.94062,
251			0.841575,
252			0.716664,
253			0.582053,
254			0.449595,
255			0.331336,
256			0.234332
257			};
258
259
260			static double weight_16[25] = {
261			4.03159e-06,
262			2.41895e-05,
263			0.000141106,
264			0.00081942,
265			0.00314565,
266			0.00990662,
267			0.026293,
268			0.0603881,
269			0.120973,
270			0.214532,
271			0.343708,
272			0.501141,
273			0.665978,
274			0.820107,
275			0.938149,
276			1,
277			0.99941,
278			0.940768,
279			0.837813,
280			0.703086,
281			0.564023,
282			0.42928,
283			0.312515,
284			0.216251,
285			0.14561
286			};
287
288
289			static double weight_15[25] = {
290			9.76084e-07,
291			5.07564e-05,
292			0.000303562,
293			0.00174036,
294			0.00617959,
295			0.0188579,
296			0.047465,
297			0.101656,
298			0.189492,
299			0.315673,
300			0.474383,
301			0.646828,
302			0.809462,
303			0.934107,
304			0.998874,
305			1,
306			0.936163,
307			0.827473,
308			0.689675,
309			0.544384,
310			0.40907,
311			0.290648,
312			0.198861,
313			0.12951,
314			0.0808051
315			};
316
317
318			static double weight_14[25] = {
319			1.13288e-05,
320			0.000124617,
321			0.000753365,
322			0.00345056,
323			0.0123909,
324			0.0352712,
325			0.0825463,
326			0.16413,
327			0.287213,
328			0.44615,
329			0.625826,
330			0.796365,
331			0.930624,
332			0.999958,
333			1,
334			0.934414,
335			0.816456,
336			0.672939,
337			0.523033,
338			0.386068,
339			0.269824,
340			0.180342,
341			0.114669,
342			0.0698288,
343			0.0406496
344			};
345
346
347			static double weight_13[25] = {
348			2.54296e-05,
349			0.000261561,
350			0.00167018,
351			0.00748083,
352			0.0241308,
353			0.0636801,
354			0.138222,
355			0.255814,
356			0.414275,
357			0.600244,
358			0.779958,
359			0.92256,
360			0.999155,
361			1,
362			0.927126,
363			0.804504,
364			0.651803,
365			0.497534,
366			0.35976,
367			0.245834,
368			0.160904,
369			0.0991589,
370			0.0585434,
371			0.0332437,
372			0.0180159
373			};
374
375			static double weight_12[25] = {
376			5.85742e-05,
377			0.000627706,
378			0.00386677,
379			0.0154068,
380			0.0465892,
381			0.111683,
382			0.222487,
383			0.381677,
384			0.5719,
385			0.765001,
386			0.915916,
387			1,
388			0.999717,
389			0.921443,
390			0.791958,
391			0.632344,
392			0.475195,
393			0.334982,
394			0.223666,
395			0.141781,
396			0.0851538,
397			0.048433,
398			0.0263287,
399			0.0133969,
400			0.00696683
401			};
402
403
404			static double weight_11[25] = {
405			0.00015238,
406			0.00156064,
407			0.00846044,
408			0.0310939,
409			0.0856225,
410			0.187589,
411			0.343579,
412			0.541892,
413			0.74224,
414			0.909269,
415			0.998711,
416			1,
417			0.916889,
418			0.77485,
419			0.608819,
420			0.447016,
421			0.307375,
422			0.198444,
423			0.121208,
424			0.070222,
425			0.0386492,
426			0.0201108,
427			0.0100922,
428			0.00484937,
429			0.00222458
430			};
431
432			static double weight_10[25] = {
433			0.000393044,
434			0.00367001,
435			0.0179474,
436			0.060389,
437			0.151477,
438			0.302077,
439			0.503113,
440			0.720373,
441			0.899568,
442			1,
443			0.997739,
444			0.909409,
445			0.75728,
446			0.582031,
447			0.415322,
448			0.277663,
449			0.174147,
450			0.102154,
451			0.0566719,
452			0.0298642,
453			0.0147751,
454			0.00710995,
455			0.00319628,
456			0.00140601,
457			0.000568796
458			};
459
460
461			static double weight_9[25] = {
462			0.00093396,
463			0.00854448,
464			0.0380306,
465			0.113181,
466			0.256614,
467			0.460894,
468			0.690242,
469			0.888781,
470			1,
471			0.998756,
472			0.899872,
473			0.735642,
474			0.552532,
475			0.382726,
476			0.246114,
477			0.147497,
478			0.0825541,
479			0.0441199,
480			0.0218157,
481			0.0103578,
482			0.00462959,
483			0.0019142,
484			0.000771598,
485			0.000295893,
486			0.000111529
487			};
488
489
490			static double weight_8[25] = {
491			0.00240233,
492			0.0192688,
493			0.0768653,
494			0.205008,
495			0.410958,
496			0.65758,
497			0.875657,
498			0.999886,
499			1,
500			0.889476,
501			0.711446,
502			0.517781,
503			0.345774,
504			0.212028,
505			0.121208,
506			0.0644629,
507			0.0324928,
508			0.0152492,
509			0.00673527,
510			0.0028547,
511			0.00117213,
512			0.000440177,
513			0.000168471,
514			5.80689e-05,
515			1.93563e-05
516			};
517
518			static double weight_7[25] = {
519			0.00617233,
520			0.0428714,
521			0.150018,
522			0.350317,
523			0.612535,
524			0.856525,
525			0.999923,
526			1,
527			0.87544,
528			0.679383,
529			0.478345,
530			0.303378,
531			0.176923,
532			0.0950103,
533			0.0476253,
534			0.0222211,
535			0.00972738,
536			0.00392962,
537			0.0015258,
538			0.000559168,
539			0.000183928,
540			6.77983e-05,
541			1.67818e-05,
542			7.38398e-06,
543			6.71271e-07
544			};
545
546			static double weight_6[25] = {
547			0.0154465,
548			0.0923472,
549			0.277322,
550			0.55552,
551			0.833099,
552			0.999035,
553			1,
554			0.855183,
555			0.641976,
556			0.428277,
557			0.256804,
558			0.139798,
559			0.0700072,
560			0.0321586,
561			0.0137971,
562			0.00544756,
563			0.00202316,
564			0.000766228,
565			0.000259348,
566			8.45836e-05,
567			1.80362e-05,
568			8.70713e-06,
569			3.73163e-06,
570			6.21938e-07,
571			0
572			};
573
574
575			static double weight_5[25] = {
576			0.0382845,
577			0.191122,
578			0.478782,
579			0.797314,
580			1,
581			0.997148,
582			0.831144,
583			0.59461,
584			0.371293,
585			0.205903,
586			0.103102,
587			0.0471424,
588			0.0194997,
589			0.00749415,
590			0.00273709,
591			0.000879189,
592			0.000286049,
593			0.000102364,
594			1.70606e-05,
595			3.98081e-06,
596			2.27475e-06,
597			0,
598			0,
599			0,
600			0
601			};
602
603
604			static double weight_4[25] = {
605			0.0941305,
606			0.373824,
607			0.750094,
608			1,
609			0.997698,
610			0.800956,
611			0.532306,
612			0.304597,
613			0.152207,
614			0.0676275,
615			0.0270646,
616			0.00975365,
617			0.00326077,
618			0.00101071,
619			0.000301781,
620			7.41664e-05,
621			1.58563e-05,
622			3.58045e-06,
623			1.02299e-06,
624			0,
625			5.11493e-07,
626			0,
627			0,
628			0,
629			0
630			};
631
632
633			static double weight_3[25] = {
634			0.222714,
635			0.667015,
636			1,
637			0.999208,
638			0.750609,
639			0.449854,
640			0.224968,
641			0.0965185,
642			0.0361225,
643			0.012084,
644			0.00359618,
645			0.000977166,
646			0.000239269,
647			6.29422e-05,
648			1.16064e-05,
649			1.78559e-06,
650			0,
651			4.46398e-07,
652			0,
653			0,
654			0,
655			0,
656			0,
657			0,
658			0
659			};
660
661			static double weight_2[25] = {
662			0.499541,
663			0.999607,
664			1,
665			0.666607,
666			0.333301,
667			0.13279,
668			0.0441871,
669			0.0127455,
670			0.00318434,
671			0.00071752,
672			0.000132204,
673			2.69578e-05,
674			5.16999e-06,
675			2.21571e-06,
676			0,
677			0,
678			0,
679			0,
680			0,
681			0,
682			0,
683			0,
684			0,
685			0,
686			0
687			};
688
689			static double weight_1[25] = {
690			0.999165,
691			1,
692			0.499996,
693			0.166868,
694			0.0414266,
695			0.00831053,
696			0.00137472,
697			0.000198911,
698			2.66302e-05,
699			2.44563e-06,
700			2.71737e-07,
701			2.71737e-07,
702			0,
703			0,
704			0,
705			0,
706			0,
707			0,
708			0,
709			0,
710			0,
711			0,
712			0,
713			0,
714			0
715			};
716
717			static double weight_0[25] = {
718			1,
719			0,
720			0,
721			0,
722			0,
723			0,
724			0,
725			0,
726			0,
727			0,
728			0,
729			0,
730			0,
731			0,
732			0,
733			0,
734			0,
735			0,
736			0,
737			0,
738			0,
739			0,
740			0,
741			0,
742			0
743			};
744
745			//WeightOOTPU_ = {0};
746
747			double* WeightPtr = 0;
748
749			for(int iint = 0; iint<25; ++iint){
750			if(iint ==0) WeightPtr = weight_0;
751			if(iint ==1) WeightPtr = weight_1;
752			if(iint ==2) WeightPtr = weight_2;
753			if(iint ==3) WeightPtr = weight_3;
754			if(iint ==4) WeightPtr = weight_4;
755			if(iint ==5) WeightPtr = weight_5;
756			if(iint ==6) WeightPtr = weight_6;
757			if(iint ==7) WeightPtr = weight_7;
758			if(iint ==8) WeightPtr = weight_8;
759			if(iint ==9) WeightPtr = weight_9;
760			if(iint ==10) WeightPtr = weight_10;
761			if(iint ==11) WeightPtr = weight_11;
762			if(iint ==12) WeightPtr = weight_12;
763			if(iint ==13) WeightPtr = weight_13;
764			if(iint ==14) WeightPtr = weight_14;
765			if(iint ==15) WeightPtr = weight_15;
766			if(iint ==16) WeightPtr = weight_16;
767			if(iint ==17) WeightPtr = weight_17;
768			if(iint ==18) WeightPtr = weight_18;
769			if(iint ==19) WeightPtr = weight_19;
770			if(iint ==20) WeightPtr = weight_20;
771			if(iint ==21) WeightPtr = weight_21;
772			if(iint ==22) WeightPtr = weight_22;
773			if(iint ==23) WeightPtr = weight_23;
774			if(iint ==24) WeightPtr = weight_24;
775
776			for(int ibin = 0; ibin<25; ++ibin){
777			WeightOOTPU_[iint][ibin] = *(WeightPtr+ibin);
778			}
779			}
780
781			}
782
783
784
785			void LumiReWeighting_setpileupWeights(TH1D hdata, TH1D hmc){
786			TH1D num = (TH1D)hdata->Clone("datawt");
787
788			cout<<"dd" <<endl;
789			weights_ = (TH1D*)hdata->Clone("normwt");
790			double total = weights_->Integral();
791			cout<<"dd1" <<endl;
792
793			num ->Scale( 1.0/total);
794			TH1D den = (TH1D)hmc->Clone("mcwt");
795			double total1 = den->Integral();
796			den ->Scale(1.0/total1);
797			////weights_ ->Divide( den);
798			int nbin = weights_->GetNbinsX();
799
800			double sum_wt = 0;
801			for(int b=1; b<= nbin; b++){
802			if( den->GetBinContent(b) > 0){
803			weights_ ->SetBinContent(b, num->GetBinContent(b)/ den->GetBinContent(b));
804			}else{
805			weights_ ->SetBinContent(b, 0);
806			}
807			sum_wt += weights_ ->GetBinContent(b) * den->GetBinContent(b);
808			cout<<"intimepileweight "<< b-1 <<" "<< weights_ ->GetBinContent(b)<<endl;
809			}
810			cout<<"sum wt "<< sum_wt <<" "<< weights_->GetMean()<<endl;
811
812			LumiReWeighting_weightOOT_init();
813
814			}
815
816
817
818			///Now use this version
819			vector<double> LumiReWeighting_getWeights(TH1D hdata, TH1D hmc){
820			TH1D num = (TH1D)hdata->Clone("datawt");
821			double total = num->Integral();
822			num ->Scale( 1.0/total);
823
824			TH1D den = (TH1D)hmc->Clone("mcwt");
825			double total1 = den->Integral();
826			den ->Scale(1.0/total1);
827
828			int nbindata = num->GetNbinsX();
829			int nbinmc = den->GetNbinsX();
830
831			if( nbinmc != nbindata ){
832			cout<<" waring! nbinmc != nbindata " << nbinmc <<" "<< nbindata <<endl;
833			}
834
835
836			vector<double> wts;
837			int nbin = nbinmc < nbindata ? nbinmc: nbindata;
838
839			for(int b=1; b<= nbin; b++){
840			if( den->GetBinContent(b) > 0){
841			wts.push_back( num->GetBinContent(b)/ den->GetBinContent(b));
842			}else{
843			wts.push_back(0);
844			}
845			cout<<" LumiReWeighting_getWeights " << wts[b-1] <<endl;
846
847			}
848			if( nbinmc > nbindata){
849			for(int b= nbindata; b < nbinmc; b++){
850			wts.push_back(0);
851			}
852			}
853
854			//final check
855			if( int(wts.size()) != nbinmc){
856			cout<<" problem of seeting pileup wts for MC " << wts.size() <<" "<< nbinmc <<endl;
857			exit(1);
858			}
859
860			return wts;
861
862
863			///LumiReWeighting_weightOOT_init();
864
865			}
866
867
868
869
870
871			double LumiReWeighting_weightINT_v1(vector<double> weight){
872
873			if( weight.size() <=0){
874			cout<<" LumiReWeighting_weightINT_v1 bad weight vector " << endl;
875			exit(1);
876			}
877
878			int npv = -1;
879			for(unsigned int j =0; j < pileupBunchX->size(); j++){
880			int BX = pileupBunchX->at(j);
881			if(BX==0){
882			npv = pileupNInteraction->at(j);
883			break;
884			}
885			}
886			if(npv < 0) {
887			cout << " no in-time beam crossing found\n! " <<endl;
888			return 0;
889			}
890			if(npv > int(weight.size()) ) {
891			cout<<"warning! LumiReWeighting_weightINT_v1 " << npv <<" "<< weight.size() <<endl;
892			return 0;
893			}
894			return weight[npv];
895
896			}
897
898
899
900
901			double LumiReWeighting_weightINT_v2(int npv, vector<double> weight){
902
903			if( weight.size() <=0){
904			cout<<" LumiReWeighting_weightINT_v2 bad weight vector " << endl;
905			exit(1);
906			}
907
908			// int npv = -1;
909			// for(unsigned int j =0; j < pileupBunchX->size(); j++){
910			// int BX = pileupBunchX->at(j);
911			// if(BX==0){
912			// npv = pileupNInteraction->at(j);
913			// break;
914			// }
915			// }
916			if(npv < 0) {
917			cout << " no in-time beam crossing found\n! " <<endl;
918			return 0;
919			}
920			if(npv > int(weight.size()) ) {
921			cout<<"warning! LumiReWeighting_weightINT_v2 " << npv <<" "<< weight.size() <<endl;
922			return 0;
923			}
924			return weight[npv];
925
926			}
927
928
929
930
931
932			double LumiReWeighting_weightINT(){
933
934			int npv = -1;
935			for(unsigned int j =0; j < pileupBunchX->size(); j++){
936			int BX = pileupBunchX->at(j);
937			if(BX==0){
938			npv = pileupNInteraction->at(j);
939			break;
940			}
941			}
942			if(npv < 0) cout << " no in-time beam crossing found\n! " <<endl;
943
944			if( weights_ == NULL){
945			cout<<" weights histogram not filled"<<endl;
946			return 1;
947			}
948			int bin = weights_->GetXaxis()->FindBin( npv );
949			double inTimeWeight = weights_->GetBinContent( bin );
950			return inTimeWeight;
951
952			}
953
954
955			double LumiReWeighting_weightOOT(string dataset) {
956
957			if( dataset.find("SIM") == string::npos){
958			cout<<" weight not on MC ? " << dataset.c_str() <<endl;
959			exit(1);
960			}
961
962
963			static double Correct_Weights2011[25] = { // residual correction to match lumi spectrum
964			5.30031,
965			2.07903,
966			1.40729,
967			1.27687,
968			1.0702,
969			0.902094,
970			0.902345,
971			0.931449,
972			0.78202,
973			0.824686,
974			0.837735,
975			0.910261,
976			1.01394,
977			1.1599,
978			1.12778,
979			1.58423,
980			1.78868,
981			1.58296,
982			2.3291,
983			3.86641,
984			0,
985			0,
986			0,
987			0,
988			0
989			};
990
991
992			int npv = -1;
993			int npv50ns = -1;
994			for(unsigned int j =0; j < pileupBunchX->size(); j++){
995			int BX = pileupBunchX->at(j);
996			if(BX==0){
997			npv = pileupNInteraction->at(j);
998			}
999			if(BX == 1) {
1000			npv50ns = pileupNInteraction->at(j);
1001			}
1002
1003			}
1004
1005			// Note: for the "uncorrelated" out-of-time pileup, reweighting is only done on the 50ns
1006			// "late" bunch (BX=+1), since that is basically the only one that matters in terms of
1007			// energy deposition.
1008
1009			if(npv < 0) {
1010			std::cerr << " no in-time beam crossing found\n! " ;
1011			std::cerr << " Returning event weight=0\n! ";
1012			return 0.;
1013			}
1014			// if(npv50ns < 0) {
1015			// std::cerr << " no out-of-time beam crossing found\n! " ;
1016			// std::cerr << " Returning event weight=0\n! ";
1017			// return 0.;
1018			// }
1019
1020			int bin = weights_->GetXaxis()->FindBin( npv );
1021
1022			double inTimeWeight = weights_->GetBinContent( bin );
1023
1024			double TotalWeight = 1.0;
1025
1026
1027			if( dataset.find("S3") != string::npos){ //if(Reweight_4_2_2p2_) {
1028			TotalWeight = inTimeWeight * WeightOOTPU_[bin-1][npv50ns] * Correct_Weights2011[bin-1];
1029			}
1030			else {
1031			TotalWeight = inTimeWeight;
1032			}
1033
1034			// if( TotalWeight < 1E-5 && nWarning < 10 ){
1035			// cout<<"warnign 0 weight " << TotalWeight <<" "<< inTimeWeight <<" "<< bin <<" "<< npv<<" "<< npv50ns <<" "<< WeightOOTPU_[bin-1][npv50ns] <<" "<<Correct_Weights2011[bin-1] <<endl;
1036			// nWarning ++;
1037			// }
1038
1039
1040			return TotalWeight;
1041
1042			}
1043
1044
1045
1046
1047			double LumiReWeighting_weightOOT_givenWeightVector(string dataset, vector<double> puwts) {
1048
1049			if( dataset.find("SIM") == string::npos){
1050			cout<<" weight not on MC ? " << dataset.c_str() <<endl;
1051			exit(1);
1052			}
1053
1054
1055			static double Correct_Weights2011[25] = { // residual correction to match lumi spectrum
1056			5.30031,
1057			2.07903,
1058			1.40729,
1059			1.27687,
1060			1.0702,
1061			0.902094,
1062			0.902345,
1063			0.931449,
1064			0.78202,
1065			0.824686,
1066			0.837735,
1067			0.910261,
1068			1.01394,
1069			1.1599,
1070			1.12778,
1071			1.58423,
1072			1.78868,
1073			1.58296,
1074			2.3291,
1075			3.86641,
1076			0,
1077			0,
1078			0,
1079			0,
1080			0
1081			};
1082
1083
1084			int npv = -1;
1085			int npv50ns = -1;
1086			for(unsigned int j =0; j < pileupBunchX->size(); j++){
1087			int BX = pileupBunchX->at(j);
1088			if(BX==0){
1089			npv = pileupNInteraction->at(j);
1090			}
1091			if(BX == 1) {
1092			npv50ns = pileupNInteraction->at(j);
1093			}
1094
1095			}
1096
1097			// Note: for the "uncorrelated" out-of-time pileup, reweighting is only done on the 50ns
1098			// "late" bunch (BX=+1), since that is basically the only one that matters in terms of
1099			// energy deposition.
1100
1101			if(npv < 0) {
1102			std::cerr << " no in-time beam crossing found\n! " ;
1103			std::cerr << " Returning event weight=0\n! ";
1104			return 0.;
1105			}
1106			// if(npv50ns < 0) {
1107			// std::cerr << " no out-of-time beam crossing found\n! " ;
1108			// std::cerr << " Returning event weight=0\n! ";
1109			// return 0.;
1110			// }
1111
1112
1113
1114			double inTimeWeight = puwts[npv];
1115
1116			double TotalWeight = 1.0;
1117
1118
1119
1120			if( dataset.find("S3") != string::npos){ //if(Reweight_4_2_2p2_) {
1121			TotalWeight = inTimeWeight * WeightOOTPU_[npv][npv50ns] * Correct_Weights2011[npv];
1122			}
1123			else {
1124			TotalWeight = inTimeWeight;
1125			}
1126
1127			// if( TotalWeight < 1E-5 && nWarning < 10 ){
1128			// cout<<"warnign 0 weight " << TotalWeight <<" "<< inTimeWeight <<" "<< bin <<" "<< npv<<" "<< npv50ns <<" "<< WeightOOTPU_[bin-1][npv50ns] <<" "<<Correct_Weights2011[bin-1] <<endl;
1129			// nWarning ++;
1130			// }
1131
1132
1133			return TotalWeight;
1134
1135			}
1136