#include <iostream>
#include <math.h>

// Local
#include "CorrectCaloIso.h"

//////////////////////////////////////////////////
// To get the corrected relative Isolation call:
// CorrectEtCone30Rel(pt, isoCone30Rel, nvtx, eta)
//    -- author: Doug Schaefer <schae@cern.ch>
//////////////////////////////////////////////////

using namespace std;

//-----------------------------------------------------------------------------
CorrectCaloIso::CorrectCaloIso() : fDebug(false)
{
  //
  // Configuring the eta dependent Isolation corrections
  //
  ConfigStrawmanIso();
}

//-----------------------------------------------------------------------------
CorrectCaloIso::~CorrectCaloIso()
{
}

//-----------------------------------------------------------------------------
void CorrectCaloIso::ConfigStrawmanIso()
{ //                   
  // Eta dependent isolation from Peter's Ztag and probe  
  //   -- 11/7/2011                                      
  //   -- units are MeV  
  //                                                                                                                                               
  fSlope[0] =81.4380;  fConstant[0] =256.798;  //fEta[0] =0.10; 
  fSlope[1] =88.7883;  fConstant[1] =222.295;  //fEta[1] =0.20; 
  fSlope[2] =90.4055;  fConstant[2] =246.810;  //fEta[2] =0.30; 
  fSlope[3] =88.5415;  fConstant[3] =256.253;  //fEta[3] =0.40; 
  fSlope[4] =89.8105;  fConstant[4] =253.627;  //fEta[4] =0.50; 
  fSlope[5] =87.2666;  fConstant[5] =269.239;  //fEta[5] =0.60; 
  fSlope[6] =89.8660;  fConstant[6] =231.288;  //fEta[6] =0.70; 
  fSlope[7] =88.4182;  fConstant[7] =219.546;  //fEta[7] =0.80; 
  fSlope[8] =83.7103;  fConstant[8] =238.926;  //fEta[8] =0.90; 
  fSlope[9] =88.3318;  fConstant[9] =193.716;  //fEta[9] =1.00; 
  fSlope[10]=85.1312;  fConstant[10]=207.372;  //fEta[10]=1.10; 
  fSlope[11]=86.0413;  fConstant[11]=306.257;  //fEta[11]=1.20; 
  fSlope[12]=81.6753;  fConstant[12]=250.002;  //fEta[12]=1.30; 
  fSlope[13]=78.6961;  fConstant[13]=377.195;  //fEta[13]=1.40; 
  fSlope[14]=66.8804;  fConstant[14]=259.919;  //fEta[14]=1.50; 
  fSlope[15]=71.9913;  fConstant[15]=187.506;  //fEta[15]=1.60; 
  fSlope[16]=73.2569;  fConstant[16]=275.101;  //fEta[16]=1.70; 
  fSlope[17]=78.2649;  fConstant[17]=165.755;  //fEta[17]=1.80; 
  fSlope[18]=77.9631;  fConstant[18]=199.264;  //fEta[18]=1.90; 
  fSlope[19]=79.0986;  fConstant[19]=198.490;  //fEta[19]=2.00; 
  fSlope[20]=85.7765;  fConstant[20]=153.297;  //fEta[20]=2.10; 
  fSlope[21]=77.4238;  fConstant[21]=177.392;  //fEta[21]=2.20; 
  fSlope[22]=75.9160;  fConstant[22]=120.396;  //fEta[22]=2.30; 
  fSlope[23]=70.5277;  fConstant[23]=123.942;  //fEta[23]=2.40; 

}

//-----------------------------------------------------------------------------    
unsigned CorrectCaloIso::GetEtaBin(const float &eta) const
{
  //                           
  // Find Eta bin
  //
  unsigned bin = static_cast<unsigned>(fabs(eta)/0.10);

  if(bin<24)  return bin;
  else{
    std::cout << "ERROR bin is bigger than known isolation corrections!" << std::endl;
    return 23;
  }
}

//-----------------------------------------------------------------------------
float CorrectCaloIso::CorrectEtCone30Rel(const float &pt, const float &isoCone30Rel, const float &nvtx, const float &eta)
{
  // 
  // Calculate new relative isolation   
  //    -- pt           is in MeV and is the combined muon pt                 
  //    -- isoCone30Rel is the EtCone30/pt
  //    -- nVtx         is number of primary vertices as calculated in the HSG3 recommendation
  //    -- eta          is the combined muon eta
  //
  if(pt==0.0) return isoCone30Rel;
  float isoCone30 = isoCone30Rel*pt;

  //
  // Calculate the cone corrections  
  //   -- (etcone - correction)/pt
  //
  float newIsoCone = (isoCone30 - GetCorrectionEtCone30(nvtx,eta))/pt; 

  if(fDebug){
    std::cout << "Uncorrected Relative Isolation Cone 30:  " << isoCone30Rel
	      << " Pileup Corrected Relative Isolation 30: " << newIsoCone 
	      << std::endl;
  }

  return newIsoCone;
}

//-----------------------------------------------------------------------------
float CorrectCaloIso::GetCorrectionEtCone30(const unsigned &nvtx,const float &eta)
{
  const unsigned bin     = GetEtaBin(eta);
  const float slope      = fSlope   [bin];
  const float constant   = fConstant[bin];

  return CaloCor(nvtx, slope, constant);
}

//-----------------------------------------------------------------------------
float CorrectCaloIso::CaloCor(const unsigned &nvtx, const float &slope, const float &constant)
{
  //
  // Returns the correction to the calorimeter isolation
  //   -- this is the amount to subtract from the corrected cone
  //
  return float(nvtx)*slope + constant;
}