JRootfitToGauss.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <type_traits>
#include <chrono>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JMath/JMathlib.hh"
 
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootTestkit.hh"
#include "JROOT/JRootfit.hh"
#include "JROOT/JGraphErrors.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     = 25;
    double background =  5;
    double center     = 0.0;
    double sigma      = 1.0;
  } parameters;
 
  string     inputFile;
  string     outputFile;
  range_type X;
  size_t     N;
  string     option;
  bool       writeFits;
  ULong_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['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 nx   = 20;
  const double xmin = parameters.center  -  5.0 * parameters.sigma;
  const double xmax = parameters.center  +  5.0 * parameters.sigma;
 
  TH1D* h1 = NULL;
 
  if (inputFile == "") {
 
    h1 = new TH1D("h1", NULL, nx, xmin, xmax);
 
    h1->Sumw2();
 
    auto f0 = JP0<1>(parameters.signal) * JGauss<1>(parameters.center, parameters.sigma)  +  JP0<2>(parameters.background);
 
    if (N == 0)
      FillRandom(h1, f0);
    else
      FillRandom(h1, f0, N);
 
  } else {
 
    TFile* in = TFile::Open(inputFile.c_str(), "exist");
 
    in->GetObject("h1", h1);
  }
 
 
  // determine start values from data
 
  const double center     = h1->GetMean();
  const double sigma      = h1->GetStdDev()  * 0.66;
  const double signal     = h1->GetMaximum();
  const double background = h1->GetMinimum() + 0.1;
 
 
  // fit
 
  auto f1 = JP0<1>(signal) * JGauss<1>(center, sigma)  +  Exp(JP0<2>(log(background)));
 
  typedef decltype(f1)  function_type;
 
  JRootfit<function_type>::debug = debug;
 
  for (size_t i = 0; i != getNumberOfParameters<function_type>(); ++i) {
    cout << setw(2) << i << ' ' << FIXED(15,9) << f1[i] << endl;
  }
 
  {
    const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
 
    const auto result = (option == count_t ?
                         (writeFits ? Fit<m_count>(h1, f1, {}, X) : Fit<m_count>(*h1, f1, {}, X)) :
                         (writeFits ? Fit<m_value>(h1, f1, {}, X) : Fit<m_value>(*h1, f1, {}, X)));
 
    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<chrono::microseconds>(t1 - t0).count() << endl;
 
    for (size_t i = 0; i != result.getNumberOfParameters(); ++i) {
      cout << setw(2)     << i                  << ' ' 
           << FIXED(15,9) << result.getValue(i) << " +/- "
           << FIXED(15,9) << result.getError(i) << endl;
    }
  }
  {
    JGraphErrors_t gs;
 
    for (Int_t ix = 1; ix <= h1->GetXaxis()->GetNbins(); ++ix) {
 
      const double x     = h1->GetXaxis()->GetBinCenter(ix);
      const double value = h1->GetBinContent(ix);
      const double error = h1->GetBinError  (ix);
 
      gs.put(x, value, 0.0, error);
    }
 
    JGraphErrors g1(gs, "");
 
    const auto result = Fit(g1, f1, {}, X);
 
    cout << "chi2/NDF " << FIXED(7,3) << result.getChi2() << "/" << result.getNDF() << endl;
    cout << "Number of iterations " << result.numberOfIterations << endl;
 
    for (size_t i = 0; i != result.getNumberOfParameters(); ++i) {
      cout << setw(2)     << i                  << ' ' 
           << FIXED(15,9) << result.getValue(i) << " +/- "
           << FIXED(15,9) << result.getError(i) << endl;
    }
  }
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << *h1;
 
  out.Write();
  out.Close();
}