JSimplexFitArray.cc File Reference

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

#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"

Go to the source code of this file.

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

Detailed Description

Program to test JFIT::JSimplex algorithm.

Author

mdejong

Definition in file JSimplexFitArray.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 160 of file JSimplexFitArray.cc.

{
  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();
  }
}