#include #include #include #include #include #include "TRandom3.h" #include "JLang/JException.hh" #include "JMath/JPolynome.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JGrid.hh" #include "JTools/JQuantile.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to test 1D interpolation of a polynome. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; unsigned int numberOfEvents; unsigned int numberOfBins; JPolynome f1; double x0; int debug; try { JParser<> zap("Example program to test 1D interpolation of a polynome."); zap['n'] = make_field(numberOfEvents) = 0; zap['N'] = make_field(numberOfBins) = 21; zap['P'] = make_field(f1); zap['x'] = make_field(x0) = 0.0; zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (f1.empty()) { FATAL("Invalid polynomial."); } const int N = 3; // degree of polynomial interpolation typedef JResultPolynome JResult_t; typedef JPolintFunction1D, JCollection, JResult_t> JFunction1D_t; JFunction1D_t polint; const double xmin = -1.0; const double xmax = +1.0; polint.configure(JGrid(numberOfBins, xmin, xmax), f1); polint.compile(); polint.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero)); vector fp(1, f1); for (int i = 1; i != (int) f1.size(); ++i) { fp.push_back(fp.rbegin()->getDerivative()); } DEBUG("polynome: "); for (int i = 0; i != (int) f1.size(); ++i) { DEBUG(' ' << fp[i](x0)); } DEBUG(endl); DEBUG("result: "); JResult_t result = polint(x0); for (int i = 0; i != N+1; ++i) { DEBUG(' ' << result.y[i]); } DEBUG(endl); if (numberOfEvents != 0) { JQuantile Q[N+1]; for (int j = 0; j != sizeof(Q)/sizeof(Q[0]); ++j) { ostringstream os; os << "f^" << j; Q[j].setTitle(os.str()); } const int M = min(fp.size(), sizeof(Q)/sizeof(Q[0])); for (unsigned int i = 0; i != numberOfEvents; ++i) { const double x = gRandom->Uniform(xmin, xmax); JResult_t result = polint(x); for (int j = 0; j != M; ++j) { Q[j].put(fp[j](x) - result.y[j]); } } for (int j = 0; j != M; ++j) { Q[j].print(cout, j == 0); } } }