#include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TF1.h" #include "TH1D.h" #include "TMath.h" #include "JMath/JModel.hh" #include "JMath/JConstants.hh" #include "JFit/JGandalf.hh" #include "JTools/JElement.hh" #include "JTools/JQuantile.hh" #include "JROOT/JRootToolkit.hh" #include "Jeep/JTimer.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using namespace JPP; using TMath::PoissonI; typedef JElement2D JElement_t; /** * Gauss. */ struct JGauss : public JMATH::JModel_t { /** * Default constructor. */ JGauss() {} /** * Constructor. * * \param mean mean * \param sigma sigma * \param signal signal * \param background background */ JGauss(const double mean, const double sigma, const double signal, const double background) { this->push_back(mean); this->push_back(sigma); this->push_back(signal); this->push_back(background); } }; typedef JGandalf::result_type result_type; /** * Fit function. * * \param gauss gauss * \param point point * \return chi2 */ inline result_type g1(const JGauss& gauss, const JElement_t& point) { result_type result; const double u = (point.getX() - gauss[0]) / gauss[1]; const double fs = gauss[2] * exp(-0.5*u*u) / (sqrt(2.0*PI) * gauss[1]); const double fb = gauss[3]; const double f1 = fs + fb; const double p = PoissonI(point.getY(), f1); result.chi2 = -log(p); result.gradient[0] = fs * u / gauss[1]; // d(f)/d(mean) result.gradient[1] = fs * u*u / gauss[1] - fs / gauss[1]; // d(f)/d(sigma) result.gradient[2] = fs / gauss[2]; // d(f)/d(signal) result.gradient[3] = 1.0; // d(f)/d(background) result.gradient *= 1.0 - point.getY()/f1; // d(chi2)/d(f) return result; } } /** * \file * * Program to test JFIT::JGandalf algorithm. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; string outputFile; int numberOfEvents; JGauss gauss; JGauss precision; int debug; try { JParser<> zap("Program to test JGandalf algorithm."); zap['o'] = make_field(outputFile) = ""; zap['n'] = make_field(numberOfEvents) = 100; zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 100.0); zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0); zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } using namespace JFIT; ASSERT(numberOfEvents > 0); JGandalf::debug = debug; TF1 fs("fs", "exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))"); TF1 fb("fb", "1.0"); fs.FixParameter(0, gauss[0]); fs.FixParameter(1, gauss[1]); const Int_t nx = 21; const Double_t xmin = -5.0; const Double_t xmax = +5.0; JQuantile Q[] = { JQuantile("mean "), JQuantile("sigma "), JQuantile("signal "), JQuantile("background") }; TH1D H[] = { TH1D("ha", "", 101, -0.1, +0.1), TH1D("hb", "", 101, -0.1, +0.1), TH1D("hc", "", 101, -100.0, +100.0), TH1D("hd", "", 101, -100.0, +100.0) }; JTimer timer; for (int i = 0; i != numberOfEvents; ++i) { STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl); TH1D h0("h0", NULL, nx, xmin, xmax); h0.Sumw2(); h0.FillRandom("fs", (Int_t) gauss[2]); h0.FillRandom("fb", (Int_t) gauss[3]); vector data; for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) { data.push_back(JElement_t(h0.GetBinCenter (i), h0.GetBinContent(i))); } JGandalf fit; fit.parameters.resize(4); for (size_t i = 0; i != 4; ++i) { fit.parameters[i] = i; } fit.value = JGauss(h0.GetMean(), h0.GetRMS(), h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(), h0.GetMinimum()); DEBUG("Start value " << fit.value << endl); timer.start(); const double chi2 = fit(g1, data.begin(), data.end()); timer.stop(); DEBUG("Final value " << fit.value << endl); DEBUG("Chi2 " << chi2 << endl); const double Y[] = { fit.value[0] - gauss[0], fit.value[1] - gauss[1], fit.value[2] * nx / (xmax - xmin) - gauss[2], fit.value[3] * nx - gauss[3] }; for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { Q[i].put (Y[i]); H[i].Fill(Y[i]); } } for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { NOTICE((i == 0 ? longprint : shortprint) << Q[i]); } if (debug >= notice_t) { timer.print(cout, true, micro_t); } if (outputFile != "") { TFile out(outputFile.c_str(), "recreate"); for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) { out << H[i]; } out.Write(); out.Close(); } for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev()); } ASSERT(Q[0].getSTDev() < precision[0]); ASSERT(Q[1].getSTDev() < precision[1]); ASSERT(Q[2].getSTDev() < precision[2]); ASSERT(Q[3].getSTDev() < precision[3]); return 0; }