LumiReWeighting.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_interface_LumiReWeighting_h
#define PhysicsTools_Utilities_interface_LumiReWeighting_h


#include "TH1.h"
#include "TFile.h"
#include <cmath>
#include <string>
#include <boost/shared_ptr.hpp>
#include <vector>
#include "FWCore/Utilities/interface/InputTag.h" 

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.

      constexpr double p0_minus[20] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
      constexpr 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
      };
      constexpr 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
      };

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

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


      constexpr double p0_plus[20] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
      constexpr 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
      };
      constexpr 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
      };

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

      constexpr 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 ) {
      return ShiftWeight(int(pvnum));
    };

  private:

    double Pweight_[25];

  };



} // end of reweight namespace


namespace edm {
  class EventBase;
  class LumiReWeighting {
  public:
    LumiReWeighting( std::string generatedFile,
             std::string dataFile,
             std::string GenHistName = "pileup",
             std::string DataHistName = "pileup",
             const edm::InputTag& PileupSumInfoInputTag = edm::InputTag( "addPileupInfo" ) );
    
    LumiReWeighting( const std::vector< float >& MC_distr, const std::vector< float >& Lumi_distr, const edm::InputTag& PileupSumInfoInputTag = edm::InputTag( "addPileupInfo" ) );

    LumiReWeighting ( ) { } ;

    double weight( int npv ) ;

    double weight( float npv ) ;

    double weight( const edm::EventBase &e ) ;

    double weightOOT( const edm::EventBase &e ) ;

  protected:

    std::string generatedFileName_;
    std::string dataFileName_;
    std::string GenHistName_;
    std::string DataHistName_;
    edm::InputTag pileupSumInfoTag_;
    boost::shared_ptr<TFile>     generatedFile_;
    boost::shared_ptr<TFile>     dataFile_;
    boost::shared_ptr<TH1>      weights_;

    //keep copies of normalized distributions:

    boost::shared_ptr<TH1>      MC_distr_;
    boost::shared_ptr<TH1>      Data_distr_;

    int  OldLumiSection_;
    bool FirstWarning_;


  };
}



#endif