text: 
  #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
  #include "CondFormats/PhysicsToolsObjects/interface/PerformancePayloadFromTFormula.h"
  
  #include "CondFormats/ESObjects/interface/ESEEIntercalibConstants.h"
  
  #include <TMath.h>
  #include <cmath>
  #include <vector>
  #include <TF1.h>
  #include <map>
  #include <algorithm>
  #include <numeric>
  
  using namespace std;
  
  PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(nullptr), esEEInterCalib_(nullptr)
  {
    initializeCalibrationFunctions();
  }
  
  PFEnergyCalibration::~PFEnergyCalibration() 
  {
  
  }
  
  void
  PFEnergyCalibration::initializeCalibrationFunctions() {
  
    threshE = 3.5;
    threshH = 2.5;
  
  
   //calibChrisClean.C calibration parameters bhumika Nov, 2018
    faBarrel = std::make_unique<TF1>("faBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
    faBarrel->SetParameter(0,-30.7141);
    faBarrel->SetParameter(1,31.7583);
    faBarrel->SetParameter(2,4.40594);
    faBarrel->SetParameter(3,1.70914);
    faBarrel->SetParameter(4,0.0613696);
    faBarrel->SetParameter(5,0.000104857);
    faBarrel->SetParameter(6,-1.38927);
    faBarrel->SetParameter(7,-0.743082);
    fbBarrel = std::make_unique<TF1>("fbBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
    fbBarrel->SetParameter(0,2.25366);	
    fbBarrel->SetParameter(1,0.537715);
    fbBarrel->SetParameter(2,-4.81375);
    fbBarrel->SetParameter(3,12.109);
    fbBarrel->SetParameter(4,1.80577);
    fbBarrel->SetParameter(5,0.187919);
    fbBarrel->SetParameter(6,-6.26234);
    fbBarrel->SetParameter(7,-0.607392);
    fcBarrel = std::make_unique<TF1>("fcBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
    fcBarrel->SetParameter(0,1.5125962);
    fcBarrel->SetParameter(1,0.855057);
    fcBarrel->SetParameter(2,-6.04199);
    fcBarrel->SetParameter(3,2.08229);
    fcBarrel->SetParameter(4,0.592266);
    fcBarrel->SetParameter(5,0.0291232);
    fcBarrel->SetParameter(6,0.364802);
    fcBarrel->SetParameter(7,-1.50142);
    faEtaBarrelEH = std::make_unique<TF1>("faEtaBarrelEH","[0]+[1]*exp(-x/[2])",0.,1000.);
    faEtaBarrelEH->SetParameter(0,0.0185555);
    faEtaBarrelEH->SetParameter(1,-0.0470674);
    faEtaBarrelEH->SetParameter(2,396.959);
    fbEtaBarrelEH = std::make_unique<TF1>("fbEtaBarrelEH","[0]+[1]*exp(-x/[2])",0.,1000.);
    fbEtaBarrelEH->SetParameter(0,0.0396458);
    fbEtaBarrelEH->SetParameter(1,0.114128);
    fbEtaBarrelEH->SetParameter(2,251.405);
    faEtaBarrelH = std::make_unique<TF1>("faEtaBarrelH","[0]+[1]*x",0.,1000.);
    faEtaBarrelH->SetParameter(0,0.00434994);
    faEtaBarrelH->SetParameter(1,-5.16564e-06);
    fbEtaBarrelH = std::make_unique<TF1>("fbEtaBarrelH","[0]+[1]*exp(-x/[2])",0.,1000.);
    fbEtaBarrelH->SetParameter(0,-0.0232604);
    fbEtaBarrelH->SetParameter(1,0.0937525);
    fbEtaBarrelH->SetParameter(2,34.9935);
  
    faEndcap = std::make_unique<TF1>("faEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
    faEndcap->SetParameter(0,1.17227);
    faEndcap->SetParameter(1,13.1489);
    faEndcap->SetParameter(2,-29.1672);
    faEndcap->SetParameter(3,0.604223);
    faEndcap->SetParameter(4,0.0426363);
    faEndcap->SetParameter(5,3.30898e-15);
    faEndcap->SetParameter(6,0.165293);
    faEndcap->SetParameter(7,-7.56786);
    fbEndcap = std::make_unique<TF1>("fbEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
    fbEndcap->SetParameter(0,-0.974251);
    fbEndcap->SetParameter(1,1.61733);
    fbEndcap->SetParameter(2,0.0629183);
    fbEndcap->SetParameter(3,7.78495);
    fbEndcap->SetParameter(4,-0.774289);
    fbEndcap->SetParameter(5,7.81399e-05);
    fbEndcap->SetParameter(6,0.139116);
    fbEndcap->SetParameter(7,-4.25551);
    fcEndcap = std::make_unique<TF1>("fcEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
    fcEndcap->SetParameter(0,1.01863);
    fcEndcap->SetParameter(1,1.29787);
    fcEndcap->SetParameter(2,-3.97293);
    fcEndcap->SetParameter(3,21.7805);
    fcEndcap->SetParameter(4,0.810195);
    fcEndcap->SetParameter(5,0.234134);
    fcEndcap->SetParameter(6,1.42226);
    fcEndcap->SetParameter(7,-0.0997326);
    faEtaEndcapEH = std::make_unique<TF1>("faEtaEndcapEH","[0]+[1]*exp(-x/[2])",0.,1000.);
    faEtaEndcapEH->SetParameter(0,0.0112692);
    faEtaEndcapEH->SetParameter(1,-2.68063);
    faEtaEndcapEH->SetParameter(2,2.90973);
    fbEtaEndcapEH = std::make_unique<TF1>("fbEtaEndcapEH","[0]+[1]*exp(-x/[2])",0.,1000.);
    fbEtaEndcapEH->SetParameter(0,-0.0192991);
    fbEtaEndcapEH->SetParameter(1,-0.265);
    fbEtaEndcapEH->SetParameter(2,80.5502);
    faEtaEndcapH = std::make_unique<TF1>("faEtaEndcapH","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",0.,1000.);
    faEtaEndcapH->SetParameter(0,-0.0106029);
    faEtaEndcapH->SetParameter(1,-0.692207);
    faEtaEndcapH->SetParameter(2,0.0542991);
    faEtaEndcapH->SetParameter(3,-0.171435);
    faEtaEndcapH->SetParameter(4,-61.2277);
    fbEtaEndcapH = std::make_unique<TF1>("fbEtaEndcapH","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",0.,1000.);
    fbEtaEndcapH->SetParameter(0,0.0214894);
    fbEtaEndcapH->SetParameter(1,-0.266704);
    fbEtaEndcapH->SetParameter(2,5.2112);
    fbEtaEndcapH->SetParameter(3,0.303578);
    fbEtaEndcapH->SetParameter(4,-104.367);
  
    //added by Bhumika on 2 august 2018
  
    fcEtaBarrelH = std::make_unique<TF1>("fcEtaBarrelH","[3]*((x-[0])^[1])+[2]",0.,1000.);
    fcEtaBarrelH->SetParameter(0,0);
    fcEtaBarrelH->SetParameter(1,2);
    fcEtaBarrelH->SetParameter(2,0);
    fcEtaBarrelH->SetParameter(3,1);
  
    fcEtaEndcapH = std::make_unique<TF1>("fcEtaEndcapH","[3]*((x-[0])^[1])+[2]",0.,1000.);
    fcEtaEndcapH->SetParameter(0,0);
    fcEtaEndcapH->SetParameter(1,0);
    fcEtaEndcapH->SetParameter(2,0.05);
    fcEtaEndcapH->SetParameter(3,0);
    
    fdEtaEndcapH = std::make_unique<TF1>("fdEtaEndcapH","[3]*((x-[0])^[1])+[2]",0.,1000.);
    fdEtaEndcapH->SetParameter(0,1.5);
    fdEtaEndcapH->SetParameter(1,4);
    fdEtaEndcapH->SetParameter(2,-1.1);
    fdEtaEndcapH->SetParameter(3,1.0);
  
    fcEtaBarrelEH = std::make_unique<TF1>("fcEtaBarrelEH","[3]*((x-[0])^[1])+[2]",0.,1000.);
    fcEtaBarrelEH->SetParameter(0,0);
    fcEtaBarrelEH->SetParameter(1,2);
    fcEtaBarrelEH->SetParameter(2,0);
    fcEtaBarrelEH->SetParameter(3,1);
  
    fcEtaEndcapEH = std::make_unique<TF1>("fcEtaEndcapEH","[3]*((x-[0])^[1])+[2]",0.,1000.);
    fcEtaEndcapEH->SetParameter(0,0);
    fcEtaEndcapEH->SetParameter(1,0);
    fcEtaEndcapEH->SetParameter(2,0);
    fcEtaEndcapEH->SetParameter(3,0);
    
    fdEtaEndcapEH = std::make_unique<TF1>("fdEtaEndcapEH","[3]*((x-[0])^[1])+[2]",0.,1000.);
    fdEtaEndcapEH->SetParameter(0,1.5);
    fdEtaEndcapEH->SetParameter(1,2.0);
    fdEtaEndcapEH->SetParameter(2,0.6);
    fdEtaEndcapEH->SetParameter(3,1.0);
    
  
  
  }
  
  void 
  PFEnergyCalibration::energyEmHad(double t, double& e, double&h, double eta, double phi) const { 
  
    // Use calorimetric energy as true energy for neutral particles
    double tt = t;
    double ee = e;
    double hh = h;
    double a = 1.;
    double b = 1.;
    double etaCorrE = 1.;
    double etaCorrH = 1.;
    auto absEta = std::abs(eta);
    t = min(999.9,max(tt,e+h));
    if ( t < 1. ) return;
  
    // Barrel calibration
    if ( absEta < 1.48 ) { 
      // The energy correction
      a = e>0. ? aBarrel(t) : 1.;
      b = e>0. ? bBarrel(t) : cBarrel(t);
      double thresh = e > 0. ? threshE : threshH;
  
      // Protection against negative calibration
      if ( a < -0.25 || b < -0.25 ) { 
        a = 1.;
        b = 1.;
        thresh = 0.;
      }
  
      // The new estimate of the true energy
      t = min(999.9,max(tt, thresh+a*e+b*h));
  
      // The angular correction 
      if ( e > 0. && thresh > 0. ) {
        etaCorrE = 1.0 + aEtaBarrelEH(t) + 1.3*bEtaBarrelEH(t)*cEtaBarrelEH(absEta);
        etaCorrH = 1.0;
      } else {
        etaCorrE = 1.0 + aEtaBarrelH(t) + 1.3*bEtaBarrelH(t)*cEtaBarrelH(absEta);
        etaCorrH = 1.0 + aEtaBarrelH(t) + bEtaBarrelH(t)*cEtaBarrelH(absEta);
  
  
      }
      if ( e > 0. && thresh > 0. ) 
        e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
      if ( h > 0. && thresh > 0. ) {
        h = threshH + etaCorrH * b * h;
      }
  
    // Endcap calibration   
    } else {
      // The energy correction
      a = e>0. ? aEndcap(t) : 1.;
      b = e>0. ? bEndcap(t) : cEndcap(t);
      double thresh = e > 0. ? threshE : threshH;
  
      // Protection against negative calibration
      if ( a < -0.25 || b < -0.25 ) { 
        a = 1.;
        b = 1.;
        thresh = 0.;
      }
  
      // The new estimate of the true energy
      t = min(999.9,max(tt, thresh+a*e+b*h));
      
      // The angular correction 
      double dEta = std::abs( absEta - 1.5 );
      double etaPow = dEta * dEta * dEta * dEta;
  
  
      if ( e > 0. && thresh > 0. ) {
        if(absEta<2.5) {
          etaCorrE = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t)*cEtaEndcapEH(absEta);
        }
        else {
  
          etaCorrE = 1. + aEtaEndcapEH(t) + 1.3*bEtaEndcapEH(t)*dEtaEndcapEH(absEta);
        }
  
          etaCorrH = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t)*(0.04 + etaPow);
      } else {
        etaCorrE = 1.;
        if(absEta<2.5) {
          etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t)*cEtaEndcapH(absEta);
        }
        else {
          etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t)*dEtaEndcapH(absEta);
        }
      }
  
      //t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
  
      if ( e > 0. && thresh > 0. ) 
        e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
      if ( h > 0. && thresh > 0. ) {
        h = threshH + etaCorrH * b * h;
      }
    }
  
    // Protection
    if ( e < 0. || h < 0. ) {
    
      // Some protection against crazy calibration
      if ( e < 0. ) e = ee;
      if ( h < 0. ) h = hh;
    }
  
    // And that's it !
  
    
  }
  
  // The calibration functions
  double 
  PFEnergyCalibration::aBarrel(double x) const { 
  
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfa_BARREL,point); 
  
    } else { 
  
      return faBarrel->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::bBarrel(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfb_BARREL,point); 
  
    } else { 
  
      return fbBarrel->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::cBarrel(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfc_BARREL,point); 
  
    } else { 
  
      return fcBarrel->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::aEtaBarrelEH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARRELEH,point); 
  
    } else { 
  
      return faEtaBarrelEH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::bEtaBarrelEH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARRELEH,point); 
  
    } else { 
  
      return fbEtaBarrelEH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::aEtaBarrelH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARRELH,point); 
  
    } else { 
  
      return faEtaBarrelH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::bEtaBarrelH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARRELH,point); 
  
    } else { 
  
      return fbEtaBarrelH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::aEndcap(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfa_ENDCAP,point); 
  
    } else { 
  
      return faEndcap->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::bEndcap(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfb_ENDCAP,point); 
  
    } else { 
  
      return fbEndcap->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::cEndcap(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfc_ENDCAP,point); 
  
    } else { 
  
      return fcEndcap->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::aEtaEndcapEH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAPEH,point); 
  
    } else { 
  
      return faEtaEndcapEH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::bEtaEndcapEH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAPEH,point); 
  
    } else { 
  
      return fbEtaEndcapEH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::aEtaEndcapH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAPH,point); 
  
    } else { 
  
      return faEtaEndcapH->Eval(x); 
  
    }
  }
  
  double 
  PFEnergyCalibration::bEtaEndcapH(double x) const { 
  
    if ( pfCalibrations ) { 
  
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAPH,point); 
    } else { 
      return fbEtaEndcapH->Eval(x); 
  
    }
  }
  
  //added by Bhumika Kansal on 3 august 2018
  
  double 
  PFEnergyCalibration::cEtaBarrelH(double x) const { 
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfcEta_BARRELH,point); 
      
    } else { 
      return fcEtaBarrelH->Eval(x); 
    }
  }
  double 
  PFEnergyCalibration::cEtaEndcapH(double x) const { 
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfcEta_ENDCAPH,point); 
      
    } else{  
  
      return fcEtaEndcapH->Eval(x); 
    }
  }
  
  double 
  PFEnergyCalibration::dEtaEndcapH(double x) const { 
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfdEta_ENDCAPH,point); 
      
    } else{ 
  
      return fdEtaEndcapH->Eval(x); 
    }
  }
  
  double 
  PFEnergyCalibration::cEtaBarrelEH(double x) const { 
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfcEta_BARRELEH,point); 
      
    } else{  
  
      return fcEtaBarrelEH->Eval(x); 
    }
  }
  
  double 
  PFEnergyCalibration::cEtaEndcapEH(double x) const { 
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfcEta_ENDCAPEH,point); 
      
    } else{  
  
      return fcEtaEndcapEH->Eval(x); 
    }
  }
  
  double 
  PFEnergyCalibration::dEtaEndcapEH(double x) const { 
    if ( pfCalibrations ) { 
      BinningPointByMap point;
      point.insert(BinningVariables::JetEt, x);
      return pfCalibrations->getResult(PerformanceResult::PFfdEta_ENDCAPEH,point); 
      
    } else{  
  
      return fdEtaEndcapEH->Eval(x); 
    }
  }
  //
  double
  PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
  			      std::vector<double> &EclustersPS1,
  			      std::vector<double> &EclustersPS2,
  			      bool crackCorrection ) const {
    double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
    double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
    return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
  }
  
  double
  PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
  			      double ePS1,
  			      double ePS2,
  			      bool crackCorrection ) const {
    double eEcal = clusterEcal.energy();
    //temporaty ugly fix
    reco::PFCluster myPFCluster=clusterEcal;
    myPFCluster.calculatePositionREP();
    double eta = myPFCluster.positionREP().eta();
    double phi = myPFCluster.positionREP().phi();
  
    double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
    // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<"  eta = "<<eta<<"  phi = "<<phi<<std::endl; 
    return calibrated; 
  }
  
  double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
  				     std::vector<double> &EclustersPS1,
  				     std::vector<double> &EclustersPS2,
  				     double& ps1,double& ps2,
  				     bool crackCorrection) const {
    double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
    double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
    return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
  }
  double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
  				     double ePS1, double ePS2,
  				     double& ps1,double& ps2,
  				     bool crackCorrection) const {
    double eEcal = clusterEcal.energy();
    //temporaty ugly fix
    reco::PFCluster myPFCluster=clusterEcal;
    myPFCluster.calculatePositionREP();
    double eta = myPFCluster.positionREP().eta();
    double phi = myPFCluster.positionREP().phi();
  
    double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
    // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<"  eta = "<<eta<<"  phi = "<<phi<<std::endl; 
    return calibrated; 
  }
  
  
  std::ostream& operator<<(std::ostream& out,
  			 const PFEnergyCalibration& calib) {
  
    if(!out ) return out;
  
    out<<"PFEnergyCalibration -- "<<endl;
  
    if ( calib.pfCalibrations ) { 
  
      static const std::map<std::string, PerformanceResult::ResultType> functType = {
        {"PFfa_BARREL", PerformanceResult::PFfa_BARREL},
        {"PFfa_ENDCAP", PerformanceResult::PFfa_ENDCAP},
        {"PFfb_BARREL", PerformanceResult::PFfb_BARREL},
        {"PFfb_ENDCAP", PerformanceResult::PFfb_ENDCAP},
        {"PFfc_BARREL", PerformanceResult::PFfc_BARREL},
        {"PFfc_ENDCAP", PerformanceResult::PFfc_ENDCAP},
        {"PFfaEta_BARRELH", PerformanceResult::PFfaEta_BARRELH},
        {"PFfaEta_ENDCAPH", PerformanceResult::PFfaEta_ENDCAPH},
        {"PFfbEta_BARRELH", PerformanceResult::PFfbEta_BARRELH},
        {"PFfbEta_ENDCAPH", PerformanceResult::PFfbEta_ENDCAPH},
        {"PFfaEta_BARRELEH", PerformanceResult::PFfaEta_BARRELEH},
        {"PFfaEta_ENDCAPEH", PerformanceResult::PFfaEta_ENDCAPEH},
        {"PFfbEta_BARRELEH", PerformanceResult::PFfbEta_BARRELEH},
        {"PFfbEta_ENDCAPEH", PerformanceResult::PFfbEta_ENDCAPEH},
        {"PFfcEta_BARRELH", PerformanceResult::PFfcEta_BARRELH},
        {"PFfcEta_ENDCAPH", PerformanceResult::PFfcEta_ENDCAPH},
        {"PFfdEta_ENDCAPH", PerformanceResult::PFfdEta_ENDCAPH},
        {"PFfcEta_BARRELEH", PerformanceResult::PFfcEta_BARRELEH},
        {"PFfcEta_ENDCAPEH", PerformanceResult::PFfcEta_ENDCAPEH},
        {"PFfdEta_ENDCAPEH", PerformanceResult::PFfdEta_ENDCAPEH}
  
      };
  
          
      for(std::map<std::string,PerformanceResult::ResultType>::const_iterator 
  	  func = functType.begin(); 
          func != functType.end(); 
          ++func) {    
        
        cout << "Function: " << func->first << endl;
        PerformanceResult::ResultType fType = func->second;
        calib.pfCalibrations->printFormula(fType);
      }
  
    } else { 
      
      std::cout << "Default calibration functions : " << std::endl;
      
      calib.faBarrel->Print();
      calib.fbBarrel->Print();
      calib.fcBarrel->Print();
      calib.faEtaBarrelEH->Print();
      calib.fbEtaBarrelEH->Print();
      calib.faEtaBarrelH->Print();
      calib.fbEtaBarrelH->Print();
      calib.faEndcap->Print();
      calib.fbEndcap->Print();
      calib.fcEndcap->Print();
      calib.faEtaEndcapEH->Print();
      calib.fbEtaEndcapEH->Print();
      calib.faEtaEndcapH->Print();
      calib.fbEtaEndcapH->Print();
      //
  
    }
      
    return out;
  }
  
  
  
  
  ///////////////////////////////////////////////////////////////
  ////                                                       ////  
  ////             CORRECTION OF PHOTONS' ENERGY             ////
  ////                                                       ////
  ////              Material effect: No tracker              ////
  ////       Tuned on CMSSW_2_1_0_pre4, Full Sim events      ////
  ////                                                       ////
  ///////////////////////////////////////////////////////////////
  ////                                                       ////
  ////            Jonathan Biteau - June 2008                ////
  ////                                                       ////
  ///////////////////////////////////////////////////////////////
  
  
  
  ///////////////////////////////////////////////////////////////
  ////                                                       ////  
  ////  USEFUL FUNCTIONS FOR THE CORRECTION IN THE BARREL    ////
  ////                                                       ////
  ///////////////////////////////////////////////////////////////
  
  
  //useful to compute the signed distance to the closest crack in the barrel
  double
  PFEnergyCalibration::minimum(double a,double b) const {
    if(TMath::Abs(b)<TMath::Abs(a)) a=b;
    return a;
  }
  
  namespace {
    constexpr double pi= M_PI;// 3.14159265358979323846;
  
  
    std::vector<double> fillcPhi() {
      std::vector<double> retValue;
      retValue.resize(18,0);
      retValue[0]=2.97025;
      for(unsigned i=1;i<=17;++i) retValue[i]=retValue[0]-2*i*pi/18;
      
      return retValue;
    }
    
    //Location of the 18 phi-cracks
    const std::vector<double> cPhi = fillcPhi();
  }
  
  //compute the unsigned distance to the closest phi-crack in the barrel
  double
  PFEnergyCalibration::dCrackPhi(double phi, double eta) const {
  
    
    //Shift of this location if eta<0
    constexpr double delta_cPhi=0.00638;
  
    double m; //the result
  
    //the location is shifted
    if(eta<0) phi +=delta_cPhi;
  
    if (phi>=-pi && phi<=pi){
  
      //the problem of the extrema
      if (phi<cPhi[17] || phi>=cPhi[0]){
        if (phi<0) phi+= 2*pi;
        m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);        	
      }
  
      //between these extrema...
      else{
        bool OK = false;
        unsigned i=16;
        while(!OK){
  	if (phi<cPhi[i]){
  	  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
  	  OK=true;
  	}
  	else i-=1;
        }
      }
    }
    else{
      m=0.;        //if there is a problem, we assum that we are in a crack
      std::cout<<"Problem in dminphi"<<std::endl;
    }
    if(eta<0) m=-m;   //because of the disymetry
    return m;
  }
  
  // corrects the effect of phi-cracks
  double
  PFEnergyCalibration::CorrPhi(double phi, double eta) const {
  
    // we use 3 gaussians to correct the phi-cracks effect
    constexpr double p1=   5.59379e-01;
    constexpr double p2=   -1.26607e-03;
    constexpr double p3=  9.61133e-04;
  
    constexpr double p4=   1.81691e-01;
    constexpr double p5=   -4.97535e-03;
    constexpr double p6=   1.31006e-03;
  
    constexpr double p7=   1.38498e-01;
    constexpr double p8=   1.18599e-04;
    constexpr double p9= 2.01858e-03;
    
  
    double dminphi = dCrackPhi(phi,eta);
    
    double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
  
    return result;
  }   
  
  
  // corrects the effect of  |eta|-cracks
  double
  PFEnergyCalibration::CorrEta(double eta) const {
    
    // we use a gaussian with a screwness for each of the 5 |eta|-cracks
    constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01}; // amplitude
    constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; // mean
    constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04}; // sigma
    constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01}; // screwness amplitude
    constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03}; // screwness sigma
    double result = 1;
  
    for(unsigned i=0;i<=4;i++) result+=a[i]*TMath::Gaus(eta,m[i],s[i])*(1+sa[i]*TMath::Sign(1.,eta-m[i])*TMath::Exp(-TMath::Abs(eta-m[i])/ss[i]));
  
    return result;
  }
  
  
  //corrects the global behaviour in the barrel
  double
  PFEnergyCalibration::CorrBarrel(double E, double eta) const {
  
    //Energy dependency
    /*
    //YM Parameters 52XX:
    constexpr double p0=1.00000e+00;
    constexpr double p1=3.27753e+01;
    constexpr double p2=2.28552e-02;
    constexpr double p3=3.06139e+00;
    constexpr double p4=2.25135e-01;
    constexpr double p5=1.47824e+00;
    constexpr double p6=1.09e-02;
    constexpr double p7=4.19343e+01;
    */
    constexpr double p0 = 0.9944;
    constexpr double p1 = 9.827;
    constexpr double p2 = 1.503;
    constexpr double p3 = 1.196;
    constexpr double p4 = 0.3349;
    constexpr double p5 = 0.89;
    constexpr double p6 = 0.004361;
    constexpr double p7 = 51.51;
    //Eta dependency
    constexpr double p8=2.705593e-03;
    
    double result = (p0+1/(p1+p2*TMath::Power(E,p3))+p4*TMath::Exp(-E/p5)+p6*TMath::Exp(-E*E/(p7*p7)))*(1+p8*eta*eta);
  
    return result;
  }
  
  
  
  ///////////////////////////////////////////////////////////////
  ////                                                       ////  
  ////  USEFUL FUNCTIONS FOR THE CORRECTION IN THE ENDCAPS   ////
  ////  Parameters tuned for:                                ////
  ////          dR(ClustersPS1,ClusterEcal) < 0.08           ////
  ////          dR(ClustersPS2,ClusterEcal) < 0.13           ////
  ////                                                       ////
  ///////////////////////////////////////////////////////////////
  
  
  //Alpha, Beta, Gamma give the weight of each sub-detector (PS layer1, PS layer2 and Ecal) in the areas of the endcaps where there is a PS
  // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2) 
  
  double
  PFEnergyCalibration::Alpha(double eta) const {
  
    //Energy dependency
    constexpr double p0 = 5.97621e-01;
  
    //Eta dependency
    constexpr double p1 =-1.86407e-01;
    constexpr double p2 = 3.85197e-01; 
  
    //so that <feta()> = 1
    constexpr double norm = (p1+p2*(2.6+1.656)/2);
  
    double result = p0*(p1+p2*eta)/norm;
  
    return result;
  }
  
  double
  PFEnergyCalibration::Beta(double E, double eta) const {
  
   //Energy dependency
    constexpr double p0 = 0.032;
    constexpr double p1 = 9.70394e-02;
    constexpr double p2 = 2.23072e+01;
    constexpr double p3 = 100;
  
    //Eta dependency
    constexpr double p4 = 1.02496e+00 ;
    constexpr double p5 = -4.40176e-03 ;
  
    //so that <feta()> = 1
    constexpr double norm = (p4+p5*(2.6+1.656)/2);
  
    double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;			  
    return result;
  }
  
  
  double
  PFEnergyCalibration::Gamma(double etaEcal) const {
  
   //Energy dependency
    constexpr double p0 = 2.49752e-02;
  
    //Eta dependency
    constexpr double p1 = 6.48816e-02;
    constexpr double p2 = -1.59517e-02; 
   
    //so that <feta()> = 1
    constexpr double norm = (p1+p2*(2.6+1.656)/2);
  
    double result = p0*(p1+p2*etaEcal)/norm;					  
  
    return result;
  }
  
  
  
  ///////////////////////////////////////////////////////////////
  ////                                                       ////  
  ////   THE CORRECTIONS IN THE BARREL AND IN THE ENDCAPS    ////
  ////                                                       ////
  ///////////////////////////////////////////////////////////////
  
  
  // returns the corrected energy in the barrel (0,1.48)
  // Global Behaviour, phi and eta cracks are taken into account
  double
  PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
  				 bool crackCorrection ) const {
  
    // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
    double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
    double result = E * CorrBarrel(E,eta) * correction;
  
    return result;
  }
  
  
  // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
  double
  PFEnergyCalibration::EcorrZoneBeforePS(double E, double eta) const {
  
   //Energy dependency
    constexpr double p0 =1; 
    constexpr double p1 =0.18;
    constexpr double p2 =8.;
  
    //Eta dependency
    constexpr double p3 =0.3;
    constexpr double p4 =1.11;
    constexpr double p5 =0.025;
    constexpr double p6 =1.49;
    constexpr double p7 =0.6;
  
    //so that <feta()> = 1
    constexpr double norm = 1.21;
  
    double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
  
    return result;
  }
  
  
  // returns the corrected energy in the PS (1.65,2.6)
  // only when (ePS1>0)||(ePS2>0)
  double
  PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) const {
  
    // gives the good weights to each subdetector
    double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
  
    //Correction of the residual energy dependency
    constexpr double p0 = 1.00;
    constexpr double p1 = 2.18;
    constexpr double p2 =1.94;
    constexpr double p3 =4.13;
    constexpr double p4 =1.127;
  
    double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
  
    return result;
  } 
  
  // returns the corrected energy in the PS (1.65,2.6)
  // only when (ePS1>0)||(ePS2>0)
  double
  PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) const {
  
    // gives the good weights to each subdetector
    double gammaprime=Gamma(etaEcal)/9e-5;
  
    if(outputPS1 == 0 && outputPS2 == 0 && esEEInterCalib_ != nullptr){
      // both ES planes working
      // scaling factor accounting for data-mc                                                                                 
      outputPS1=gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
      outputPS2=gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
    }
    else if(outputPS1 == 0 && outputPS2 == -1 && esEEInterCalib_ != nullptr){
      // ESP1 only working
      double corrTotES = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0() * esEEInterCalib_->getGammaLow1();
      outputPS1 = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
      outputPS2 = corrTotES - outputPS1;
    }
    else if(outputPS1 == -1 && outputPS2 == 0 && esEEInterCalib_ != nullptr){
      // ESP2 only working
      double corrTotES = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3() * esEEInterCalib_->getGammaLow2();
      outputPS2 = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
      outputPS1 = corrTotES - outputPS2;
    }
    else{
      // none working
      outputPS1 = gammaprime*ePS1;
      outputPS2 = gammaprime*Alpha(etaEcal)*ePS2;
    }
  
    double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
  
    //Correction of the residual energy dependency
    constexpr double p0 = 1.00;
    constexpr double p1 = 2.18;
    constexpr double p2 =1.94;
    constexpr double p3 =4.13;
    constexpr double p4 =1.127;
    
    double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
    outputPS1*=corrfac;
    outputPS2*=corrfac;
    double result = E*corrfac;
  
    return result;
  } 
  
  
  // returns the corrected energy in the PS (1.65,2.6)
  // only when (ePS1=0)&&(ePS2=0)
  double 
  PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta) const {
  
    //Energy dependency
    constexpr double p0= 1.02;
    constexpr double p1= 0.165;
    constexpr double p2= 6.5 ;
    constexpr double p3=  2.1 ;
  
    //Eta dependency
    constexpr double p4 = 1.02496e+00 ;
    constexpr double p5 = -4.40176e-03 ;
  
    //so that <feta()> = 1
    constexpr double norm = (p4+p5*(2.6+1.656)/2);
  
    double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
  		  
    return result;
  }
  
  
  // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
  double
  PFEnergyCalibration::EcorrZoneAfterPS(double E, double eta) const {
  
    //Energy dependency
    constexpr double p0 =1; 
    constexpr double p1 = 0.058;
    constexpr double p2 =12.5;
    constexpr double p3 =-1.05444e+00;
    constexpr double p4 =-5.39557e+00;
    constexpr double p5 =8.38444e+00;
    constexpr double p6 = 6.10998e-01  ;
  
    //Eta dependency
    constexpr double p7 =1.06161e+00;
    constexpr double p8 = 0.41;
    constexpr double p9 =2.918;
    constexpr double p10 =0.0181;
    constexpr double p11= 2.05;
    constexpr double p12 =2.99;
    constexpr double p13=0.0287;
  
    //so that <feta()> = 1
    constexpr double norm=1.045;
  
    double result = E*(p0+p1*TMath::Exp(-(E-p3)/p2)+1/(p4+p5*TMath::Power(E,p6)))*(p7+p8*TMath::Gaus(eta,p9,p10)+p11*TMath::Gaus(eta,p12,p13))/norm;
    return result;
  }
  
  
  
  
  // returns the corrected energy everywhere
  // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
  double
  PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
  			   double eta,double phi,
  			   bool crackCorrection ) const {
  
    constexpr double endBarrel=1.48;
    constexpr double beginingPS=1.65;
    constexpr double endPS=2.6;
    constexpr double endEndCap=2.98;
   
    double result=0;
  
    eta=TMath::Abs(eta);
  
    if(eEcal>0){
      if(eta <= endBarrel)                         result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
      else if(eta <= beginingPS)                   result = EcorrZoneBeforePS(eEcal,eta);
      else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
      else if(eta < endPS)                         result = EcorrPS(eEcal,ePS1,ePS2,eta);
      else if(eta < endEndCap)                     result = EcorrZoneAfterPS(eEcal,eta); 
      else result =eEcal;
    }
    else result = eEcal;// useful if eEcal=0 or eta>2.98
    //protection
    if(result<eEcal) result=eEcal;
    return result;
  }
  
  // returns the corrected energy everywhere
  // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
  double
  PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection)  const {
  
    constexpr double endBarrel=1.48;
    constexpr double beginingPS=1.65;
    constexpr double endPS=2.6;
    constexpr double endEndCap=2.98;
   
    double result=0;
  
    eta=TMath::Abs(eta);
  
    if(eEcal>0){
      if(eta <= endBarrel)                         result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
      else if(eta <= beginingPS)                   result = EcorrZoneBeforePS(eEcal,eta);
      else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
      else if(eta < endPS)                         result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
      else if(eta < endEndCap)                     result = EcorrZoneAfterPS(eEcal,eta); 
      else result =eEcal;
    }
    else result = eEcal;// useful if eEcal=0 or eta>2.98
    // protection
    if(result<eEcal) result=eEcal;
    return result;
  }
  