JPolynome1P.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
 
#include "TRandom3.h"
 
#include "JLang/JException.hh"
#include "JMath/JPolynome.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantile.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to test 1D 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        x0;
  int           debug;
 
  try {
 
    JParser<> zap("Example program to test 1D interpolation of a polynome.");
 
    zap['n'] = make_field(numberOfEvents)   =  0;
    zap['N'] = make_field(numberOfBins)     = 21;
    zap['P'] = make_field(f1);
    zap['x'] = make_field(x0)               =  0.0;
    zap['d'] = make_field(debug)            =  1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  if (f1.empty()) {
    FATAL("Invalid polynomial.");
  }
 
 
  const int N = 3;                      // degree of polynomial interpolation
  
  typedef JResultPolynome<N, double>                                                JResult_t;
  typedef JPolintFunction1D<N, JElement2D<double, double>, JCollection, JResult_t>  JFunction1D_t;
 
 
  JFunction1D_t  polint;
 
 
  const double xmin = -1.0;
  const double xmax = +1.0;  
 
  polint.configure(JGrid<double>(numberOfBins, xmin, xmax), f1);
  
  polint.compile();
  
  polint.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
    
 
  vector<JPolynome> fp(1, f1);
 
  for (int i = 1; i != (int) f1.size(); ++i) {
    fp.push_back(fp.rbegin()->getDerivative());
  }
 
  DEBUG("polynome: ");
 
  for (int i = 0; i != (int) f1.size(); ++i) {
    DEBUG(' ' << fp[i](x0));
  }
  DEBUG(endl);
 
  DEBUG("result:   ");
 
  JResult_t result = polint(x0);
 
  for (int i = 0; i != N+1; ++i) {
    DEBUG(' ' << result.y[i]);
  }
  DEBUG(endl);
 
 
  if (numberOfEvents != 0) {
 
    JQuantile Q[N+1];
 
    for (int j = 0; j != sizeof(Q)/sizeof(Q[0]); ++j) {
 
      ostringstream os;
 
      os << "f^" << j;
 
      Q[j].setTitle(os.str());
    }
 
    const int M = min(fp.size(), sizeof(Q)/sizeof(Q[0]));
 
    for (unsigned int i = 0; i != numberOfEvents; ++i) {
 
      const double x   = gRandom->Uniform(xmin, xmax);
      
      JResult_t result = polint(x);
 
      for (int j = 0; j != M; ++j) {
        Q[j].put(fp[j](x) - result.y[j]);
      }
    }
 
    for (int j = 0; j != M; ++j) {
      Q[j].print(cout, j == 0);
    }
  }
}