JRootfitToGauss3D.cc File Reference

Program to test JRootfit algorithm. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <type_traits>
#include <chrono>
#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Program to test JRootfit algorithm.

Author

mdejong

Definition in file JRootfitToGauss3D.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 30 of file JRootfitToGauss3D.cc.

{
  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<function_type>::debug = debug;

  for (size_t i = 0; i != getNumberOfParameters<function_type>(); ++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<m_count>(h3, f3, {}, X, Y, Z) : Fit<m_count>(*h3, f3, {}, X, Y, Z)) :
                         (writeFits ? Fit<m_value>(h3, f3, {}, X, Y, Z) : Fit<m_value>(*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<chrono::microseconds>(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();
}