#include "ausa/eloss/RangeInterpolator.h"
#include "ausa/eloss/RangeTabulation.h"
#include "ausa/util/memory"

#include <Math/Interpolator.h>

using namespace std;
using namespace AUSA::EnergyLoss;
using namespace ROOT::Math;

struct RangeInterpolator::Pimpl {
    Pimpl(const RangeTabulation& tabulation) {
        // get range table
        const auto& R = tabulation.getRange();
        // create interpolation objects, feed data to them, and assign pointer
        p = unique_ptr<Interpolator>( new Interpolator( R.first, R.second, Interpolation::kCSPLINE) );     //::kAKIMA) );
        pinv = unique_ptr<Interpolator>( new Interpolator( R.second, R.first, Interpolation::kCSPLINE) );  //::kAKIMA) );
        // limits of interpolation
        low.first = R.first.front();
        high.first = R.first.back();
        low.second = R.second.front();
        high.second = R.second.back();
    }

    double range(double energy) {
        double result;
        if (energy<=low.first) result = low.second;
        else if (energy>=high.first) result = high.second;
        else result = p -> Eval(energy);
        return result;
    }

    double energy(double range) {
        double result;
        if (range<=low.second) result = low.first;
        else if (range>=high.second) result = high.first;
        else result = pinv -> Eval(range);
        return result;
    }

    std::unique_ptr<ROOT::Math::Interpolator> p, pinv;
    using limit = std::pair<double,double>;
    limit low, high;
};

RangeInterpolator::RangeInterpolator( const RangeTabulation& tabulation )
    : pimpl(std::make_unique<Pimpl>(tabulation))
{
}

RangeInterpolator::~RangeInterpolator() {
    // Nothing
}

double RangeInterpolator::getRange( double energy ) const {
    return pimpl -> range(energy);
}   

double RangeInterpolator::getEnergy( double range ) const {
    return pimpl -> energy(range);
}

double RangeInterpolator::getMinimumEnergy() const {
    return pimpl -> low.first;
}