JSchwarzenegger.cc File Reference

Program to test JFIT::JGandalf algorithm. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "TRandom3.h"
#include "JFit/JGandalf.hh"
#include "JTools/JElement.hh"
#include "JTools/JStats.hh"
#include "JROOT/JRootToolkit.hh"
#include "JMath/JGauss.hh"
#include "JMath/JConstants.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 JFIT::JGandalf algorithm.

Author

mdejong

Definition in file JSchwarzenegger.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 74 of file JSchwarzenegger.cc.

{
  using namespace std;
  using namespace JPP;
 
  string   outputFile;
  JGauss   gauss;
  Int_t    seed;
  int      debug;
 
  try {
 
    JParser<> zap("Program to test JGandalf algorithm.");
 
    zap['o'] = make_field(outputFile)          = "";
    zap['@'] = make_field(gauss)               = JGauss(0.0,  1.0,  1000.0, 100.0);
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 3;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  gRandom->SetSeed(seed);
  
  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.mean);
  fs.FixParameter(1, gauss.sigma);
 
  const Int_t    nx   = 21;
  const Double_t xmin = -5.25;
  const Double_t xmax = +5.25;
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D h0("h0", NULL, nx,     xmin, xmax);
  TH1D h1("h1", NULL, nx*100, xmin, xmax);
 
  h0.Sumw2();
  h1.Sumw2();
  /*
  for (size_t i = 0; i != (size_t) (gauss.signal);          ++i) { h0.Fill(fs.GetRandom(gRandom)); }
  for (size_t i = 0; i != (size_t) (gauss.background * nx); ++i) { h0.Fill(fb.GetRandom(gRandom)); }
  */
  h0.FillRandom("fs", (Int_t) gauss.signal);
  h0.FillRandom("fb", (Int_t) gauss.background * nx);
  
  const Int_t    ny   = 5000;
  const Double_t ymin =  0.0;
  const Double_t ymax = gauss.signal + 3.0 * sqrt(gauss.signal) * (xmax - xmin) / nx + gauss.background;
 
  TH2D h2("h2", NULL,
          nx, xmin, xmax,
          ny, ymin, ymax);
 
 
  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;
 
  JGandalf<JGauss>::EPSILON = 1.0e-4;
 
  fit.parameters.resize(4);
 
  fit.parameters[0] = &JGauss::mean;
  fit.parameters[1] = &JGauss::sigma;
  fit.parameters[2] = &JGauss::signal;
  fit.parameters[3] = &JGauss::background;
 
  fit.value = JGauss(h0.GetMean(),
                     h0.GetRMS(),
                     h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
                     h0.GetMinimum());
    
  const double chi2 = fit(g1, data.begin(), data.end());
 
  DEBUG("Final value " << fit.value << endl);
  DEBUG("Chi2 " << chi2 << endl);
 
  for (Int_t ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {
 
    const double x = h1.GetXaxis()->GetBinCenter(ix);
    const double y = fit.value(x);
 
    h1.SetBinContent(ix, y);
  }
    
  for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
    for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
 
      const double x = h2.GetXaxis()->GetBinCenter(ix);
      const double y = h2.GetYaxis()->GetBinCenter(iy);
      const double z = TMath::Poisson(y, fit.value(x));
 
      h2.SetBinContent(ix, iy, z);
    }
  }
 
  out.Write();
  out.Close();
}