JMakePDF.cc File Reference

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

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.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 muon.

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

• R is the minimal distance of approach of the muon to the PMT;
• (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 muon travels along the z-axis and the PMT is located in the x-z plane.

Author

mdejong

Definition in file JMakePDF.cc.

Function Documentation

double getAbsorptionLength ( const double  lambda )

inline

Definition at line 26 of file JMakePDF.cc.

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

double getScatteringLength ( const double  lambda )

inline

Definition at line 32 of file JMakePDF.cc.

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

int main	(	int	argc,
		char **	argv
	)

Definition at line 52 of file JMakePDF.cc.

{
  using namespace std;
  using namespace JPP;

  string              outputFile;
  int                 numberOfPoints;
  double              epsilon;
  int                 function;
  set<double>         R;                          // 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 muon.");

    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_MUON,
      DIRECT_LIGHT_FROM_EMSHOWERS,
      DIRECT_LIGHT_FROM_DELTARAYS,
      SCATTERED_LIGHT_FROM_MUON,
      SCATTERED_LIGHT_FROM_EMSHOWERS,
      SCATTERED_LIGHT_FROM_DELTARAYS;
    zap['R'] = make_field(R,                      "distance of approach [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;


  // 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,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  typedef JPDFTransformer<3, JFunction1D_t::argument_type>        JFunction3DTransformer_t;
  typedef JArray<3, JFunction1D_t::argument_type>                 JArray_t;

  JPDF_t pdf;


  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
  

  const double kmin = pdf_c.getKappa(wmax);
  const double kmax = pdf_c.getKappa(wmin);
  const double cmin = pdf_c.getKmin (wmax);

  map<int, pair<fcn,JFunction3DTransformer_t> > zmap;

  zmap[DIRECT_LIGHT_FROM_MUON]         = make_pair((fcn) &JPDF::getDirectLightFromMuon,         JFunction3DTransformer_t(21.5, 2, kmin, kmax, NAMESPACE::getAngularAcceptance, 0.001));
  zmap[SCATTERED_LIGHT_FROM_MUON]      = make_pair((fcn) &JPDF::getScatteredLightFromMuon,      JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  zmap[DIRECT_LIGHT_FROM_EMSHOWERS]    = make_pair((fcn) &JPDF::getDirectLightFromEMshowers,    JFunction3DTransformer_t(21.5, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  zmap[SCATTERED_LIGHT_FROM_EMSHOWERS] = make_pair((fcn) &JPDF::getScatteredLightFromEMshowers, JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  zmap[DIRECT_LIGHT_FROM_DELTARAYS]    = make_pair((fcn) &JPDF::getDirectLightFromDeltaRays,    JFunction3DTransformer_t(21.5, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  zmap[SCATTERED_LIGHT_FROM_DELTARAYS] = make_pair((fcn) &JPDF::getScatteredLightFromDeltaRays, JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));

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

  fcn                      f           = zmap[function].first;    // PDF
  JFunction3DTransformer_t transformer = zmap[function].second;   // transformer


  if (R.empty()) {
    R.insert(  0.10);
    R.insert(  0.30);
    R.insert(  0.50);
    R.insert(  1.00);
    R.insert(  2.00);
    R.insert(  3.00);
    R.insert(  4.00);
    R.insert(  5.00);
    R.insert(  6.00);
    R.insert(  7.00);
    R.insert(  8.00);
    R.insert(  9.00);
    R.insert( 10.00);
    R.insert( 11.00);
    R.insert( 12.00);
    R.insert( 13.00);
    R.insert( 14.00);
    R.insert( 15.00);
    R.insert( 16.00);
    R.insert( 17.00);
    R.insert( 18.00);
    R.insert( 19.00);
    R.insert( 20.00);
    R.insert( 22.00);
    R.insert( 24.00);
    R.insert( 26.00);
    R.insert( 28.00);
    R.insert( 30.00);
    R.insert( 32.00);
    R.insert( 34.00);
    R.insert( 36.00);
    R.insert( 38.00);
    R.insert( 40.00);
    R.insert( 42.00);
    R.insert( 44.00);
    R.insert( 46.00);
    R.insert( 48.00);
    R.insert( 50.00);
    R.insert( 55.00);
    R.insert( 60.00);
    R.insert( 65.00);
    R.insert( 70.00);
    R.insert( 75.00);
    R.insert( 80.00);
    R.insert( 85.00);
    R.insert( 90.00);
    R.insert( 95.00);
    R.insert(100.00);
    R.insert(110.00);
    R.insert(120.00);
    R.insert(130.00);
    R.insert(140.00);
    R.insert(150.00);
    R.insert(170.00);
    R.insert(190.00);
    R.insert(210.00);
    R.insert(250.00);
  }

  set<double> X;

  if (function == DIRECT_LIGHT_FROM_MUON) {
 
    for (double buffer[] = { -0.01, -0.005, 0.0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, -1.0 }, *x = buffer; *x != -1.0; ++x) {
      X.insert(0.0 + *x);
      X.insert(1.0 - *x);
    }

    for (double x = 0.01; x < 0.1; x += 0.0025) {
      X.insert(0.0 + x);
      X.insert(1.0 - x);
    }

    for (double x = 0.10; x < 0.5; x += 0.010) {
      X.insert(0.0 + x);
      X.insert(1.0 - x);
    }
     
  } else {

    X.insert(  0.00);
    X.insert(  0.01);
    X.insert(  0.02);
    X.insert(  0.03);
    X.insert(  0.04);
    X.insert(  0.05);
    X.insert(  0.06);
    X.insert(  0.07);
    X.insert(  0.08);
    X.insert(  0.09);
    X.insert(  0.10);
    X.insert(  0.12);
    X.insert(  0.15);
    X.insert(  0.20);
    X.insert(  0.25);
    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.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( 12.0);
    X.insert( 14.0);
    X.insert( 16.0);
    X.insert( 18.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 r = R.begin(); r != R.end(); ++r) {

    const double R_m = *r;

    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[R_m][theta][phi];
        
        const JArray_t array(R_m, 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)(R_m, theta, phi, t);

          if (y != 0.0) {

            if (*x < 0.0) {
              WARNING("dt < 0 " << *x << ' ' << R_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 22 of file JMakePDF.cc.

double scatteringLengthFactor

Definition at line 23 of file JMakePDF.cc.