// ReconCalibrator.cc
// Contact person: Logan Lebanowski <loganleb@hep.upenn.edu>
// See ReconCalibrator.hh for more details
//———————————————————————//

#include <RAT/DU/ReconCalibrator.hh>
#include <RAT/DB.hh>
#include <RAT/Log.hh>

namespace RAT {
  namespace DU {

    ReconCalibrator::ReconCalibrator() : fInitDone(false) {

    }

    void ReconCalibrator::Init()
    {

      if (fInitDone) return;
      warn << "\n\n\n" << BYELLOW
           << "This calibration should be applied on corrected energies reconstructed with RAT release 6.5.X. \nThe calibration is less accurate for other versions.\n\n\n"
           << CLR << newline;

      if (!DB::Get()->AreDefaultsLoaded()) {
        DB::Get()->LoadDefaults();
      }
      DBLinkPtr dbRCalibLink;
      try {
        dbRCalibLink = DB::Get()->GetLink("ReconCalibrator");
        fCoefficientsData = dbRCalibLink->GetDArray("CoefficientsData");
        fCoefficientsSim  = dbRCalibLink->GetDArray("CoefficientsSim");
      }
      catch(RAT::DBNotFoundError &e) {
        warn << "ReconCalibrator::Init : Failed to load calibration table. Message: " << e.what() << newline;
        fInitDone = false;
        return;
      }
      fInitDone = true;

    }


    double ReconCalibrator::CalibrateEnergyRSP( bool isData, double energy, double rho, double z ) const
    {
      if (!fInitDone) {
        warn << "ReconCalibrator::CalibrateEnergyRSP : Calibration being called without tables being loaded. Returning input, but this is likely wrong." << newline;
        return energy;
      }

      double scale, A, B, C, D, E;
      if( isData ) {
        A = fCoefficientsData[0];
        B = fCoefficientsData[1];
        C = fCoefficientsData[2];
        D = fCoefficientsData[3];
        E = fCoefficientsData[4];
      } else {
        A = fCoefficientsSim[0];
        B = fCoefficientsSim[1];
        C = fCoefficientsSim[2];
        D = fCoefficientsSim[3];
        E = fCoefficientsSim[4];
      }

      scale = A + ( (1. + B*rho*rho)*(1. + C*z + D*z*z + E*z*z*z) - 1. );
      return (1.+scale)*energy;

    }

  }
}