JDiffPDF.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTable.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \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;

  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-4;
    zap['d'] = make_field(debug)            = 0;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace JPP;

  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 JEntry_t {

    const char* file_name;
    JPDF_t&     pdf;

  } buffer[] = { 
    { inputFileA.c_str(), pdfA },
    { inputFileB.c_str(), pdfB }
  };


  for (JEntry_t* p = buffer; p != buffer + sizeof(buffer)/sizeof(buffer[0]); ++p) {

    try {
      
      NOTICE("loading input from file " << p->file_name << "... " << flush);
      
      p->pdf.load(p->file_name);
      p->pdf.compile();
      p->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.f[] " << p->getX() << endl);
        DEBUG("b.f[] " << q->getX() << endl);

        ok = false;
      }

      if (fabs(p->getY() - q->getY()) > precision) {

        DEBUG("a.f() " << p->getX() << ' ' << p->getY() << endl);
        DEBUG("b.f() " << 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.f() " << p->getX() << endl);

      ok = false;
    }

    for ( ; q != fb.end(); ++q) {

      DEBUG("b.f() " << 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);
  }
}