JMultiPDF.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TMath.h"
#include "TRandom.h"

#include "JTools/JHistogram1D_t.hh"
#include "JTools/JHistogramMap_t.hh"
#include "JTools/JMultiHistogram.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JMultiPDF.hh"
#include "JTools/JToolsToolkit.hh"
#include "JTools/JQuadrature.hh"
#include "JTools/JQuantile.hh"
#include "JIO/JFileStreamIO.hh"
#include "JLang/JObjectIO.hh"

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


namespace {

  /**
   * Normalised function in 2D.
   *
   * \param  x       x value
   * \param  y       y value
   * \return         function value
   */
  inline Double_t g2(const Double_t x,
                     const Double_t y)
  {
    return (TMath::Gaus(x, 0.0, 1.0, kTRUE) *
            TMath::Gaus(y, 0.0, 1.0, kTRUE));
  }


  /**
   * Function in 4D.
   *
   * \param  x0      x0 value
   * \param  x1      x1 value
   * \param  x2      x2 value
   * \param  x3      x3 value
   * \return         function value
   */
  inline Double_t g4(const Double_t x0,
                     const Double_t x1,
                     const Double_t x2,
                     const Double_t x3)
  {
    return g2(x2 - x0, x3 - x1);
  }
}

  
/**
 * \file
 *
 * Example program to test conditional probability density function \f$P(x_0,x_1\|x_2,x_3)\f$ using JTOOLS::JMultiPDF and JTOOLS::JMultiFunction.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  string  inputFile;
  string  outputFile;
  int     numberOfEvents;
  int     numberOfBins;
  int     debug;

  try {

    JParser<> zap("Example program to test conditional probability density function.");

    zap['f'] = make_field(inputFile)      = "";
    zap['o'] = make_field(outputFile)     = "";
    zap['n'] = make_field(numberOfEvents);
    zap['N'] = make_field(numberOfBins)   = 15;
    zap['d'] = make_field(debug)          = 2;

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


  if (numberOfEvents <= 0) {
    FATAL("No events." << endl);
  }

  using namespace JPP;


  const double xmin = -10.0;
  const double xmax = +10.0;
  const double dx   = (xmax - xmin) / (numberOfBins - 1);


  if (outputFile != "") {

    typedef JMultiHistogram<JHistogram1D_t, 
                            JMAPLIST<JHistogramGridMap_t,
                                     JHistogramGridMap_t,
                                     JHistogramMap_t>::maplist> JMultiHistogram_t;

    JMultiHistogram_t histogram;

    const JGaussHermite bounds(numberOfBins);

    for (double x0 = xmin; x0 <= xmax + 0.5*dx; x0 += dx) {
      for (double x1 = xmin; x1 <= xmax + 0.5*dx; x1 += dx) {
        for (JGaussHermite::const_iterator x2 = bounds.begin(); x2 != bounds.end(); ++x2) {
          for (JGaussHermite::const_iterator x3 = bounds.begin(); x3 != bounds.end(); ++x3) {
            histogram[x0][x1][x0 + x2->getX()][x1 + x3->getX()] = 0.0;
          }
        }
      }
    }    

    // fill

    JTimer timer("histogram");

    for (int i = 0; i != numberOfEvents; ++i) {

      if (i%1000 == 0) {
        STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);
      }

      const double x0 = gRandom->Uniform(xmin, xmax);
      const double x1 = gRandom->Uniform(xmin, xmax);
      const double x2 = x0 + gRandom->Gaus(0.0, 1.0);
      const double x3 = x1 + gRandom->Gaus(0.0, 1.0);
      const double w  = 1.0;

      timer.start();

      histogram.fill(x0, x1, x2, x3, w);

      timer.stop();
    }
    STATUS(endl);

    if (debug >= status_t) {
      timer.print(cout, true, micro_t);
    }

    try {

      NOTICE("storing output to file " << outputFile << "... " << flush);      

      JLANG::store<JIO::JFileStreamWriter>(outputFile.c_str(), histogram);

      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }


  if (inputFile != "") {

    typedef JMultiHistogram<JHistogram1D_t, 
                            JMAPLIST<JHistogramMap_t>::maplist>          JHistogram2D_t;
  
    typedef JMultiHistogram<JHistogram2D_t, 
                            JMAPLIST<JHistogramGridMap_t,
                                     JHistogramGridMap_t>::maplist>      JMultiHistogram_t;

    JMultiHistogram_t histogram;

    try {

      NOTICE("loading input to file " << inputFile << "... " << flush);      

      JLANG::load<JIO::JFileStreamReader>(inputFile.c_str(), histogram);

      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }


    typedef JMultiPDF<JPolint1Function1D_t,
                      JMAPLIST<JPolint1FunctionalMap>::maplist>          JFunction2D_t;     // PDF

    typedef JMultiFunction<JFunction2D_t,
                           JMAPLIST<JPolint1FunctionalGridMap,
                                    JPolint1FunctionalGridMap>::maplist> JMultiFunction_t;  // interpolation


    JMultiFunction_t pdf(histogram);


    // test

    JQuantile Q;

    for (int i = 0; i != numberOfEvents; ++i) {
      
      const double x0 = gRandom->Uniform(xmin, xmax);
      const double x1 = gRandom->Uniform(xmin, xmax);
      const double x2 = x0 + gRandom->Gaus(0.0, 1.0);
      const double x3 = x1 + gRandom->Gaus(0.0, 1.0);
      
      try {
        
        const double u = g4 (x0, x1, x2, x3);
        const double v = pdf(x0, x1, x2, x3);
        
        Q.put(u - v);
      }
      catch(const std::exception& error) {}
    }
    

    JQuantile V;

    for (JMultiFunction_t::super_const_iterator i = pdf.super_begin(); i != pdf.super_end(); ++i) {
      V.put(getIntegral((*i).getValue().getMultiFunction()));
    }

    cout << "normalisation " 
         << FIXED(6,4) << V.getMean()  << " +/- " 
         << FIXED(6,4) << V.getSTDev() << endl;
    cout << "efficiency    " 
         << FIXED(5,3) << (double) Q.getCount() / (double) numberOfEvents << endl;
    
    cout << "mean " << SCIENTIFIC(10,2) << Q.getMean()  << endl;
    cout << "RMS  " << SCIENTIFIC(10,2) << Q.getSTDev() << endl;
  }
}