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