JMakePD0.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>
#include <cmath>
 
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuadrature.hh"
 
#include "JPhysics/JPDF.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JPDFSupportkit.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a bright point.
 *
 * The PDFs are tabulated as a function of <tt>(D, cos(theta), t)</tt>, where:
 * - <tt>D</tt> is the distance between the bright point and the PMT;
 * - <tt>cos(theta)</tt> the cosine of the PMT angle;
 * - <tt>t</tt> the arrival time of the light with respect to the Cherenkov hypothesis.
 *
 * The orientation of the PMT is defined in the coordinate system in which
 * the position of the bright point and that of the PMT are along the x-axis and the PMT is oriented in the x-y plane.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string                  outputFile;
  int                     numberOfPoints;
  double                  epsilon;
  JAbsorptionLength       absorptionLength;
  JScatteringLength       scatteringLength;
  JScatteringProbability  scatteringProbability;
  int                     function;
  set<double>             D;                          // distance [m]
  int                     debug;
  
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(absorptionLength));
    properties.insert(gmake_property(scatteringLength));
    properties.insert(gmake_property(scatteringProbability));
 
    JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a bright point.");
 
    zap['@'] = make_field(properties, endl
                          << "possible options absorptionLength: "      << get_keys(absorptionLength)      << endl
                          << "possible options scatteringLength: "      << get_keys(scatteringLength)      << endl
                          << "possible options scatteringProbability: " << get_keys(scatteringProbability) << endl)  = JPARSER::initialised();
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfPoints,         "points for integration")     = 25;
    zap['e'] = make_field(epsilon,                "precision for integration")  = 1.0e-10;
    zap['F'] = make_field(function,               "PDF type")                   =
      DIRECT_LIGHT_FROM_BRIGHT_POINT,
      SCATTERED_LIGHT_FROM_BRIGHT_POINT;
    zap['D'] = make_field(D,                      "distance [m]")               = JPARSER::initialised();
    zap['d'] = make_field(debug)                  = 0;
 
    zap['F'] = JPARSER::not_initialised();
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  typedef double (JPDF::*fcn)(const double,
                              const double,
                              const double) const;
 
 
  // set global parameters
 
  const double P_atm = NAMESPACE::getAmbientPressure();
  const double wmin  = getMinimalWavelength();
  const double wmax  = getMaximalWavelength();
  
 
  const JPDF_C
    pdf_c(NAMESPACE::getPhotocathodeArea(),
          NAMESPACE::getQE,
          NAMESPACE::getAngularAcceptance,
          JAbsorptionLength::getAbsorptionLength,
          JScatteringLength::getScatteringLength,
          JScatteringProbability::getScatteringProbability,
          P_atm,
          wmin,
          wmax,
          numberOfPoints,
          epsilon);
 
 
  typedef JSplineFunction1D_t                                     JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;
 
  typedef JPDFTransformer<2, JFunction1D_t::argument_type>        JFunction2DTransformer_t;
  typedef JArray<2, JFunction1D_t::argument_type>                 JArray_t;
 
  JPDF_t pdf;
 
 
  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
  
  const double ng[] = {
    pdf_c.getIndexOfRefractionGroup(wmax),
    pdf_c.getIndexOfRefractionGroup(wmin)
  };
  
  map<int, pair<fcn,JFunction2DTransformer_t> > zmap;
 
  zmap[DIRECT_LIGHT_FROM_BRIGHT_POINT]    = make_pair((fcn) &JPDF::getDirectLightFromBrightPoint,    JFunction2DTransformer_t(21.5, 2, ng[0], ng[1]));
  zmap[SCATTERED_LIGHT_FROM_BRIGHT_POINT] = make_pair((fcn) &JPDF::getScatteredLightFromBrightPoint, JFunction2DTransformer_t(21.5, 2, ng[0], 0.0));
 
  if (zmap.find(function) == zmap.end()) {
    FATAL("illegal function specifier" << endl);
  }
 
  fcn                      f           = zmap[function].first;    // PDF
  JFunction2DTransformer_t transformer = zmap[function].second;   // transformer
 
 
  if (D.empty()) {
    D.insert(  0.10);
    D.insert(  0.50);
    D.insert(  1.00);
    D.insert(  5.00);
    D.insert( 10.00);
    D.insert( 20.00);
    D.insert( 30.00);
    D.insert( 40.00);
    D.insert( 50.00);
    D.insert( 60.00);
    D.insert( 70.00);
    D.insert( 80.00);
    D.insert( 90.00);
    D.insert(100.00);
  }
 
  set<double> X;        // time [ns]
 
  if (function == DIRECT_LIGHT_FROM_BRIGHT_POINT) {
 
    for (double buffer[] = { 0.0, 0.005, 0.01, 0.015, -1 }, *x = buffer; *x >= 0; ++x) {
      X.insert(0.0 + *x);
      X.insert(1.0 - *x);
    }
 
    for (double x = 0.02; x < 0.99; x += 0.01)
      X.insert(x);
 
  } else {
 
    X.insert(  0.00);
    X.insert(  0.10);
    X.insert(  0.20);
    X.insert(  0.30);
    X.insert(  0.40);
    X.insert(  0.50);
    X.insert(  0.60);
    X.insert(  0.70);
    X.insert(  0.80);
    X.insert(  0.90);
    X.insert(  1.00);
    X.insert(  1.00);
    X.insert(  1.10);
    X.insert(  1.20);
    X.insert(  1.30);
    X.insert(  1.40);
    X.insert(  1.50);
    X.insert(  1.60);
    X.insert(  1.70);
    X.insert(  1.80);
    X.insert(  1.90);
    X.insert(  2.00);
    X.insert(  2.20);
    X.insert(  2.40);
    X.insert(  2.60);
    X.insert(  2.80);
    X.insert(  3.00);
    X.insert(  3.25);
    X.insert(  3.50);
    X.insert(  3.75);
    X.insert(  4.00);
    X.insert(  4.25);
    X.insert(  4.50);
    X.insert(  4.75);
    X.insert(  5.0);
    X.insert(  6.0);
    X.insert(  7.0);
    X.insert(  8.0);
    X.insert(  9.0);
    X.insert( 10.0);
    X.insert( 15.0);
    X.insert( 20.0);
    X.insert( 25.0);
    X.insert( 30.0);
    X.insert( 40.0);
    X.insert( 50.0);
    X.insert( 60.0);
    X.insert( 70.0);
    X.insert( 80.0);
    X.insert( 90.0);
    X.insert(100.0);
    X.insert(120.0);
    X.insert(140.0);
    X.insert(160.0);
    X.insert(180.0);
    X.insert(200.0);
  }
    
 
  for (set<double>::const_iterator d = D.begin(); d != D.end(); ++d) {
 
    const double D_m = *d;
 
    for (double dc = 0.1, ct = -1.0; ct < +1.0 + 0.5*dc; ct += dc) {
 
      JFunction1D_t& f1 = pdf[D_m][ct];
 
      const JArray_t array(D_m, ct);
 
      double t_old = transformer.getXn(array, *X.begin());
      double y_old = 0.0;
 
      for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
          
        const double t = transformer.getXn(array, *x);
        const double y = (pdf_c.*f)(D_m, ct, t);
 
        if (y != 0.0) {
 
          if (*x < 0.0) {
            WARNING("dt < 0 " << *x << ' ' << D_m << ' ' << t << ' ' << y << endl);
          }
 
          if (y_old == 0.0) {
            f1[t_old] = y_old;
          }
 
          f1[t] = y;
              
        } else {
              
          if (y_old != 0.0) {
            f1[t] = y;
          }
        }
            
        t_old = t;
        y_old = y;
      }
    }
  }
 
  pdf.transform(transformer);
  pdf.compile();
 
  NOTICE("OK" << endl);
 
  try {
 
    NOTICE("storing output to file " << outputFile << "... " << flush);
 
    pdf.store(outputFile.c_str());
 
    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
}