JGandalfFitAlaMinuitToGauss.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
 
#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<double, double>     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<JGauss>::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<JGauss>::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<JElement_t> data;
 
    for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
      data.push_back(JElement_t(h0.GetBinCenter (i),
                                h0.GetBinContent(i)));
    }
 
    JGandalf<JGauss> 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;
}