JPlotNPE1D.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 "TH1D.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/JPDFTypes.hh"
#include "JPhysics/JNPETable.hh"
#include "JPhysics/JGeanz.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 jseneca
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<double> JHistogram_t;

  vector<string> inputFile;
  string         outputFile;
  JAngle3D       dir;
  double         D;
  JHistogram_t   x;
  double         E;
  int            numberOfPoints;
  int            debug;

  try { 

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

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)                              = "";
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['R'] = make_field(D,         "distance from vertex to PMT [m]");
    zap['x'] = make_field(x,         "histogram x-binning")        = JHistogram_t();
    zap['E'] = make_field(E,         "Energy [GeV]")               = 0.0;
    zap['N'] = make_field(numberOfPoints)                          = 10;
    zap['d'] = make_field(debug)                                   = 2;

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


  vector<double> Z;

  if (E > 0.0) {
    for (int i = 0; i != numberOfPoints; ++i) {
      Z.push_back(geanz.getLength(E, (i + 0.5) / (double) numberOfPoints));  
    }
  }

 
  typedef JPolint0Function1D_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;

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

  const int N = inputFile.size();

  JPDF_t pdf[N];
  JNPE_t npe[N];

  try {

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

      NOTICE("loading input from file " << inputFile[i] << "... " << flush);
      
      pdf[i].load(inputFile[i].c_str());

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

      npe[i] = JNPE_t(pdf[i]);

      npe[i].setExceptionHandler(new JFunction<double, double>::JDefaultResult(JMATH::zero));

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


  if (outputFile != "") {

    TFile out(outputFile.c_str(), "recreate");
    
    if (!x.is_valid()) { x = JHistogram_t(100000, -1.0, +1.0); }
    
    TH1D h0("h0", "PDF Projection; D [m]; P [npe]", x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit());


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

      const double cd = h0.GetBinCenter(ix);
    
      double Y = 0.0;
    
      if (!Z.empty()) {
      
        const double W = 1.0 / (double) Z.size();
      
        for (vector<double>::const_iterator z = Z.begin(); z != Z.end(); ++z) {
        
          const double __D  = sqrt(D*D - 2.0*(D*cd)*(*z) + (*z)*(*z));
          const double __cd = (D * cd - (*z)) / __D;
        
          for (int i = 0; i != N; ++i) {
            try {
              Y += W * npe[i](__D, __cd, dir.getTheta(), dir.getPhi());
            }
            catch(const exception& error) {}
          }
        }
      
      } else {
      
        for (int i = 0; i != N; ++i) {
          Y += npe[i](D, cd, dir.getTheta(), dir.getPhi());
        }
      }
    
      h0.SetBinContent(ix, Y);
    }

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

  } else {

    for (double cd; ; ) {

      cout << "> " << flush;
      cin  >> D >> cd;

      for (int i = 0; i != N; ++i) {
        cout << i << ' ' << FIXED(9,5) << npe[i](D, cd, dir.getTheta(), dir.getPhi()) << endl;
      }
    }
  }
}