JPDF.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JPDF__
#define __JPHYSICS__JPDF__

#include <vector>
#include <algorithm>
#include <cmath>
#include <limits>

#include "JLang/JCC.hh"
#include "JLang/JException.hh"
#include "JTools/JConstants.hh"
#include "JTools/JQuadrature.hh"
#include "JPhysics/JDispersionInterface.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JAbstractPMT.hh"
#include "JPhysics/JAbstractMedium.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JGeometry3D/JVector3D.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  using JTOOLS::PI;
  using JLANG::JException;

  /**
   * Radius of optical module [m].
   *
   * This parameter is used to implement shadowing of the PMT by the optical module.
   */
  static double MODULE_RADIUS_M = 0.25;


  /**
   * Low level interface for the calculation of the Probability Density Functions (PDFs) of
   * the arrival time of Cherenkov light from a muon or an EM-shower on a photo-multiplier tube (PMT).
   * The calculation of the PDF covers direct light light and single-scattered light.
   * The direct light is attenuated by absorption and any form of light scattering.\n
   * The attenuation length is defined as:
   \f[
   \frac{1}{\lambda_{att}} ~\equiv~ \frac{1}{\lambda_{abs}} + \frac{1}{\lambda_{s}} 
   \f] 
   * where
   \f$\lambda_{abs}\f$ and 
   \f$\lambda_{s}\f$ 
   refer to the absorption and scattering length, respectively.\n
   * The scattered light is attenuated by absorption and light scattering weighed with the scattering angle \f$\theta_{s}\f$.
   * The attenuation length is in this case defined as:
   \f[
   \frac{1}{\lambda_{att}} ~\equiv~ \frac{1}{\lambda_{abs}} + \frac{1-\left<\cos\theta_s\right>}{\lambda_{s}} 
   \f] 
   *
   * The PDFs of a muon are evaluated as a function of:
   *   - distance of closest approach of the muon to the PMT, and
   *   - orientation of the PMT (zenith and azimuth angle).
   *
   * The PDFs of an EM-shower are evaluated as a function of:
   *   - distance between the EM-shower and the PMT, 
   *   - cosine angle EM-shower direction and EM-shower - PMT position, and
   *   - orientation of the PMT (zenith and azimuth angle).
   *
   * The intensity of light from an EM-shower is assumed to be proportional to the energy of the shower.\n
   * For a given energy of the EM-shower, the longitudinal profile of the EM-shower is taken into account. 
   *
   * The coordinate system is defined such that z-axis corresponds to the muon trajectory (EM-shower).\n
   * The orientation of the PMT (i.e. zenith and azimuth angle) are defined following a rotation around the z-axis,
   * such that the PMT is located in the x-z plane.
   *
   * Schematics of the PDF of direct light from a muon:
   *
   \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)

   \put( 1.0, 0.5){\vector(0,1){9}}
   \put( 1.0,10.0){\makebox(0,0){$z$}}
   \put( 1.0, 0.0){\makebox(0,0){muon}}

   \put( 1.0, 8.0){\line(1,0){6}}
   \put( 4.0, 8.5){\makebox(0,0)[c]{$R$}}

   \multiput( 1.0, 2.0)(0.5, 0.5){12}{\qbezier(0.0,0.0)(-0.1,0.35)(0.25,0.25)\qbezier(0.25,0.25)(0.6,0.15)(0.5,0.5)}
   \put( 4.5, 4.5){\makebox(0,0)[l]{photon}}

   \put( 1.0, 2.0){\circle*{0.2}}
   \put( 0.5, 2.0){\makebox(0,0)[r]{$(0,0,z,t_{0})$}}

   \put( 1.0, 8.0){\circle*{0.2}}
   \put( 0.5, 8.0){\makebox(0,0)[r]{$(0,0,0)$}}

   \put( 7.0, 8.0){\circle*{0.2}}
   \put( 7.0, 9.0){\makebox(0,0)[c]{PMT}}
   \put( 7.5, 8.0){\makebox(0,0)[l]{$(R,0,0,t_{1})$}}

   \put( 1.1, 2.1){
   \put(0.0,1.00){\vector(-1,0){0.1}}
   \qbezier(0.0,1.0)(0.25,1.0)(0.5,0.75)
   \put(0.5,0.75){\vector(1,-1){0.1}}
   \put(0.4,1.5){\makebox(0,0){$\theta_{c}$}}
   }

   \end{picture}
   \f}
   *
   * For the PDF of a muon, the time is defined such that \f$t ~=~ 0\f$
   * refers to the time when the muon passes through the point at minimal approach to the PMT 
   * (i.e. position \f$(0,0,0)\f$ in the above figure).
   *
   * The delay for the Cherenkov cone to actually pass through the PMT is defined as 
   \f[
   \Delta t ~\equiv~ R \tan(\theta_{c}) / c
   \f] 
   * where \f$\tan(\theta_{c})\f$ is given by method JTOOLS::getTanThetaC.
   *
   * Given a measured hit time, \f$t_{hit}\f$, the PDF should be probed at \f$(t_{hit} - t_1)\f$, 
   * where \f$t_1\f$ is given by:
   *
   \f[
   ct_1 =  ct_0 - z + R \tan(\theta_{c}) 
   \f]
   *
   * In this, \f$t_0\f$ refers to the time the muon passes through point \f$(0,0,z)\f$.
   *
   * Schematics of single scattered light from a muon:
   *
   \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)

   \put( 1.0, 0.5){\vector(0,1){9}}
   \put( 1.0,10.0){\makebox(0,0){$z$}}
   \put( 1.0, 0.0){\makebox(0,0){muon}}

   \put( 1.0, 8.0){\line(1,0){2}}
   \put( 2.0, 8.5){\makebox(0,0)[c]{$R$}}

   \multiput( 1.0, 2.0)( 0.5, 0.5){8}{\qbezier(0.0,0.0)(-0.1,0.35)( 0.25,0.25)\qbezier( 0.25,0.25)( 0.6,0.15)( 0.5,0.5)}
   \multiput( 5.0, 6.0)(-0.5, 0.5){4}{\qbezier(0.0,0.0)(+0.1,0.35)(-0.25,0.25)\qbezier(-0.25,0.25)(-0.6,0.15)(-0.5,0.5)}

   \put( 5.0, 6.0){\circle*{0.2}}

   \put( 3.5, 3.5){\makebox(0,0)[c]{$\bar{u}$}}
   \put( 5.0, 7.0){\makebox(0,0)[c]{$\bar{v}$}}

   \put( 1.0, 2.0){\circle*{0.2}}
   \put( 0.5, 2.0){\makebox(0,0)[r]{$(0,0,z,t_0)$}}

   \put( 1.0, 8.0){\circle*{0.2}}
   \put( 0.5, 8.0){\makebox(0,0)[r]{$(0,0,0)$}}

   \put( 3.0, 8.0){\circle*{0.2}}
   \put( 3.0, 9.0){\makebox(0,0)[c]{PMT}}
   \put( 3.5, 8.0){\makebox(0,0)[l]{$(R,0,0,t_1)$}}

   \end{picture}
   \f}
   *
   * In this,
   *
   \f{math}{
   \bar{u} ~\equiv~ u ~ \hat{u} 
   ~\equiv~  u \left(\begin{array}{c} \sin(\theta_{0}) \cos(\phi_{0}) \\ \sin(\theta_{0}) \sin(\phi_{0}) \\ \cos(\theta_{0})\end{array}\right) 
   \textrm{     and     }
   \bar{v} ~\equiv~  v ~ \hat{v} 
   ~\equiv~  v \left(\begin{array}{c} \sin(\theta_{1}) \cos(\phi_{1}) \\ \sin(\theta_{1}) \sin(\phi_{1}) \\ \cos(\theta_{1})\end{array}\right) 
   \f}
   *
   * The vectors \f$\bar{u}\f$ and \f$\bar{v}\f$ are subject to the following constraint:
   \f{math}{
   \begin{array}{ccccccc}

   \begin{array}{c}
   \mathrm{start}\\
   \mathrm{point}
   \end{array}

   &&

   \begin{array}{c}
   \mathrm{before}\\
   \mathrm{scattering}
   \end{array}

   &&

   \begin{array}{c}
   \mathrm{after}\\
   \mathrm{scattering}
   \end{array}

   &&

   \mathrm{PMT}

   \\
   \\

   \left(\begin{array}{c} 0 \\ 0 \\ z \end{array}\right) 
   & + &  
   u \left(\begin{array}{c} \sin(\theta_{0}) \cos(\phi_{0}) \\ \sin(\theta_{0}) \sin(\phi_{0}) \\ \cos(\theta_{0})\end{array}\right) 
   & + &
   v \left(\begin{array}{c} \sin(\theta_{1}) \cos(\phi_{1}) \\ \sin(\theta_{1}) \sin(\phi_{1}) \\ \cos(\theta_{1})\end{array}\right) 
   & = &
   \left(\begin{array}{c} R \\ 0 \\ 0 \end{array}\right)
   \end{array}
   * \f}
   *
   * Note that an expression for the cosine of the scattering angle, \f$\cos(\theta_s) ~\equiv~ \hat{u} \cdot \hat{v}\f$, 
   * can directly be obtained by multiplying the left hand side and the right hand side with \f$\hat{u}\f$:
   *
   \f[
   \cos\theta_s = \frac{R\sin\theta_0\cos\phi_0 - z\cos\theta_0 - u}{v}
   \f]
   *
   *
   * The arrival time of the single scattered light can be expressed as:
   *
   \f[
   ct_1  =  ct_0 + n (u + v)
   \f]
   *
   * where \f$n\f$ is the index of refraction.
   *
   * Schematics of direct light from a EM-shower:
   *
   \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)

   \put( 1.0, 2.0){\vector(0,1){1}}
   \put( 1.0, 3.5){\multiput(0,0)(0,0.5){12}{\line(0,1){0.2}}}
   \put( 1.0,10.0){\makebox(0,0){$z$}}
   \put( 1.0, 1.0){\makebox(0,0){EM-shower}}

   \put( 1.0, 8.0){\line(1,0){6}}
   \put( 4.0, 8.5){\makebox(0,0)[c]{$R$}}

   \multiput( 1.0, 2.0)(0.5, 0.5){12}{\qbezier(0.0,0.0)(-0.1,0.35)(0.25,0.25)\qbezier(0.25,0.25)(0.6,0.15)(0.5,0.5)}
   \put( 4.5, 4.5){\makebox(0,0)[l]{photon}}

   \put( 1.0, 2.0){\circle*{0.2}}
   \put( 0.5, 2.0){\makebox(0,0)[r]{$(0,0,z,t_0)$}}

   \put( 1.0, 8.0){\circle*{0.2}}
   \put( 0.5, 8.0){\makebox(0,0)[r]{$(0,0,0)$}}

   \put( 7.0, 8.0){\circle*{0.2}}
   \put( 7.0, 9.0){\makebox(0,0)[c]{PMT}}
   \put( 7.5, 8.0){\makebox(0,0)[l]{$(R,0,0,t_1)$}}

   \put( 1.1, 2.1){
   \put(0.0,1.00){\vector(-1,0){0.1}}
   \qbezier(0.0,1.0)(0.25,1.0)(0.5,0.75)
   \put(0.5,0.75){\vector(1,-1){0.1}}
   \put(0.4,1.5){\makebox(0,0){$\theta_0$}}
   }

   \end{picture}
   \f}
   *
   * Note that the emission angle \f$\theta_0\f$ is in general not equal to the Cherenkov angle \f$\theta_c\f$.\n
   * For an EM-shower, the arrival time of the light can be expressed as:
   *
   \f[
   ct_1  =  ct_0 + n D
   \f]
   *
   * where \f$D\f$ is the total distance traveled by the photon from its emission point to the PMT,
   * including possible scattering.
   *
   * For the computation of the various integrals, it is assumed that the index of refraction
   * decreases monotonously as a function of the wavelength of the light.
   */
  class JPDF : 
    public JTOOLS::JGaussLegendre, 
    public virtual JDispersionInterface, 
    public virtual JAbstractPMT,
    public virtual JAbstractMedium
  {

    typedef JTOOLS::JElement2D<double, double>               element_type;

  public:
    /**
     * Constructor.
     *
     * \param  Wmin            minimal wavelength for integration [nm]
     * \param  Wmax            maximal wavelength for integration [nm]
     * \param  numberOfPoints  number    of points for integration
     * \param  epsilon         precision of points for integration
     */
    JPDF(const double Wmin,
         const double Wmax,
         const int    numberOfPoints = 20,
         const double epsilon        = 1e-12) :
      JTOOLS::JGaussLegendre(numberOfPoints,epsilon),
      wmin(Wmin),
      wmax(Wmax)
    {
      // phi integration points
      
      const double dx = PI / numberOfPoints;

      for (double x = 0.5*dx; x < PI; x += dx) {
        phd.push_back(element_type(cos(x),sin(x)));
      }
    }


    /**
     * Virtual destructor.
     */
    virtual ~JPDF()
    {}

   
    /**
     * Number of Cherenkov photons per unit track length
     *
     * \return            number of photons per unit track length [m^-1]
     */
    double getNumberOfPhotons() const
    {
      double value = 0.0;

      const double xmin = 1.0 / wmax;
      const double xmax = 1.0 / wmin;

      for (const_iterator i = begin(); i != end(); ++i) {

        const double x  = 0.5 * (xmax + xmin)  +  i->getX() * 0.5 * (xmax - xmin);
        const double dx = i->getY() * 0.5 * (xmax - xmin);

        const double w  = 1.0 / x;
        const double dw = dx * w*w;

        const double n  = getIndexOfRefractionPhase(w);
        
        value += cherenkov(w,n) * dw;
      }

      return value;
    }


    /**
     * Number of photo-electrons from direct Cherenkov light from muon.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \return            P [npe]
     */
    double getDirectLightFromMuon(const double R_m,
                                  const double theta,
                                  const double phi) const
    {
      using namespace JTOOLS;

      double value = 0;

      const double R  =  std::max(R_m, getRmin());
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double pz =  cos(theta);

      for (const_iterator m = begin(); m != end(); ++m) {

        const double w  = 0.5 * (wmax + wmin)  +  m->getX() * 0.5 * (wmax - wmin);
        const double dw = m->getY() * 0.5 * (wmax - wmin);

        const double n     = getIndexOfRefractionPhase(w);

        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);

        const double npe   = cherenkov(w,n) * dw * getQE(w);
        
        const double ct0 = 1.0 / n;
        const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
        
        const double d  = R / st0;                                 // distance traveled by photon
        const double ct = st0*px + ct0*pz;                         // cosine angle of incidence on PMT
        
        const double U  = getAngularAcceptance(ct);                // PMT angular acceptance
        const double V  = exp(-d/l_abs) * exp(-d/ls);              // absorption & scattering
        const double W  = A / (2.0*PI*R*st0);                      // solid angle
        
        value += npe * U * V * W;
      }

      return value;
    }

 
    /**
     * Probability density function for direct light from muon.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getDirectLightFromMuon(const double R_m,
                                  const double theta,
                                  const double phi,
                                  const double t_ns) const
    {
      using namespace JTOOLS;

      static const int    N   = 100;                               // maximal number of iterations
      static const double eps = 1.0e-6;                            // precision index of refraction

      const double R  =  std::max(R_m, getRmin());
      const double t  =  R * getTanThetaC() / C  +  t_ns;          // time [ns]
      const double a  =  C*t/R;                                    // target value
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double pz =  cos(theta);

      // check validity range for index of refraction

      for (double buffer[] = { wmin, wmax, 0.0 }, *p = buffer; *p != 0.0; ++p) {

        const double n  = getIndexOfRefractionPhase(*p);
        const double ng = getIndexOfRefractionGroup(*p);    
        
        const double ct0 = 1.0 / n;
        const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
        
        const double b  = (ng - ct0) / st0;                        // running value

        if (*p == wmin && b < a) { return 0.0; }
        if (*p == wmax && b > a) { return 0.0; }
      }


      double umin = wmin;
      double umax = wmax;

      for (int i = 0; i != N; ++i) {                               // binary search for wavelength
        
        const double w  = 0.5 * (umin + umax);
        
        const double n  = getIndexOfRefractionPhase(w);
        const double ng = getIndexOfRefractionGroup(w);    
        
        const double ct0 = 1.0 / n;
        const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
        
        const double b  = (ng - ct0) / st0;                        // running value

        if (fabs(b-a) < a*eps) {

          const double np    = getDispersionPhase(w);
          const double ngp   = getDispersionGroup(w);

          const double l_abs = getAbsorptionLength(w);
          const double ls    = getScatteringLength(w);
      
          const double d  = R / st0;                               // distance traveled by photon
          const double ct = st0*px + ct0*pz;                       // cosine angle of incidence on PMT
          
          const double i3 = ct0*ct0*ct0 / (st0*st0*st0);

          const double U  = getAngularAcceptance(ct);              // PMT angular acceptance
          const double V  = exp(-d/l_abs - d/ls);                  // absorption & scattering
          const double W  = A / (2.0*PI*R*st0);                    // solid angle
          
          const double Ja = R * (ngp/st0 + np*(n-ng)*i3) / C;      // dt/dlambda

          return cherenkov(w,n) * getQE(w) * U * V * W / fabs(Ja);
        }       

        if (b > a)
          umin = w;
        else
          umax = w;
      }
        
      return 0.0;
    }


    /**
     * Probability density function for scattered light from muon.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getScatteredLightFromMuon(const double R_m,
                                     const double theta,
                                     const double phi,
                                     const double t_ns) const
    {
      using namespace JTOOLS;

      static const double eps = 1.0e-10;

      double value = 0;

      const double R  =  std::max(R_m, getRmin());
      const double t  =  R * getTanThetaC() / C  +  t_ns;          // time [ns]
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);
      const double ni = sqrt(R*R + C*t*C*t) / R;                   // maximal index of refraction
    
      if (n0 >= ni) return value;

      const double nj = std::min(ni,n1);
      
      double w = wmax;

      for (const_iterator i = begin(); i != end(); ++i) {
        
        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
        
        w  = getWavelength(ng, w, 1.0e-5);

        const double dw = dn / fabs(getDispersionGroup(w));
        
        const double n  = getIndexOfRefractionPhase(w);   
        
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);
          
        const double npe   = cherenkov(w,n) * dw * getQE(w);

        if (npe <= 0) { continue; }

        const double Jc  = 1.0 / ls;                               // dN/dx

        const double ct0 = 1.0 / n;                                // photon direction before scattering
        const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));

        JRoot rz(R, ng, t);                                        // square root

        if (!rz.is_valid) { continue; }

        const double zmin = rz.first;
        const double zmax = rz.second;

        const double zap  = 1.0 / l_abs;
        
        const double xmin = exp(zap*zmax);
        const double xmax = exp(zap*zmin);
        
        for (const_iterator j = begin(); j != end(); ++j) {
          
          const double x  = 0.5 * (xmax + xmin)  +  j->getX() * 0.5 * (xmax - xmin);
          const double dx = j->getY() * 0.5 * (xmax - xmin);
          
          const double z  = log(x) / zap;
          const double dz = -dx / (zap*x);
        
          const double D  = sqrt(z*z + R*R);
          const double cd = -z / D;
          const double sd =  R / D;
                
          const double d  = (C*t - z) / ng;                        // photon path

          //const double V  = exp(-d/l_abs);                         // absorption

          const double cta  = cd*ct0 + sd*st0;
          const double dca  = d - 0.5*(d+D)*(d-D) / (d - D*cta);
          const double tip  = -log(D*D/(dca*dca) + eps) / PI;
          
          const double ymin = exp(tip*PI);
          const double ymax = 1.0;
          
          for (const_iterator k = begin(); k != end(); ++k) {
            
            const double y  = 0.5 * (ymax + ymin)  +  k->getX() * 0.5 * (ymax - ymin);
            const double dy = k->getY() * 0.5 * (ymax - ymin);
            
            const double phi = log(y) / tip;
            const double dp  = -dy / (tip*y);
            
            const double cp0 = cos(phi);
            const double sp0 = sin(phi);
            
            const double u  = (R*R + z*z - d*d) / (2*R*st0*cp0 - 2*z*ct0 - 2*d);
            const double v  = d - u;

            if (u <= 0) { continue; }
            if (v <= 0) { continue; }
                  
            const double vi  = 1.0 / v;
            const double cts = (R*st0*cp0 - z*ct0 - u)*vi;         // cosine scattering angle

            const double V  = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));

            if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
            
            const double vx  =   R - u*st0*cp0;
            const double vy  =     - u*st0*sp0;
            const double vz  =  -z - u*ct0;
              
            const double ct[]  = {                                 // cosine angle of incidence on PMT
              (vx*px + vy*py + vz*pz) * vi,
              (vx*px - vy*py + vz*pz) * vi
            };

            const double U  = 
              getAngularAcceptance(ct[0]) + 
              getAngularAcceptance(ct[1]);                         // PMT angular acceptance
            const double W  = std::min(A*vi*vi, 2.0*PI);           // solid angle

            const double Ja = getScatteringProbability(cts);       // d^2P/dcos/dphi
            const double Jd = ng * (1.0 - cts) / C;                // dt/du

            value += (npe * dz * dp / (2*PI)) * U * V * W * Ja * Jc / fabs(Jd);
          }
        }
      }

      return value;
    }


    /**
     * Probability density function for direct light from EM-showers.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            dP/dt [npe/ns/GeV]
     */
    double getDirectLightFromEMshowers(const double R_m,
                                       const double theta,
                                       const double phi,
                                       const double t_ns) const
    {
      using namespace JTOOLS;

      double value = 0;

      const double R  =  std::max(R_m, getRmin());
      const double t  =  R * getTanThetaC() / C  +  t_ns;          // time [ns]
      const double A  =  getPhotocathodeArea();
      const double D  =  2.0*sqrt(A/PI);

      const double px =  sin(theta)*cos(phi);
      //const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);
      const double ni = sqrt(R*R + C*t*C*t) / R;                   // maximal index of refraction
    
      if (n0 >= ni) return value;

      const double nj = std::min(n1, ni);

      double w = wmax;

      for (const_iterator i = begin(); i != end(); ++i) {
            
        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
        
        w  = getWavelength(ng, w, 1.0e-5);

        const double dw = dn / fabs(getDispersionGroup(w));

        const double n  = getIndexOfRefractionPhase(w);   
          
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);

        const double npe   = cherenkov(w,n) * dw * getQE(w);

        JRoot rz(R, ng, t);                                        // square root

        if (!rz.is_valid) { continue; }

        for (int j = 0; j != 2; ++j) {

          const double z = rz[j];
          
          if (C*t <= z) continue;
          
          const double d  = sqrt(z*z + R*R);                       // distance traveled by photon

          const double ct0 = -z / d;
          const double st0 =  R / d;
          
          const double dz = D / st0;                               // average track length
          
          const double ct = st0*px + ct0*pz;                       // cosine angle of incidence on PMT
          
          const double U  = getAngularAcceptance(ct);              // PMT angular acceptance
          const double V  = exp(-d/l_abs - d/ls);                  // absorption & scattering
          const double W  = A/(d*d);                               // solid angle
          
          const double Ja = geant(n,ct0);                          // d^2N/dcos/dphi
          const double Jd = (1.0 - ng*ct0) / C;                    // d  t/ dz
          const double Je = ng*st0*st0*st0 / (R*C);                // d^2t/(dz)^2
          
          value += gWater() * npe * U * V * W * Ja / (fabs(Jd) + 0.5*Je*dz);
        }
      }

      return value;
    }


    /**
     * Probability density function for scattered light from EM-showers.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            dP/dt [npe/ns/GeV]
     */
    double getScatteredLightFromEMshowers(const double R_m,
                                          const double theta,
                                          const double phi,
                                          const double t_ns) const
    {
      using namespace JTOOLS;

      double value = 0;

      const double R  =  std::max(R_m, getRmin());
      const double t  =  R * getTanThetaC() / C  +  t_ns;          // time [ns]
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);
      const double ni = sqrt(R*R + C*t*C*t) / R;                   // maximal index of refraction
    
      if (n0 >= ni) return value;

      const double nj = std::min(ni,n1);

      double w = wmax;

      for (const_iterator i = begin(); i != end(); ++i) {

        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
        
        w  = getWavelength(ng, w, 1.0e-5);

        const double dw = dn / fabs(getDispersionGroup(w));
        
        const double n  = getIndexOfRefractionPhase(w);   
        
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);
          
        const double npe   = cherenkov(w,n) * dw * getQE(w);

        if (npe <= 0) { continue; }

        const double Jc = 1.0 / ls;                                // dN/dx
        
        JRoot rz(R, ng, t);                                        // square root

        if (!rz.is_valid) { continue; }

        const double zmin = rz.first;
        const double zmax = rz.second;

        const double zap  = 1.0 / l_abs;
        
        const double xmin = exp(zap*zmax);
        const double xmax = exp(zap*zmin);
        
        for (const_iterator j = begin(); j != end(); ++j) {
          
          const double x  = 0.5 * (xmax + xmin)  +  j->getX() * 0.5 * (xmax - xmin);
          const double dx = j->getY() * 0.5 * (xmax - xmin);
          
          const double z  = log(x) / zap;
          const double dz = -dx / (zap*x);
          
          const double D  = sqrt(z*z + R*R);
          const double cd = -z / D;
          const double sd =  R / D;
          
          const double qx = cd*px +  0  - sd*pz;
          const double qy =         1*py;
          const double qz = sd*px +  0  + cd*pz;

          const double d = (C*t - z) / ng;                         // photon path
            
          //const double V  = exp(-d/l_abs);                         // absorption

          const double ds = 2.0 / (size() + 1);
              
          for (double sb = 0.5*ds; sb < 1.0 - 0.25*ds; sb += ds) {
                
            for (int buffer[] = { -1, +1, 0}, *k = buffer; *k != 0; ++k) {
                  
              const double cb  = (*k) * sqrt((1.0 + sb)*(1.0 - sb));
              const double dcb = (*k) * ds*sb/cb;
                  
              const double v  = 0.5 * (d + D) * (d - D) / (d - D*cb);
              const double u  = d - v;
                  
              if (u <= 0) { continue; }
              if (v <= 0) { continue; }
                    
              const double cts = (D*cb - v) / u;                   // cosine scattering angle
                    
              const double V  = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
              
              if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }

              const double W  = std::min(A/(v*v), 2.0*PI);         // solid angle
              const double Ja = getScatteringProbability(cts);     // d^2P/dcos/dphi
              const double Jd = ng * (1.0 - cts) / C;              // dt/du
                    
              const double ca = (D - v*cb) / u;
              const double sa =      v*sb  / u;

              const double dp  = PI / phd.size();
              const double dom = dcb*dp * v*v / (u*u);

              for (const_iterator l = phd.begin(); l != phd.end(); ++l) {

                const double cp  = l->getX();
                const double sp  = l->getY();

                const double ct0 = cd*ca - sd*sa*cp;
                    
                const double vx  = sb*cp * qx;
                const double vy  = sb*sp * qy;
                const double vz  = cb    * qz;

                const double U  = 
                  getAngularAcceptance(vx + vy + vz) +
                  getAngularAcceptance(vx - vy + vz);              // PMT angular acceptance

                const double Jb = geant(n,ct0);                    // d^2N/dcos/dphi

                value += dom * gWater() * npe * dz * U * V * W * Ja * Jb * Jc / fabs(Jd);
              }
            }
          }
        }
      }

      return value;
    }


    /**
     * Probability density function for scattered light from muon.
     *
     * \param  D_m        distance between track segment and PMT [m]
     * \param  cd         cosine angle muon direction and track segment - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dx [npe/ns/m]
     */
    double getScatteredLightFromMuon(const double D_m,
                                     const double cd,
                                     const double theta,
                                     const double phi,
                                     const double t_ns) const
    {
      using namespace JTOOLS;

      static const double eps = 1.0e-10;

      double value = 0;

      const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
      const double D  =  std::max(D_m, getRmin());
      const double R  =  sd * D;                                   // minimal distance of approach [m]
      const double Z  = -cd * D;                                   // photon emission point
      const double L  =  D;
      const double t  =  D * getIndexOfRefraction() / C  +  t_ns;  // time [ns]
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);

      const double ni = C*t / L;                                   // maximal index of refraction
    
      if (n0 >= ni) {
        return value;
      }

      const double nj = std::min(ni,n1);

      double w = wmax;

      for (const_iterator i = begin(); i != end(); ++i) {
        
        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
          
        w  = getWavelength(ng, w, 1.0e-5);
          
        const double dw = dn / fabs(getDispersionGroup(w));
        
        const double n  = getIndexOfRefractionPhase(w);   
        
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);

        const double npe   = cherenkov(w,n) * dw * getQE(w);

        if (npe <= 0) { continue; }

        const double Jc    = 1.0 / ls;                             // dN/dx

        const double ct0 = 1.0 / n;                                // photon direction before scattering
        const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));

        const double d =  C*t / ng;                                // photon path
        
        //const double V  = exp(-d/l_abs);                           // absorption

        const double cta  = cd*ct0 + sd*st0;
        const double dca  = d - 0.5*(d+L)*(d-L) / (d - L*cta);
        const double tip  = -log(L*L/(dca*dca) + eps) / PI;
        
        const double ymin = exp(tip*PI);
        const double ymax = 1.0;
          
        for (const_iterator j = begin(); j != end(); ++j) {
          
          const double y  = 0.5 * (ymax + ymin)  +  j->getX() * 0.5 * (ymax - ymin);
          const double dy = j->getY() * 0.5 * (ymax - ymin);
            
          const double phi = log(y) / tip;
          const double dp  = -dy / (tip*y);

          const double cp0 = cos(phi);
          const double sp0 = sin(phi);
          
          const double u  = (R*R + Z*Z - d*d) / (2*R*st0*cp0 - 2*Z*ct0 - 2*d);
          const double v  = d - u;
          
          if (u <= 0) { continue; }
          if (v <= 0) { continue; }
          
          const double vi  = 1.0 / v;
          const double cts = (R*st0*cp0 - Z*ct0 - u)*vi;         // cosine scattering angle
          
          const double V  = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
          
          if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
          
          const double vx  =   R - u*st0*cp0;
          const double vy  =     - u*st0*sp0;
          const double vz  =  -Z - u*ct0;
          
          const double ct[]  = {                                 // cosine angle of incidence on PMT
            (vx*px + vy*py + vz*pz) * vi,
            (vx*px - vy*py + vz*pz) * vi
          };

          const double U  = 
            getAngularAcceptance(ct[0]) + 
            getAngularAcceptance(ct[1]);                         // PMT angular acceptance
          const double W  = std::min(A*vi*vi, 2.0*PI);           // solid angle

          const double Ja = getScatteringProbability(cts);       // d^2P/dcos/dphi
          const double Jd = ng * (1.0 - cts) / C;                // dt/du
          
          value += (npe * dp / (2*PI)) * U * V * W * Ja * Jc / fabs(Jd);
        }
      }
    
      return value;
    }


    /**
     * Probability density function for direct light from EM-shower.
     *
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getDirectLightFromEMshower(const double D_m,
                                      const double cd,
                                      const double theta,
                                      const double phi,
                                      const double t_ns) const
    {
      using namespace JTOOLS;

      const double ct0 = cd;
      const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));

      const double D  =  std::max(D_m, getRmin());
      const double t  =  D * getIndexOfRefraction() / C  +  t_ns;  // time [ns]
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);
      const double ng = C*t/D;                                     // index of refraction

      if (n0 >= ng) { return 0.0; }
      if (n1 <= ng) { return 0.0; }

      const double w  = getWavelength(ng);
      const double n  = getIndexOfRefractionPhase(w);

      const double l_abs = getAbsorptionLength(w);
      const double ls    = getScatteringLength(w);

      const double npe   = cherenkov(w,n) * getQE(w);

      const double ct = st0*px + ct0*pz;                           // cosine angle of incidence on PMT

      const double U  = getAngularAcceptance(ct);                  // PMT angular acceptance
      const double V  = exp(-D/l_abs - D/ls);                      // absorption & scattering
      const double W  = A/(D*D);                                   // solid angle

      const double ngp = getDispersionGroup(w);

      const double Ja = D * ngp / C;                               // dt/dlambda
      const double Jb = geant(n,ct0);                              // d^2N/dcos/dphi

      return npe * geanc() * U * V * W * Jb / fabs(Ja);
    }


    /**
     * Probability density function for scattered light from EM-shower.
     *
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getScatteredLightFromEMshower(const double D_m,
                                         const double cd,
                                         const double theta,
                                         const double phi,
                                         const double t_ns) const
    {
      using namespace JTOOLS;
    
      double value = 0;

      const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
      const double D  =  std::max(D_m, getRmin());
      const double L  =  D;
      const double t  =  D * getIndexOfRefraction() / C  +  t_ns;  // time [ns]

      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double qx = cd*px +  0  - sd*pz;
      const double qy =         1*py;
      const double qz = sd*px +  0  + cd*pz;

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);

      const double ni = C*t / L;                                   // maximal index of refraction
    
      if (n0 >= ni) {
        return value;
      }

      const double nj = std::min(ni,n1);

      double w = wmax;
      
      for (const_iterator i = begin(); i != end(); ++i) {
        
        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
          
        w  = getWavelength(ng, w, 1.0e-5);
          
        const double dw = dn / fabs(getDispersionGroup(w));
        
        const double n  = getIndexOfRefractionPhase(w);   
        
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);
        
        const double npe   = cherenkov(w,n) * dw * getQE(w);

        if (npe <= 0) { continue; }

        const double Jc    = 1.0 / ls;                             // dN/dx
        
        const double d =  C*t / ng;                                // photon path
        
        //const double V  = exp(-d/l_abs);                           // absorption
          
        const double ds = 2.0 / (size() + 1);

        for (double sb = 0.5*ds; sb < 1.0 - 0.25*ds; sb += ds) {
          
          for (int buffer[] = { -1, +1, 0}, *k = buffer; *k != 0; ++k) {
            
            const double cb  = (*k) * sqrt((1.0 + sb)*(1.0 - sb));
            const double dcb = (*k) * ds*sb/cb;
            
            const double v  = 0.5 * (d + L) * (d - L) / (d - L*cb);
            const double u  = d - v;
            
            if (u <= 0) { continue; }
            if (v <= 0) { continue; }
            
            const double cts = (L*cb - v) / u;                     // cosine scattering angle
            
            const double V  = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));

            if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }

            const double W  = std::min(A/(v*v), 2.0*PI);           // solid angle
            const double Ja = getScatteringProbability(cts);       // d^2P/dcos/dphi
            const double Jd = ng * (1.0 - cts) / C;                // dt/du
            
            const double ca = (L - v*cb) / u;
            const double sa =      v*sb  / u;
            
            const double dp  = PI / phd.size();
            const double dom = dcb*dp * v*v / (u*u);
            const double dos = sqrt(dom);
           
            for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
              
              const double cp = l->getX();
              const double sp = l->getY();
              
              const double ct0 = cd*ca - sd*sa*cp;
              
              const double vx  = -sb*cp * qx;
              const double vy  = -sb*sp * qy;
              const double vz  =  cb    * qz;
              
              const double U  = 
                getAngularAcceptance(vx + vy + vz) +
                getAngularAcceptance(vx - vy + vz);                // PMT angular acceptance
              
              //const double Jb = geant(n,ct0);                      // d^2N/dcos/dphi
              
              //value += npe * geanc() * dom * U * V * W * Ja * Jb * Jc / fabs(Jd);
             
              const double Jb = geant(n, 
                                      ct0 - 0.5*dos, 
                                      ct0 + 0.5*dos);              // dN/dphi
             
              value += npe * geanc() * dos * U * V * W * Ja * Jb * Jc / fabs(Jd);
            }
          }
        }
      }
    
      return value;
    }


    /**
     * Probability density function for direct light from EM-shower.
     *
     * \param  E          EM-shower energy [GeV]
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getDirectLightFromEMshower(const double E,
                                      const double D_m,
                                      const double cd,
                                      const double theta,
                                      const double phi,
                                      const double t_ns) const
    {
      using namespace JTOOLS;

      double value = 0;

      const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
      const double D  =  std::max(D_m, getRmin());
      const double R  =  D * sd;                                   // minimal distance of approach [m]
      const double Z  = -D * cd;
      const double t  =  D * getIndexOfRefraction() / C  +  t_ns;  // time [ns]

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);

      double zmin  =  0.0;                                         // minimal shower length [m]
      double zmax  =  geanz.getLength(E, 1.0);                     // maximal shower length [m]

      const double d  = sqrt((Z+zmax)*(Z+zmax) + R*R);

      if (C*t > std::max(n1*D, zmax + n1*d)) {
        return value;
      }

      if (C*t < n0*D) {
        
        JRoot rz(R, n0, t + Z/C);                                  // square root
        
        if (!rz.is_valid) {
          return value;
        }

        if (rz.second > Z) { zmin = rz.second - Z; }
        if (rz.first  > Z) { zmin = rz.first  - Z; }
      }
        
      if (C*t > n1*D) {
        
        JRoot rz(R, n1, t + Z/C);                                  // square root
        
        if (!rz.is_valid) {
          return value;
        }

        if (rz.second > Z) { zmin = rz.second - Z; }
        if (rz.first  > Z) { zmin = rz.first  - Z; }
      }
      
      if (C*t < zmax + n0*d) {
        
        JRoot rz(R, n0, t + Z/C);                                  // square root
        
        if (!rz.is_valid) {
          return value;
        }

        if (rz.first  > Z) { zmax = rz.first  - Z; }
        if (rz.second > Z) { zmax = rz.second - Z; }
      }

      if (zmin < 0.0) {
        zmin = 0.0;
      }
      
      if (zmax > zmin) {

        const double ymin  =  geanz.getIntegral(E, zmin);
        const double ymax  =  geanz.getIntegral(E, zmax);
        const double dy    =  (ymax - ymin) / size();
        
        if (dy > 2*std::numeric_limits<double>::epsilon()) {
         
          for (double y = ymin + 0.5*dy; y < ymax; y += dy) {

            const double z   =  Z  +  geanz.getLength(E, y);
            const double d   =  sqrt(R*R + z*z);
            const double t1  =  t  +  (Z - z) / C  -  d * getIndexOfRefraction() / C;
            
            value += dy * E * getDirectLightFromEMshower(d, -z/d, theta, phi, t1);
          }
         
        }
      }

      return value;
    }


    /**
     * Probability density function for scattered light from EM-shower.
     *
     * \param  E          EM-shower energy [GeV]
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getScatteredLightFromEMshower(const double E,
                                         const double D_m,
                                         const double cd,
                                         const double theta,
                                         const double phi,
                                         const double t_ns) const
    {
      using namespace JTOOLS;
  
      double value = 0;

      const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
      const double D  =  std::max(D_m, getRmin());
      const double R  =  D * sd;                                   // minimal distance of approach [m]
      const double Z  = -D * cd;
      const double t  =  D * getIndexOfRefraction() / C  +  t_ns;  // time [ns]

      const double n0 = getIndexOfRefractionGroup(wmax);

      double zmin  =  0.0;                                         // minimal shower length [m]
      double zmax  =  geanz.getLength(E, 1.0);                     // maximal shower length [m]

      const double d  =  sqrt((Z+zmax)*(Z+zmax) + R*R);
        
      if (C*t < n0*D) {
                     
        JRoot rz(R, n0, t + Z/C);                                  // square root
          
        if (!rz.is_valid) {
          return value;
        }

        if (rz.second > Z) { zmin = rz.second - Z; }
        if (rz.first  > Z) { zmin = rz.first  - Z; }
      }
        
      if (C*t < zmax + n0*d) {
          
        JRoot rz(R, n0, t + Z/C);                                  // square root
          
        if (!rz.is_valid) {
          return value;
        }

        if (rz.first  > Z) { zmax = rz.first  - Z; }
        if (rz.second > Z) { zmax = rz.second - Z; }
      }
        
      const double ymin  =  geanz.getIntegral(E, zmin);
      const double ymax  =  geanz.getIntegral(E, zmax);
      const double dy    =  (ymax - ymin) / size();

      if (dy > 2*std::numeric_limits<double>::epsilon()) {

        for (double y = ymin + 0.5*dy; y < ymax; y += dy) {
         
          const double z  =  Z  +  geanz.getLength(E, y);
          const double d  =  sqrt(R*R + z*z);
          const double t1 =  t  +  (Z - z) / C  -  d * getIndexOfRefraction() / C;
          
          value += dy * E * getScatteredLightFromEMshower(d, -z/d, theta, phi, t1);
        }
      }

      return value;
    }


    /**
     * Probability density function for direct light from delta-rays.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns x m/GeV]
     */
    double getDirectLightFromDeltaRays(const double R_m,
                                       const double theta,
                                       const double phi,
                                       const double t_ns) const
    {
      using namespace JTOOLS;

      double value = 0;

      const double R  =  std::max(R_m, getRmin());
      const double t  =  R * getTanThetaC() / C  +  t_ns;          // time [ns]
      const double A  =  getPhotocathodeArea();
      const double D  =  2.0*sqrt(A/PI);

      const double px =  sin(theta)*cos(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);
      const double ni = sqrt(R*R + C*t*C*t) / R;                   // maximal index of refraction
    
      if (n0 >= ni) return value;

      const double nj = std::min(n1, ni);

      double w = wmax;

      for (const_iterator i = begin(); i != end(); ++i) {
            
        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
        
        w  = getWavelength(ng, w, 1.0e-5);

        const double dw = dn / fabs(getDispersionGroup(w));

        const double n  = getIndexOfRefractionPhase(w);   
          
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);

        const double npe   = cherenkov(w,n) * dw * getQE(w);

        JRoot rz(R, ng, t);                                        // square root

        if (!rz.is_valid) { continue; }

        for (int j = 0; j != 2; ++j) {

          const double z = rz[j];
          
          if (C*t <= z) continue;

          const double d  = sqrt(z*z + R*R);                       // distance traveled by photon

          const double ct0 = -z / d;
          const double st0 =  R / d;
          
          const double dz = D / st0;                               // average track length
          const double ct = st0*px + ct0*pz;                       // cosine angle of incidence on PMT
          
          const double U  = getAngularAcceptance(ct);              // PMT angular acceptance
          const double V  = exp(-d/l_abs - d/ls);                  // absorption & scattering
          const double W  = A/(d*d);                               // solid angle

          const double Ja = 1.0/(4*PI);                            // d^2N/dcos/dphi
          const double Jd = (1.0 - ng*ct0) / C;                    // d  t/ dz
          const double Je = ng*st0*st0*st0 / (R*C);                // d^2t/(dz)^2

          value += npe * geanc() * U * V * W * Ja / (fabs(Jd) + 0.5*Je*dz);
        }
      }

      return value;
    }

    
    /**
     * Probability density function for direct light from delta-rays.
     *
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dx [npe/ns x m/GeV]
     */
    double getScatteredLightFromDeltaRays(const double R_m,
                                          const double theta,
                                          const double phi,
                                          const double t_ns) const
    {
     using namespace JTOOLS;

      double value = 0;

      const double R  =  std::max(R_m, getRmin());
      const double t  =  R * getTanThetaC() / C  +  t_ns;          // time [ns]
      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double n0 = getIndexOfRefractionGroup(wmax);
      const double n1 = getIndexOfRefractionGroup(wmin);
      const double ni = sqrt(R*R + C*t*C*t) / R;                   // maximal index of refraction
    
      if (n0 >= ni) return value;

      const double nj = std::min(ni,n1);

      double w = wmax;

      for (const_iterator i = begin(); i != end(); ++i) {

        const double ng = 0.5 * (nj + n0)  +  i->getX() * 0.5 * (nj - n0);
        const double dn = i->getY() * 0.5 * (nj - n0);
        
        w  = getWavelength(ng, w, 1.0e-5);

        const double dw = dn / fabs(getDispersionGroup(w));
        
        const double n  = getIndexOfRefractionPhase(w);   
        
        const double l_abs = getAbsorptionLength(w);
        const double ls    = getScatteringLength(w);
          
        const double npe   = cherenkov(w,n) * dw * getQE(w);

        if (npe <= 0) { continue; }

        const double Jc = 1.0 / ls;                                // dN/dx
        
        JRoot rz(R, ng, t);                                        // square root

        if (!rz.is_valid) { continue; }

        const double zmin = rz.first;
        const double zmax = rz.second;

        const double zap  = 1.0 / l_abs;
        
        const double xmin = exp(zap*zmax);
        const double xmax = exp(zap*zmin);
        
        for (const_iterator j = begin(); j != end(); ++j) {
          
          const double x  = 0.5 * (xmax + xmin)  +  j->getX() * 0.5 * (xmax - xmin);
          const double dx = j->getY() * 0.5 * (xmax - xmin);
          
          const double z  = log(x) / zap;
          const double dz = -dx / (zap*x);
          
          const double D  = sqrt(z*z + R*R);
          const double cd = -z / D;
          const double sd =  R / D;
          
          const double qx = cd*px +  0  - sd*pz;
          const double qy =         1*py;
          const double qz = sd*px +  0  + cd*pz;

          const double d = (C*t - z) / ng;                         // photon path
            
          //const double V  = exp(-d/l_abs);                         // absorption

          const double ds = 2.0 / (size() + 1);
              
          for (double sb = 0.5*ds; sb < 1.0 - 0.25*ds; sb += ds) {
                
            for (int buffer[] = { -1, +1, 0}, *k = buffer; *k != 0; ++k) {
                  
              const double cb  = (*k) * sqrt((1.0 + sb)*(1.0 - sb));
              const double dcb = (*k) * ds*sb/cb;
                  
              const double v  = 0.5 * (d + D) * (d - D) / (d - D*cb);
              const double u  = d - v;
                  
              if (u <= 0) { continue; }
              if (v <= 0) { continue; }
                    
              const double cts = (D*cb - v) / u;                   // cosine scattering angle
                    
              const double V  = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
              
              if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
              
              const double W  = std::min(A/(v*v), 2.0*PI);         // solid angle
              const double Ja = getScatteringProbability(cts);     // d^2P/dcos/dphi
              const double Jd = ng * (1.0 - cts) / C;              // dt/du

              const double dp  = PI / phd.size();
              const double dom = dcb*dp * v*v / (u*u);

              for (const_iterator l = phd.begin(); l != phd.end(); ++l) {

                const double cp  = l->getX();
                const double sp  = l->getY();

                const double vx  = sb*cp * qx;
                const double vy  = sb*sp * qy;
                const double vz  = cb    * qz;

                const double U  = 
                  getAngularAcceptance(vx + vy + vz) +
                  getAngularAcceptance(vx - vy + vz);              // PMT angular acceptance

                const double Jb = 1.0/(4*PI);                      // d^2N/dcos/dphi

                value += dom * npe * geanc() * dz * U * V * W * Ja * Jb * Jc / fabs(Jd);
              }
            }
          }
        }
      }

      return value;
    }


    /**
     * Probability density function for light from muon.
     *
     * \param  type       PDF type
     * \param  E_GeV      muon energy [GeV]
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dx [npe/ns]
     */
    double getLightFromMuon(const int    type,
                            const double E_GeV,
                            const double R_m,
                            const double theta,
                            const double phi,
                            const double t_ns) const
    {
      switch (type) {

      case DIRECT_LIGHT_FROM_MUON:
        return getDirectLightFromMuon        (R_m, theta, phi, t_ns);

      case SCATTERED_LIGHT_FROM_MUON:
        return getScatteredLightFromMuon     (R_m, theta, phi, t_ns);

      case DIRECT_LIGHT_FROM_EMSHOWERS:
        return getDirectLightFromEMshowers   (R_m, theta, phi, t_ns) * E_GeV;
        
      case SCATTERED_LIGHT_FROM_EMSHOWERS:
        return getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV;

      case DIRECT_LIGHT_FROM_DELTARAYS:
        return getDirectLightFromDeltaRays   (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);
        
      case SCATTERED_LIGHT_FROM_DELTARAYS:
        return getScatteredLightFromDeltaRays(R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);

      default:
        return 0.0;
      }      
    }


    /**
     * Probability density function for light from muon.
     *
     * \param  E_GeV      muon energy [GeV]
     * \param  R_m        distance between muon and PMT [m]
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dx [npe/ns]
     */
    double getLightFromMuon(const double E_GeV,
                            const double R_m,
                            const double theta,
                            const double phi,
                            const double t_ns) const
    {
      return (getDirectLightFromMuon        (R_m, theta, phi, t_ns)                                +
              getScatteredLightFromMuon     (R_m, theta, phi, t_ns)                                +
              getDirectLightFromEMshowers   (R_m, theta, phi, t_ns) * E_GeV                        +
              getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV                        +
              getDirectLightFromDeltaRays   (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV)  + 
              getScatteredLightFromDeltaRays(R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV));
    }

    
    /**
     * Probability density function for light from EM-shower.
     *
     * \param  type       PDF type
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getLightFromEMshower(const int    type,
                                const double D_m,
                                const double cd,
                                const double theta,
                                const double phi,
                                const double t_ns) const
    {
      switch (type) {

      case DIRECT_LIGHT_FROM_EMSHOWER:
        return getDirectLightFromEMshower   (D_m, cd, theta, phi, t_ns);

      case SCATTERED_LIGHT_FROM_EMSHOWER:
        return getScatteredLightFromEMshower(D_m, cd, theta, phi, t_ns);

      default:
        return 0.0;
      }      
    }


    /**
     * Probability density function for light from EM-shower.
     *
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getLightFromEMshower(const double D_m,
                                const double cd,
                                const double theta,
                                const double phi,
                                const double t_ns) const
    {
      return (getDirectLightFromEMshower   (D_m, cd, theta, phi, t_ns) +
              getScatteredLightFromEMshower(D_m, cd, theta, phi, t_ns));
    }

    
    /**
     * Probability density function for light from EM-shower.
     *
     * \param  E_GeV      EM-shower energy [GeV]
     * \param  type       PDF type
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getLightFromEMshower(const int    type,
                                const double E_GeV,
                                const double D_m,
                                const double cd,
                                const double theta,
                                const double phi,
                                const double t_ns) const
    {
      switch (type) {

      case DIRECT_LIGHT_FROM_EMSHOWER:
        return getDirectLightFromEMshower   (E_GeV, D_m, cd, theta, phi, t_ns);

      case SCATTERED_LIGHT_FROM_EMSHOWER:
        return getScatteredLightFromEMshower(E_GeV, D_m, cd, theta, phi, t_ns);

      default:
        return 0.0;
      }      
    }


    /**
     * Probability density function for light from EM-shower.
     *
     * \param  E_GeV      EM-shower energy [GeV]
     * \param  D_m        distance between EM-shower and PMT [m]
     * \param  cd         cosine angle EM-shower direction and EM-shower - PMT position
     * \param  theta      zenith  angle orientation PMT [rad]
     * \param  phi        azimuth angle orientation PMT [rad]
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getLightFromEMshower(const double E_GeV,
                                const double D_m,
                                const double cd,
                                const double theta,
                                const double phi,
                                const double t_ns) const
    {
      return (getDirectLightFromEMshower   (E_GeV, D_m, cd, theta, phi, t_ns) +
              getScatteredLightFromEMshower(E_GeV, D_m, cd, theta, phi, t_ns));
    }

    
    /**
     * Auxiliary class to find solution(s) to \f$z\f$ of the square root expression:
     \f{eqnarray*}
                    ct(z=0) & = & z + n \sqrt(z^2 + R^2)
     \f}
     * where \f$n = 1/\cos(\theta_{c})\f$ is the index of refraction.
     */
    class JRoot {
    public:
      /**
       * Determine solution(s) of quadratic equation.
       *
       * \param  R          minimal distance of approach [m]
       * \param  n          index of refraction
       * \param  t          time at z = 0 [ns]
       */
      JRoot(const double R,
            const double n,
            const double t) :
        is_valid(false),
        first (0.0),
        second(0.0)
      {
        using JTOOLS::C;

        const double a  =  n*n        -  1.0;
        const double b  =  2*C*t;
        const double c  =  R*n * R*n  -  C*t * C*t;
        
        const double q  =  b*b - 4*a*c;
                
        if (q >= 0.0) {

          first    = (-b - sqrt(q)) / (2*a);
          second   = (-b + sqrt(q)) / (2*a);

          is_valid =  C*t > second;
        }
      }


      /**
       * Get result by index.
       *
       * \param  i       index (0 or 1)
       * \return         i == 0, first; i == 1, second; else throws exception
       */
      double operator[](const int i) const
      {
        switch (i) {

        case 0:
          return first;
          
        case 1:
          return second;
          
        default:
          throw JException("JRoot::operator[] invalid index");
        }
      }

      bool   is_valid;    //!< validity of solutions
      double first;       //!< most upstream   solution
      double second;      //!< most downstream solution
    };


  protected:    
    /**
     * Determine wavelength for a given index of refraction corresponding to the group velocity.
     *
     * \param  n          index of refraction
     * \param  eps        precision index of refraction
     * \return            wavelength [nm]
     */
    virtual double getWavelength(const double n,
                                 const double eps = 1.0e-10) const
    {
      //return getWavelength(n, 0.5*(wmin + wmax), eps);

      double vmin = wmin;
      double vmax = wmax;

      for (int i = 0; i != 1000; ++i) {

        const double v = 0.5 * (vmin + vmax); 
        const double y = getIndexOfRefractionGroup(v);
        
        if (fabs(y - n) < eps) {
          return v;
        }

        if (y < n) 
          vmax = v;
        else
          vmin = v;
      }

      return 0.5 * (vmin + vmax);
    }


    /**
     * Determine wavelength for a given index of refraction corresponding to the group velocity.
     * The estimate of the wavelength is made by successive linear extrapolations.
     * The procedure starts from the given wavelength and terminates if the index of refraction 
     * is equal to the target value within the given precision.
     *
     * \param  n          index of refraction
     * \param  w          start  value wavelength [nm]
     * \param  eps        precision index of refraction
     * \return            return value wavelength [nm]
     */
    virtual double getWavelength(const double n,
                                 const double w,
                                 const double eps) const
    {
      double v = w;

      // determine wavelength by linear extrapolation
      
      for ( ; ; ) {
        
        const double y = getIndexOfRefractionGroup(v);
        
        if (fabs(y - n) < eps) {
          break;
        }

        v += (n - y) / getDispersionGroup(v);   
      }

      return v;
    }


    static double getRmin() { return 1.0e-1; }       //!< minimal distance of approach of muon to PMT [m]


    /**
     * Get the inverse of the attenuation length.
     *
     * \param  l_abs      absorption length [m]
     * \param  ls         scattering length [m]
     * \param  cts        cosine scattering angle
     * \return            inverse attenuation length [m^-1]
     */
    virtual double getInverseAttenuationLength(const double l_abs, const double ls, const double cts) const
    {
      static JTOOLS::JGridPolint1Function1D_t f1;
      
      if (f1.empty()) {
        
        const int    nx   = 100000;
        const double xmin =   -1.0;
        const double xmax =   +1.0;
        
        const double dx   = (xmax - xmin) / (nx - 1);
        
        for (double x = xmin, W = 0.0; x < xmax; x += dx) {

          f1[x] = W;

          W += 2*PI * dx * getScatteringProbability(x + 0.5*dx);
        }

        f1[xmin] = 0.0;
        f1[xmax] = 1.0;

        f1.compile();
      }
        
      return 1.0/l_abs  +  f1(cts)/ls;
    }

    
  protected:
    /**
     * Integration limits
     */
    const double wmin;               //!<  minimal wavelength for integration [nm]
    const double wmax;               //!<  maximal wavelength for integration [nm]

    std::vector<element_type> phd;   //!<  fast evaluation of phi integral
  };


  /**
   * Probability Density Functions of the time response of a PMT
   * with an implementation for the JDispersionInterface interface.
   */
  class JAbstractPDF : 
    public JPDF, 
    public JDispersion 
  {
  public:
    /**
     * Constructor.
     *
     * \param  P_atm           ambient pressure [atm]
     * \param  Wmin            minimal wavelength for integration [nm]
     * \param  Wmax            maximal wavelength for integration [nm]
     * \param  numberOfPoints  number    of points for integration
     * \param  epsilon         precision of points for integration
     */
    JAbstractPDF(const double P_atm,
                 const double Wmin,
                 const double Wmax,
                 const int    numberOfPoints = 20,
                 const double epsilon        = 1e-12) :
      JPDF(Wmin, Wmax, numberOfPoints, epsilon),
      JDispersion(P_atm)
    {}
  };


  /**
   * Probability Density Functions of the time response of a PMT 
   * with an implementation of the JAbstractPMT and JAbstractMedium interfaces via C-like methods.
   */
  class JPDF_C :
    public JAbstractPDF
  {
  public:
    /**
     * Constructor.
     *
     * \param  Apmt            photo-cathode area [m^2]
     * \param  pF_qe           pointer to function for quantum efficiency of PMT
     * \param  pF_pmt          pointer to function for angular acceptence of PMT
     * \param  pF_l_abs        pointer to function for absorption length [m]
     * \param  pF_ls           pointer to function for scattering length [m]
     * \param  pF_ps           pointer to model specific function to describe light scattering in water
     * \param  P_atm           ambient pressure [atm]
     * \param  Wmin            minimal wavelength for integration [nm]
     * \param  Wmax            maximal wavelength for integration [nm]
     * \param  numberOfPoints  number    of points for integration
     * \param  epsilon         precision of points for integration
     */
    JPDF_C(const double Apmt,
           //
           // pointers to static functions
           //
           double (*pF_qe)   (const double),
           double (*pF_pmt)  (const double),
           double (*pF_l_abs)(const double),
           double (*pF_ls)   (const double),
           double (*pF_ps)   (const double),
           //
           // parameters for base class
           //
           const double P_atm,
           const double Wmin,
           const double Wmax,
           const int    numberOfPoints = 20,
           const double epsilon        = 1e-12) :
      JAbstractPDF(P_atm, Wmin, Wmax, numberOfPoints, epsilon),
      A    (Apmt),
      qe   (pF_qe),
      l_abs(pF_l_abs),
      ls   (pF_ls),
      pmt  (pF_pmt),
      ps   (pF_ps)
    {}


    /**
     * Photo-cathode area of PMT.
     *
     * \return            photo-cathode area [m^2]
     */
    virtual double getPhotocathodeArea() const 
    { 
      return A;
    }


    /**
     * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.).
     *
     * \param  lambda     wavelenth [nm]
     * \return            QE
     */
    virtual double getQE(const double lambda) const 
    { 
      return qe(lambda);
    }

 
    /**
     * Angular acceptence of PMT (normalised to one at cos() = -1).
     *
     * \param  ct         cosine angle of incidence
     * \return            relative efficiency
     */
    virtual double getAngularAcceptance(const double ct) const
    { 
      return pmt(ct);
    }


    /**
     * Absorption length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            absorption length [m]
     */
    virtual double getAbsorptionLength(const double lambda) const
    { 
      return l_abs(lambda);
    }


    /**
     * Scattering length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            scattering length [m]
     */
    virtual double getScatteringLength(const double lambda) const
    { 
      return ls(lambda);
    }


    /**
     * Model specific function to describe light scattering in water
     * (integral over full solid angle normalised to one).
     *
     * \param  ct         cosine scattering angle
     * \return            probability
     */
    virtual double getScatteringProbability(const double ct) const
    { 
      return ps(ct);
    }

  protected:
    /**
     * photo-cathode area [m2]
     */
    const double A;

    /**
     * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.)
     *
     * \param  lambda     wavelenth [nm]
     * \return            QE 
     */ 
    double (*qe)(const double lambda);

    /**
     * Absorption length
     *
     * \param  lambda     wavelenth [nm]
     * \return            absorption length [m]
     */ 
    double (*l_abs)(const double lambda);

    /**
     * Scattering length
     *
     * \param  lambda     wavelenth [nm]
     * \return            scattering length [m]
     */ 
    double (*ls)(const double lambda);

    /**
     * Angular acceptence of PMT (normalised to one at cos() = -1)
     *
     * \param  ct         cosine angle of incidence
     * \return            relative efficiency
     */
    double (*pmt)(const double ct);

    /**
     * Model specific function to describe light scattering in water
     *
     * \param  ct         cosine scattering angle
     * \return            probability
     */
    double (*ps)(const double ct);
  };
}


#endif