#include #include #include #include "TMath.h" #include "TRandom.h" #include "JTools/JHistogram1D_t.hh" #include "JTools/JHistogramMap_t.hh" #include "JTools/JMultiHistogram.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JMultiPDF.hh" #include "JTools/JToolsToolkit.hh" #include "JTools/JQuadrature.hh" #include "JTools/JQuantile.hh" #include "JIO/JFileStreamIO.hh" #include "JLang/JObjectIO.hh" #include "Jeep/JPrint.hh" #include "Jeep/JTimer.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * Normalised function in 2D. * * \param x x value * \param y y value * \return function value */ inline Double_t g2(const Double_t x, const Double_t y) { return (TMath::Gaus(x, 0.0, 1.0, kTRUE) * TMath::Gaus(y, 0.0, 1.0, kTRUE)); } /** * Function in 4D. * * \param x0 x0 value * \param x1 x1 value * \param x2 x2 value * \param x3 x3 value * \return function value */ inline Double_t g4(const Double_t x0, const Double_t x1, const Double_t x2, const Double_t x3) { return g2(x2 - x0, x3 - x1); } } /** * \file * * Example program to test conditional probability density function \f$P(x_0,x_1\|x_2,x_3)\f$ using JTOOLS::JMultiPDF and JTOOLS::JMultiFunction. * \author mdejong */ int main(int argc, char **argv) { using namespace std; string inputFile; string outputFile; int numberOfEvents; int numberOfBins; int debug; try { JParser<> zap("Example program to test conditional probability density function."); zap['f'] = make_field(inputFile) = ""; zap['o'] = make_field(outputFile) = ""; zap['n'] = make_field(numberOfEvents); zap['N'] = make_field(numberOfBins) = 15; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (numberOfEvents <= 0) { FATAL("No events." << endl); } using namespace JPP; const double xmin = -10.0; const double xmax = +10.0; const double dx = (xmax - xmin) / (numberOfBins - 1); if (outputFile != "") { typedef JMultiHistogram::maplist> JMultiHistogram_t; JMultiHistogram_t histogram; const JGaussHermite bounds(numberOfBins); for (double x0 = xmin; x0 <= xmax + 0.5*dx; x0 += dx) { for (double x1 = xmin; x1 <= xmax + 0.5*dx; x1 += dx) { for (JGaussHermite::const_iterator x2 = bounds.begin(); x2 != bounds.end(); ++x2) { for (JGaussHermite::const_iterator x3 = bounds.begin(); x3 != bounds.end(); ++x3) { histogram[x0][x1][x0 + x2->getX()][x1 + x3->getX()] = 0.0; } } } } // fill JTimer timer("histogram"); for (int i = 0; i != numberOfEvents; ++i) { if (i%1000 == 0) { STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl); } const double x0 = gRandom->Uniform(xmin, xmax); const double x1 = gRandom->Uniform(xmin, xmax); const double x2 = x0 + gRandom->Gaus(0.0, 1.0); const double x3 = x1 + gRandom->Gaus(0.0, 1.0); const double w = 1.0; timer.start(); histogram.fill(x0, x1, x2, x3, w); timer.stop(); } STATUS(endl); if (debug >= status_t) { timer.print(cout, true, micro_t); } try { NOTICE("storing output to file " << outputFile << "... " << flush); JLANG::store(outputFile.c_str(), histogram); NOTICE("OK" << endl); } catch(const JException& error) { FATAL(error.what() << endl); } } if (inputFile != "") { typedef JMultiHistogram::maplist> JHistogram2D_t; typedef JMultiHistogram::maplist> JMultiHistogram_t; JMultiHistogram_t histogram; try { NOTICE("loading input to file " << inputFile << "... " << flush); JLANG::load(inputFile.c_str(), histogram); NOTICE("OK" << endl); } catch(const JException& error) { FATAL(error.what() << endl); } typedef JMultiPDF::maplist> JFunction2D_t; // PDF typedef JMultiFunction::maplist> JMultiFunction_t; // interpolation JMultiFunction_t pdf(histogram); // test JQuantile Q; for (int i = 0; i != numberOfEvents; ++i) { const double x0 = gRandom->Uniform(xmin, xmax); const double x1 = gRandom->Uniform(xmin, xmax); const double x2 = x0 + gRandom->Gaus(0.0, 1.0); const double x3 = x1 + gRandom->Gaus(0.0, 1.0); try { const double u = g4 (x0, x1, x2, x3); const double v = pdf(x0, x1, x2, x3); Q.put(u - v); } catch(const std::exception& error) {} } JQuantile V; for (JMultiFunction_t::super_const_iterator i = pdf.super_begin(); i != pdf.super_end(); ++i) { V.put(getIntegral((*i).getValue().getMultiFunction())); } cout << "normalisation " << FIXED(6,4) << V.getMean() << " +/- " << FIXED(6,4) << V.getSTDev() << endl; cout << "efficiency " << FIXED(5,3) << (double) Q.getCount() / (double) numberOfEvents << endl; cout << "mean " << SCIENTIFIC(10,2) << Q.getMean() << endl; cout << "RMS " << SCIENTIFIC(10,2) << Q.getSTDev() << endl; } }