JMakePDG.cc File Reference

Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower. 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 shower.

The PDFs are tabulated as a function of (D, cos(theta), theta, phi, t), where:

• D is the distance between the vertex and the PMT;
• cos(theta) the cosine of the photon emission angle;
• (theta, phi) the orientation of the PMT; and
• 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 shower developes along the z-axis and the PMT is located in the x-z plane.

Author

mdejong

Definition in file JMakePDG.cc.

Function Documentation

double getAbsorptionLength ( const double  lambda )

inline

Definition at line 28 of file JMakePDG.cc.

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

double getScatteringLength ( const double  lambda )

inline

Definition at line 34 of file JMakePDG.cc.

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

int main	(	int	argc,
		char **	argv
	)

Definition at line 55 of file JMakePDG.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 shower.");

    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")                   =
      SCATTERED_LIGHT_FROM_MUON_5D,
      DIRECT_LIGHT_FROM_EMSHOWER,
      SCATTERED_LIGHT_FROM_EMSHOWER;
    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 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,
                   JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  typedef JPDFTransformer<4, JFunction1D_t::argument_type>        JFunction4DTransformer_t;
  typedef JArray<4, 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,JFunction4DTransformer_t> > zmap;

  zmap[SCATTERED_LIGHT_FROM_MUON_5D]  = make_pair((fcn) &JPDF::getScatteredLightFromMuon,     JFunction4DTransformer_t(21.5, 2, ng[0], 0.0,   JGeant(JGeanx(0.33, -9.5)), 6e-4, NAMESPACE::getAngularAcceptance, 0.06));
  zmap[DIRECT_LIGHT_FROM_EMSHOWER]    = make_pair((fcn) &JPDF::getDirectLightFromEMshower,    JFunction4DTransformer_t(21.5, 2, ng[0], ng[1], JGeant(JGeanx(0.35, -5.4)), 1e-5, NAMESPACE::getAngularAcceptance, 0.001));
  zmap[SCATTERED_LIGHT_FROM_EMSHOWER] = make_pair((fcn) &JPDF::getScatteredLightFromEMshower, JFunction4DTransformer_t(21.5, 2, ng[0], 0.0,   JGeant(JGeanx(0.55, -4.5)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));

  if (zmap.find(function) == zmap.end()) {
    FATAL("illegal function specifier" << endl);
  }

  fcn                      f           = zmap[function].first;    // PDF
  JFunction4DTransformer_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);
    D.insert(120.00);
    D.insert(150.00);
    D.insert(170.00);
    D.insert(190.00);
    D.insert(210.00);
    D.insert(230.00);
    D.insert(250.00);
    D.insert(270.00);
    D.insert(290.00);
    D.insert(310.00);
  }

  set<double> C;        // cosine emission angle

  JQuadrature qeant(-1.0, +1.0, 30, geanx);

  for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i)
    C.insert(i->getX());

  C.insert(-1.00);
  C.insert(+1.00);


  set<double> X;        // time [ns]

  if (function == DIRECT_LIGHT_FROM_EMSHOWER) {

    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);
    X.insert(250.0);
    X.insert(300.0);
    X.insert(350.0);
    X.insert(400.0);
    X.insert(450.0);
    X.insert(500.0);
    X.insert(600.0);
    X.insert(700.0);
    X.insert(800.0);
    X.insert(900.0);
    X.insert(1200.0);
    X.insert(1500.0);
  }
    
  const double grid  =  5.0;                                      // [deg]
      
  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0));  // azimuth angle unit step size


  for (set<double>::const_iterator d = D.begin(); d != D.end(); ++d) {

    const double D_m = *d;

    for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {

      const double cd  = *c;

      const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));

      for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {

        const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
      
        for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {

          JFunction1D_t& f1 = pdf[D_m][cd][theta][phi];

          const JArray_t array(D_m, cd, theta, phi);

          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, cd, theta, phi, 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 JMakePDG.cc.

double scatteringLengthFactor

Definition at line 25 of file JMakePDG.cc.