JPlotPDG.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
 
#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/JGeanz.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Auxiliary program to plot PDF of Cherenkov light from EM-shower using interpolation tables.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JAbstractHistogram<double> JHistogram_t;
 
  vector<string> inputFile;
  string         outputFile;
  double         E;
  double         D;
  double         cd;
  JAngle3D       dir;
  double         TTS_ns;
  size_t         profile;
  bool           cms;
  JHistogram_t   histogram;
  int            debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to plot PDF of Cherenkov light from EM-shower using interpolation tables.");
 
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)                              = "";
    zap['E'] = make_field(E,         "shower energy [GeV]")        = 1.0;
    zap['R'] = make_field(D,         "distance [m]");
    zap['c'] = make_field(cd,        "cosine emission angle");
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['T'] = make_field(TTS_ns,    "PMT time smearing [ns]")     = 0.0;   // [ns]
    zap['N'] = make_field(profile,   "steps shower profile")       = 0;
    zap['C'] = make_field(cms,       "center shower profile");
    zap['H'] = make_field(histogram, "histogram binning")          = JHistogram_t();
    zap['d'] = make_field(debug)                                   = 0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  typedef JSplineFunction1S_t                                     JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;
 
  const int N = inputFile.size();
 
  int    type[N];
  JPDF_t pdf [N];
 
  try {
 
    for (int i = 0; i != N; ++i) {
 
      NOTICE("loading input from file " << inputFile[i] << "... " << flush);
      
      type[i] = getPDFType(inputFile[i]);
 
      pdf [i].load(inputFile[i].c_str());
 
      pdf [i].setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
      
      if (TTS_ns > 0.0) {
        pdf[i].blur(TTS_ns);
      }
      
      NOTICE("OK" << endl);
    }
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
 
 
  if (outputFile == "") {
 
    for (double dt; cin >> dt; ) {
 
      for (int i = 0; i != N; ++i) {
 
        JFunction1D_t::result_type y = pdf[i](D, cd, dir.getTheta(), dir.getPhi(), dt);
      
        cout << setw(2)         << type[i]               << ' '
             << SCIENTIFIC(7,1) << E                     << ' '
             << FIXED(5,1)      << D                     << ' '
             << FIXED(5,2)      << cd                    << ' '
             << FIXED(5,2)      << dir.getTheta()        << ' '
             << FIXED(5,2)      << dir.getPhi()          << ' '
             << FIXED(5,1)      << dt                    << ' '
             << SCIENTIFIC(9,3) << get_value(y)          << ' '
             << SCIENTIFIC(9,3) << get_derivative(y)     << ' '
             << SCIENTIFIC(9,3) << get_integral(y)       << ' '
             << SCIENTIFIC(9,3) << get_total_integral(y) << endl;
      }
    }
 
    return 0;
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
    
  int function = 0;
    
  if (inputFile.size() == 1 &&
      inputFile.begin()->find(getLabel(SCATTERED_LIGHT_FROM_EMSHOWER)) == string::npos) {
    function = 1;
  }
    
  //const double t0   =  D * getIndexOfRefraction() / C;        // time [ns]
  const double t0   =  0.0;                                   // time [ns]
 
  if (!histogram.is_valid()) {
 
    if (function == 1) {
 
      histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);
 
      histogram.setBinWidth(0.1);
 
    } else {
 
      histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);
 
      histogram.setBinWidth(0.5);
    }
  }
 
  TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
 
  const double D0 = D;
  const double z0 = D0 * cd;
  const double R  = D0 * sqrt((1.0 + cd) * (1.0 - cd));
 
  double W = 1.0;
 
  vector<double> Z;
 
  if (profile >= 2) {
    
    const double z1 = (cms ? geanz.getMaximum(E) :0.0);
 
    W = 1.0 / (double) profile;
 
    for (double x = 0.5*W; x < 1.0; x += W) {
      Z.push_back(geanz.getLength(E, x) - z1);
    }
 
  } else {
 
    W = 1.0;
    Z = { 0.0 };
  }
 
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
      
    const double dt = h0.GetBinCenter(i) - t0;
 
    JFunction1D_t::result_type Y = JMATH::zero;
 
    for (const double z1 : Z) {
 
      const double z = z0 - z1;
      
      D  = sqrt(R*R + z*z);
      cd = z / D;
 
      const double t1 = dt  -  (z1 + (D - D0) * getIndexOfRefraction()) / C;
 
      for (int j = 0; j != N; ++j) {
        Y += pdf[j](D, cd, dir.getTheta(), dir.getPhi(), t1) * E * W;
      }
 
      DEBUG("sum: "
            << FIXED(7,2) << dt << ' '
            << FIXED(7,2) << z1 << ' '
            << FIXED(7,2) << t1 << ' '
            << SCIENTIFIC(12,3) << get_value(Y) << endl);
    }
    DEBUG(endl);
 
    h0.SetBinContent(i, get_value(Y));
  }
 
  out.Write();
  out.Close();
}