JRootFitToGauss.cc File Reference

Program to test ROOT fit. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1D.h"
#include "JROOT/JMinimizer.hh"
#include "JTools/JStats.hh"
#include "JROOT/JRootToolkit.hh"
#include "JMath/JGauss.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

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

Detailed Description

Program to test ROOT fit.

Author

mdejong

Definition in file JRootFitToGauss.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 30 of file JRootFitToGauss.cc.

{
  using namespace std;
  using namespace JPP;
 
  string   outputFile;
  int      numberOfEvents;
  JGauss   gauss;
  JGauss   precision;
  int      debug;
 
  try {
 
    JParser<> zap("Program to test ROOT fit.");
 
    zap['o'] = make_field(outputFile)          = "";
    zap['n'] = make_field(numberOfEvents)      = 1000;
    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);
  }
 
 
  ASSERT(numberOfEvents > 0);
 
  TF1 fs("fs", "exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
  TF1 fb("fb", "1.0");
 
  fs.FixParameter(0, gauss.mean);
  fs.FixParameter(1, gauss.sigma);
 
  const Int_t    nx   = 21;
  const Double_t xmin = -5.0;
  const Double_t xmax = +5.0;
 
  JStats Q[] = {
    JStats("mean      "),
    JStats("sigma     "),
    JStats("signal    "),
    JStats("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.signal);
    h0.FillRandom("fb", (Int_t) gauss.background);
 
    TF1 f1("f1", "[2]*exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1])) / (TMath::Sqrt(2.0*TMath::Pi())*[1]) + [3]");
 
    f1.SetParameter(0, h0.GetMean());
    f1.SetParameter(1, h0.GetRMS());
    f1.SetParameter(2, h0.GetEntries() - h0.GetMinimum() * h0.GetNbinsX());
    f1.SetParameter(3, h0.GetMinimum());
 
    f1.SetParError(0, 1.0e-3);
    f1.SetParError(1, 1.0e-3);
    f1.SetParError(2, 0.5);
    f1.SetParError(3, 0.5);
 
    string option = "NWL";
 
    if (debug < debug_t && option.find('Q') == string::npos) {
      option += "Q";
    }
 
    timer.start();
 
    h0.Fit(&f1, option.c_str());
 
    timer.stop();
 
    const double Y[] = { f1.GetParameter(0)                       - gauss.mean,
                         f1.GetParameter(1)                       - gauss.sigma,
                         f1.GetParameter(2) * nx / (xmax - xmin)  - gauss.signal,
                         f1.GetParameter(3) * nx                  - gauss.background };
 
    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.mean);
  ASSERT(Q[1].getSTDev() < precision.sigma);
  ASSERT(Q[2].getSTDev() < precision.signal);
  ASSERT(Q[3].getSTDev() < precision.background);
 
  return 0;
}