JGradientFitToGauss.cc File Reference

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

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TRandom3.h"
#include "JFit/JGradient.hh"
#include "JTools/JElement.hh"
#include "JTools/JQuantile.hh"
#include "JMath/JGauss.hh"
#include "JMath/JConstants.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::JGradient algorithm.

Author

mdejong

Definition in file JGradientFitToGauss.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 121 of file JGradientFitToGauss.cc.

{
  using namespace std;
  using namespace JPP;

  int      numberOfEvents;
  JGauss   gauss;
  JGauss   precision;
  UInt_t   seed;
  int      debug;

  try {

    JParser<> zap("Program to test JGradient algorithm.");

    zap['n'] = make_field(numberOfEvents)      = 0;
    zap['@'] = make_field(gauss)               = JGauss(0.0,  1.0,  1000.0, 10.0);
    zap['e'] = make_field(precision)           = JGauss(0.05, 0.05,   25.0,  25.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);

  if (numberOfEvents == 0) {

    buffer.push_back(element_type(-3.00000,     16.00000));
    buffer.push_back(element_type(-2.50000,     22.00000));
    buffer.push_back(element_type(-2.00000,     70.00000));
    buffer.push_back(element_type(-1.50000,    152.00000));
    buffer.push_back(element_type(-1.00000,    261.00000));
    buffer.push_back(element_type(-0.50000,    337.00000));
    buffer.push_back(element_type( 0.00000,    378.00000));
    buffer.push_back(element_type( 0.50000,    360.00000));
    buffer.push_back(element_type( 1.00000,    253.00000));
    buffer.push_back(element_type( 1.50000,    129.00000));
    buffer.push_back(element_type( 2.00000,     72.00000));
    buffer.push_back(element_type( 2.50000,     22.00000));
    buffer.push_back(element_type( 3.00000,     11.00000));

    JGradient gradient(1000000, 0, 1.0e-4, debug);

    // start values
    
    fit = JGauss(0.5, 0.5, 700.0, 0.0);

    gradient.push_back(JModifier_t("mean",       new JGaussEditor(fit, JGauss(1.0, 0.0, 0.0, 0.0)), 1.0e-2));
    gradient.push_back(JModifier_t("sigma",      new JGaussEditor(fit, JGauss(0.0, 1.0, 0.0, 0.0)), 1.0e-2));
    gradient.push_back(JModifier_t("signal",     new JGaussEditor(fit, JGauss(0.0, 0.0, 1.0, 0.0)), 5.0e-0));
    gradient.push_back(JModifier_t("background", new JGaussEditor(fit, JGauss(0.0, 0.0, 0.0, 1.0)), 5.0e-1));
    
    gradient(getChi2);

  } else {

    JQuantile Q[] = { JQuantile("mean      "),
                      JQuantile("sigma     "),
                      JQuantile("signal    "),
                      JQuantile("background") };

    const size_t nx   = 21;
    const double xmin = -5.0;
    const double xmax = +5.0;
  
    for (int i = 0; i != numberOfEvents; ++i) {
    
      STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);

      buffer.clear();
    
      for (double x = xmin, dx = (xmax - xmin) / (nx - 1); x < xmax + 0.5*dx; x += dx) {

        const double value = gauss(x);
      
        buffer.push_back(element_type(x, gRandom->Poisson(value)));
      }

      JGradient gradient(1000000, 0, 1.0e-4, debug);

      // start values
    
      fit = JGauss(0.5, 0.5, 700.0, 0.0);

      gradient.push_back(JModifier_t("mean",       new JGaussEditor(fit, JGauss(1.0, 0.0, 0.0, 0.0)), 1.0e-2));
      gradient.push_back(JModifier_t("sigma",      new JGaussEditor(fit, JGauss(0.0, 1.0, 0.0, 0.0)), 1.0e-2));
      gradient.push_back(JModifier_t("signal",     new JGaussEditor(fit, JGauss(0.0, 0.0, 1.0, 0.0)), 5.0e-0));
      gradient.push_back(JModifier_t("background", new JGaussEditor(fit, JGauss(0.0, 0.0, 0.0, 1.0)), 5.0e-1));

      const double chi2 = gradient(getChi2);

      DEBUG("Final value " << fit << endl);
      DEBUG("Chi2 " << chi2 << endl);

      Q[0].put(fit.mean       - gauss.mean);
      Q[1].put(fit.sigma      - gauss.sigma);
      Q[2].put(fit.signal     - gauss.signal);
      Q[3].put(fit.background - gauss.background);
    }

    for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
      NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
    }

    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;
}