#include #include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH3D.h" #include "JROOT/JRootToolkit.hh" #include "JROOT/JRootTestkit.hh" #include "JROOT/JRootfit.hh" #include "JMath/JMathlib3D.hh" #include "JLang/JManip.hh" #include "Jeep/JProperties.hh" #include "Jeep/JParser.hh" /** * \file * * Program to test JRootfit algorithm. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; const char* const count_t = "count"; const char* const value_t = "value"; struct { double signal = 100; double background = 5; double center = 0.0; double sigma = 1.0; } parameters; string inputFile; string outputFile; range_type X; range_type Y; range_type Z; size_t N; string option; bool writeFits; UInt_t seed; int debug; try { JProperties properties; properties.insert(gmake_property(parameters.signal)); properties.insert(gmake_property(parameters.background)); properties.insert(gmake_property(parameters.center)); properties.insert(gmake_property(parameters.sigma)); JParser<> zap("Program to test JRootfit algorithm."); zap['f'] = make_field(inputFile) = ""; zap['o'] = make_field(outputFile); zap['@'] = make_field(properties) = JPARSER::initialised(); zap['x'] = make_field(X) = range_type(); zap['y'] = make_field(Y) = range_type(); zap['z'] = make_field(Z) = range_type(); zap['n'] = make_field(N) = 0; zap['O'] = make_field(option) = count_t, value_t; zap['w'] = make_field(writeFits); zap['S'] = make_field(seed) = 0; zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } gRandom->SetSeed(seed); const size_t ls = 2; const size_t nx = 20 * ls; const double xmin = -5.0 * ls; const double xmax = +5.0 * ls; const size_t ny = 20 * ls; const double ymin = -5.0 * ls; const double ymax = +5.0 * ls; const size_t nz = 20 * ls; const double zmin = -5.0 * ls; const double zmax = +5.0 * ls; TH3D* h3 = NULL; if (inputFile == "") { h3 = new TH3D("h3", NULL, nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax); h3->Sumw2(); auto f3 = (JGauss3X<1>(parameters.center, parameters.sigma) * JGauss3Y<2>(parameters.center, parameters.sigma) * JGauss3Z<3>(parameters.center, parameters.sigma) * JP0<1>(parameters.signal) + JP0<2>(parameters.background)); if (N == 0) FillRandom(h3, f3); else FillRandom(h3, f3, N); } else { TFile* in = TFile::Open(inputFile.c_str(), "exist"); in->GetObject("h3", h3); } // determine start values from data const double x0 = h3->GetMean(1); const double y0 = h3->GetMean(2); const double z0 = h3->GetMean(3); const double xs = h3->GetStdDev(1) * 0.17; const double ys = h3->GetStdDev(2) * 0.17; const double zs = h3->GetStdDev(3) * 0.17; const double signal = h3->GetMaximum(); const double background = h3->GetMinimum() + 0.10; // fit auto f3 = (JGauss3X<1>(x0, xs) * JGauss3Y<2>(y0, ys) * JGauss3Z<3>(z0, zs) * JP0<1>(signal)) + JP0<2>(background); typedef decltype(f3) function_type; JRootfit::debug = debug; for (size_t i = 0; i != getNumberOfParameters(); ++i) { cout << setw(2) << i << ' ' << FIXED(12,5) << getParameter(f3, i) << endl; } { const chrono::steady_clock::time_point t0 = chrono::steady_clock::now(); const auto result = (option == count_t ? (writeFits ? Fit(h3, f3, {}, X, Y, Z) : Fit(*h3, f3, {}, X, Y, Z)) : (writeFits ? Fit(h3, f3, {}, X, Y, Z) : Fit(*h3, f3, {}, X, Y, Z))); const chrono::steady_clock::time_point t1 = chrono::steady_clock::now(); cout << "chi2/NDF " << FIXED(7,3) << result.getChi2() << "/" << result.getNDF() << endl; cout << "Number of iterations " << result.numberOfIterations << endl; cout << "Elapsed time [us] " << setw(8) << chrono::duration_cast(t1 - t0).count() << endl; for (size_t i = 0; i != result.getNumberOfParameters(); ++i) { cout << setw(2) << i << ' ' << FIXED(12,5) << result.getValue(i) << " +/- " << FIXED(12,5) << result.getError(i) << endl; } } TFile out(outputFile.c_str(), "recreate"); out << *h3; out.Write(); out.Close(); }