interpolate.h

Go to the documentation of this file.

#ifndef NPSTAT_INTERPOLATE_HH_
#define NPSTAT_INTERPOLATE_HH_

#include "JetMETCorrections/InterpolationTables/interface/ProperDblFromCmpl.h"

namespace npstat {
    template<typename T>
    inline T interpolate_linear(const double x, const T& f0, const T& f1)
    {
        const typename ProperDblFromCmpl<T>::type dx = 1.0 - x;
        return f0*dx + f1*static_cast<typename ProperDblFromCmpl<T>::type>(x);
    }

    template<typename T>
    inline T interpolate_quadratic(const double x, const T& f0,
                                   const T& f1, const T& f2)
    {
        static const typename ProperDblFromCmpl<T>::type two = 2.0;
        const typename ProperDblFromCmpl<T>::type dx = x - 1.0;
        return f1 + ((f2 - f0)/two + ((f2 - f1) + (f0 - f1))*(dx/two))*dx;
    }

    template<typename T>
    inline T interpolate_cubic(const double x, const T& f0, const T& f1,
                               const T& f2, const T& f3)
    {
        return interpolate_linear(x*(3.0 - x)/2.0,
                                  interpolate_linear(x/3.0, f0, f3),
                                  interpolate_linear(x - 1.0, f1, f2));
    }


    template<typename T>
    unsigned interpolation_coefficients(T* buffer, unsigned bufLen,
                                        const T& f0, const T& f1);
    template<typename T>
    unsigned interpolation_coefficients(T* buffer, unsigned bufLen,
                                        const T& f0, const T& f1, const T& f2);
    template<typename T>
    unsigned interpolation_coefficients(T* buffer, unsigned bufLen,
                                        const T& f0, const T& f1,
                                        const T& f2, const T& f3);
}

#include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"

namespace npstat {
    template<typename T>
    unsigned interpolation_coefficients(T* buffer, const unsigned bufLen,
                                        const T& f0, const T& f1)
    {
        if (bufLen <= 1U) throw npstat::NpstatInvalidArgument(
            "In npstat::interpolation_coefficients: "
            "insufficient length of the output buffer");
        buffer[0] = f0;
        buffer[1] = f1 - f0;
        return 2U;
    }

    template<typename T>
    unsigned interpolation_coefficients(T* buffer, const unsigned bufLen,
                                        const T& f0, const T& f1, const T& f2)
    {
        if (bufLen <= 2U) throw npstat::NpstatInvalidArgument(
            "In npstat::interpolation_coefficients: "
            "insufficient length of the output buffer");
        buffer[0] = f0;
        buffer[1] = static_cast<T>((f1 - f2 + 3*(f1 - f0))/2.0);
        buffer[2] = static_cast<T>(((f0 - f1) + (f2 - f1))/2.0);
        return 3U;
    }

    template<typename T>
    unsigned interpolation_coefficients(T* buffer, const unsigned bufLen,
                                        const T& f0, const T& f1,
                                        const T& f2, const T& f3)
    {
        if (bufLen <= 3U) throw npstat::NpstatInvalidArgument(
            "In npstat::interpolation_coefficients: "
            "insufficient length of the output buffer");
        buffer[0] = f0;
        buffer[1] = static_cast<T>((11*(f1 - f0)+7*(f1 - f2)+2*(f3 - f2))/6.0);
        buffer[2] = static_cast<T>((2*(f0 - f1)+3*(f2 - f1)+(f2 - f3))/2.0);
        buffer[3] = static_cast<T>(((f3 - f0) + 3*(f1 - f2))/6.0);
        return 4U;
    }
}


#endif // NPSTAT_INTERPOLATE_HH_