JPolynome2Dx1D.cc

Go to the documentation of this file.

 
#include <iostream>
#include <iomanip>
#include <cmath>
 
#include "TRandom3.h"
 
#include "JTools/JTable1D.hh"
#include "JTools/JGrid.hh"
#include "JTools/JElement.hh"
#include "JTools/JPolint.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JMapList.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JStats.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * 3D function.
   * smooth in x and y, arbitrary in z
   * 
   * \param  x       x value
   * \param  y       y value
   * \param  z       z value
   * \return         function value
   */
  inline double f3(const double x,
                   const double y,
                   const double z)
  { 
    return sin(z/2)*exp(-fabs(z)) * (x*x - x*y - y*y + x -2*y + 3) / 25;
  }
}
 
 
/**
 * \file
 *
 * Example program to test interpolation on 2D grid between 1D tables.
 * \author hbaetsen
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  int     numberOfEvents;
  int     debug;
  double  precision;
 
  try {
 
    JParser<> zap("Example program to test interpolation between 1D tables on 2D grid.");
 
    zap['n'] = make_field(numberOfEvents);
    zap['d'] = make_field(debug)            = 2;
    zap['p'] = make_field(precision)        = 1e-6;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const int NX = 19;    // 11
  const int NY = 35;    // 11
  const int NZ = 16384; // 21
 
  const JGrid<double> X(NX,  -5.0,  +5.0);
  const JGrid<double> Y(NY,  -5.0,  +5.0);
  const JGrid<double> Z(NZ,  -10.0,  +10.0);
 
 
  // Define 2D 2nd order interpolants
  typedef JTable1D<NZ>                                             JTable1D_t;
  typedef JElement2D<double, JTable1D_t>                           JElement2D_t;
  typedef JPolintFunction1D<2, JElement2D_t, JGridCollection>      JFunction1D_t;
  typedef JMapList<JPolint2FunctionalGridMap>                      JMaplist_t;
  typedef JMultiFunction<JFunction1D_t, JMaplist_t>                JMultiFunction_t;
 
  JMultiFunction_t g3;
 
  for (int ix = 0; ix != X.getSize(); ++ix) {
    for (int iy = 0; iy != Y.getSize(); ++iy) {
 
      DEBUG("table[" << ix << ", " << iy << "]" << endl);
 
      const double x = X.getX(ix);
      const double y = Y.getX(iy);
 
      JTable1D_t& table = g3[x][y];
 
      for (int iz = 0; iz != Z.getSize(); ++iz) {
 
        const double z = Z.getX(iz);  
        table[iz] = f3(x,y,z);
 
        DEBUG(' ' << FIXED(8,3) << table[iz]);
        
      }
      DEBUG(endl);
    }
  }
 
  g3.compile();
 
  // Random events
  JTimer    t1("JTable1D interpolant");
 
  for (int i = 0; i != numberOfEvents; ++i) {
 
    const double x = gRandom->Uniform(X.getXmin(), X.getXmax());
    const double y = gRandom->Uniform(Y.getXmin(), Y.getXmax());
    DEBUG("x = " << x << " y = " << y << endl);
 
    t1.start();
    JTable1D_t res = g3(x,y);
    t1.stop();
 
    int valid = 1;
 
    DEBUG(FIXED(9,3) << "f(x,y,z)" << FIXED(9,3) << "interp" << endl);
 
    for (int iz = 0; iz != Z.getSize(); ++iz) {
      const double z = Z.getX(iz);  
      double val_f3 = f3(x,y,z);
 
      DEBUG(FIXED(9,3) << val_f3 << FIXED(9,3) << res[iz] << endl);
      if (fabs((val_f3 - res[iz]) / res[iz]) > precision) {valid = 0;}
 
    }
 
    ASSERT(valid, FIXED(5,0) << i << ": [x, y] = " << FIXED(9,3) << x << FIXED(9,3) << y);
 
  }
 
  for (JTimer buffer[] = { t1, JTimer() }, *i = buffer; i->getTitle() != ""; ++i) {
    i->print(cout, true, micro_t);
  }
 
}