JPolfit1D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
 
#include "TRandom3.h"
 
#include "JLang/JException.hh"
#include "JMath/JPolynome.hh"
#include "JTools/JPolfit.hh"
#include "JTools/JPolint.hh"
#include "JTools/JElement.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JToolsToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to test 1D Legendre interpolation of a polynome.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  unsigned int  numberOfEvents;
  unsigned int  numberOfBins;
  JPolynome     f1;
  double        sigma;
  double        P;
  int           debug;
 
  try {
 
    JParser<> zap("Example program to test 1D Legendre interpolation of a polynome.");
 
    zap['n'] = make_field(numberOfEvents)   =  0;
    zap['N'] = make_field(numberOfBins)     = 21;
    zap['P'] = make_field(f1);
    zap['s'] = make_field(sigma)            = 0.0;
    zap['p'] = make_field(P)                = 0.0;
    zap['d'] = make_field(debug)            =  3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const int N = 7;                        // number of points for Legendre interpolation
  const int M = 3;                        // degree of Legendre polynome
 
  typedef JPolfitFunction1D<N, 
                            M, 
                            JElement2D<double, double>, 
                            JGridCollection, 
                            double>       JFunction1D_t;
 
  JFunction1D_t g1;
 
  JPolintFunction1D<M, 
                    JElement2D<double, double>, 
                    JGridCollection, 
                    double> h1;
 
  const double xmin = -10.0;
  const double xmax = +10.0;  
 
  const JGrid<double> grid(numberOfBins, xmin, xmax);
 
  for (int i = 0; i != grid.getSize(); ++i) {
 
    const double x = grid.getX(i);
 
    g1[x]= f1(x) + (sigma > 0.0 && gRandom->Rndm() < P ? gRandom->Gaus(0.0, sigma) : 0.0);
  }
 
  copy(g1, h1);
 
  g1.compile();
  h1.compile();
 
  for (JFunction1D_t::const_iterator i = g1.begin(); i != g1.end(); ++i) {
    DEBUG(FIXED(12,5) << i->getX() << ' ' << FIXED(12,5) << i->getY() << endl);
  }
  
    
  if (numberOfEvents != 0) {
 
    JQuantile Q[2];
    
    for (unsigned int i = 0; i != numberOfEvents; ++i) {
      
      const double x = gRandom->Uniform(xmin, xmax);
      const double y = f1(x);
      const double u = g1(x);
      const double v = h1(x);
 
      DEBUG(FIXED(12,5) << x << ' ' << FIXED(12,5) << y << ' ' << FIXED(12,5) << u << ' ' << FIXED(12,5) << v << endl);
 
      Q[0].put(u - y);
      Q[1].put(v - y);
    }
 
    Q[0].print(cout);
    Q[1].print(cout, false);
  }
 
  return 0;
}