JHistogram3D.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/JMultiSet.hh"
#include "JTools/JQuadrature.hh"
#include "JTools/JToolsToolkit.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JMultiPDF.hh"
#include "JIO/JFileStreamIO.hh"
#include "JLang/JObjectIO.hh"

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


namespace {

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

  
/**
 * \file
 *
 * Example program to test JTOOLS::JMultiHistogram.
 * \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 3D histogram.");

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

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


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

  using namespace JPP;


  typedef JMultiPDF<JPolint1Function1D_t,
                    JMAPLIST<JPolint1FunctionalMap,
                             JPolint1FunctionalMap>::maplist>         JMultiPDF_t;


  if (outputFile != "") {

    typedef JMAPLIST<JHistogramMap_t,
                     JHistogramMap_t>::maplist                        JMapList_t;

    typedef JMultiHistogram<JHistogram1D_t, JMapList_t>               JMultiHistogram_t;


    JMultiHistogram_t histogram;

    configure(histogram, make_set(JGaussHermite(numberOfBins)));


    // fill

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

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

      const double x = gRandom->Gaus(0.0, 1.0);
      const double y = gRandom->Gaus(0.0, 1.0);
      const double z = gRandom->Gaus(0.0, 1.0);
      const double w = 1.0;

      histogram.fill(x, y, z, w);
    }
    STATUS(endl);

    histogram.div((double) numberOfEvents);


    try {

      JMultiPDF_t pdf(histogram);

      JLANG::store<JIO::JFileStreamWriter>(outputFile.c_str(), pdf);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }


  if (inputFile != "") {

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

    JQuantile Q;

    for (int i = 0; i != numberOfEvents; ++i) {
      
      const double x = gRandom->Gaus(0.0, 1.0);
      const double y = gRandom->Gaus(0.0, 1.0);
      const double z = gRandom->Gaus(0.0, 1.0);
      
      try {
        
        const double u = g3 (x, y, z);
        const double v = pdf(x, y, z);
        
        Q.put(u - v);
      }
      catch(const std::exception& error) {}
    }
    
    typedef JMultiFunction<JPolint1Function1D_t,
      JMapList<JPolint1FunctionalMap,
      JMapList<JPolint1FunctionalMap> > >                               JMultiFunction_t;

    cout << "normalisation " << FIXED(5,3) << getIntegral(static_cast<const JMultiFunction_t&>(pdf)) << endl;
    cout << "efficiency    " << FIXED(5,3) << (double) Q.getCount() / (double) numberOfEvents        << endl;

    Q.print(cout);
  }
}