#include <RAT/DU/EnergySystematics.hh>
#include <TMath.h>
#include <CLHEP/Random/RandGauss.h>

ClassImp(RAT::DU::EnergySystematics)

namespace RAT {
namespace DU {
    EnergySystematics::EnergySystematics(double _energy_scale, double _resolution_uncertainty) {
        energy_scale_systematic = _energy_scale;
        resolution_uncertainty = _resolution_uncertainty;
        reconstruction_smearing_factor = sqrt(pow(1 + resolution_uncertainty, 2) -1);
    }

    double EnergySystematics::ApplyEnergyScale(const double& energy, const bool& negate) const {
        const int factor = negate ? -1 : 1;
        double new_energy = (1 + factor*energy_scale_systematic)*energy;
        return new_energy;
    }

    double EnergySystematics::ApplyEnergyResolution(const double& energy) const {
        // The smearing sigma scales with sqrt energy
        double sigma = sqrt(energy)*reconstruction_smearing_factor;
        // Return a new energy that's drawn from a guassian centered at the
        // energy, with the above sigma
        return  CLHEP::RandGauss::shoot(energy, sigma);
    }
}
}