AkimaFunction.h

Go to the documentation of this file.

#ifndef AKIMAFUNCTION_H
#define AKIMAFUNCTION_H
 
#include "Math/Interpolator.h"
#include "TNamed.h"
 
#include <algorithm>
#include <iostream>
#include <memory>
#include <vector>
 
class AkimaFunction : public TNamed
{
  public:
    std::vector<double> x;
    std::vector<double> y;
 
    std::unique_ptr<ROOT::Math::Interpolator> interpolator;   //!
 
    AkimaFunction()
        : TNamed()
    {}
 
    AkimaFunction(const char* name, const std::vector<double>& xVals, const std::vector<double>& yVals)
        : TNamed(name, name)
        , x(xVals)
        , y(yVals)
    {
        SortPoints();
        BuildInterpolator();
    }
 
    virtual ~AkimaFunction() {}
 
    void SortPoints()
    {
        if (x.size() != y.size())
            return;
 
        std::vector<std::pair<double, double>> pairs;
        pairs.reserve(x.size());
 
        for (size_t i = 0; i < x.size(); ++i) {
            pairs.emplace_back(x[i], y[i]);
        }
 
        std::sort(pairs.begin(), pairs.end());
 
        for (size_t i = 0; i < pairs.size(); ++i) {
            x[i] = pairs[i].first;
            y[i] = pairs[i].second;
        }
    }
 
    void BuildInterpolator()
    {
        if (x.size() < 2)
            return;
 
        try {
            auto type = (x.size() < 5) ? ROOT::Math::Interpolation::kLINEAR : ROOT::Math::Interpolation::kAKIMA;
 
            interpolator = std::make_unique<ROOT::Math::Interpolator>(x, y, type);
        } catch (...) {
            std::cerr << "Interpolator build failed for " << GetName() << std::endl;
            interpolator.reset();
        }
    }
 
    void Init()
    {
        if (!interpolator) {
            BuildInterpolator();
        }
    }
 
    double Eval(double xx)
    {
        if (x.empty())
            return 0.0;
 
        if (!interpolator) {
            Init();
        }
 
        if (!interpolator) {
            return 0.0;
        }
 
        if (xx <= x.front())
            return y.front();
        if (xx >= x.back())
            return y.back();
 
        return interpolator->Eval(xx);
    }
 
    void Print(Option_t* option = "") const override
    {
        std::cout << "AkimaFunction: " << GetName() << " | N points = " << x.size() << std::endl;
    }
 
    ClassDefOverride(AkimaFunction, 1);
};
 
#endif   // AKIMAFUNCTION_H