JDiffPDF.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <algorithm>
 
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTable.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
  /**
   * Compare two function values.
   * If average function value less than one, absolute comparison; else relative comparison.
   *
   * \param  u             first  value
   * \param  v             second value
   * \param  precision     precision
   * \return               true if comparable; else false
   */
  inline bool compare(const double u, const double v, const double precision)
  {
    using namespace std;
 
    return fabs(u - v)  <=  precision * max(1.0, 0.5*(u+v));
  }
}
 
/**
 * \file
 *
 * Auxiliary program to compare PDF tables of the arrival time of the Cherenkov light from a muon.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string      inputFileA;
  string      inputFileB;
  double      precision;
  int         debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to compare PDF tables of the arrival time of the Cherenkov light from a muon.");
 
    zap['a'] = make_field(inputFileA);
    zap['b'] = make_field(inputFileB);
    zap['p'] = make_field(precision)        = 1.0e-6;
    zap['d'] = make_field(debug)            = 0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  typedef JSplineFunction1D_t                                     JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;
 
  typedef JFunction1D_t::argument_type                            argument_type;
  typedef JArray<3, argument_type>                                JArray_t;
 
  JPDF_t pdfA;
  JPDF_t pdfB;
 
  const JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(0.0));
 
  struct pair_type {
    const string file_name;
    JPDF_t&      pdf;
  } buffer[] = { 
    { inputFileA, pdfA },
    { inputFileB, pdfB }
  };
 
  for (const pair_type& i : buffer) {
 
    try {
 
      NOTICE("loading input from file " << i.file_name << "... " << flush);
 
      i.pdf.load(i.file_name.c_str());
      i.pdf.compile();
      i.pdf.setExceptionHandler(supervisor);
 
      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }    
 
  bool ok = true;
 
  JPDF_t::super_const_iterator i = pdfA.super_begin();
  JPDF_t::super_const_iterator j = pdfB.super_begin();
 
  for ( ; i != pdfA.super_end() && 
          j != pdfB.super_end(); ++i, ++j) {
 
    const double Wa = pdfA.transformer->getWeight(JArray_t((*i).getKey()));
    const double Wb = pdfB.transformer->getWeight(JArray_t((*j).getKey()));
 
    if (fabs(i->first                 - j->first)                  >  precision ||
        fabs(i->second->first         - j->second->first)          >  precision ||
        fabs(i->second->second->first - j->second->second->first)  >  precision ||
        fabs(Wa                       - Wb)                        >  precision) {
 
      DEBUG("a[] " 
            << i->first                 << ' '
            << i->second->first         << ' '
            << i->second->second->first << ' '
            << Wa                       << endl);
 
      DEBUG("b[] " 
            << j->first                 << ' '
            << j->second->first         << ' '
            << j->second->second->first << ' '
            << Wb                       << endl);
 
      ok = false;
    }
 
    const JFunction1D_t& fa = (*i).getValue();
    const JFunction1D_t& fb = (*j).getValue();
 
    JFunction1D_t::const_iterator p = fa.begin();
    JFunction1D_t::const_iterator q = fb.begin();
 
    while (p != fa.end() && 
           q != fb.end()) {
      
      if (fabs(p->getX() - q->getX()) > precision) {
 
        DEBUG("a.x " << p->getX() << endl);
        DEBUG("b.x " << q->getX() << endl);
 
        ok = false;
      }
 
      if (!compare(p->getY(), q->getY(), precision)) {
 
        DEBUG("a.y " << p->getX() << ' ' << p->getY() << endl);
        DEBUG("b.y " << q->getX() << ' ' << q->getY() << endl);
 
        ok = false;
      }
 
      if      (p->getX() < q->getX())
        ++p;
      else if (q->getX() < p->getX())
        ++q;
      else {
        ++p;
        ++q;
      }
    }
 
    for ( ; p != fa.end(); ++p) {
 
      DEBUG("a.x " << p->getX() << endl);
 
      ok = false;
    }
 
    for ( ; q != fb.end(); ++q) {
 
      DEBUG("b.x " << q->getX() << endl);
 
      ok = false;
    }
  }
 
  for ( ; i != pdfA.super_end(); ++i) {
 
    DEBUG("a[] " 
          << i->first                 << ' '
          << i->second->first         << ' '
          << i->second->second->first << endl); 
    
    ok = false;
  }
 
  for ( ; j != pdfB.super_end(); ++j) {
 
    DEBUG("b[] " 
          << j->first                 << ' '
          << j->second->first         << ' '
          << j->second->second->first << endl);
 
    ok = false;
  }
 
  if (!ok) {
    ERROR(inputFileA << " and " << inputFileB << " differ." << endl);
  }
 
  return (ok ? 0 : 1);
}