JPhysics/JShowerNPE.cc

Go to the documentation of this file.

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

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

#include "JTools/JAbstractHistogram.hh"
#include "JPhysics/JNPE_t.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Program to plot 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> JHistogram_t;

  string         pdfFile;
  string         outputFile;
  JAngle3D       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['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]")               = 0.0;
    zap['N'] = make_field(numberOfPoints)                          = 0;
    zap['d'] = make_field(debug)                                   = 0;

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


  const JShowerNPE_t npe(pdfFile, numberOfPoints);


  TFile out(outputFile.c_str(), "recreate");
    
  if (!x.is_valid()) { x = JHistogram_t(150,  0.0, 150.0); }
  if (!y.is_valid()) { y = JHistogram_t(200, -1.0,  +1.0); }
    
  TH2D h0("h0", "PDF Projection; D [m]; cos #theta_{0}",
          x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit(),
          y.getNumberOfBins(), y.getLowerLimit(), y.getUpperLimit());

  for (int ix = 1; ix <= h0.GetNbinsX(); ++ix) {
    for (int iy = 1; iy <= h0.GetNbinsY(); ++iy) {

      const double cd = h0.GetYaxis()->GetBinCenter(iy);
      const double D  = h0.GetXaxis()->GetBinCenter(ix);

      const double Y  = npe.calculate(E, D, cd, dir.getTheta(), dir.getPhi());

      h0.SetBinContent(ix, iy, Y/E);
    }
  }

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