/data/refman/pasoursint/CMSSW_6_1_1/src/PhysicsTools/Utilities/interface/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>

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,
                     std::string DataHistName);
    
    LumiReWeighting( std::vector< float > MC_distr, std::vector< float > Lumi_distr);

    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_;
    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