JPlotNPE-PDG.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JNPETable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JROOT/JRootToolkit.hh"

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


typedef std::pair<double,double> orientation;

std::istream& operator>>(std::istream& in, orientation& x) { return in >> x.first >> x.second; }
std::ostream& operator<<(std::ostream& out, const orientation& x) { return out << x.first << ' ' << x.second; }


/**
 * \file
 *
 * Program to compare integrals of PDF of Cherenkov light from EM-shower using interpolation tables.
 * \author lquinn
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<Double_t> JHistogram_t;

  string         pdfFile;
  string         cdfFile;
  string         outputFile;
  orientation    dir;
  JHistogram_t   x;
  JHistogram_t   y;
  double         E;
  int            numberOfPoints;
  int            debug;

  try { 

    JParser<> zap("Program to plot PDF of Cherenkov light from EM-shower using interpolation tables.");

    zap['P'] = make_field(pdfFile);
    zap['C'] = make_field(cdfFile);
    zap['o'] = make_field(outputFile);
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['x'] = make_field(x,         "histogram x-binning")        = JHistogram_t();
    zap['y'] = make_field(y,         "histogram y-binning")        = JHistogram_t();
    zap['E'] = make_field(E,         "Energy [GeV]")               = 1.0;
    zap['N'] = make_field(numberOfPoints)                          = 10;
    zap['d'] = make_field(debug)                                   = 0;

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


  const double theta = dir.first;
  const double phi   = dir.second;
  
  typedef JSplineFunction1S_t                                     JFunction1D_t;
  typedef
    JMapList<JPolint1FunctionalMap,
    JMapList<JPolint1FunctionalMap,
    JMapList<JPolint1FunctionalGridMap,
    JMapList<JPolint1FunctionalGridMap> > > >                     JMapList_t;
  typedef JPDFTable<JFunction1D_t,  JMapList_t>                   JPDF_t;
  typedef JNPETable<double, double, JMapList_t>                   JNPE_t;
  typedef JCDFTable<JFunction1D_t,  JMapList_t>                   JCDF_t;

  JDistance<double>::precision = 1.0e-10;

  JPDF_t pdf;
  JNPE_t npe;
  JCDF_t cdf;

  try {

    NOTICE("loading input from file " << pdfFile << "... " << flush);
      
    pdf.load(pdfFile.c_str());

    pdf.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));

    npe = JNPE_t(pdf);

    NOTICE("OK" << endl);

    NOTICE("loading input from file " << cdfFile << "... " << flush);
      
    cdf.load(cdfFile.c_str());

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

  if (!x.is_valid()) { x = JHistogram_t(250,  0.0, 250.0); }
  if (!y.is_valid()) { y = JHistogram_t(200, -1.0,  +1.0); }

  TH2D h0("pdf", NULL, x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit(), y.getNumberOfBins(), y.getLowerLimit(), y.getUpperLimit());
  TH2D h1("npe", NULL, x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit(), y.getNumberOfBins(), y.getLowerLimit(), y.getUpperLimit());
  TH2D h2("cdf", NULL, x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit(), y.getNumberOfBins(), y.getLowerLimit(), y.getUpperLimit());

  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
    for (int j = 1; j <= h0.GetNbinsY(); ++j) {

      const double D  = h0.GetXaxis()->GetBinCenter(i);
      const double cd = h0.GetYaxis()->GetBinCenter(j);

      try {

        const double y0 = get_integral(pdf(D, cd, theta, phi, 1.0e3));
        const double y1 = npe             (D, cd, theta, phi);
        const double y2 = cdf.getNPE      (D, cd, theta, phi);

        DEBUG("npe "
              << FIXED(5,1)       << D  << ' ' 
              << FIXED(5,1)       << cd << ' '
              << SCIENTIFIC(12,3) << y0 << ' '
              << SCIENTIFIC(12,3) << y1 << ' '
              << SCIENTIFIC(12,3) << y2 << endl);
        
        h0.SetBinContent(i, j, E * y0);
        h1.SetBinContent(i, j, E * y1);
        h2.SetBinContent(i, j, E * y2);
      }
      catch(const exception& error) {
        ERROR(error.what());
      }
    }
  }
  
  if (outputFile != "") {

    TFile out(outputFile.c_str(), "recreate");

    out << h0 << h1 << h2;

    out.Write();
    out.Close();
  }
}