#include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TMath.h" #include "TString.h" #include "JTools/JFunction1D_t.hh" #include "JTools/JQuadrature.hh" #include "JTools/JGrid.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" // command line options int numberOfPoints; bool quadrature; int debug; namespace { /** * Test function f1(). * * \param x abscissa value * \return ordinate value */ inline Double_t f1(const Double_t x) { return TMath::Gaus(x, 0.0, 1.0, kTRUE); } /** * Derivative of f1(). * * \param x abscissa value * \return derivative value */ inline Double_t g1(const Double_t x) { return f1(x) * -x; } /** * Integral of f1(). * * \param x abscissa value * \return integral value */ inline Double_t G1(const Double_t x) { return 0.5 * (1.0 + TMath::Erf(sqrt(0.5)*x)); } /** * Method to test interpolations for PDF. */ template void test(const char* title) { using namespace JPP; JFunction1D_t buffer; const int nx = 1000; const double xmin = -3.5; const double xmax = +3.5; if (!quadrature) buffer.configure(make_grid(numberOfPoints, xmin, xmax), f1); else buffer.configure(JGaussHermite(numberOfPoints), f1); buffer.compile(); TH1D f0(TString("f0[") + title + "]", NULL, nx, xmin, xmax); TH1D g0(TString("g0[") + title + "]", NULL, nx, xmin, xmax); TH1D G0(TString("G0[") + title + "]", NULL, nx, xmin, xmax); for (int i = 1; i <= f0.GetNbinsX(); ++i) { const Double_t x = f0.GetBinCenter(i); try { const typename JFunction1D_t::result_type result = buffer(x); f0.SetBinContent(i, result.f); g0.SetBinContent(i, result.fp); G0.SetBinContent(i, result.v); } catch(const std::exception& error) { //ERROR(error.what() << endl); } } f0.Write(); g0.Write(); G0.Write(); } } /** * \file * * Example program to test interpolations for PDF. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; string outputFile; try { JParser<> zap("Example program to test interpolations for PDF."); zap['o'] = make_field(outputFile); zap['N'] = make_field(numberOfPoints) = 15; zap['Q'] = make_field(quadrature); zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } TFile out(outputFile.c_str(), "recreate"); const int nx = 1000; const double xmin = -3.5; const double xmax = +3.5; TH1D f0(TString("f0[") + "true" + "]", NULL, nx, xmin, xmax); TH1D g0(TString("g0[") + "true" + "]", NULL, nx, xmin, xmax); TH1D G0(TString("G0[") + "true" + "]", NULL, nx, xmin, xmax); for (int i = 1; i <= f0.GetNbinsX(); ++i) { const Double_t x = f0.GetBinCenter(i); f0.SetBinContent(i, f1(x)); g0.SetBinContent(i, g1(x)); G0.SetBinContent(i, G1(x)); } { const int N = 4; typedef JPolintFunction1D, JCollection, JResultPDF > JFunction1D_t; test("Polint"); } { typedef JSplineFunction1S_t JFunction1D_t; test("Spline"); } out.Write(); out.Close(); }