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