Newer
Older
#ifndef CH_INTERPOLATOR_H
#define CH_INTERPOLATOR_H
#include "ConfigParser.h"
#include "Math/Interpolator.h"
#include <string>
#include <vector>
#include <iostream>
ChInterpolator(ROOT::Math::Interpolation::Type type = ROOT::Math::Interpolation::kLINEAR)
: interp(0, type) {}
ChInterpolator(const char *path, double xf = 1.0, double yf = 1.0,
ROOT::Math::Interpolation::Type type = ROOT::Math::Interpolation::kLINEAR)
: interp(0, type)
{
LoadData(path, xf, yf);
}
bool LoadData(const char *path, double x_fac = 1.0, double y_fac = 1.0)
{
return LoadData(path,
[x_fac] (double xv) { return x_fac*xv; },
[y_fac] (double yv) { return y_fac*yv; });
}
template<class OperatorX, class OperatorY>
bool LoadData(const char *path, OperatorX &&xop)
{
return LoadData(path, xop, [] (double yv) { return yv; });
}
template<class OperatorX, class OperatorY>
bool LoadData(const char *path, OperatorX &&xop, OperatorY &&yop)
{
std::string pathstr = path;
if (pathstr.empty()) { return false; }
index.clear();
vals.clear();
ConfigParser parser;
if (!parser.ReadFile(pathstr)) {
std::cerr << "ChInterpolator Error: Cannot read file \"" << path << "\"" << std::endl;
interp.SetData(index, vals);
return false;
}
if (parser.NbofElements() < 2) { continue; }
parser >> x >> y;
// strictly increasing
if (index.empty() || index.back() != x) {
index.push_back(xop(x));
vals.push_back(yop(y));
} else if (index.size()) {
std::cout << "ChInterpolator Warning: Skip the same index "
<< x << ", " << y << std::endl;
}
}
interp.SetData(index, vals);
return true;
}
void Print()
{
for (size_t i = 0; i < index.size(); ++i) {
std::cout << std::setw(15) << index[i]
<< std::setw(15) << vals[i]
<< std::endl;
std::cout << index.size() << " entries." << std::endl;
}
inline double Eval(double x)
{
if (index.empty()) {
std::cout << "Trying to interpolate from an empty ChInterpolator!" << std::endl;
return 0.;
if (x < index.front()) { return vals.front(); }
if (x > index.back()) { return vals.back(); }
return interp.Eval(x);
}
double *GetX() { return &index[0]; }
double *GetY() { return &vals[0]; }
int GetEntries() { return index.size(); }
private:
ROOT::Math::Interpolator interp;
std::vector<double> index, vals;
};
#endif // CH_INTERPOLATOR_H