JMakePD0.cc File Reference

Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a bright point. More...

#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 "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
double	getAbsorptionLength (const double lambda)
double	getScatteringLength (const double lambda)
int	main (int argc, char **argv)

Variables
double	absorptionLengthFactor
	Scaling of absorption and scattering length. More...
double	scatteringLengthFactor

Detailed Description

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 (D, cos(theta), t), where:

• D is the distance between the bright point and the PMT;
• cos(theta) the cosine of the PMT angle;
• t 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

Definition in file JMakePD0.cc.

Function Documentation

double getAbsorptionLength ( const double  lambda )

inline

Definition at line 28 of file JMakePD0.cc.

{
  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
}

double getScatteringLength ( const double  lambda )

inline

Definition at line 34 of file JMakePD0.cc.

{
  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 54 of file JMakePD0.cc.

{
  using namespace std;
  using namespace JPP;

  string              outputFile;
  int                 numberOfPoints;
  double              epsilon;
  int                 function;
  set<double>         D;                          // distance [m]
  int                 debug;
  
  try {

    JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a bright point.");

    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['A'] = make_field(absorptionLengthFactor, "scaling factor")             = 1.0;
    zap['S'] = make_field(scatteringLengthFactor, "scaling factor")             = 1.0;
    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,
          getAbsorptionLength,
          getScatteringLength,
          NAMESPACE::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);
  }
}

Variable Documentation

double absorptionLengthFactor

Scaling of absorption and scattering length.

Definition at line 24 of file JMakePD0.cc.

double scatteringLengthFactor

Definition at line 25 of file JMakePD0.cc.