JMakePDE.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <cmath>
 
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiMapTransformer.hh"
#include "JTools/JQuadrature.hh"
 
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JGeanx.hh"
#include "JPhysics/JGeanz.hh"
 
#include "JSystem/JSystemToolkit.hh"
 
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JConstants.hh" 
 
/**
 * \file
 *
 * 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 <tt>(E, D, cos(theta), theta, phi, t)</tt>, where:
 * - <tt>E</tt> is the logarithm of the energy;
 * - <tt>D</tt> is the distance between the vertex and the PMT;
 * - <tt>cos(theta)</tt> the cosine of the photon emission angle;
 * - <tt>(theta, phi)</tt> the orientation of the PMT; and
 * - <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 shower developes along the z-axis and the PMT is located in the x-z plane.
 * \author jseneca, mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string              fileDescriptor;
  string              outputFile;
  int                 numberOfPoints;
  size_t              elong_sample_count;
  size_t              cosine_angle_count; // Cosine angle is not as sharp when elongated
  double              epsilon;
  double              TTS;
  set<double>         Ds;                          // distance [m]
  set<double>         Es;                         // logarithm of energy (GeV)
  bool                option;
  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['F'] = make_field(fileDescriptor,         "file name descriptor for PDF tables");
    zap['o'] = make_field(outputFile,             "file name output PDF");
    zap['n'] = make_field(numberOfPoints,         "points for integration")     = 25;
    zap['e'] = make_field(epsilon,                "precision for integration")  = 1.0e-10;
    zap['D'] = make_field(Ds,                     "distance [m]")               = JPARSER::initialised();
    zap['E'] = make_field(Es,                     "energy [log10(GeV)]")        = JPARSER::initialised();
    zap['C'] = make_field(cosine_angle_count,     "points for cosine angle")    = 60;
    zap['N'] = make_field(elong_sample_count,     "points for elongation")      = 20;
    zap['T'] = make_field(TTS,                    "sigma of time bluring [ns]") = 2.0;
    zap['O'] = make_field(option,                 "optional transformation");
    zap['d'] = make_field(debug)                  = 0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  if (cosine_angle_count == 0) {
    FATAL("Invalid option -C <" << cosine_angle_count << ">" << endl);
  }
 
  typedef JPolint1Function1D_t                                    JFunction1D_t;
  
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            J4DMapList_t;
 
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalMap,
                   JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            J5DMapList_t;
 
  typedef JPDFTable<JFunction1D_t, J4DMapList_t>                  JPDF4d_t;
  typedef JPDFTable<JFunction1D_t, J5DMapList_t>                  JPDF5d_t;
 
  typedef JPDFTransformer<4, JFunction1D_t::argument_type>        JFunction4DTransformer_t;
  typedef JPDFTransformer<5, JFunction1D_t::argument_type>        JFunction5DTransformer_t;
  typedef JArray<5, JFunction1D_t::argument_type>                 JArray_t;
 
  JPDF4d_t i_pdf;
 
  DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
  {
    const string file_name = getFilename(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER);
 
    NOTICE("loading input from file " << file_name << "... " << flush);
 
    i_pdf.load(file_name.c_str());
 
    NOTICE("OK" << endl);
  }
  DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
  {
    JPDF4d_t pdf;
 
    const string file_name = getFilename(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER);
 
    NOTICE("loading input from file " << file_name << "... " << flush);
 
    pdf.load(file_name.c_str());
 
    NOTICE("OK" << endl);
    DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
 
    pdf.setExceptionHandler(new JPDF4d_t::JDefaultResult(JMATH::zero));
 
    NOTICE("Add PDF... " << flush);
 
    JTimer timer;
 
    timer.start();
 
    i_pdf.add(pdf);
 
    timer.stop();
 
    if (debug >= debug_t) {
      timer.print(cout, unit_t);
    }
 
    NOTICE("OK" << endl);
  }
  DEBUG("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
 
  i_pdf.setExceptionHandler(new JPDF4d_t::JDefaultResult(JMATH::zero));
 
  if ( TTS > 0.0 ){
 
    NOTICE("Smearing combined table with sigma [ns] = " << TTS << " ..." << flush);
 
    i_pdf.blur( TTS );
 
    NOTICE("OK" << endl);
  }
 
  JPDF5d_t o_pdf;
  
  const JFunction4DTransformer_t* i_transformer = dynamic_cast<JFunction4DTransformer_t*>(i_pdf.transformer->clone());
  const JFunction5DTransformer_t* o_transformer = (i_transformer != NULL ? new JFunction5DTransformer_t(*i_transformer) : NULL);
 
  if (o_transformer != NULL && debug >= debug_t) {
    o_transformer->print(cout);
  }
 
  if (Es.empty()) {
    Es.insert( 1.0 );
    Es.insert( 2.0 );
    Es.insert( 3.0 );
    Es.insert( 3.5 );
    Es.insert( 4.0 );
    Es.insert( 4.5 );
    Es.insert( 5.0 );
    Es.insert( 5.5 );
    Es.insert( 6.0 );
    Es.insert( 6.5 );
    Es.insert( 7.0 );
    Es.insert( 7.5 );
    Es.insert( 8.0 );
  };
 
  if (Ds.empty()) {
    Ds.insert(  0.10);
    Ds.insert(  0.50);
    Ds.insert(  1.00);
    Ds.insert(  2.00);
    Ds.insert(  5.00);
    Ds.insert(  8.00);
    Ds.insert( 10.00);
    Ds.insert( 20.00);
    Ds.insert( 30.00);
    Ds.insert( 40.00);
    Ds.insert( 50.00);
    Ds.insert( 60.00);
    Ds.insert( 70.00);
    Ds.insert( 80.00);
    Ds.insert( 90.00);
    Ds.insert(100.00);
    Ds.insert(120.00);
    Ds.insert(150.00);
    Ds.insert(170.00);
    Ds.insert(190.00);
    Ds.insert(210.00);
    Ds.insert(230.00);
    Ds.insert(250.00);
    Ds.insert(270.00);
    Ds.insert(290.00);
    Ds.insert(310.00);
    Ds.insert(400.00);
    Ds.insert(800.00);   // ARCA is big, we go to 800 m
  }
 
  set<double> C;        // cosine emission angle
 
  JQuadrature qeant(-1.0, +1.0, cosine_angle_count, 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]
  
  X.insert(  0.00);
  X.insert(  0.50);
  X.insert(  1.00);
  X.insert(  1.50);
  X.insert(  2.00);
  X.insert(  2.50);
  X.insert(  3.00);
  X.insert(  3.50);
  X.insert(  4.00);
  X.insert(  5.0);
  X.insert(  6.0);
  X.insert(  8.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( 80.0);
  X.insert(100.0);
  X.insert(120.0);
  X.insert(150.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  = 7.0;                                       // [deg]  
  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0));  // azimuth angle unit step size
  
  NOTICE("building multi-dimensional function object" << endl);
      
  for (const double E : Es) {
 
    DEBUG("E: " << E << endl);
        
    vector< double > elong_lengths;
 
    for (size_t i = 0; i != elong_sample_count ; ++i ){
 
      const double fraction = (double) (i + 0.5) / (double) elong_sample_count;
      const double len      = geanz.getLength(pow(10.0, E), fraction, 1e-6);
 
      DEBUG("fraction: " << fraction << ", length [m]: "<< len << endl);
 
      elong_lengths.push_back(len);
    }
 
    if (elong_lengths.empty()) {
      elong_lengths.push_back(0.0);
    }
 
      
    for (const double D_m : Ds) {
 
      DEBUG("D_m: " << D_m << endl);
      
      for (const double cd : C) {
 
        NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
        
        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 = o_pdf[E][D_m][cd][theta][phi];
 
            f1.reserve(X.size());
            
            const JArray_t array(E, D_m, cd, theta, phi);
 
            double t_old = (o_transformer != NULL ? o_transformer->getXn(array, *X.begin()) : *X.begin());
            double y_old = 0.0;
 
            for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
          
              const double t = (o_transformer != NULL ? o_transformer->getXn(array, *x) : *x);
 
              double y = 0.0;
              
              for (const double Z : elong_lengths) {
        
                const double Z0 = cd * D_m;
                const double Z1 = Z0 - Z;
                
                const double D_m_prime = sqrt(D_m * D_m + Z1 * Z1 - Z0 * Z0);
                const double cd_prime  = Z1 / D_m_prime;
                const double t_prime   = t - (Z) / getSpeedOfLight() - (D_m_prime - D_m) * getIndexOfRefraction() / getSpeedOfLight();
                
                y += i_pdf(D_m_prime, cd_prime, theta, phi, t_prime);
              }
 
              y /= (double) elong_lengths.size();
              
              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;
            }
          }
        }
      }
    }
  }
  NOTICE("OK" << endl);
  NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
 
  o_pdf.compile();
 
  if (o_transformer != NULL && option) {
 
    NOTICE("transform... " << flush);
 
    o_pdf.transform(*o_transformer);
 
    NOTICE("OK" << endl);
  }
 
  try {
 
    NOTICE("storing output to file " << outputFile << "... " << flush);
 
    o_pdf.store(outputFile.c_str());
 
    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
}