#include #include #include #include #include #include #include "TRandom.h" #include "JTools/JQuantile.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to test JTOOLS::JQuantile calculation. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; int numberOfEvents; double x; double sigma; UInt_t seed; double precision; int debug; try { JParser<> zap("Example program to test quantile calculation."); zap['n'] = make_field(numberOfEvents); zap['x'] = make_field(x) = 0.0; zap['s'] = make_field(sigma) = 1.0; zap['S'] = make_field(seed) = 0; zap['e'] = make_field(precision) = 1.0e-2; zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (numberOfEvents < 2) { FATAL("Fatal error: number of events " << numberOfEvents << endl); } gRandom->SetSeed(seed); { vector buffer; JQuantile Q; for (int i = 0; i != numberOfEvents; ++i) { double x = gRandom->Uniform(-1.0e6, +1.0e6); buffer.push_back(x); Q.put(x); } random_device rd; mt19937 g(rd()); for (int i = 0; i != 100; ++i) { std::shuffle(buffer.begin(), buffer.end(), g); DEBUG("buffer[0] " << FIXED(15,6) << buffer[0] << endl); JQuantile R; for (const double x : buffer) { R.put(x); } ASSERT( Q.getCount() == R.getCount(), "Test shuffle counts."); ASSERT(fabs(Q.getMean() - R.getMean()) <= precision, "Test shuffle means."); ASSERT(fabs(Q.getSTDev() - R.getSTDev()) <= precision, "Test shuffle sigmas."); } } { JQuantile Q; for (int i = 0; i != numberOfEvents; ++i) { Q.put(gRandom->Gaus(x, sigma)); } NOTICE("quantity " << CENTER(10) << "calculated" << " | " << CENTER(10) << "true" << endl); NOTICE("mean " << FIXED(10,3) << Q.getMean() << " | " << FIXED(10,3) << x << endl); NOTICE("RMS " << FIXED(10,3) << Q.getSTDev() << " | " << FIXED(10,3) << sigma << endl); ASSERT(numberOfEvents > 0, "Test number of events."); ASSERT(fabs(Q.getMean() - x) <= precision, "Test mean."); ASSERT(fabs(Q.getSTDev() - sigma) <= precision, "Test sigma."); } { JQuantile Q("", true); for (int i = 0; i != numberOfEvents; ++i) { Q.put(gRandom->Uniform(0.0, 1.0)); } for (const double x : { 0.1, 0.5, 0.9}) { NOTICE("forward quantile " << FIXED(6,3) << x << ' ' << FIXED(6,3) << Q.getQuantile(x, JQuantile::forward_t) << endl); NOTICE("symmetric quantile " << FIXED(6,3) << x << ' ' << FIXED(6,3) << Q.getQuantile(x, JQuantile::symmetric_t) << endl); NOTICE("backward quantile " << FIXED(6,3) << x << ' ' << FIXED(6,3) << Q.getQuantile(x, JQuantile::backward_t) << endl); ASSERT(fabs(Q.getQuantile(x, JQuantile::forward_t) - ( x )) <= precision, "Test forward quantile "); ASSERT(fabs(Q.getQuantile(x, JQuantile::symmetric_t) - ( x )) <= precision, "Test symmetric quantile "); ASSERT(fabs(Q.getQuantile(x, JQuantile::backward_t) - (1.0 - x)) <= precision, "Test backward quantile "); } } return 0; }