JSimplexFitArray.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TRandom3.h"
#include "TMath.h"
 
#include "JFit/JSimplex.hh"
#include "JTools/JArray.hh"
#include "JTools/JArrayIterator.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JTOOLS::JArray;
  using TMath::PoissonI;
 
  static const int NUMBER_OF_DIMENSIONS  =  3;
 
 
  /**
   * Type definition of a point in ND-space.
   */
  typedef JArray<NUMBER_OF_DIMENSIONS, double>  point_type;
 
 
  /**
   * Model.
   */
  struct model_type :
    public JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>  
  {
    /**
     * Default constructor.
     */
    model_type() :
      JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  value            first value
     * \param  args             remaining values
     */
    template<class ...Args>
    model_type(argument_type value, const Args& ...args) :
      JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>(value, args...)
    {}
 
 
    /**
     * Model function.
     *
     * \param  point              point
     * \return                    value
     */
    inline double operator()(const point_type& point) const
    {
      double value = (*this)[0];
 
      for (int i = 0; i != NUMBER_OF_DIMENSIONS;++i) {
        value *= gauss(point[i] - (*this)[2*i + 1], (*this)[2*i + 2]);
      }
      
      return value;
    }
 
 
    /**
     * Gauss function (normalised to 1 at x = 0).
     *
     * \param  x                    x
     * \param  sigma                sigma
     * \return                      function value
     */
    static inline double gauss(const double x, const double sigma)
    {
      const double u = x / sigma;
 
      if (fabs(u) < 10.0)
        return exp(-0.5*u*u);
      else
        return 0.0;
    }
  };
 
 
  /**
   * Hit.
   */
  struct hit_type :
    public point_type
  {
    /**
     * Constructor.
     *
     * \param  point              point
     * \param  value              value
     * \param  sigma              sigma
     */
    hit_type(const point_type& point, const double value, const double sigma) :
      point_type(point),
      value(value),
      sigma(sigma)
    {}
 
    double value;
    double sigma;
  };
 
 
  /**
   * Fit function based on classical chi2.
   *
   * \param  model              model
   * \param  hit                hit
   * \return                    chi2
   */
  inline double chi2(const model_type& model, const hit_type& hit) 
  {
    const double u = (model(hit) - hit.value) / hit.sigma;
 
    return u*u;
  }
 
 
  /**
   * Fit function based on logarithm of likelihood.
   *
   * \param  model              model
   * \param  hit                hit
   * \return                    chi2
   */
  inline double likelihood(const model_type& model, const hit_type& hit) 
  {
    return -log(PoissonI(hit.value, model(hit)));
  }
}
 
 
/**
 * \file
 *
 * Program to test JFIT::JSimplex algorithm.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JArray<2, double> array_type;
 
  string      outputFile;
  int         numberOfEvents;
  array_type  parameters;
  ULong_t     seed;
  int         debug;
 
  try {
 
    JParser<> zap("Program to test JSimplex algorithm.");
 
    zap['o'] = make_field(outputFile)          = "simplex.root";
    zap['n'] = make_field(numberOfEvents);
    zap['M'] = make_field(parameters, "model parameters") = array_type(100.0, 1.0);
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  gRandom->SetSeed(seed);
 
 
  const double top   = parameters[0];
  const double sigma = parameters[1];
 
  model_type model;
  
  model[0] = top;
 
  for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
    model[i*2 + 1] = gRandom->Uniform(-sigma, +sigma);
    model[i*2 + 2] = gRandom->Gaus(sigma, 0.2*sigma);
  }
 
  const JGrid<double> grid(11, -3.0*sigma, +3.0*sigma);
 
 
  JSimplex<model_type> simplex;
 
  JSimplex<model_type>::debug               =  debug;
  JSimplex<model_type>::MAXIMUM_ITERATIONS  =  10000;
 
  double (*fit)(const model_type&, const hit_type&) = likelihood;  // fit function
 
 
  vector<TH1D*> H1;
 
  TH1D h0("chi2/NDF", NULL, 200, 0.0, 10.0);
 
  H1.push_back(new TH1D(MAKE_CSTRING("p" << 0), NULL, 101, -20.0, +20.0));
 
  for (size_t i =1; i != model.size(); ++i) {
    H1.push_back(new TH1D(MAKE_CSTRING("p" << i), NULL, 101, -0.1, +0.1));
  }
 
 
  for (int counter = 0; counter != numberOfEvents; ++counter) {
 
    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
 
 
    // generate data
    
    vector<hit_type> data;
 
    for (JArrayIterator<NUMBER_OF_DIMENSIONS, double> g1(grid); g1; ++g1) {
 
      const JArray<NUMBER_OF_DIMENSIONS, double> p1 = *g1;
 
      const double value = gRandom->Poisson(model(p1));
      const double sigma = (value > 1.0 ? sqrt(value) : 0.7);
 
      data.push_back(hit_type(p1, value, sigma));         
    }
 
 
    // start values
 
    JQuantile Q[NUMBER_OF_DIMENSIONS];
 
    double value = 0.0;
 
    for (const auto& hit : data) {
 
      for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
        Q[i].put(hit[i], hit.value);
      }
 
      if (hit.value > value) {
 
        simplex.value[0] = hit.value;
 
        for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
          simplex.value[2*i + 1] = hit[i];
        }
 
        value = hit.value;
      }
    }
 
    for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
      simplex.value[2*i + 2] = 1.0 * Q[i].getSTDev();
    }
 
 
    // steps
 
    simplex.step.resize(2*NUMBER_OF_DIMENSIONS + 1);
 
    model_type step;
 
    simplex.step[0]    = model_type();
    simplex.step[0][0] = 1.0;
 
    for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
 
      simplex.step[2*i + 1]          = model_type();
      simplex.step[2*i + 1][2*i + 1] = 0.1;
 
      simplex.step[2*i + 2]          = model_type();
      simplex.step[2*i + 2][2*i + 2] = 0.1;
    }
 
    DEBUG("start values " << simplex.value << endl);
 
    for (size_t i = 0; i != simplex.step.size(); ++i) {
      DEBUG("step: " << setw(2) << i << ' ' << simplex.step[i] << endl);
    }
 
 
    // fit
 
    const double chi2 = simplex(fit, data.begin(), data.end());
    const int    ndf  = data.size()  -  simplex.step.size();
 
 
    h0.Fill(chi2 / ndf);
 
    for (size_t i = 0; i != model.size(); ++i) {
      H1[i]->Fill(model[i] - simplex.value[i]);
    }
 
    DEBUG("chi2 / NDF" << FIXED(12,3) << chi2 << " / " << ndf << endl);
 
    DEBUG("final values " << simplex.value << endl);
    DEBUG("model values " << model         << endl);
 
    if (debug >= debug_t) {
 
      for (const auto& hit : data) {
 
        cout << "hit: ";
 
        for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
          cout << FIXED(7,3) << hit[i] << ' ';
        }
 
        const double value = simplex.value(hit);
 
        cout << FIXED(7,3) << hit.value << ' '
             << FIXED(7,3) << value  << ' '
             << FIXED(7,3) << (hit.value - value) / hit.sigma << endl;
      }
    }
 
  }
  STATUS(endl);
  
 
  if (outputFile != "") {
 
    TFile out(outputFile.c_str(), "recreate");
 
    out << h0;
 
    for (size_t i = 0; i != sizeof(H1)/sizeof(H1[0]); ++i) {
      out << *(H1[i]);
    }
 
    out.Write();
    out.Close();
  }
}