#include #include #include #include "TMath.h" #include "TRandom.h" #include "JTools/JFunction1D_t.hh" #include "JTools/JTransformableMultiFunction.hh" #include "JTools/JFunctionalMap_t.hh" #include "JTools/JQuantile.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using namespace JPP; /** * Test function. * * \param x x value * \param y y value * \return function value */ inline Double_t f2(const Double_t x, const Double_t y) { return TMath::Gaus(y, 0.0, x, kTRUE); } /** * Derivative of test function. * * \param x x value * \param y y value * \return derivative value */ inline Double_t g2(const Double_t x, const Double_t y) { return TMath::Gaus(y, 0.0, x, kTRUE) * -y/x; } /** * Integral of test function. * * \param x x value * \param y y value * \return integral value */ inline Double_t G2(const Double_t x, const Double_t y) { return 0.5 * (1.0 + TMath::Erf(sqrt(0.5)*y/x)); } typedef JTransformableMultiFunction::maplist> function_type; typedef typename function_type::transformer_type transformer_type; typedef typename function_type::result_type result_type; /** * Transformer. */ struct JTransformer_t : public transformer_type { typedef typename transformer_type::clone_type clone_type; typedef typename transformer_type::argument_type argument_type; typedef typename transformer_type::const_array_type const_array_type; /** * Evaluate y as a function of x. * * \param buffer x value * \param xn y value * \return y value */ virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const { return xn * buffer[0]; } /** * Evaluate y as a function of x. * * \param buffer x value * \param xn y value * \return y value */ virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const { return xn / buffer[0]; } /** * Weight function. * * \param buffer x0 - xn-1 values * \return weight */ virtual double getWeight(const_array_type& buffer) const { return 1.0; } /** * Read from input. * * \param in reader * \return reader */ virtual JReader& read(JReader& in) { return in; } /** * Write to output. * * \param out writer * \return writer */ virtual JWriter& write(JWriter& out) const { return out; } /** * Get clone of this object. * * \return pointer to newly created object */ virtual clone_type clone() const { return new JTransformer_t(*this); } }; } /** * \file * * Example program to test transformable function. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; int numberOfEvents; int numberOfBins; bool transform; int debug; try { JParser<> zap("Example program to test transformable function."); zap['n'] = make_field(numberOfEvents); zap['N'] = make_field(numberOfBins) = 35; zap['U'] = make_field(transform); zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (numberOfEvents <= 0) { FATAL("No events." << endl); } const double xmin = 0.1; const double xmax = 1.0; const double dx = (xmax - xmin) / 3; function_type h2; for (double x = xmin; x < xmax + 0.5*dx; x += dx) { const double ymin = -5.0 * xmax; const double ymax = +5.0 * xmax; const double dy = (ymax - ymin) / (numberOfBins - 1); for (double y = ymin; y < ymax + 0.5*dy; y += dy) { h2[x][y] = f2(x,y); } } h2.setExceptionHandler(new typename function_type::JDefaultResult(JMATH::zero)); if (transform) { h2.transform(JTransformer_t()); } h2.compile(); JQuantile Q[] = { JQuantile("function"), JQuantile("derivative"), JQuantile("integral") }; const double ymin = -5.0 * xmax; const double ymax = +5.0 * xmax; for (int i = 0; i != numberOfEvents; ++i) { const double x = gRandom->Uniform(xmin, xmax); const double y = gRandom->Uniform(ymin, ymax); try { const result_type value = h2(x,y); const double f = f2(x,y); const double fp = g2(x,y); const double v = G2(x,y); Q[0].put(value.f - f); Q[1].put(value.fp - fp); Q[2].put(value.v - v); } catch(const std::exception& error) {} } if (debug >= notice_t) { cout << " "; for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { cout << " " << setw(10) << Q[i].getTitle(); } cout << endl; cout << "mean "; for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { cout << " " << SCIENTIFIC(10,2) << Q[i].getMean(); } cout << endl; cout << "RMS "; for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { cout << " " << SCIENTIFIC(10,2) << Q[i].getSTDev(); } cout << endl; } }