#ifndef __JPHYSICS__JPDF__
#define __JPHYSICS__JPDF__

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

#include "JLang/JCC.hh"
#include "JLang/JException.hh"
#include "JPhysics/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 "JPhysics/JGeane.hh"
#include "JPhysics/JGeant.hh"
#include "JPhysics/JGeanz.hh"


/**
 * \author mdejong
 */

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

namespace JPHYSICS {

  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 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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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
    {
      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 = getDeltaRayProbability(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 += 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
    {
      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 = getDeltaRayProbability(ct0);     // d^2N/dcos/dphi

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

      return value;
    }


    /**
     * Probability density function for direct light from isotropic light source.
     *
     * \param  D_m        distance between light source and PMT [m]
     * \param  ct         cosine angle PMT
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getDirectLightFromBrightPoint(const double D_m,
					 const double ct,
					 const double t_ns) const
    {
      const double D  =  std::max(D_m, getRmin());
      const double t  =  D * getIndexOfRefraction() / C  +  t_ns;  // time [ns]
      const double A  =  getPhotocathodeArea();

      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 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 = 1.0 / (4.0*PI);                            // d^2N/dOmega

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


    /**
     * Probability density function for scattered light from isotropic light source.
     *
     * \param  D_m        distance between light source and PMT [m]
     * \param  ct         cosine angle PMT
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getScatteredLightFromBrightPoint(const double D_m,
					    const double ct,
					    const double t_ns) const
    {
      double value = 0;

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

      const double A  =  getPhotocathodeArea();

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

      const double ni = C*t / D;                                   // 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 npe   = cherenkov(w,n) * dw * getQE(w);

	if (npe <= 0) { continue; }

	const double l_abs = getAbsorptionLength(w);
	const double ls    = getScatteringLength(w);
	
	const double Jc    = 1.0 / ls;                             // dN/dx
	const double Jb    = 1.0 / (4.0*PI);                       // d^2N/dcos/dphi
	
	const double d =  C*t / ng;                                // photon path
	
	const double dcb = 2.0 / (size() + 1);

	for (double cb = -1.0 + 0.5*dcb; cb < +1.0; cb += dcb) {

	  const double sb = sqrt((1.0 + cb)*(1.0 - 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);               // dOmega

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

	    const double cp  = l->getX();
	    const double dot = cb*ct + sb*cp*st;

	    const double U  = 2 * getAngularAcceptance(dot);     // PMT angular acceptance

	    value += npe * geanc() * dom * 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  type       PDF type
     * \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 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));
    }

    
    /**
     * Probability density function for direct light from isotropic light source.
     *
     * \param  type       PDF type
     * \param  D_m        distance between light source and PMT [m]
     * \param  ct         cosine angle PMT
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getLightFromBrightPoint(const int    type,
				   const double D_m,
				   const double ct,
				   const double t_ns) const
    {
      switch (type) {

      case DIRECT_LIGHT_FROM_BRIGHT_POINT:
	return getDirectLightFromBrightPoint   (D_m, ct, t_ns);

      case SCATTERED_LIGHT_FROM_BRIGHT_POINT:
	return getScatteredLightFromBrightPoint(D_m, ct, t_ns); 

      default:
	return 0.0;
      }      
    }

    
    /**
     * Probability density function for direct light from isotropic light source.
     *
     * \param  D_m        distance between light source and PMT [m]
     * \param  ct         cosine angle PMT
     * \param  t_ns       time difference relative to direct Cherenkov light [ns]
     * \return            d^2P/dt/dE [npe/ns/GeV]
     */
    double getLightFromBrightPoint(const double D_m,
				   const double ct,
				   const double t_ns) const
    {
      return (getDirectLightFromBrightPoint   (D_m, ct, t_ns)  +
	      getScatteredLightFromBrightPoint(D_m, ct, 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)
      {
	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 acceptance 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 override 
    { 
      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 override 
    { 
      return qe(lambda);
    }

 
    /**
     * Angular acceptance of PMT.
     *
     * \param  ct         cosine angle of incidence
     * \return            relative efficiency
     */
    virtual double getAngularAcceptance(const double ct) const override 
    { 
      return pmt(ct);
    }


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


    /**
     * Scattering length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            scattering length [m]
     */
    virtual double getScatteringLength(const double lambda) const override 
    { 
      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 override 
    { 
      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 acceptance of PMT
     *
     * \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