src/jpp/JMath/JMathSupportkit.hh

+#ifndef __JMATHSUPPORTKIT__
+#define __JMATHSUPPORTKIT__
+
+#include <limits>
+#include <cmath>
+
+#include "JMath/JConstants.hh"
+#include "JLang/JException.hh"
+
+
+/**
+ * \file
+ *
+ * Auxiliary methods for mathematics.
+ * \author mdejong
+ */
+
+namespace JMATH {}
+namespace JPP { using namespace JMATH; }
+
+namespace JMATH {
+
+  using JLANG::JValueOutOfRange;
+
+  
+  /**
+   * Gauss function (normalised to 1 at x = 0).
+   *
+   * \param  x                    x
+   * \param  sigma                sigma
+   * \return                      function value
+   */
+  inline double gauss(const double x, const double sigma)
+  {
+    const double u = x / sigma;
+
+    if (fabs(u) < 20.0)
+      return exp(-0.5*u*u);
+    else
+      return 0.0;
+  }
+
+
+  /**
+   * Gauss function (normalised to 1 at x = x0).
+   *
+   * \param  x                    x
+   * \param  x0                   central value
+   * \param  sigma                sigma
+   * \return                      function value
+   */
+  inline double gauss(const double x, const double x0, const double sigma)
+  {
+    return gauss(x - x0, sigma);
+  }
+
+
+  /**
+   * Normalised Gauss function.
+   *
+   * \param  x                    x
+   * \param  sigma                sigma
+   * \return                      function value
+   */
+  inline double Gauss(const double x, const double sigma)
+  {
+    return gauss(x, sigma) / sqrt(2.0*PI) / sigma;
+  }
+
+
+  /**
+   * Normalised Gauss function.
+   *
+   * \param  x                    x
+   * \param  x0                   central value
+   * \param  sigma                sigma
+   * \return                      function value
+   */
+  inline double Gauss(const double x, const double x0, const double sigma)
+  {
+    return Gauss(x - x0, sigma);
+  }
+
+
+  /**
+   * Incomplete gamma function.
+   *
+   * Source code is taken from reference:
+   * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
+   * Cambridge University Press.
+   *
+   * \param  a                    a
+   * \param  x                    x
+   * \return                      function value
+   */
+  inline double Gamma(const double a, const double x)
+  {
+    using namespace std;
+
+    const int max = 1000000;
+
+    if (x <  0.0) { THROW(JValueOutOfRange, "x <  0 " << x); }
+    if (a <= 0.0) { THROW(JValueOutOfRange, "a <= 0 " << a); }
+ 
+    const double gln = lgamma(a);
+
+    if (x < a + 1.0) {
+
+      if (x < 0.0) { 
+	THROW(JValueOutOfRange, "x <= 0 " << x);
+      }
+
+      if (x == 0.0) {
+	return 0.0;
+      }
+
+      double ap  = a;
+      double sum = 1.0 / a;
+      double del = sum;
+
+      for (int i = 0; i != max; ++i) {
+
+	ap  += 1.0;
+	del *= x/ap;
+	sum += del;
+
+	if (fabs(del) < fabs(sum)*numeric_limits<double>::epsilon()) {
+	  return sum*exp(-x + a*log(x) - gln);
+	}
+      }
+
+      THROW(JValueOutOfRange, "i == " << max);
+
+    } else {
+
+      const double FPMIN = numeric_limits<double>::min() / numeric_limits<double>::epsilon();
+
+      double b = x + 1.0 - a;
+      double c = 1.0 / FPMIN;
+      double d = 1.0 / b;
+      double h = d;
+
+      for (int i = 1; i != max; ++i) {
+
+	const double an = -i * (i-a);
+
+	b += 2.0;
+	d  = an*d + b;
+
+	if (fabs(d) < FPMIN) {
+	  d = FPMIN;
+	}
+
+	c  = b + an/c;
+
+	if (fabs(c) < FPMIN) {
+	  c = FPMIN;
+	}
+
+	d  = 1.0/d;
+
+	const double del = d*c;
+	
+	h *= del;
+	
+	if (fabs(del - 1.0) < numeric_limits<double>::epsilon()) {
+	  return 1.0 - exp(-x + a*log(x) - gln) * h;
+	}
+      }
+
+      THROW(JValueOutOfRange, "i == " << max);
+    }
+  }
+
+
+  /**
+   * Legendre polynome.
+   *
+   * \param  n                    degree 
+   * \param  x                    x
+   * \return                      function value
+   */
+  inline double legendre(const size_t n, const double x)
+  {
+    switch (n) {
+
+    case 0:
+      return 1.0;
+
+    case 1:
+      return x;
+
+    default:
+      {
+	double p0;
+	double p1 = 1.0;
+	double p2 = x;
+
+	for (size_t i = 2; i <= n; ++i) {
+	  p0 = p1;
+	  p1 = p2;
+	  p2 = ((2*i-1) * x*p1  -  (i-1) * p0) / i;
+	}
+
+	return p2;
+      }
+    }
+  }
+
+
+  /**
+   * Binomial function.
+   *
+   * \param  n                    n
+   * \param  k                    k
+   * \return                      function value
+   */
+  inline double binomial(const size_t n, const size_t k)
+  {
+    if (n == 0 || n < k) {
+      return 0.0;
+    }
+
+    if (k == 0 || n == k) {
+      return 1.0;
+    }
+
+    const int k1 = std::min(k, n - k);
+    const int k2 = n - k1;
+
+    double value = k2 + 1;
+
+    for (int i = k1; i != 1; --i) {
+      value *= (double) (k2 + i) / (double) i;
+    }
+
+    return value;
+  }
+
+
+  /**
+   * Poisson probability density distribition.
+   *
+   * \param  n                    number of occurences
+   * \param  mu                   expectation value
+   * \return                      probability
+   */
+  inline double poisson(const size_t n, const double mu)
+  {
+    using namespace std;
+
+    if (mu > 0.0) {
+
+      if (n > 0)
+	return exp(n*log(mu) - lgamma(n+1) - mu);
+      else
+	return exp(-mu);
+    } else if (mu == 0.0) {
+
+      return (n == 0 ? 1.0 : 0.0);
+    }
+ 
+    THROW(JValueOutOfRange, "mu <= 0 " << mu);
+  }
+
+
+  /**
+   * Poisson cumulative density distribition.
+   *
+   * \param  n                    number of occurences
+   * \param  mu                   expectation value
+   * \return                      probability
+   */
+  inline double Poisson(const size_t n, const double mu)
+  {
+    return 1.0 - Gamma(n + 1, mu);
+  }
+}
+
+#endif

src/jpp/JMath/JZero.hh

+#ifndef __JMATH__JZERO__
+#define __JMATH__JZERO__
+
+/**
+ * \file
+ *
+ * Definition of zero value for any class.
+ * \author mdejong
+ */
+namespace JMATH {}
+namespace JPP { using namespace JMATH; }
+
+namespace JMATH {
+
+  /**
+   * Get zero value for a given data type.
+   *
+   * The default implementation of this method returns an object which 
+   * is created with the default constructor.
+   * This method should be specialised if this value does not correspond
+   * to the equivalent of a zero result.
+   *
+   * \return                     zero
+   */
+  template<class T>
+  inline T getZero() 
+  {
+    return T();
+  }
+
+
+  /**
+   * Get zero value for <tt>bool</tt>.
+   *
+   * \return                     false
+   */
+  template<> inline bool getZero<bool>()
+  {
+    return false;
+  }
+
+
+  /**
+   * Get zero value for <tt>float</tt>.
+   *
+   * \return                     zero
+   */
+  template<> inline float getZero<float>()
+  { 
+    return float(0.0);
+  }
+
+
+  /**
+   * Get zero value for <tt>double</tt>.
+   *
+   * \return                     zero
+   */
+  template<>
+  inline double getZero<double>() 
+  {
+    return double(0.0);
+  }
+
+
+  /**
+   * Get zero value for <tt>long double</tt>.
+   *
+   * \return                     zero
+   */
+  template<>
+  inline long double getZero<long double>() 
+  {
+    return (long double)(0.0);
+  }
+
+
+  /**
+   * Auxiliary class to assign zero value.
+   */
+  struct JZero {
+    /** 
+     * Default constructor.
+     */
+    JZero()
+    {}
+
+
+    /**
+     * Type conversion operator.
+     *
+     * \return                     zero
+     */
+    template<class T>
+    operator T() const
+    {
+      return getZero<T>();
+    }
+  };
+
+
+  /**
+   * Function object to assign zero value.
+   */
+  static const JZero zero;
+}
+
+#endif

src/jpp/JOscProb/JBaselineCalculator.hh

+#ifndef __JOSCPROB__JBASELINECALCULATOR__
+#define __JOSCPROB__JBASELINECALCULATOR__
+
+#include "JIO/JSerialisable.hh"
+
+
+/**
+ * \author bjung
+ * Auxiliary data structure for storing and computing oscillation baselines.
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JIO::JReader;
+  using JIO::JWriter;
+  using JIO::JSerialisable;
+
+
+  /**
+   * Auxiliary data structure for storing and calculating baselines.
+   */
+  struct JBaselineCalculator :
+    public JSerialisable
+  {
+    /**
+     * Default constructor.
+     */
+    JBaselineCalculator() :
+      Lmin(0.0),
+      Lmax(0.0)
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  Lmin     Minimum baseline [km]
+     * \param  Lmax     Maximum baseline [km]
+     */
+    JBaselineCalculator(const double Lmin,
+			const double Lmax) :
+      Lmin(Lmin),
+      Lmax(Lmax)
+    {}
+
+
+    /**
+     * Get minimum baseline.
+     *
+     * \return                 maximum baseline [km]
+     */
+    double getMinimumBaseline() const
+    {
+      return Lmin;
+    }
+
+    
+    /**
+     * Get maximum baseline.
+     *
+     * \return                 maximum baseline [km]
+     */
+    double getMaximumBaseline() const
+    {
+      return Lmax;
+    }
+
+
+    /**
+     * Get inner radius.
+     *
+     * \return                inner radius [km]
+     */
+    double getInnerRadius() const
+    {
+      return 0.5 * (Lmax - Lmin);
+    }
+
+
+    /**
+     * Get outer radius.
+     *
+     * \return                outer radius [km]
+     */
+    double getOuterRadius() const
+    {
+      return 0.5 * (Lmax + Lmin);;
+    }
+    
+
+    /**
+     * Get cosine zenith angle for a given baseline.
+     *
+     * \param  L              baseline [km]
+     * \return                cosine zenith angle
+     */
+    double getCosth(const double L) const
+    {
+      static const double r = getInnerRadius();
+      static const double R = getOuterRadius();
+      
+      return (R*R - r*r - L*L) / (2*L*r);
+    }
+
+    
+    /**
+     * Get baseline for a given cosine zenith angle.
+     *
+     * \param  costh           cosine zenith angle
+     * \return                 baseline [km]
+     */
+    double getBaseline(const double costh) const
+    {
+      static const double r = getInnerRadius();
+      static const double R = getOuterRadius();
+
+      const double ct = (fabs(costh) < 1.0 ? costh : (costh < 0 ? -1.0 : 1.0));
+    
+      return (-r * ct + sqrt(R*R - r*r * (1 - ct) * (1 + ct)));
+    }
+
+
+    /**
+     * Get baseline for a given cosine zenith angle.
+     *
+     * \param  costh           cosine zenith angle
+     * \return                 baseline [km]
+     */
+    double operator()(const double costh) const
+    {
+      return getBaseline(costh);
+    }
+
+
+    /**
+     * Binary stream input of baseline extrema.
+     *
+     * \param  in              input stream
+     * \return                 input stream
+     */
+    JReader& read(JReader& in) override
+    {
+      return in >> Lmin >> Lmax;
+    }
+
+    
+    /**
+     * Binary stream output of oscillation parameters.
+     *
+     * \param  out             output stream
+     * \return                 output stream
+     */
+    JWriter& write(JWriter& out) const override
+    {
+      return out << Lmin << Lmax;
+    }
+
+
+    /**
+     * Stream input of baseline calculator.
+     *
+     * \param  in              input stream
+     * \param  object          object
+     * \return                 input stream
+     */
+    friend inline std::istream& operator>>(std::istream& in, JBaselineCalculator& object)
+    {
+      return in >> object.Lmin >> object.Lmax;
+    }
+
+
+    /**
+     * Stream output of baseline calculator.
+     *
+     * \param  out             output stream
+     * \param  object          object
+     * \return                 output stream
+     */
+    friend inline std::ostream& operator<<(std::ostream& out, const JBaselineCalculator& object)
+    {
+      return out << FIXED(15,5) << object.Lmin << FIXED(15,5) << object.Lmax;
+    }
+
+    
+  protected:
+    
+    double Lmin; //!< Minimum baseline [km]
+    double Lmax; //!< Maximum baseline [km]
+  };
+}
+
+#endif

src/jpp/JOscProb/JOscChannel.hh

+#ifndef __JOSCPROB__JOSCCHANNEL__
+#define __JOSCPROB__JOSCCHANNEL__
+
+#include <iostream>
+
+#include "JLang/JComparable.hh"
+
+#include "Jeep/JProperties.hh"
+
+
+/**
+ * \author bjung
+ * \file
+ * Oscillation channels and auxiliary methods.
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JLANG::JEquationParameters;
+
+
+  /**
+   * Neutrino flavours.
+   */
+  enum class JFlavour_t      { ELECTRON           = 12,
+			       MUON               = 14,
+			       TAU                = 16,
+			       FLAVOUR_UNDEFINED  =  0 };
+
+
+  /**
+   * Charge parities.
+   */
+  enum class JChargeParity_t { ANTIPARTICLE       = -1,
+			       PARTICLE           = +1,
+			       CPARITY_UNDEFINED  =  0 };
+  
+  
+  /**
+   * Auxiliary function for retrieving the flavour corresponding to a given PDG identifier.
+   *
+   * \param  pdgType             PDG particle identifier
+   */
+  inline JFlavour_t getFlavour(const int pdgType)
+  {
+    const int type = abs(pdgType);
+      
+    switch (type) {
+	
+    case (int) JFlavour_t::ELECTRON:
+    case (int) JFlavour_t::MUON:
+    case (int) JFlavour_t::TAU:
+      return static_cast<JFlavour_t>(type);
+    default:
+      return JFlavour_t::FLAVOUR_UNDEFINED;
+    }
+  }
+
+
+  /**
+   * Auxiliary function for retrieving the charge-parity of a given PDG type.
+   *
+   * \param  pdgType             PDG particle identifier
+   * \return                     charge-parity (1 for neutrinos; -1 for anti-neutrinos)
+   */
+  inline JChargeParity_t getChargeParity(const int pdgType)
+  {
+    if        (pdgType < 0) {
+      return JChargeParity_t::ANTIPARTICLE;
+    } else if (pdgType > 0) {
+      return JChargeParity_t::PARTICLE;
+    } else {
+      return JChargeParity_t::CPARITY_UNDEFINED;
+    }
+  }
+  
+  
+  /**
+   * Neutrino oscillation channel.
+   */  
+  struct JOscChannel :
+    public JLANG::JComparable<JOscChannel>
+  {
+    /**
+     * Default constructor.
+     */
+    JOscChannel() :
+      in     (JFlavour_t::FLAVOUR_UNDEFINED),
+      out    (JFlavour_t::FLAVOUR_UNDEFINED),
+      Cparity(JChargeParity_t::CPARITY_UNDEFINED)
+    {}
+
+
+    /**
+     * Constructor.				
+     *
+     * \param  in              input  flavour
+     * \param  out             output flavour
+     * \param  Cparity         charge parity
+     */
+    JOscChannel(const JFlavour_t      in,
+		const JFlavour_t      out,
+		const JChargeParity_t Cparity) :
+      in (in),
+      out(out),
+      Cparity(Cparity)
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  in              input  flavour
+     * \param  out             output flavour
+     * \param  Cparity         charge parity
+     */
+    JOscChannel(const int in,
+		const int out,
+		const int Cparity) :
+      in (getFlavour(in)),
+      out(getFlavour(out)),
+      Cparity(getChargeParity(Cparity))
+    {}
+
+
+    /**
+     * Check validity of this oscillation channel.
+     *
+     * \return                 true if this oscillation channel is valid; else false.
+     */
+    inline bool is_valid() const
+    {
+      return (in      != JFlavour_t::FLAVOUR_UNDEFINED &&
+	      out     != JFlavour_t::FLAVOUR_UNDEFINED &&
+	      Cparity != JChargeParity_t::CPARITY_UNDEFINED);
+    }
+
+
+    /**
+     * Less-than method
+     *
+     * \param  channel         channel
+     * \return                 true this channel less than given channel; else false 
+     */
+    inline bool less(const JOscChannel& channel) const
+    {
+      if (this->Cparity == channel.Cparity) {
+	
+	if (this->in == channel.in) {
+	  
+	  return this->out < channel.out;
+	  
+	} else {
+	  
+	  return this->in < channel.in;
+	}
+	
+      } else {
+	
+	return this->Cparity < channel.Cparity;
+      }
+    }    
+
+
+    /**
+     * Write channel to output.
+     *
+     * \param  out           output stream
+     * \param  object        oscillation channel
+     * \return               output stream
+     */
+    friend inline std::ostream& operator<<(std::ostream& out, const JOscChannel& object)
+    {
+      return out << object.getProperties();
+    }
+
+
+    /** 
+     * Read channel from input.
+     *
+     * \param  in            input stream
+     * \param  object        oscillation channel
+     * \return               input stream
+     */
+    friend inline std::istream& operator>>(std::istream& in, JOscChannel& object)
+    {
+      JProperties properties(object.getProperties());
+
+      in >> properties;
+
+      object.setProperties(properties);
+      
+      return in;
+    }
+
+
+    /**
+     * Get equation parameters.
+     *
+     * \return                 equation parameters
+     */
+    static inline JEquationParameters& getEquationParameters()
+    {
+      static JEquationParameters equation("=", "\n\r;,", "./", "#");
+
+      return equation;
+    }
+      
+
+    /**
+     * Set equation parameters.
+     *
+     * \param  equation        equation parameters
+     */
+    static inline void setEquationParameters(const JEquationParameters& equation)
+    {
+      getEquationParameters() = equation;
+    }
+
+    
+    /**
+     * Get properties of this class.
+     *
+     * \param  equation        equation parameters
+     */
+    JProperties getProperties(const JEquationParameters& equation = JOscChannel::getEquationParameters())
+    {
+      return JOscChannelHelper(*this, equation);
+    }
+
+
+    /**
+     * Get properties of this class.
+     *
+     * \param  equation        equation parameters
+     */
+    JProperties getProperties(const JEquationParameters& equation = JOscChannel::getEquationParameters()) const
+    {
+      return JOscChannelHelper(*this, equation);
+    }
+
+
+    /**
+     * Set properties of this class
+     *
+     * \param  properties      properties
+     */
+    void setProperties(const JProperties& properties)
+    {
+      this->in      = getFlavour     (properties.getValue<const int>("in"));
+      this->out     = getFlavour     (properties.getValue<const int>("out"));
+      this->Cparity = getChargeParity(properties.getValue<const int>("Cparity"));
+    }
+
+
+    JFlavour_t      in;      //!< Incoming  flavour
+    JFlavour_t      out;     //!< Outcoming flavour
+    JChargeParity_t Cparity; //!< Charge-parity
+
+    
+  private:
+    /**
+     * Auxiliary class for I/O of oscillation channel.
+     */
+    struct JOscChannelHelper :
+      public JProperties
+    {
+      /**
+       * Constructor.
+       *
+       * \param  object        oscillation channel
+       * \param  equation      equation    parameters 
+       */
+      template<class JOscChannel_t>
+      JOscChannelHelper(JOscChannel_t&             object,
+			const JEquationParameters& equation) :
+	JProperties(equation, 1),
+	in         ((int) object.in),
+	out        ((int) object.out),
+	Cparity    ((int) object.Cparity)
+      {
+	this->insert(gmake_property(in));
+	this->insert(gmake_property(out));
+	this->insert(gmake_property(Cparity));
+      }
+
+      int in;
+      int out;
+      int Cparity;
+    };
+  };
+
+
+  /**
+   * Declare group of neutrino oscillation channels.
+   */
+  static const JOscChannel getOscChannel[]  =  {
+    JOscChannel(JFlavour_t::ELECTRON,  JFlavour_t::ELECTRON,  JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::ELECTRON,  JFlavour_t::MUON,      JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::ELECTRON,  JFlavour_t::TAU,       JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::MUON,      JFlavour_t::ELECTRON,  JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::MUON,      JFlavour_t::MUON,      JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::MUON,      JFlavour_t::TAU,       JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::TAU,       JFlavour_t::ELECTRON,  JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::TAU,       JFlavour_t::MUON,      JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::TAU,       JFlavour_t::TAU,       JChargeParity_t::PARTICLE),
+    JOscChannel(JFlavour_t::ELECTRON,  JFlavour_t::ELECTRON,  JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::ELECTRON,  JFlavour_t::MUON,      JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::ELECTRON,  JFlavour_t::TAU,       JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::MUON,      JFlavour_t::ELECTRON,  JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::MUON,      JFlavour_t::MUON,      JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::MUON,      JFlavour_t::TAU,       JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::TAU,       JFlavour_t::ELECTRON,  JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::TAU,       JFlavour_t::MUON,      JChargeParity_t::ANTIPARTICLE),
+    JOscChannel(JFlavour_t::TAU,       JFlavour_t::TAU,       JChargeParity_t::ANTIPARTICLE)
+  };
+
+
+  /**
+   * Number of neutrino oscillation channels.
+   */
+  static const unsigned int NUMBER_OF_OSCCHANNELS = sizeof(getOscChannel) / sizeof(JOscChannel);
+}
+
+#endif

src/jpp/JOscProb/JOscParameters.hh

+#ifndef __JOSCPROB__JOSCPARAMETERS__
+#define __JOSCPROB__JOSCPARAMETERS__
+
+#include "JLang/JException.hh"
+
+#include "JOscProb/JOscParametersInterface.hh"
+
+
+/**
+ * \author bjung, mdejong
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JLANG::JEquationParameters;
+
+  using JIO::JReader;
+  using JIO::JWriter;
+  
+
+  /**
+   * Data structure for single set of oscillation parameters.
+   */
+  struct JOscParameters :
+    public JOscParametersInterface<double>
+  {
+    typedef JOscParametersInterface<double>                           JOscParameters_t;
+    typedef typename JOscParameters_t::JParameter_t                       JParameter_t;
+    typedef typename JOscParameters_t::argument_type                     argument_type;
+
+    
+    /**
+     * Default constructor.
+     */
+    JOscParameters() :
+      JOscParameters_t()
+    {}
+
+
+    /**
+     * Constructor.
+     * 
+     * \param  dM21sq          Squared mass difference between the first and second neutrino mass eigenstates               [eV2]
+     * \param  dM31sq          Squared mass difference between the first and third neutrino mass eigenstates                [eV2]
+     * \param  deltaCP         PMNS phase angle                                                                             [pi rad]
+     * \param  sinsqTh12       Squared sine of the PMNS mixing angle between the first and second neutrino mass eigenstates [-]
+     * \param  sinsqTh13       Squared sine of the PMNS mixing angle between the first and third neutrino mass eigenstates  [-]
+     * \param  sinsqTh23       Squared sine of the PMNS mixing angle between the second and third neutrino mass eigenstates [-]
+     */
+    JOscParameters(const double dM21sq,
+		   const double dM31sq,
+		   const double deltaCP,
+		   const double sinsqTh12,
+		   const double sinsqTh13,
+		   const double sinsqTh23) :
+      JOscParameters_t(dM21sq,
+		       dM31sq,
+		       deltaCP,
+		       sinsqTh12,
+		       sinsqTh13,
+		       sinsqTh23)
+    {
+      if (!is_valid()) {
+	THROW(JLANG::JValueOutOfRange, "JOscParameters::JOscParameters(...): Invalid parameters " << *this);
+      }
+    }
+
+
+    /**
+     * Constructor.
+     * 
+     * \param  name            parameter name
+     * \param  value           parameter value
+     * \param  args            remaining pairs of parameter names and values
+     */
+    template<class ...Args>
+    JOscParameters(const std::string& name,
+		   const double       value,
+		   const Args&     ...args) :
+      JOscParameters_t(name, value, args...)
+    {
+      if (!is_valid()) {
+	THROW(JLANG::JValueOutOfRange, "JOscParameters::JOscParameters(...): Invalid parameters " << *this);
+      }
+    }
+    
+    
+    /**
+     * Constructor.
+     * 
+     * Values taken from the NuFIT 5.1 three-flavour global analysis best fit values in:\n
+     * https://arxiv.org/abs/2111.03086?context=hep-ex\n
+     * including the Super-Kamiokande atmospheric data.
+     *
+     * \param  useIO           toggle inverted ordering
+     */
+    JOscParameters(const bool useIO) :
+      JOscParameters_t(          7.42e-5,
+			useIO ? -2.490e-3 + 7.42e-5 : 2.510e-3,
+			useIO ?  1.544              : 1.278,
+				 0.304,
+			useIO ?  0.02241            : 0.02246,
+			useIO ?  0.570              : 0.450    )
+    {}
+
+
+    /**
+     * Check validity of oscillation parameters.
+     *
+     * \return                 true if valid; else false
+     */
+    bool is_valid() const override
+    {
+      if ((this->sinsqTh12.isDefined() && this->sinsqTh12.getValue() < 0.0) ||
+	  (this->sinsqTh13.isDefined() && this->sinsqTh13.getValue() < 0.0) ||
+	  (this->sinsqTh23.isDefined() && this->sinsqTh23.getValue() < 0.0)) {
+	return false;
+      }
+      
+      return true;
+    }
+  };
+}
+
+#endif

src/jpp/JOscProb/JOscParametersGrid.hh

+#ifndef __JOSCPROB__JOSCPARAMETERSGRID__
+#define __JOSCPROB__JOSCPARAMETERSGRID__
+
+#include "JLang/JException.hh"
+
+#include "JTools/JGrid.hh"
+
+#include "JOscProb/JOscParametersInterface.hh"
+
+
+/**
+ * \author bjung, mdejong
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JTOOLS::JGrid;
+  using JTOOLS::make_grid;
+  
+
+  /**
+   * Data structure for oscillation parameter grids.
+   */
+  struct JOscParametersGrid :
+    public JOscParametersInterface<JGrid<double> >
+  {
+    typedef JGrid<double>                                                      JGrid_t;
+    typedef JOscParametersInterface<JGrid_t>                          JOscParameters_t;
+    typedef typename JOscParameters_t::JParameter_t                       JParameter_t;
+    typedef typename JOscParameters_t::argument_type                     argument_type;
+
+    
+    /**
+     * Default constructor.
+     */
+    JOscParametersGrid() :
+      JOscParameters_t()
+    {}
+
+
+    /**
+     * Constructor.
+     * 
+     * \param  dM21sq          Squared mass difference between the first and second neutrino mass eigenstates               [eV2]
+     * \param  dM31sq          Squared mass difference between the first and third neutrino mass eigenstates                [eV2]
+     * \param  deltaCP         PMNS phase angle                                                                             [pi rad]
+     * \param  sinsqTh12       Squared sine of the PMNS mixing angle between the first and second neutrino mass eigenstates [-]
+     * \param  sinsqTh13       Squared sine of the PMNS mixing angle between the first and third neutrino mass eigenstates  [-]
+     * \param  sinsqTh23       Squared sine of the PMNS mixing angle between the second and third neutrino mass eigenstates [-]
+     */
+    JOscParametersGrid(const JGrid_t& dM21sq,
+		       const JGrid_t& dM31sq,
+		       const JGrid_t& deltaCP,
+		       const JGrid_t& sinsqTh12,
+		       const JGrid_t& sinsqTh13,
+		       const JGrid_t& sinsqTh23) :
+      JOscParameters_t(dM21sq,
+		       dM31sq,
+		       deltaCP,
+		       sinsqTh12,
+		       sinsqTh13,
+		       sinsqTh23)
+    {
+      if (!is_valid()) {
+	THROW(JLANG::JValueOutOfRange, "JOscParametersGrid::JOscParametersGrid(...): Invalid parameters " << *this);
+      }
+    }
+
+
+    /**
+     * Constructor.
+     * 
+     * \param  name            parameter name
+     * \param  value           parameter value
+     * \param  args            remaining pairs of parameter names and values
+     */
+    template<class ...Args>
+    JOscParametersGrid(const std::string& name,
+		       const JGrid_t&     value,
+		       const Args&     ...args) :
+      JOscParameters_t(name, value, args...)      
+    {
+      if (!is_valid()) {
+	THROW(JLANG::JValueOutOfRange, "JOscParametersGrid::JOscParametersGrid(...): Invalid parameters " << *this);
+      }
+    }
+
+
+    /**
+     * Constructor.
+     * 
+     * Values taken from the NuFIT 5.1 three-flavour global analysis best fit values in:\n
+     * https://arxiv.org/abs/2111.03086?context=hep-ex\n
+     * including the Super-Kamiokande atmospheric data.
+     *
+     * \param  useIO           toggle inverted ordering
+     */
+    JOscParametersGrid(const bool useIO) :
+      JOscParameters_t( make_grid(         7.42e-5                       ),
+			make_grid(useIO ? -2.490e-3 + 7.42e-5 : 2.510e-3 ),
+			make_grid(useIO ?  1.544              : 1.278    ),
+			make_grid(         0.304                         ),
+			make_grid(useIO ?  0.02241            : 0.02246  ),
+			make_grid(useIO ?  0.570              : 0.450)     )
+    {}
+   
+
+    /**
+     * Check validity of oscillation parameter grids.
+     *
+     * \return                 true if valid; else false
+     */
+    bool is_valid() const override
+    {
+      if ((this->sinsqTh12.isDefined() && this->sinsqTh12.getValue().getXmin() < 0.0) ||
+	  (this->sinsqTh13.isDefined() && this->sinsqTh13.getValue().getXmin() < 0.0) ||
+	  (this->sinsqTh23.isDefined() && this->sinsqTh23.getValue().getXmin() < 0.0)) {
+	return false;
+      }
+      
+      return true;
+    }
+  };
+}
+
+#endif

src/jpp/JOscProb/JOscParametersInterface.hh

+#ifndef __JOSCPROB__JOSCPARAMETERSINTERFACE__
+#define __JOSCPROB__JOSCPARAMETERSINTERFACE__
+
+#include <iostream>
+#include <iomanip>
+#include <string>
+
+#include "JSystem/JStat.hh"
+
+#include "Jeep/JPrint.hh"
+#include "Jeep/JProperties.hh"
+
+#include "JIO/JSerialisable.hh"
+
+#include "JLang/JEquals.hh"
+#include "JLang/JParameter.hh"
+#include "JLang/JVectorize.hh"
+#include "JLang/JStringStream.hh"
+#include "JLang/JObjectStreamIO.hh"
+#include "JLang/JEquationParameters.hh"
+
+
+/**
+ * \author bjung, mdejong
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JIO::JReader;
+  using JIO::JWriter;
+  using JIO::JSerialisable;
+
+  using JLANG::JEquals;
+  using JLANG::JParameter;
+  using JLANG::JObjectStreamIO;
+  using JLANG::JValueOutOfRange;
+  using JLANG::JEquationParameters;
+
+  
+  /**
+   * Abstract base class for sets of oscillation parameters.
+   */
+  template<class T>
+  struct JOscParametersInterface :
+    public JSerialisable,
+    public JObjectStreamIO<JOscParametersInterface<T> >,
+    public JEquals        <JOscParametersInterface<T> >
+  {
+    typedef JOscParametersInterface<T>                      JOscParameters_t;
+    typedef JParameter<T>                                       JParameter_t;
+    typedef typename JParameter_t::argument_type               argument_type;
+
+
+    /**
+     * Default constructor.
+     */
+    JOscParametersInterface() :
+      dM21sq   (),
+      dM31sq   (),
+      deltaCP  (),
+      sinsqTh12(),
+      sinsqTh13(),
+      sinsqTh23()
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  dM21sq          Squared mass difference between the first and second neutrino mass eigenstates               [eV2]
+     * \param  dM31sq          Squared mass difference between the first and third neutrino mass eigenstates                [eV2]
+     * \param  deltaCP         PMNS phase angle                                                                             [pi rad]
+     * \param  sinsqTh12       Squared sine of the PMNS mixing angle between the first and second neutrino mass eigenstates [-]
+     * \param  sinsqTh13       Squared sine of the PMNS mixing angle between the first and third neutrino mass eigenstates  [-]
+     * \param  sinsqTh23       Squared sine of the PMNS mixing angle between the second and third neutrino mass eigenstates [-]
+     */
+    JOscParametersInterface(const T& dM21sq,
+			    const T& dM31sq,
+			    const T& deltaCP,
+			    const T& sinsqTh12,
+			    const T& sinsqTh13,
+			    const T& sinsqTh23) :
+      dM21sq   (dM21sq),
+      dM31sq   (dM31sq),
+      deltaCP  (deltaCP),
+      sinsqTh12(sinsqTh12),
+      sinsqTh13(sinsqTh13),
+      sinsqTh23(sinsqTh23)
+    {}
+
+
+    /**
+     * Constructor.
+     * 
+     * \param  name            parameter name
+     * \param  value           parameter value
+     * \param  args            remaining pairs of parameter names and values
+     */
+    template<class ...Args>
+    JOscParametersInterface(const std::string& name,
+			    const T&           value,
+			    const Args&     ...args)
+    {
+      set(name, value, args...);
+    }
+
+
+    /**
+     * Set value for a given oscillation parameter.
+     *
+     * \param  name            parameter name
+     * \param  value           parameter value
+     */
+    void set(const std::string& name,
+	     const T&           value)
+    {
+      JProperties properties = this->getProperties();
+
+      JProperties::iterator i = properties.find(name);
+
+      if (i != properties.end()) {
+
+	i->second.setValue(JParameter_t(value));
+	
+      } else {
+	
+        THROW(JValueOutOfRange,
+	      "template<class T> JOscParametersInterface<T>::set(const std::string&, const T&): " <<
+	      "Invalid oscillation parameter name " << name << "; Valid options:\n"  << JLANG::get_keys(properties));
+      }
+    }
+
+
+    /**
+     * Set value for given list of oscillation parameters.
+     *
+     * \param  name            parameter name
+     * \param  value           parameter value
+     * \param  args            remaining pairs of parameter names and values
+     */
+    template<class ...Args>
+    void set(const std::string& name,
+	     const T&           value,
+	     const Args&     ...args)
+    {
+      set(name, value);
+      set(args...);
+    }
+
+
+    /**
+     * Join the given oscillation parameters with these oscillation parameters.
+     *
+     * \param  parameters      oscillation parameters
+     */
+    JOscParameters_t& join(const JOscParameters_t& parameters)
+    {
+      if (parameters.dM21sq.isDefined())    { this->dM21sq    = parameters.dM21sq;    }
+      if (parameters.dM31sq.isDefined())    { this->dM31sq    = parameters.dM31sq;    }
+      if (parameters.deltaCP.isDefined())   { this->deltaCP   = parameters.deltaCP;   }
+      if (parameters.sinsqTh12.isDefined()) { this->sinsqTh12 = parameters.sinsqTh12; }
+      if (parameters.sinsqTh13.isDefined()) { this->sinsqTh13 = parameters.sinsqTh13; }
+      if (parameters.sinsqTh23.isDefined()) { this->sinsqTh23 = parameters.sinsqTh23; }
+      
+      return *this;
+    }
+   
+
+    /**
+     * Get oscillation parameters.
+     *
+     * \return                 oscillation parameters
+     */
+    const JOscParameters_t& getOscParameters() const
+    {
+      return static_cast<const JOscParameters_t&>(*this);
+    }
+
+
+    /**
+     * Set oscillation parameters.
+     *
+     * \param  parameters      oscillation parameters
+     */
+    void setOscParameters(const JOscParameters_t& parameters)
+    {
+      static_cast<JOscParameters_t&>(*this) = parameters;
+    }
+
+
+    /**
+     * Check validity of oscillation parameters.
+     *
+     * \return                 true if valid; else false
+     */
+    virtual bool is_valid() const = 0;
+
+
+    /**
+     * Get size of this oscillation parameters set.
+     *
+     * \return                 size (= number of defined parameters)
+     */
+    virtual unsigned int size() const
+    {
+      return ((unsigned int) this->dM21sq.isDefined()    +
+	      (unsigned int) this->dM31sq.isDefined()    +
+	      (unsigned int) this->deltaCP.isDefined()   +
+	      (unsigned int) this->sinsqTh12.isDefined() +
+	      (unsigned int) this->sinsqTh13.isDefined() +
+	      (unsigned int) this->sinsqTh23.isDefined());
+    }
+
+
+    /**
+     * Check if this oscillations parameter set contains the given oscillation parameters.
+     *
+     * \param  parameters      oscillation parameters
+     * \return                 true if all given oscillation parameters\n
+     *                         are also defined in this oscillation parameters set
+     */
+    virtual bool contains(const JOscParameters_t& parameters) const
+    {
+      if ( (parameters.dM21sq.isDefined()    && !this->dM21sq.isDefined())    ||
+	   (parameters.dM31sq.isDefined()    && !this->dM31sq.isDefined())    ||
+	   (parameters.deltaCP.isDefined()   && !this->deltaCP.isDefined())   ||
+	   (parameters.sinsqTh12.isDefined() && !this->sinsqTh12.isDefined()) ||
+	   (parameters.sinsqTh13.isDefined() && !this->sinsqTh13.isDefined()) ||
+	   (parameters.sinsqTh23.isDefined() && !this->sinsqTh23.isDefined()) ) {
+	return false;
+      } else {
+	return true;
+      }
+    }
+
+
+    bool equals(const JOscParameters_t& parameters) const
+    {
+      return (this->dM21sq    == parameters.dM21sq    &&
+	      this->dM31sq    == parameters.dM31sq    &&
+	      this->deltaCP   == parameters.deltaCP   &&
+	      this->sinsqTh12 == parameters.sinsqTh12 &&
+	      this->sinsqTh13 == parameters.sinsqTh13 &&
+	      this->sinsqTh23 == parameters.sinsqTh23);
+    }
+
+
+    /**
+     * Stream input of oscillation parameters.
+     *
+     * \param  in              input stream
+     * \param  parameters      oscillation parameters
+     * \return                 input stream
+     */
+    friend inline std::istream& operator>>(std::istream& in, JOscParameters_t& parameters)
+    {
+      using namespace std;
+      using namespace JPP;
+      
+      JStringStream is(in);
+
+      if (getFileStatus(is.str().c_str())) {
+	is.load();
+      }
+
+      JProperties properties(parameters.getProperties());
+
+      is >> properties;
+
+      parameters.setProperties(properties);
+
+      return in;
+    }
+
+
+    /**
+     * Stream output of oscillation parameters.
+     *
+     * \param  out             output stream
+     * \param  parameters      oscillation parameters
+     * \return                 output stream
+     */
+    friend inline std::ostream& operator<<(std::ostream& out, const JOscParameters_t& parameters)
+    {
+      return out << parameters.getProperties();
+    }
+
+
+    /**
+     * Binary stream input of oscillation parameters.
+     *
+     * \param  in              input stream
+     * \return                 input stream
+     */
+    JReader& read(JReader& in) override
+    {
+      JProperties properties = getProperties();
+
+      for (JProperties::iterator i = properties.begin(); i != properties.end(); ++i) {
+	
+	bool is_defined;
+	T    value;
+	
+	if ((in >> is_defined >> value) && is_defined) {
+
+	  JParameter_t& parameter = i->second.getValue<JParameter_t>();
+	  
+	  parameter.setValue(value);
+	}
+      }
+
+      return in;
+    }
+    
+
+    /**
+     * Binary stream output of oscillation parameters.
+     *
+     * \param  out             output stream
+     * \return                 output stream
+     */
+    JWriter& write(JWriter& out) const override
+    {
+      const JProperties properties = getProperties();
+
+      for (JProperties::const_iterator i = properties.cbegin(); i != properties.cend(); ++i) {
+
+	const JParameter_t& parameter = i->second.getValue<const JParameter_t>();	
+
+	out << parameter.isDefined() << parameter.getValue();
+      }
+
+      return out;
+    }
+    
+
+    /**
+     * Get equation parameters.
+     *
+     * \return                 equation parameters
+     */
+    static inline JEquationParameters& getEquationParameters()
+    {
+      static JEquationParameters equation("=", "\n\r;,", "./", "#");
+
+      return equation;
+    }
+      
+
+    /**
+     * Set equation parameters.
+     *
+     * \param  equation        equation parameters
+     */
+    static inline void setEquationParameters(const JEquationParameters& equation)
+    {
+      getEquationParameters() = equation;
+    }
+
+
+    /**
+     * Get properties of this class.
+     *
+     * \param  equation        equation parameters
+     */
+    JProperties getProperties(const JEquationParameters& equation = JOscParameters_t::getEquationParameters())
+    {
+      return JOscParametersHelper(*this, equation);
+    }
+
+
+    /**
+     * Get properties of this class.
+     *
+     * \param  equation        equation parameters
+     */
+    JProperties getProperties(const JEquationParameters& equation = JOscParameters_t::getEquationParameters()) const
+    {
+      return JOscParametersHelper(*this, equation);
+    }
+
+
+    /**
+     * Set properties of this class
+     *
+     * \param  properties      properties
+     */
+    void setProperties(const JProperties& properties)
+    {
+      this->dM21sq    = properties.getValue<JParameter_t>("dM21sq");
+      this->dM31sq    = properties.getValue<JParameter_t>("dM31sq");
+      this->deltaCP   = properties.getValue<JParameter_t>("deltaCP");
+      this->sinsqTh12 = properties.getValue<JParameter_t>("sinsqTh12");
+      this->sinsqTh13 = properties.getValue<JParameter_t>("sinsqTh13");
+      this->sinsqTh23 = properties.getValue<JParameter_t>("sinsqTh23");
+    }
+
+
+    JParameter_t dM21sq;     //!< Squared mass difference between the first and second neutrino mass eigenstates               [eV2]
+    JParameter_t dM31sq;     //!< Squared mass difference between the first and third neutrino mass eigenstates                [eV2]
+    JParameter_t deltaCP;    //!< PMNS phase angle                                                                             [pi * rad]
+    JParameter_t sinsqTh12;  //!< Squared sine of the PMNS mixing angle between the first and second neutrino mass eigenstates [-]
+    JParameter_t sinsqTh13;  //!< Squared sine of the PMNS mixing angle between the first and third neutrino mass eigenstates  [-]
+    JParameter_t sinsqTh23;  //!< Squared sine of the PMNS mixing angle between the second and third neutrino mass eigenstates [-]
+
+    
+  private:
+    
+    /**
+     * Auxiliary class for I/O of oscillation parameters.
+     */
+    struct JOscParametersHelper :
+      public JProperties
+    {
+      /**
+       * Constructor.
+       *
+       * \param  parameters    oscillation parameters
+       * \param  equation      equation    parameters 
+       */
+      template<class JOscParameters_t>
+      JOscParametersHelper(JOscParameters_t&          parameters,
+			   const JEquationParameters& equation) :
+	JProperties(equation, 1)
+      {
+	this->insert(gmake_property(parameters.dM21sq));
+	this->insert(gmake_property(parameters.dM31sq));
+	this->insert(gmake_property(parameters.deltaCP));	
+	this->insert(gmake_property(parameters.sinsqTh12));
+	this->insert(gmake_property(parameters.sinsqTh13));
+	this->insert(gmake_property(parameters.sinsqTh23));
+      }
+    };
+  };
+}
+
+#endif

src/jpp/JOscProb/JOscProb.hh

+#ifndef __JOSCPROB__JOSCPROB__
+#define __JOSCPROB__JOSCPROB__
+
+#include "JLang/JClonable.hh"
+
+#include "JOscProb/JOscChannel.hh"
+
+
+/**
+ * \author bjung
+ */
+
+namespace JOSCPROB {
+
+  using JLANG::JClonable;
+  
+
+  /**
+   * Low-level interface for retrieving the oscillation probability\n
+   * corresponding to a given oscillation channel, neutrino energy and zenith angle.
+   */
+  struct JOscProb :
+    public JClonable<JOscProb>
+  {
+    /**
+     * Virtual destructor.
+     */
+    virtual ~JOscProb()
+    {}
+
+    
+    /**
+     * Get oscillation probability for given oscillation channel.
+     *
+     * \param  channel            oscillation channel
+     * \param  energy             neutrino energy     [GeV]
+     * \param  costh              cosine zenith angle
+     * \return                    oscillation probability
+     */
+    virtual double getOscProb(const JOscChannel& channel,
+			      const double       energy,
+			      const double       costh) const = 0;
+  };
+}
+
+#endif

src/jpp/JOscProb/JOscProbInterpolator.hh

+#ifndef __JOSCPROB__JOSCPROBINTERPOLATOR__
+#define __JOSCPROB__JOSCPROBINTERPOLATOR__
+
+#include "Jeep/JMessage.hh"
+#include "Jeep/JProperties.hh"
+
+#include "JIO/JSerialisable.hh"
+#include "JIO/JFileStreamIO.hh"
+
+#include "JLang/JClonable.hh"
+#include "JLang/JObjectIO.hh"
+#include "JLang/JException.hh"
+
+#include "JTools/JPolint.hh"
+#include "JTools/JMapList.hh"
+#include "JTools/JCollection.hh"
+#include "JTools/JFunction1D_t.hh"
+#include "JTools/JMultiFunction.hh"
+#include "JTools/JFunctionalMap_t.hh"
+
+#include "JOscProb/JOscChannel.hh"
+#include "JOscProb/JOscParameters.hh"
+#include "JOscProb/JOscProbToolkit.hh"
+#include "JOscProb/JBaselineCalculator.hh"
+#include "JOscProb/JOscProbInterpolatorInterface.hh"
+
+
+/**
+ * \author bjung, mdejong
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JEEP::JMessage;
+
+  using JLANG::JClonable;
+
+  using JTOOLS::JMultiFunction;
+  
+
+  /** 
+   * Template definition of a multi-dimensional oscillation probability interpolation table.
+   */
+  template<template<class, class> class JCollection_t = JTOOLS::JCollection,
+	   class JFunction1D_t                        = JTOOLS::JPolintFunction1D <1,
+										   JTOOLS::JElement2D<double, JTOOLS::JArray<NUMBER_OF_OSCCHANNELS, double> >,
+										   JCollection_t,
+										   JTOOLS::JArray<NUMBER_OF_OSCCHANNELS, double> >,
+	   class JFunctionalMaplist_t                 = JTOOLS::JMAPLIST          <JTOOLS::JPolint1FunctionalMap,
+										   JTOOLS::JPolint1FunctionalMap,
+										   JTOOLS::JPolint1FunctionalMap,
+										   JTOOLS::JPolint1FunctionalMap,
+										   JTOOLS::JPolint1FunctionalMap,
+										   JTOOLS::JPolint1FunctionalMap,
+										   JTOOLS::JPolint2FunctionalMap>::maplist >
+  class JOscProbInterpolator :
+    public JMultiFunction                                                               <JFunction1D_t, JFunctionalMaplist_t>,
+    public JClonable<JOscProbInterpolatorInterface, JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >,
+    public JMessage                                <JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >
+  {
+  public:
+
+    typedef JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t>                interpolator_type;
+    typedef JMultiFunction<JFunction1D_t, JFunctionalMaplist_t>                                    multifunction_type;
+
+    enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS };
+
+    typedef typename multifunction_type::abscissa_type                                                  abscissa_type;    
+    typedef typename multifunction_type::argument_type                                                  argument_type;
+    typedef typename multifunction_type::result_type                                                      result_type;
+    typedef typename multifunction_type::value_type                                                        value_type;
+
+    typedef typename multifunction_type::multimap_type                                                  multimap_type;
+
+    typedef typename multifunction_type::super_const_iterator                                    super_const_iterator;
+    typedef typename multifunction_type::super_iterator                                                super_iterator;
+    typedef typename multifunction_type::function_type                                                  function_type;
+
+    using JMessage<interpolator_type>::debug;
+
+    
+    /**
+     * Default constructor.
+     */
+    JOscProbInterpolator() :
+      multifunction_type(),
+      parameters(),
+      getBaseline()
+    {
+      this->set(JOscParameters(false)); // Initialize buffer with NuFIT NO best fit parameters
+    }
+
+
+    /**
+     * Constructor.
+     *
+     * \param  fileName             oscillation probability table filename
+     */
+    JOscProbInterpolator(const char* fileName) :
+      multifunction_type(),
+      parameters(),
+      getBaseline()
+    {
+      this->load(fileName);
+      this->set(JOscParameters(false)); // Initialize buffer with NuFIT NO best fit parameters      
+    }
+
+
+    /**
+     * Constructor.
+     *
+     * \param  fileName             oscillation probability table filename
+     * \param  parameters           oscillation parameters
+     */
+    JOscProbInterpolator(const char*           fileName,
+			 const JOscParameters& parameters) :
+      multifunction_type(),
+      parameters(),
+      getBaseline()
+    {
+      this->load(fileName);
+      this->set(parameters);
+    }
+
+
+    /**
+     * Load oscillation probability table from file.
+     *
+     * \param  fileName             oscillation probability table filename
+     */
+    void load(const char* fileName) override
+    {
+      using namespace std;
+      using namespace JPP;
+      
+      try {
+
+	NOTICE("loading oscillation probability table from file " << fileName << "... " << flush);
+	
+	JLANG::load<JIO::JFileStreamReader>(fileName, *this);
+
+	NOTICE("OK" << endl);
+      }
+      catch(const JException& error) {
+	THROW(JFileReadException, "JOscProbInterpolator::load(): Error reading file " << fileName);
+      }
+    }    
+
+
+    /**
+     * Get fixed oscillation parameters associated with this interpolation table.
+     *
+     * \return                    oscillation parameters
+     */
+    const JOscParameters& getTableParameters() const override
+    {
+      return parameters;
+    }
+
+
+    /**
+     * Get baseline calculator associated with this interpolation table.
+     *
+     * \return                    baseline calculator
+     */
+    const JBaselineCalculator& getBaselineCalculator() const override
+    {
+      return getBaseline;
+    }
+
+
+    /**
+     * Set oscillation parameters for interpolation.
+     *
+     * \return                      oscillation parameters
+     */
+    void set(JOscParameters parameters) override
+    {
+      using namespace JPP;
+      
+      parameters.join(this->parameters);
+      
+      const JProperties properties = parameters.getProperties();
+
+      for (JProperties::const_iterator i = properties.cbegin(); i != properties.cend(); ++i) {
+	
+	const JOscParameters::JParameter_t& parameter = i->second.getValue<JOscParameters::JParameter_t>();
+	
+	if (parameter.isDefined()) {
+	    
+	  const int index = std::distance(properties.cbegin(), i);
+
+	  this->buffer[index] = parameter.getValue();
+	  
+	} else {
+
+	  THROW(JNoValue,
+		"JOscProbInterpolator<...>::set(JOscParameters): " <<
+		"No value for parameter "                          << i->first);
+	}
+      }
+    }
+
+
+    /**
+     * Get oscillation probability for a given oscillation channel.
+     *
+     * \param  channel            oscillation channel
+     * \param  E                  neutrino energy [GeV]
+     * \param  costh              cosine zenith angle
+     * \return                    oscillation probability
+     */
+    double operator()(const JOscChannel&    channel,
+		      const double          E,
+		      const double          costh) const override
+    {
+      using namespace std;
+      using namespace JPP;
+
+      const JOscChannel* p = find(getOscChannel, getOscChannel + NUMBER_OF_OSCCHANNELS, channel);
+	
+      if (p != end(getOscChannel)) {
+	
+	const double L = getBaseline(costh);
+
+	this->buffer[NUMBER_OF_DIMENSIONS-2] = L/E;
+	this->buffer[NUMBER_OF_DIMENSIONS-1] = costh;
+
+	const argument_type* arguments = this->buffer.data();
+	
+	const size_t       index         = distance(getOscChannel, p);
+	const result_type& probabilities = this->evaluate(arguments);
+
+	return probabilities[index];
+	
+      } else {
+	
+	THROW(JValueOutOfRange, "JOscProbInterpolator<...>::operator(): Invalid oscillation channel " << channel << endl);
+      }
+    }
+
+
+    /**
+     * Get oscillation probability for a given set of oscillation parameters\n
+     * and a given oscillation channel.
+     *
+     * \param  channel            oscillation channel
+     * \param  parameters         oscillation parameters
+     * \param  E                  neutrino energy [GeV]
+     * \param  costh              cosine zenith angle
+     * \return                    oscillation probability
+     */
+    double operator()(const JOscParameters& parameters,
+		      const JOscChannel&    channel,
+		      const double          E,
+		      const double          costh) override
+    {
+      set(parameters);
+
+      return (*this)(channel, E, costh);      
+    }
+    
+    
+    /**
+     * Read from input.
+     *
+     * \param  in                   reader
+     * \return                      reader
+     */
+    JReader& read(JReader& in) override 
+    {
+      parameters.read(in);
+      getBaseline.read(in);
+      
+      in >> static_cast<multifunction_type&>(*this);
+      
+      this->compile();
+      
+      return in;
+    }
+    
+    
+    /**
+     * Write from input.
+     *
+     * \param  out                  writer
+     * \return                      writer
+     */
+    JWriter& write(JWriter& out) const override 
+    {
+      parameters.write(out);
+      getBaseline.write(out);
+      
+      out << static_cast<const multifunction_type&>(*this);
+
+      return out;
+    }
+
+    
+  private:
+
+    JOscParameters      parameters;  //!< Fixed oscillation parameters corresponding to the oscillation probability table
+    JBaselineCalculator getBaseline; //!< Baseline functor
+  };
+}
+
+
+namespace JEEP {
+
+  /**
+   * JMessage template specialization for oscillation probability interpolators.
+   */
+  template<template<class, class> class JCollection_t, class JFunction1D_t, class JFunctionalMaplist_t>
+  struct JMessage<JOSCPROB::JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >
+  {
+    static int debug;
+  };
+
+
+  /**
+   * Default verbosity for oscillation probability interpolators.
+   */
+  template<template<class, class> class JCollection_t, class JFunction1D_t, class JFunctionalMaplist_t>
+  int JMessage<JOSCPROB::JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >::debug = (int) notice_t;
+  
+}
+
+#endif

src/jpp/JOscProb/JOscProbInterpolatorInterface.hh

+#ifndef __JOSCPROB__JOSCPROBINTERPOLATORINTERFACE__
+#define __JOSCPROB__JOSCPROBINTERPOLATORINTERFACE__
+
+
+#include "JLang/JClonable.hh"
+
+#include "JIO/JSerialisable.hh"
+
+#include "JOscProb/JOscChannel.hh"
+#include "JOscProb/JOscParameters.hh"
+#include "JOscProb/JBaselineCalculator.hh"
+
+
+/**
+ * \author bjung, mdejong
+ */
+
+namespace JOSCPROB {
+
+  using JLANG::JClonable;
+
+  using JIO::JSerialisable;
+
+
+  /**
+   * Low-level interface for oscillation probability tables.
+   */
+  class JOscProbInterpolatorInterface :
+    public JSerialisable,
+    public JClonable<JOscProbInterpolatorInterface>
+  {
+  public:
+
+    /**
+     * Default constructor.
+     */
+    JOscProbInterpolatorInterface()
+    {}
+
+
+    /**
+     * Virtual destructor.
+     */
+    virtual ~JOscProbInterpolatorInterface()
+    {}
+    
+
+    /**
+     * Load oscillation probability table from file.
+     *
+     * \param  fileName        oscillation probability table fileName
+     */
+    virtual void load(const char* fileName) = 0;
+    
+    
+    /**
+     * Get oscillation parameters.
+     *
+     * \return                 oscillation parameters
+     */
+    virtual const JOscParameters& getTableParameters() const = 0;
+
+
+    /**
+     * Get baseline calculator associated with this interpolation table.
+     *
+     * \return                    baseline calculator
+     */
+    virtual const JBaselineCalculator& getBaselineCalculator() const = 0;
+    
+
+    /**
+     * Set oscillation parameters.
+     *
+     * \param  parameters      oscillation parameters
+     */
+    virtual void set(JOscParameters parameters) = 0;
+
+
+    /**
+     * Get oscillation probability for a given oscillation channel.
+     *
+     * \param  channel            oscillation channel
+     * \param  E                  neutrino energy [GeV]
+     * \param  costh              cosine zenith angle
+     * \return                    oscillation probability
+     */
+    virtual double operator()(const JOscChannel&    channel,
+			      const double          E,
+			      const double          costh) const = 0;
+    
+    
+    /**
+     * Get oscillation probability for a given set of oscillation parameters\n
+     * and a given oscillation channel.
+     *
+     * \param  channel            oscillation channel
+     * \param  parameters         oscillation parameters
+     * \param  E                  neutrino energy [GeV]
+     * \param  costh              cosine zenith angle
+     * \return                    oscillation probability
+     */
+    virtual double operator()(const JOscParameters& parameters,
+			      const JOscChannel&    channel,
+			      const double          E,
+			      const double          costh)
+    {
+      set(parameters);
+
+      return (*this)(channel, E, costh);      
+    }
+  };
+}
+
+#endif

src/jpp/JOscProb/JOscProbToolkit.hh

+#ifndef __JOSCPROB__JOSCPROBTOOLKIT__
+#define __JOSCPROB__JOSCPROBTOOLKIT__
+
+#include <string>
+
+#include "JLang/JException.hh"
+
+#include "JOscProb/JOscChannel.hh"
+
+
+/**
+ * \author bjung
+ * Auxiliary methods for oscillation probabilities.
+ */
+
+namespace JOSCPROB {}
+namespace JPP { using namespace JOSCPROB; }
+
+namespace JOSCPROB {
+
+  using JLANG::JValueOutOfRange;
+  
+
+  /**
+   * OscProb neutrino flavour identifiers.
+   */
+  enum class OscProbFlavour_t { ELECTRON,
+				MUON,
+				TAU };
+  
+  
+  /**
+   * Auxiliary function for retrieving the OscProb flavour identifier corresponding to a JOscProb flavour identifier.
+   *
+   * \param  flavour             flavour identifier
+   * \return                     OscProb flavour identifier
+   */
+  inline OscProbFlavour_t getOscProbFlavour(const JFlavour_t flavour)
+  {
+    switch(flavour) {
+    case JFlavour_t::ELECTRON:
+      return OscProbFlavour_t::ELECTRON;
+    case JFlavour_t::MUON:
+      return OscProbFlavour_t::MUON;
+    case JFlavour_t::TAU:
+      return OscProbFlavour_t::TAU;
+    default:
+      THROW(JLANG::JValueOutOfRange, "getOscProbFlavour(...): Invalid flavour " << (int) flavour);
+    }
+  }
+
+
+  /**
+   * Auxiliary function for retrieving the OscProb flavour identifier corresponding to a JOscProb flavour identifier.
+   *
+   * \param  flavour             flavour identifier
+   * \return                     OscProb flavour identifier
+   */
+  inline OscProbFlavour_t getOscProbFlavour(const int pdgType)
+  {
+    JFlavour_t flavour = getFlavour(pdgType);
+    
+    return getOscProbFlavour(flavour);
+  }
+
+
+  /**
+   * Auxiliary data structure to hold oscillation variable names.
+   */
+  struct JOscVars
+  {
+    /**
+     * Oscillation variable types.
+     */
+    enum type {
+      COSTH,
+      SINTH,
+      ENERGY,
+      LOG10E,
+      LOE,
+      BASELINE,
+      UNDEFINED
+    };
+
+    
+    static const char* const energy() { return "energy"; }   //!< energy [GeV]
+    static const char* const log10E() { return "log10E"; }   //!< logarithmic energy [GeV]
+    static const char* const LoE()    { return "LoE"; }      //!< L/E [km GeV-1]
+    
+    static const char* const costh()  { return "costh"; }    //!< cosine of zenith-angle
+    static const char* const sinth()  { return "sinth"; }    //!< sine of zenith-angle
+    static const char* const L()      { return "L"; }        //!< sine of zenith-angle
+
+    
+    /**
+     * Get oscillation variable type.
+     *
+     * \param  name                              oscillation variable name
+     * \return                                   oscillation variable type
+     */
+    static inline type getType(const std::string& name)
+    {
+      if        (name == energy()) {
+	return ENERGY;
+      } else if (name == costh()) {
+	return COSTH;
+      } else if (name == sinth()) {
+	return SINTH;
+      } else if (name == log10E()) {
+	return LOG10E;	
+      } else if (name == LoE()) {
+	return LOE;
+      } else if (name == L()) {
+	return BASELINE;	
+      } else {
+	THROW(JValueOutOfRange, "JOscVars::getType(): Invalid oscillation variable " << name);
+      }
+    }
+  };
+}
+
+#endif

src/jpp/JPhysics/JConstants.hh

+#ifndef __JPHYSICS__JCONSTANTS__
+#define __JPHYSICS__JCONSTANTS__
+
+#include <math.h>
+
+#include "JMath/JConstants.hh"
+
+
+/**
+ * \file
+ * Physics constants. 
+ * \author mdejong
+ */
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  using JMATH::PI;
+  using JMATH::EULER;
+
+  /**
+   * Physics constants.
+   */
+  static const double C                         = 0.299792458;          //!< Speed of light in vacuum [m/ns]
+  static const double C_INVERSE                 = 1.0/C;                //!< Inverse speed of light in vacuum [ns/m]
+  static const double AVOGADRO                  = 6.0221415e23;         //!< Avogadro's number [gr^-1]
+  static const double NUCLEON_MOLAR_MASS        = 1.0;                  //!< nucleon molar mass [g/mol]
+  static const double H                         = 4.13566733e-15;       //!< Planck constant [eV s]
+  static const double HBAR                      = H/(2*PI);             //!< Planck constant [eV s]
+  static const double HBARC                     = HBAR*C*1.0e9;         //!< Planck constant [eV m]
+  static const double ALPHA_ELECTRO_MAGNETIC    = 1.0/137.036;          //!< Electro-Magnetic coupling constant
+  static const double THETA_MCS                 = 13.6e-3;              //!< Multiple Coulomb scattering constant [GeV]
+
+  /**
+   * Fixed environment values.
+   */
+  static const double DENSITY_SEA_WATER         = 1.038;                //!< Density  of sea water [g/cm^3]
+  static const double DENSITY_ROCK              = 2.65;                 //!< Density  of rock      [g/cm^3]
+  static const double SALINITY_SEA_WATER        = 0.035;                //!< Salinity of sea water
+  static const double INDEX_OF_REFRACTION_WATER = 1.3800851282;         //!< Average index of refraction of water corresponding to the group velocity 
+  static const double X0_WATER_M                = 0.36;                 //!< Radiation length pure water [m]
+
+  /**
+   * Derived quantities of optical medium.
+   */
+  static const double TAN_THETA_C_WATER         = sqrt((INDEX_OF_REFRACTION_WATER - 1.0) * (INDEX_OF_REFRACTION_WATER + 1.0)); //!< Average tangent corresponding to the group velocity 
+  static const double COS_THETA_C_WATER         = 1.0 / INDEX_OF_REFRACTION_WATER;                                             //!< Average cosine  corresponding to the group velocity 
+  static const double SIN_THETA_C_WATER         = TAN_THETA_C_WATER * COS_THETA_C_WATER;                                       //!< Average sine    corresponding to the group velocity 
+  static const double KAPPA_WATER               = 0.96;                                                                        //!< Average R-dependence of arrival time of Cherenkov light
+
+  /**
+   * Particle masses.
+   * Note that the neutrino masses are set to zero.
+   */
+  static const double MASS_PHOTON               = 0.0;                 //!< photon     mass [GeV]
+
+  static const double MASS_ELECTRON_NEUTRINO    = 0.0;                 //!< electron neutrino mass [GeV]
+  static const double MASS_MUON_NEUTRINO        = 0.0;                 //!< muon     neutrino mass [GeV]
+  static const double MASS_TAU_NEUTRINO         = 0.0;                 //!< tau      neutrino mass [GeV]
+  
+  static const double MASS_ELECTRON             = 0.510998946e-3;      //!< electron       mass [GeV]
+  static const double MASS_MUON                 = 0.1056583745;        //!< muon           mass [GeV]
+  static const double MASS_TAU                  = 1.77682;             //!< tau            mass [GeV]
+
+  static const double MASS_NEUTRAL_PION         = 0.1349766;           //!< pi_0           mass [GeV]
+  static const double MASS_CHARGED_PION         = 0.13957018;          //!< pi^+/-         mass [GeV]
+  static const double MASS_NEUTRAL_KAON         = 0.497614;            //!< K_0            mass [GeV]
+  static const double MASS_CHARGED_KAON         = 0.493677;            //!< K^+/-          mass [GeV]
+  static const double MASS_NEUTRAL_RHO          = 0.77526;             //!< rho_0          mass [GeV]
+  static const double MASS_CHARGED_RHO          = 0.77511;             //!< rho^+/-        mass [GeV]
+  static const double MASS_NEUTRAL_D            = 1.86483;             //!< D_0            mass [GeV]
+  static const double MASS_CHARGED_D            = 1.86965;             //!< D^+/-          mass [GeV]
+  static const double MASS_CHARGED_D_S          = 1.96834;             //!< D_s^+/-        mass [GeV]
+
+  static const double MASS_PROTON               = 0.9382720813;        //!< proton         mass [GeV]
+  static const double MASS_NEUTRON              = 0.9395654133;        //!< neutron        mass [GeV]
+  static const double MASS_DELTA_1232           = 1.232;               //!< Delta (1232)   mass [GeV]
+  static const double MASS_LAMBDA               = 1.115683;            //!< Lambda         mass [GeV]
+  static const double MASS_NEUTRAL_SIGMA        = 1.192642;            //!< Sigma_0        mass [GeV]
+  static const double MASS_CHARGED_SIGMA        = 1.18937;             //!< Sigma^+/-      mass [GeV]
+  static const double MASS_NEUTRAL_XI           = 1.31486;             //!< Xi_0           mass [GeV]
+  static const double MASS_CHARGED_XI           = 1.32171;             //!< Xi^+/-         mass [GeV]
+  static const double MASS_CHARGED_OMEGA        = 1.67245;             //!< Omega^+/-      mass [GeV]
+  static const double MASS_CHARGED_LAMBDA_C     = 2.28646;             //!< Lambda_c^+/-   mass [GeV]
+  static const double MASS_DOUBLYCHARGED_SIGMA_C = 2.45397;            //!< Sigma_c^++/--  mass [GeV]  
+  static const double MASS_CHARGED_SIGMA_C      = 2.4529;              //!< Sigma_c^+/-    mass [GeV]
+  static const double MASS_NEUTRAL_SIGMA_C      = 2.45375;             //!< Sigma_c_0      mass [GeV]
+  static const double MASS_CHARGED_XI_C         = 2.46793;             //!< Xi_c^+/-       mass [GeV]
+  static const double MASS_NEUTRAL_XI_C         = 2.47091;             //!< Xi_c_0         mass [GeV]
+  static const double MASS_NEUTRAL_OMEGA_C      = 2.6952;              //!< Omega_c_0      mass [GeV]
+  static const double MASS_NEUTRAL_B            = 5.27958;             //!< B_0            mass [GeV]
+  static const double MASS_CHARGED_B            = 5.27926;             //!< B^+/-          mass [GeV]
+  static const double MASS_NEUTRAL_B_S          = 5.36677;             //!< B_s^0          mass [GeV]
+  static const double MASS_NEUTRAL_LAMBDA_B     = 5.6194;              //!< Lambda_b^0     mass [GeV]
+  static const double MASS_NEUTRAL_XI_B         = 5.7878;              //!< Xi_b^0         mass [GeV]
+  static const double MASS_CHARGED_XI_B         = 5.7911;              //!< Xi_b^+/-       mass [GeV]
+  static const double MASS_CHARGED_OMEGA_B      = 6.071;               //!< Omega_b^+/-    mass [GeV]
+  static const double MASS_CHARGED_B_C          = 6.2756;              //!< B_c^+/-        mass [GeV]
+  
+  /** 
+   * Get speed of light.
+   * 
+   * return           speed of light [m/ns]
+   */
+  inline const double getSpeedOfLight()
+  {
+    return C;
+  }
+
+
+  /**
+   * Get inverse speed of light.
+   *
+   * return           inverse speed of light [ns/m]
+   */
+  inline const double getInverseSpeedOfLight()
+  { 
+    return C_INVERSE; 
+  }
+
+
+  /**
+   * Get average index of refraction of water corresponding to group velocity.
+   *
+   * \return          index of refraction
+   */
+  inline double getIndexOfRefraction() 
+  {
+    return INDEX_OF_REFRACTION_WATER;
+  }
+
+
+  /**
+   * Get average index of refraction of water corresponding to phase velocity.
+   *
+   * \return          index of refraction
+   */
+  inline double getIndexOfRefractionPhase() 
+  {
+    return 1.35;
+  }
+
+
+  /**
+   * Get average tangent of Cherenkov angle of water corresponding to group velocity.
+   *
+   * \return          tan(theta_C)
+   */
+  inline double getTanThetaC() 
+  {
+    return TAN_THETA_C_WATER;
+  }
+
+
+  /**
+   * Get average cosine of Cherenkov angle of water corresponding to group velocity.
+   *
+   * \return          cos(theta_C)
+   */
+  inline double getCosThetaC() 
+  {
+    return COS_THETA_C_WATER;
+  }
+
+
+  /**
+   * Get average sine of Cherenkov angle of water corresponding to group velocity.
+   *
+   * \return          sin(theta_C)
+   */
+  inline double getSinThetaC() 
+  {
+    return SIN_THETA_C_WATER;
+  }
+
+
+  /**
+   * Get average R-dependence of arrival time of Cherenkov light (a.k.a. kappa).
+   *
+   * \return          kappa
+   */
+  inline double getKappaC() 
+  {
+    return KAPPA_WATER;
+  }
+}
+
+#endif

src/jpp/JPhysics/JGeane.hh

+#ifndef __JPHYSICS__JGEANE__
+#define __JPHYSICS__JGEANE__
+
+#include <cmath>
+#include <map>
+
+#include "JPhysics/JConstants.hh"
+
+
+/**
+ * \file
+ * Energy loss of muon.
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  /**
+   * Equivalent muon track length per unit shower energy.
+   *
+   * See ANTARES internal note ANTARES-SOFT-2002-015, J.\ Brunner.
+   *
+   * \return        equivalent muon track length [m/Gev]
+   */
+  inline double geanc()
+  {
+    return 4.7319;                   // dx/dE [m/GeV]
+  }
+
+
+  /**
+   * Interface for muon energy loss.
+   *
+   * This interface provides for the various function object operators.
+   */
+  class JGeane {
+  public:
+   /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to ionisation [GeV/m]
+     */
+    virtual double getA() const = 0;
+
+
+    /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to pair production and bremsstrahlung [m^-1]
+     */
+    virtual double getB() const = 0;
+
+
+    /**
+     * Get energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    virtual double getE(const double E, const double dx) const = 0;
+
+
+    /**
+     * Get distance traveled by muon.
+     *
+     * \param  E0      Energy of muon at start [GeV]
+     * \param  E1      Energy of muon at end   [GeV]
+     * \return         distance traveled [m]
+     */
+    virtual double getX(const double E0, 
+			const double E1) const = 0;
+
+
+    /**
+     * Energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    double operator()(const double E, const double dx) const
+    {
+      return this->getE(E, dx);
+    }
+    
+
+    /**
+     * Range of muon.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \return         range [m]
+     */
+    double operator()(const double E) const
+    {
+      return this->getX(E, 0.0);
+    }
+
+
+    /**
+     * Equivalent unit track length per unit shower energy and per unit track length.
+     *
+     * \return        equivalent unit track length [Gev^-1]
+     */
+    double operator()() const
+    {
+      return this->getB() * geanc();
+    }
+  };
+
+
+  /**
+   * Function object for the energy loss of the muon.\n
+   * The energy loss can be formulated as:
+   *
+   *     \f[ -\frac{dE}{dx}  =  a + bE \f]
+   *
+   * N.B:
+   * \f$ a \f$ and \f$ b \f$ are assumed constant (internal units m and GeV, respectively).
+   */
+  class JGeane_t :
+    public JGeane
+  {
+  public:
+    /**
+     * constructor
+     * \param  __a     Energy loss due to ionisation [GeV/m]
+     * \param  __b     Energy loss due to pair production and bremsstrahlung [m^-1]
+     */ 
+    JGeane_t(const double __a,
+	     const double __b) :
+      a(__a),
+      b(__b)
+    {}
+    
+
+    /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to ionisation [GeV/m]
+     */
+    virtual double getA() const override 
+    {
+      return a;
+    }
+
+
+    /**
+     * Get energy loss constant.
+     *
+     * \return         Energy loss due to pair production and bremsstrahlung [m^-1]
+     */
+    virtual double getB() const override 
+    {
+      return b;
+    }
+
+
+    /**
+     * Get energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    virtual double getE(const double E, const double dx) const override 
+    {
+      const double y = (a/b + E) * exp(-b*dx)  -  a/b;
+
+      if (y > 0.0)
+	return y;
+      else
+	return 0.0;
+    }
+
+
+    /**
+     * Get distance traveled by muon.
+     *
+     * \param  E0      Energy of muon at start [GeV]
+     * \param  E1      Energy of muon at end   [GeV]
+     * \return         distance traveled [m]
+     */
+    virtual double getX(const double E0, 
+			const double E1) const override
+    {
+      return -log((a + b*E1) / (a+b*E0)) / b;
+    }
+
+  protected:
+    const double a;
+    const double b;
+  };
+
+
+  /**
+   * Function object for energy dependent energy loss of the muon.
+   *
+   * Approximate values of energy loss parameters taken from reference:
+   * Proceedings of ICRC 2001, "Precise parametrizations of muon energy losses in water",
+   * S. Klimushin, E. Bugaev and I. Sokalski.
+   */
+  class JGeaneWater : 
+    public JGeane,
+    protected std::map<double, JGeane_t>
+  {    
+  public:    
+    /**
+     * Default constructor.
+     */
+    JGeaneWater()
+    {
+      using namespace std;
+
+      this->insert(make_pair( 0.0e0,  JGeane_t( 2.30e-1 * DENSITY_SEA_WATER, 15.50e-4 * DENSITY_SEA_WATER)));
+      this->insert(make_pair(30.0e0,  JGeane_t( 2.67e-1 * DENSITY_SEA_WATER,  3.40e-4 * DENSITY_SEA_WATER)));
+      this->insert(make_pair(35.3e3,  JGeane_t(-6.50e-1 * DENSITY_SEA_WATER,  3.66e-4 * DENSITY_SEA_WATER)));
+    }
+    
+
+    /**
+     * Get energy loss constant.
+     *
+     * N.B. The return value corresponds to the low-energy regime.
+     *
+     * \return         Energy loss due to ionisation [GeV/m]
+     */
+    virtual double getA() const override 
+    {
+      return 2.30e-1 * DENSITY_SEA_WATER; //This is the value for low energy (<30 GeV), the value used for step(ds) and getRange())
+    }
+
+
+    /**
+     * Get energy loss constant.
+     *
+     * N.B. The return value corresponds to the medium-energy regime.
+     *
+     * \return         Energy loss due to pair production and bremsstrahlung [m^-1]
+     */
+    virtual double getB() const override 
+    {
+      return 3.40e-4 * DENSITY_SEA_WATER;
+    }
+
+
+    /**
+     * Get energy of muon after specified distance.
+     *
+     * \param  E       Energy of muon [GeV]
+     * \param  dx      distance traveled [m]
+     * \return         Energy of muon [GeV]
+     */
+    virtual double getE(const double E, const double dx) const override 
+    {
+      double E1 = E;
+      double x1 = dx;
+
+      if (E1 > MASS_MUON / getSinThetaC()) {
+
+	const_iterator p = this->lower_bound(E1);
+
+	do {
+
+	  --p;
+
+	  const double x2 = p->second.getX(E1, p->first);
+
+	  if (x2 > x1) {
+	    return p->second.getE(E1, x1);
+	  }
+
+	  E1  = p->first;
+	  x1 -= x2;
+
+	} while (p != this->begin());
+      }
+
+      return E1;
+    }
+
+
+
+
+    /**
+     * Get energy loss due to ionisation.
+     *
+     * \param  E           initial energy                [GeV]
+     * \param  dx          distance traveled             [m]
+     * \return             energy loss due to ionisation [GeV]
+     */
+    double getEa(const double E, const double dx) const
+    {
+      using namespace std;
+      using namespace JPP;
+      
+      double Ea = 0.0;
+      
+      double E1 = E;
+      double x1 = dx;
+
+      if (E1 > MASS_MUON / getSinThetaC()) {
+
+	map<double, JGeane_t>::const_iterator p = this->lower_bound(E1);
+
+	do {
+
+	  --p;
+
+	  const double x2 = p->second.getX(E1, p->first);
+
+	  Ea += (x2 > x1 ? x1 : x2) * p->second.getA();
+	  E1  = p->first;
+	  
+	  x1 -= x2;
+	  
+	} while (p != this->cbegin() && x1 > 0.0);
+      }
+
+      return Ea;
+    }
+
+
+    /**
+     * Get energy loss due to pair production and bremsstrahlung.
+     *
+     * \param  E           initial energy                                        [GeV]
+     * \param  dx          distance traveled                                     [m]
+     * \return             energy loss due to pair production and bremsstrahlung [GeV]
+     */
+    double getEb(const double E, const double dx) const
+    {
+      const double dE = E - getE(E, dx);
+      
+      return dE - getEa(E, dx);
+    }
+
+
+    /**
+     * Get distance traveled by muon.
+     *
+     * \param  E0      Energy of muon at start [GeV]
+     * \param  E1      Energy of muon at end   [GeV]
+     * \return         distance traveled       [m]
+     */
+    virtual double getX(const double E0, 
+			const double E1) const override
+    {
+      double E  = E0;
+      double dx = 0.0;
+
+      if (E > MASS_MUON / getSinThetaC()) {
+
+	const_iterator p = this->lower_bound(E);
+
+	do {
+
+	  --p;
+
+	  if (E1 > p->first) {
+	    return dx + p->second.getX(E, E1);
+	  }
+
+	  dx += p->second.getX(E, p->first);
+	  E   = p->first;
+
+	} while (p != this->begin());
+      }
+
+      return dx;
+    }
+  };
+
+
+  /**
+   * Function object for energy loss of muon in sea water.
+   */
+  static const JGeaneWater gWater;
+
+
+  /**
+   * Function object for energy loss of muon in rock.
+   */
+  static const JGeane_t    gRock(2.67e-1 * 0.9 * DENSITY_ROCK, 
+				 3.40e-4 * 1.2 * DENSITY_ROCK);
+
+}
+
+#endif

src/jpp/JPhysics/JGeant_t.hh

+#ifndef __JPHYSICS__JGEANT_T__
+#define __JPHYSICS__JGEANT_T__
+
+#include "JTools/JFunction1D_t.hh"
+#include "JIO/JSerialisable.hh"
+
+/**
+ * \file
+ * Base class for photon emission profile EM-shower.
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  using JIO::JReader;
+  using JIO::JWriter;
+
+  typedef JTOOLS::JGridPolint1Function1D_t  JGeantFunction1D_t;
+
+
+  /**
+   * Base class for the probability density function of photon emission from EM-shower 
+   * as a function of the index of refraction and the cosine of the emission angle.
+   *
+   * The implementation of this function is based on a linear interpolation of tabulated values.
+   * In this, a linear approximation of the dependence of the normalisation constant on 
+   * the index of refraction is assumed. This assumption is valid to within 10^-3.
+   */
+  class JGeant_t : 
+    public JGeantFunction1D_t
+  {
+  public:
+    /**
+     * Default constructor.
+     */
+    JGeant_t() 
+    {}
+ 
+
+    /**
+     * Number of photons from EM-shower as a function of emission angle.
+     * The integral over full solid angle is normalised to one.
+     *
+     * \param  n       index of refraction
+     * \param  ct      cosine angle of emmision
+     * \return         d^2P/dcos()dphi
+     */
+    double operator()(const double n,
+		      const double ct) const
+    {
+      const double y = JGeantFunction1D_t::operator()(ct - 1.0/n);
+
+      return y * (a0 - a1*n);
+    }
+
+
+    /**
+     * Integral number of photons from EM-shower between two emission angles.
+     * The integral over full solid angle is normalised to one.
+     *
+     * \param  n       index of refraction
+     * \param  xmin    minimal cosine angle of emmision
+     * \param  xmax    maximal cosine angle of emmision
+     * \return         dnpe/dphi 
+     */
+    double operator()(const double n, 
+		      const double xmin, 
+		      const double xmax) const
+    {
+      const double x_min = std::max(xmin - 1.0/n, buffer. begin()->getX());
+      const double x_max = std::min(xmax - 1.0/n, buffer.rbegin()->getX());
+
+      const double y = buffer(x_max) - buffer(x_min);
+
+      return y * (a0 - a1*n);
+    }
+
+
+    /**
+     * Read geant from input.
+     *
+     * \param  in       reader
+     * \param  geant    geant
+     * \return          reader
+     */
+    friend inline JReader& operator>>(JReader& in, JGeant_t& geant)
+    {
+      in >> geant.a0;
+      in >> geant.a1;
+      in >> static_cast<JGeantFunction1D_t&>(geant);
+
+      geant.compile();
+
+      return in;
+    }
+
+
+    /**
+     * Write geant to output.
+     *
+     * \param  out      writer
+     * \param  geant    geant
+     * \return          writer
+     */
+    friend inline JWriter& operator<<(JWriter& out, const JGeant_t& geant)
+    {
+      out << geant.a0;
+      out << geant.a1;
+      out << static_cast<const JGeantFunction1D_t&>(geant);
+
+      return out;
+    }
+
+
+  protected:
+    /**
+     * Function compilation.
+     */
+    virtual void do_compile() override 
+    {
+      JGeantFunction1D_t::do_compile();
+
+      buffer.clear();
+
+      JTOOLS::integrate(*this, buffer);
+
+      buffer.compile();
+    }
+ 
+
+    double a0;     //!< offset of the normalisation dependence
+    double a1;     //!< slope  of the normalisation dependence
+    JGeantFunction1D_t buffer;
+  };
+}
+
+#endif

src/jpp/JPhysics/JGeanz.hh

+#ifndef __JPHYSICS__JGEANZ__
+#define __JPHYSICS__JGEANZ__
+
+#include <cmath>
+
+#include "JMath/JMathSupportkit.hh"
+
+/**
+ * \file
+ * Longitudinal emission profile EM-shower.
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  /**
+   * Function object for longitudinal profile of EM-shower.
+   *
+   *
+   *         \f[P(z)  \propto  z^{a-1} \times e^{-z/b}\f]
+   *
+   * where: \f[a = a_{0} + a_{1} \times \ln(E)\f]
+   *
+   * The parametrisation is taken from reference:
+   * C. Kopper, "Performance Studies for the KM3NeT Neutrino Telescope.",
+   * PhD thesis, University of Erlangen.
+   */
+  class JGeanz {
+  public:
+    /**
+     * constructor
+     * \param  __a0    power term (constant)
+     * \param  __a1    power term (E dependence)
+     * \param  __b     expontial slope
+     */ 
+    JGeanz(const double __a0,
+	   const double __a1,
+	   const double __b) :
+      a0(__a0),
+      a1(__a1),
+      b (__b),
+      Emin(exp(-a0/a1))
+    {}
+
+
+    /**
+     * Probability Density Function
+     *
+     * \param  E       EM-shower energy [GeV]
+     * \param  z       z position of light emission point relative to vertex location (z >= 0) [m]
+     * \return         dP/dz
+     */
+    double getProbability(const double E, 
+			  const double z) const
+    {
+      if (E > Emin) {
+
+	const double a = a0 + a1 * log(E);
+	const double y = pow(z,a-1.0) * exp(-z/b) / (pow(b,a) * std::tgamma(a));
+	
+	return y;
+      } 
+	  
+      if (z <= getMinimalShowerSize()) 
+	return 1.0 / getMinimalShowerSize();
+      else
+	return 0.0;
+    }
+    
+
+    /**
+     * Probability Density Function
+     *
+     * \param  E       EM-shower energy [GeV]
+     * \param  z       z position of light emission point relative to vertex location (z >= 0) [m]
+     * \return         dP/dz
+     */
+    double operator()(const double E, 
+		      const double z) const
+    {
+      return getProbability(E, z);
+    }
+   
+
+    /**
+     * Integral of PDF (starting from 0).
+     *
+     * \param  E       EM-shower energy [GeV]
+     * \param  z       z position [m] (>= 0)
+     * \return         dP
+     */
+    double getIntegral(const double E, 
+		       const double z) const
+    {
+      if (E > Emin) {
+
+	const double a = a0 + a1 * log(E);
+	const double x = z / b;
+	const double y = JMATH::Gamma(a,x);
+	
+	return y;
+      }
+      
+      if (z <= getMinimalShowerSize())
+	return z / getMinimalShowerSize();
+      else
+	return 1.0;
+    }
+    
+
+    /**
+     * Get shower length for a given integrated probability.
+     *
+     * \param  E       EM-shower energy [GeV]
+     * \param  P       integrated probability [0,1]
+     * \param  eps     relative precision
+     * \return         shower length [m]
+     */
+    double getLength(const double E, 
+		     const double P,
+		     const double eps = 1.0e-3) const
+    {
+      double zmin  =   0.0;     // [m]
+      double zmax  =  30.0;     // [m]
+
+      if (E > Emin) {
+
+	const double Q = P * (1.0 - eps);
+	
+	for (int i = 100; i != 0; --i) {
+	  
+	  const double z = 0.5 * (zmin + zmax);
+	  const double p = getIntegral(E, z);
+	  
+	  if (fabs(p-Q) < p*eps) {
+	    return z;
+	  }
+	  
+	  if (p > P)
+	    zmax = z;
+	  else
+	    zmin = z;
+	}
+	
+	return 0.5 * (zmin + zmax);
+
+      } else
+
+	return 0.0;
+    }
+
+    /**
+     * Get depth of shower maximum
+     *
+     * \param  E       EM-shower energy[GeV]
+     * \return         depth of maximum [m]
+     */
+    
+    double getMaximum(const double E) const
+    {
+      const double a = a0 + a1 * log(E); 
+
+      return (a-1)*b;
+    }
+
+
+    /**
+     * Get minimal shower size.
+     *
+     * \return         size [m]
+     */
+    static double getMinimalShowerSize()
+    {
+      return 1e-6;
+    }
+
+  protected:
+    const double a0;
+    const double a1;
+    const double b;
+    const double Emin;
+  };
+
+
+  /**
+   * Function object for longitudinal EM-shower profile
+   */
+  static const JGeanz geanz(1.85, 0.62, 0.54);
+}
+
+#endif

src/jpp/JPhysics/JNPETable.hh

+#ifndef __JPHYSICS__JNPETABLE__
+#define __JPHYSICS__JNPETABLE__
+
+#include "JLang/JSharedPointer.hh"
+#include "JLang/JException.hh"
+
+#include "JTools/JConstantFunction1D.hh"
+#include "JTools/JMultiFunction.hh"
+#include "JTools/JTransformableMultiFunction.hh"
+#include "JTools/JToolsToolkit.hh"
+
+
+/**
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  using JTOOLS::JConstantFunction1D;
+  using JTOOLS::JMapList;
+  using JTOOLS::JMultiFunction;
+  using JTOOLS::JMultiMapTransformer;
+  using JTOOLS::JTransformableMultiFunction;
+
+
+  /**
+   * Custom class for integrated values of the PDF of the arrival time of Cherenkov light.
+   *
+   * This class provides for the number of photo-electrons as a function
+   * of the leading <tt>(n - 1)</tt> parameter values. 
+   */
+  template<class JArgument_t, 
+	   class JResult_t,
+	   class JMaplist_t,
+	   class JDistance_t = JTOOLS::JDistance<JArgument_t> >
+  class JNPETable :
+    public JMultiFunction<JConstantFunction1D<JArgument_t, JResult_t>,
+			  JMaplist_t, 
+			  JDistance_t>
+  {
+  public:
+
+    typedef JMultiFunction<JConstantFunction1D<JArgument_t, JResult_t>,
+			   JMaplist_t, 
+			   JDistance_t>                                   multifunction_t;
+
+    using multifunction_t::NUMBER_OF_DIMENSIONS;
+
+    typedef JConstantFunction1D<JArgument_t, JResult_t>                   function_type;
+    typedef typename multifunction_t::map_type                            map_type;
+
+    typedef typename multifunction_t::value_type                          value_type;
+    typedef typename multifunction_t::argument_type                       argument_type;
+    typedef typename multifunction_t::supervisor_type                     supervisor_type;
+
+    typedef typename multifunction_t::abscissa_type                       abscissa_type;
+    typedef typename multifunction_t::ordinate_type                       ordinate_type;
+    typedef typename multifunction_t::result_type                         result_type;
+
+    typedef typename multifunction_t::const_iterator                      const_iterator;
+    typedef typename multifunction_t::const_reverse_iterator              const_reverse_iterator;
+    typedef typename multifunction_t::iterator                            iterator;
+    typedef typename multifunction_t::reverse_iterator                    reverse_iterator;
+
+    typedef typename multifunction_t::super_iterator                      super_iterator;
+    typedef typename multifunction_t::super_const_iterator                super_const_iterator;
+
+    typedef JMultiMapTransformer<NUMBER_OF_DIMENSIONS, argument_type>     transformer_type;
+
+
+    /**
+     * Default constructor.
+     */
+    JNPETable() :
+      transformer(transformer_type::getClone())
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  input               multi-dimensional PDF
+     */
+    template<class JPDF_t, class JPDFMaplist_t, class JPDFDistance_t>
+    JNPETable(const JTransformableMultiFunction<JPDF_t, JPDFMaplist_t, JPDFDistance_t>& input) :
+      transformer(transformer_type::getClone())
+    {
+      using namespace JTOOLS;
+
+      typedef JTransformableMultiFunction<JPDF_t, JPDFMaplist_t, JPDFDistance_t>                       JTransformableMultiFunction_t;
+      typedef JMultiKey<JTransformableMultiFunction_t::NUMBER_OF_DIMENSIONS - 1, argument_type>        JMultiKey_t;
+      typedef typename JTransformableMultiFunction_t::transformer_type                                 transformer_type;
+
+
+      this->transformer.reset(input.transformer->clone());
+
+      for (typename JTransformableMultiFunction_t::super_const_iterator i = input.super_begin(); i != input.super_end(); ++i) {
+
+	const JMultiKey_t& key   = (*i).getKey();
+	const JPDF_t&      value = (*i).getValue();
+
+        const typename transformer_type::array_type array(key);
+
+        const double        V = getIntegral(value);
+        const argument_type z = input.transformer->getXn(array, 1.0) - input.transformer->getXn(array, 0.0);
+
+        this->insert(key, function_type(z*V));
+      }
+
+      this->compile();
+    }
+
+
+    /**
+     * Add NPE table.
+     *
+     * Note that the summation is made via iteration of the elements in this multidimensional table.
+     *
+     * \param  input               NPE table
+     */
+    void add(const JNPETable& input)
+    {
+      using namespace JTOOLS;
+
+      for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
+
+        map_type& f1 = (*i).getValue();
+
+	for (typename map_type::iterator j = f1.begin(); j != f1.end(); ++j) {
+
+	  try {
+
+	    const JArray<NUMBER_OF_DIMENSIONS, argument_type> buffer((*i).getKey(), j->getX());
+
+	    const double npe = get_value(input.evaluate(buffer.data()));
+	    const double W   = this->transformer->getWeight(buffer);
+	  
+	    j->getY() += npe/W;
+	  }
+	  catch(JLANG::JException& error) {}
+	}
+      }
+    }
+
+
+    /**
+     * Get number of photo-electrons.
+     *
+     * \param  args               comma separated argument list
+     * \return                    number of photo-electrons
+     */
+    template<class ...Args>
+    result_type operator()(const Args& ...args) const
+    {
+      this->buffer.set(args...);
+
+      return this->evaluate(this->buffer.data());
+    }
+
+
+    /**
+     * Recursive function value evaluation.
+     *
+     * \param  pX              pointer to abscissa values
+     * \return                 function value
+     */
+    virtual result_type evaluate(const argument_type* pX) const override 
+    {
+      for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
+        this->buffer[i] = pX[i];
+      }
+
+      const double      W   = transformer->getWeight(buffer);
+      const result_type npe = multifunction_t::evaluate(buffer.data());
+
+      return W * npe;
+    }
+
+
+    /**
+     * Application of weight function.
+     *
+     * \param  transformer     function transformer
+     */
+    void transform(const transformer_type& transformer)
+    {
+      using namespace JTOOLS;
+
+      for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
+
+        map_type& f1 = (*i).getValue();
+
+	for (typename map_type::iterator j = f1.begin(); j != f1.end(); ++j) {
+
+	  const JArray<NUMBER_OF_DIMENSIONS, argument_type> array((*i).getKey(), j->getX());
+
+	  j->getY() *= this->transformer->getWeight(array) / transformer.getWeight(array);
+	}
+      }
+
+      this->transformer.reset(transformer.clone());
+      this->compile();
+    }
+
+
+    JLANG::JSharedPointer<transformer_type> transformer;
+
+  protected:
+    mutable JTOOLS::JArray<NUMBER_OF_DIMENSIONS, argument_type> buffer;
+  };
+}
+
+#endif

src/jpp/JPhysics/JNPE_t.hh

+
+#include "JLang/JException.hh"
+#include "JTools/JCollection.hh"
+#include "JTools/JMap.hh"
+#include "JTools/JGridMap.hh"
+#include "JTools/JMapList.hh"
+#include "JTools/JSpline.hh"
+#include "JTools/JPolint.hh"
+#include "JTools/JElement.hh"
+#include "JPhysics/JNPETable.hh"
+#include "JPhysics/JPDFTable.hh"
+#include "JPhysics/JPDFToolkit.hh"
+#include "JPhysics/JPDFTypes.hh"
+#include "JPhysics/JGeanz.hh"
+#include "JMath/JZero.hh"
+
+
+/**
+ * \file
+ *
+ * Auxiliary data structure for muon PDF.
+ * \author mdejong
+ */
+struct JMuonNPE_t {
+
+  typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
+			JPP::JPolint1FunctionalGridMap,
+			JPP::JPolint1FunctionalGridMap>::maplist   JNPEMaplist_t;
+  typedef JPP::JNPETable<double, double, JNPEMaplist_t>            JNPE_t;
+
+
+  /**
+   * Constructor.
+   *
+   * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\n
+   *
+   * \param  fileDescriptor     PDF file descriptor
+   */
+  JMuonNPE_t(const std::string& fileDescriptor)
+  {
+    using namespace std;
+    using namespace JPP;
+
+    const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
+				 SCATTERED_LIGHT_FROM_MUON,
+				 DIRECT_LIGHT_FROM_DELTARAYS,
+				 SCATTERED_LIGHT_FROM_DELTARAYS,
+ 				 DIRECT_LIGHT_FROM_EMSHOWERS,
+				 SCATTERED_LIGHT_FROM_EMSHOWERS };
+
+    const  int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
+    
+    typedef JPP::JSplineFunction1D<JSplineElement2D<double, double>,
+                                   JCollection,
+                                   double>                          JFunction1D_t;
+    typedef JPDFTable<JFunction1D_t, JNPEMaplist_t>                 JPDF_t;
+
+
+    const JNPE_t::JSupervisor supervisor(new JNPE_t::JDefaultResult(zero));
+
+    for (int i = 0; i != N; ++i) {
+
+      JPDF_t pdf;
+
+      const JPDFType_t type      = pdf_t[i];
+      const string     file_name = getFilename(fileDescriptor, type);
+
+      cout << "loading PDF from file " << file_name << "... " << flush;
+
+      pdf.load(file_name.c_str());
+
+      cout << "OK" << endl;
+
+      pdf.setExceptionHandler(supervisor);
+
+      if      (is_bremsstrahlung(type))
+	YB.push_back(JNPE_t(pdf));
+      else if (is_deltarays(type))
+	YA.push_back(JNPE_t(pdf));
+      else
+	Y1.push_back(JNPE_t(pdf));
+    }
+
+    // Add PDFs
+
+    cout << "adding PDFs... " << flush;
+
+    Y1[1].add(Y1[0]); Y1.erase(Y1.begin());
+    YA[1].add(YA[0]); YA.erase(YA.begin());
+    YB[1].add(YB[0]); YB.erase(YB.begin());
+  
+    cout << "OK" << endl;
+  }
+
+
+  /**
+   * Get PDF.
+   *
+   * The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
+   * In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
+   *
+   * \param  E                  muon energy at minimum distance of approach [GeV]
+   * \param  R                  minimum distance of approach [m]
+   * \param  theta              PMT zenith  angle [rad]
+   * \param  phi                PMT azimuth angle [rad]
+   * \return                    number of photo-electrons
+   */
+  double calculate(const double E,
+		   const double R,
+		   const double theta,
+		   const double phi) const
+  {
+    using namespace JPP;
+
+    const double y1 = getNPE(Y1, R, theta, phi);
+    const double yA = getNPE(YA, R, theta, phi);
+    const double yB = getNPE(YB, R, theta, phi);
+
+    if (E >= MASS_MUON * INDEX_OF_REFRACTION_WATER)
+      return y1  +  getDeltaRaysFromMuon(E) * yA  +  E * yB;
+    else
+      return 0.0;
+  }
+
+private:
+  std::vector<JNPE_t> Y1;     //!< light from muon
+  std::vector<JNPE_t> YA;     //!< light from delta-rays
+  std::vector<JNPE_t> YB;     //!< light from EM showers
+
+  /**
+   * Get number of photo-electrons.
+   *
+   * \param  NPE                NPE tables
+   * \param  R                  distance between muon and PMT [m]
+   * \param  theta              zenith  angle orientation PMT [rad]
+   * \param  phi                azimuth angle orientation PMT [rad]
+   * \return                    number of photo-electrons
+   */
+  static inline double getNPE(const std::vector<JNPE_t>& NPE,
+			      const double R,
+			      const double theta,
+			      const double phi)
+  {
+    using namespace std;
+    using namespace JPP;
+
+    double npe = 0.0;
+
+    for (vector<JNPE_t>::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
+
+      if (R <= i->getXmax()) {
+
+	try {
+
+	  const double y = get_value((*i)(std::max(R, i->getXmin()), theta, phi));
+
+	  if (y > 0.0) {
+	    npe += y;
+	  }
+	}
+	catch(const exception& error) {
+	  cerr << error.what() << endl;
+	}
+      }
+    }
+
+    return npe;
+  }
+};
+
+
+/**
+ * Auxiliary data structure for shower PDF.
+ */
+struct JShowerNPE_t {
+
+  typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
+			JPP::JPolint1FunctionalMap,
+			JPP::JPolint1FunctionalGridMap,
+			JPP::JPolint1FunctionalGridMap>::maplist   JNPEMaplist_t;
+  typedef JPP::JNPETable<double, double, JNPEMaplist_t>            JNPE_t;
+
+
+  /**
+   * Constructor.
+   *
+   * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\n
+   *
+   * \param  fileDescriptor     PDF file descriptor
+   * \param  numberOfPoints     number of points for shower elongation
+   */
+  JShowerNPE_t(const std::string& fileDescriptor,
+	       const int          numberOfPoints = 0) :
+    numberOfPoints(numberOfPoints)
+  {
+    using namespace std;
+    using namespace JPP;
+    
+    const JPDFType_t pdf_t[] = { SCATTERED_LIGHT_FROM_EMSHOWER,
+				 DIRECT_LIGHT_FROM_EMSHOWER };
+
+    const  int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
+
+    typedef JPP::JSplineFunction1D<JSplineElement2D<double, double>,
+                                   JCollection,
+                                   double>                          JFunction1D_t;
+    typedef JPDFTable<JFunction1D_t, JNPEMaplist_t>                 JPDF_t;
+
+
+    const JNPE_t::JSupervisor supervisor(new JNPE_t::JDefaultResult(zero));
+
+    for (int i = 0; i != N; ++i) {
+
+      const string file_name = getFilename(fileDescriptor, pdf_t[i]);
+
+      cout << "loading input from file " << file_name << "... " << flush;
+
+      JPDF_t pdf;
+
+      pdf.load(file_name.c_str());
+
+      pdf.setExceptionHandler(supervisor);
+
+      if (npe.empty())
+	npe = JNPE_t(pdf);
+      else
+	npe.add(JNPE_t(pdf));
+
+      F[i] = JNPE_t(pdf);
+
+      cout << "OK" << endl;
+    }
+  }
+
+
+  /**
+   * Get PDF.
+   *
+   * The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
+   * In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
+   *
+   * \param  E                  shower energy at minimum distance of approach [GeV]
+   * \param  D                  distance [m]
+   * \param  cd                 cosine emission angle
+   * \param  theta              PMT zenith  angle [rad]
+   * \param  phi                PMT azimuth angle [rad]
+   * \return                    hypothesis value
+   */
+  double calculate(const double E,
+		   const double D,
+		   const double cd,
+		   const double theta,
+		   const double phi) const
+  {
+    using namespace std;
+    using namespace JPP;
+
+    double Y = 0.0;
+  
+    if (numberOfPoints > 0) {
+
+      const double W = 1.0 / (double) numberOfPoints;
+
+      for (int i = 0; i != numberOfPoints; ++i) {
+
+	const double z = geanz.getLength(E, (i + 0.5) / (double) numberOfPoints);
+
+	const double __D  = sqrt(D*D - 2.0*(D*cd)*z + z*z);
+	const double __cd = (D * cd - z) / __D;
+
+	try {
+	  Y += W * npe (__D, __cd, theta, phi);
+	}
+	catch(const exception& error) {
+	  //cerr << error.what() << endl;
+	}
+      }
+
+    } else {
+
+      try {
+	Y = npe(D, cd, theta, phi);
+      }
+      catch(const exception& error) {
+	//cerr << error.what() << endl;
+      }
+    }
+
+    return E * Y;
+  }
+
+private:
+  int            numberOfPoints;
+  JNPE_t npe;    //!< PDF for shower
+  JNPE_t F[2];   //!< PDF for shower
+};

src/jpp/JPhysics/JPDFTable.hh

+#ifndef __JPHYSICS__JPDFTABLE__
+#define __JPHYSICS__JPDFTABLE__
+
+#include <cmath>
+
+#include "JIO/JObjectBinaryIO.hh"
+#include "JTools/JTransformableMultiFunction.hh"
+#include "JTools/JQuantiles.hh"
+#include "JTools/JSet.hh"
+#include "JTools/JRange.hh"
+#include "JMath/JMathSupportkit.hh"
+#include "JPhysics/JConstants.hh"
+#include "JPhysics/JPDFTransformer.hh"
+
+
+/**
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  using JIO::JSerialisable;
+  using JIO::JReader;
+  using JIO::JWriter;
+  using JIO::JObjectBinaryIO;
+  using JTOOLS::JTransformableMultiFunction;
+  using JTOOLS::JTransformableMultiHistogram;
+  using JTOOLS::JRange;
+  
+
+  /**
+   * Multi-dimensional PDF table for arrival time of Cherenkov light.
+   */
+  template<class JFunction1D_t,
+	   class JMaplist_t,
+           class JDistance_t = JTOOLS::JDistance<typename JFunction1D_t::argument_type> >
+  class JPDFTable :
+    public JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>,
+    public JSerialisable,
+    public JObjectBinaryIO< JPDFTable<JFunction1D_t, JMaplist_t, JDistance_t> >
+  {
+  public:
+
+    typedef JTransformableMultiFunction<JFunction1D_t, JMaplist_t, JDistance_t>      transformablemultifunction_type;
+
+    typedef typename transformablemultifunction_type::argument_type                  argument_type;
+    typedef typename transformablemultifunction_type::result_type                    result_type;
+    typedef typename transformablemultifunction_type::value_type                     value_type;
+
+    typedef typename transformablemultifunction_type::multimap_type                  multimap_type;
+    typedef typename transformablemultifunction_type::transformer_type               transformer_type;
+
+    enum { NUMBER_OF_DIMENSIONS = transformablemultifunction_type::NUMBER_OF_DIMENSIONS };
+
+    typedef typename transformablemultifunction_type::super_const_iterator           super_const_iterator;
+    typedef typename transformablemultifunction_type::super_iterator                 super_iterator;
+    typedef typename transformablemultifunction_type::function_type                  function_type;
+
+    using transformablemultifunction_type::insert;
+
+
+    /**
+     * Default constructor.
+     */
+    JPDFTable() :
+      transformablemultifunction_type()
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  input                multi-dimensional function
+     */
+    template<class __JFunction_t, class __JMaplist_t, class __JDistance_t>    
+    JPDFTable(const JTransformableMultiFunction<__JFunction_t, __JMaplist_t, __JDistance_t>& input) :
+      transformablemultifunction_type(input)
+    {}
+
+    
+    /**
+     * Constructor.
+     *
+     * \param  input                multi-dimensional histogram
+     */
+    template<class JHistogram_t, class __JMaplist_t, class __JDistance_t>
+    JPDFTable(const JTransformableMultiHistogram<JHistogram_t, __JMaplist_t, __JDistance_t>& input) :
+      transformablemultifunction_type(input)
+    {}
+
+
+    /**
+     * Blur PDF.
+     *
+     * The arrival times of Cherenkov light are smeared according to a Gaussian distribution
+     * with the specified width (i.e.\ TTS) using Gauss-Hermite integration.
+     * An exception is made when the time range according the specified quantile is 
+     * smaller than the specified width (TTS) of the Gaussian distribution.
+     * In that case, the resulting PDF is a Gaussian distribution with the specified width (TTS)
+     * and normalisation according to the integral value of the input PDF.
+     * A smooth transition is imposed between the normal regime and this exeption.
+     *
+     * \param  TTS                  TTS [ns]
+     * \param  numberOfPoints       number of points for Gauss-Hermite integration
+     * \param  epsilon              precision
+     * \param  quantile             quantile
+     */
+    void blur(const double TTS, 
+	      const int    numberOfPoints = 25,
+	      const double epsilon        = 1.0e-10,
+	      const double quantile       = 0.99)
+    {
+      using namespace std;
+      using namespace JPP;
+
+      typedef typename transformer_type::array_type  array_type;
+
+      const JGaussHermite engine(numberOfPoints, epsilon);
+
+      for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
+	
+	const array_type array = (*i).getKey();
+	function_type&   f1    = (*i).getValue();
+
+	if (!f1.empty()) {
+
+	  const typename function_type::supervisor_type& supervisor = f1.getSupervisor();
+
+	  const JMultiMapGetTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> get(*(this->transformer), array);
+	  const JMultiMapPutTransformer<NUMBER_OF_DIMENSIONS - 1, value_type> put(*(this->transformer), array);
+
+	  f1.transform(get);
+	  f1.compile();
+
+	  const JQuantiles Q(f1, quantile);
+
+	  // abscissa 
+
+	  JSet<double> X;
+
+	  for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
+	    X.insert(Q.getX() + TTS*sqrt(2.0)*j->getX());
+	  }
+
+	  for (typename function_type::const_iterator j = f1.begin(); j != f1.end(); ++j) {
+
+	    if (j->getX() - TTS < X.getXmin()) {
+	      X.insert(j->getX() - TTS);
+	    }
+	    
+	    if (j->getX() + TTS > X.getXmax()) {
+	      X.insert(j->getX() + TTS);
+	    }
+	  }
+
+
+	  const double W = gauss(Q.getUpperLimit() - Q.getLowerLimit(), TTS);
+
+	  function_type buffer;
+
+	  for (JSet<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
+
+	    double y = 0.0;
+
+	    for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
+
+	      const double u = j->getX();
+	      const double v = j->getY() / sqrt(PI);
+	      const double w = get_value(f1(*x + u*TTS*sqrt(2.0)));
+	
+	      y += v * w;
+	    }
+
+	    buffer[*x]  =  W * Q.getIntegral() * Gauss(*x - Q.getX(), TTS)  +  (1.0 - W) * y;
+	  }
+
+	  buffer.transform(put);
+	  buffer.compile();
+
+	  f1 = buffer;
+
+	  f1.setExceptionHandler(supervisor);
+	}
+      }
+    }
+
+
+    /**
+     * Compresses PDF to given abscissa range.
+     *
+     * \param   range               abscissa range
+     */
+    void compress(const JRange<typename function_type::abscissa_type>& range)
+    {        
+      for (super_iterator i = this->super_begin(); i != this->super_end(); ++i) {
+
+	function_type& f1 = i.getValue();
+	
+	typename function_type::iterator p = f1.lower_bound(range.getLowerLimit());
+
+	f1.function_type::container_type::erase(f1.begin(), p);
+
+	typename function_type::iterator q = f1.lower_bound(range.getUpperLimit());
+
+	f1.function_type::container_type::erase(++q, f1.end());
+      }
+
+      this->compile();
+    }
+
+
+    /**
+     * Read from input.
+     *
+     * \param  in                   reader
+     * \return                      reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      if (in >> static_cast<transformablemultifunction_type&>(*this)) {
+
+	// read optional transformer
+
+	JPDFTransformer<NUMBER_OF_DIMENSIONS - 1, argument_type> buffer;
+
+	if (buffer.read(in)) {
+
+	  this->transformer.reset(buffer.clone());
+
+	} else {
+
+	  in.clear();
+	  
+	  this->transformer.reset(transformer_type::getClone());
+	}
+      }
+      
+      this->compile();
+      
+      return in;
+    }
+    
+    
+    /**
+     * Write from input.
+     *
+     * \param  out                  writer
+     * \return                      writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << static_cast<const transformablemultifunction_type&>(*this);
+
+      this->transformer->write(out);
+
+      return out;
+    }
+  };
+}
+
+#endif

src/jpp/JPhysics/JPDFToolkit.hh

+#ifndef __JPHYSICS__JPDFTOOLKIT__
+#define __JPHYSICS__JPDFTOOLKIT__
+
+#include <cmath>
+
+#include "JPhysics/JConstants.hh"
+
+
+/**
+ * \file
+ * Auxiliary methods for PDF calculations.
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+
+  /**
+   * Get minimal wavelength for PDF evaluations.
+   *
+   * \return            wavelength of light [nm]
+   */
+  inline double getMinimalWavelength()
+  {
+    return 300.0;
+  }
+
+
+  /**
+   * Get maximal wavelength for PDF evaluations.
+   *
+   * \return            wavelength of light [nm]
+   */
+  inline double getMaximalWavelength()
+  {
+    return 700.0;
+  }
+  
+
+  /**
+   * Number of Cherenkov photons per unit track length and per unit wavelength.
+   *
+   * \param  lambda     wavelength of light [nm]
+   * \param  n          index of refraction
+   * \return            number of photons per unit track length and per unit wavelength [m^-1 nm^-1]
+   */
+  inline double cherenkov(const double lambda,
+			  const double n)
+  {
+    const double x = n*lambda;
+    
+    return 1.0e9 * 2 * PI * ALPHA_ELECTRO_MAGNETIC * (n*n - 1.0) / (x*x);
+  }
+
+  
+  /**
+   * Equivalent EM-shower energy due to delta-rays per unit muon track length.
+   *
+   * Internal parameters are obtained with application [script] JDeltaRays[.sh].
+   * 
+   * \param  E      muon energy [GeV]
+   * \return        equivalent energy loss [GeV/m]
+   */  
+  inline double getDeltaRaysFromMuon(const double E)
+  {
+    static const double a =  3.186e-01;
+    static const double b =  3.384e-01;
+    static const double c = -2.759e-02;
+    static const double d =  1.630e-03;
+    static const double Emin = 0.13078; // [GeV]
+
+    if (E > Emin) {
+
+      const double x = log10(E);                      //
+      const double y = a + x*(b + x*(c + x*(d)));     // [MeV g^-1 cm^2]
+
+      return y * DENSITY_SEA_WATER * 1.0e-1;          // [GeV/m]
+    }
+
+    return 0.0;
+  }
+
+  
+  /**
+   * Equivalent EM-shower energy due to delta-rays per unit tau track length.
+   *
+   * Internal parameters are obtained with application [script] JDeltaRays[.sh].
+   * 
+   * \param  E      tau energy [GeV]
+   * \return        equivalent energy loss [GeV/m]
+   */  
+  inline double getDeltaRaysFromTau(const double E)
+  {
+    static const double a = -2.374e-01;
+    static const double b =  5.143e-01;
+    static const double c = -4.213e-02;
+    static const double d =  1.804e-03;
+    static const double Emin = 2.19500; // [GeV]
+    
+    if (E > Emin) {
+
+      const double x = log10(E);                      //
+      const double y = a + x*(b + x*(c + x*(d)));     // [MeV g^-1 cm^2]
+
+      return y * DENSITY_SEA_WATER * 1.0e-1;          // [GeV/m]
+    }
+
+    return 0.0;
+  }
+
+
+  /**
+   * Emission profile of photons from delta-rays.
+   *
+   * Profile is taken from reference ANTARES-SOFT-2002-015, J.\ Brunner (fig.\ 3).
+   *
+   * \param  x      cosine emission angle
+   * \return        probability
+   */
+  inline double getDeltaRayProbability(const double x)
+  {
+    //return 1.0 / (4.0 * PI);
+    return 0.188 * exp(-1.25 * pow(fabs(x - 0.90), 1.30));
+  }
+
+
+  /**
+   * Rayleigh cross section.
+   *
+   * \param  n          index of refraction
+   * \param  lambda     wavelength of light [nm]
+   * \return            cross section       [cm^2]
+   */
+  inline const double getRayleighCrossSection(const double n,
+					      const double lambda)
+  {
+    static const double d = 0.36;                    // size of H2O molecule [nm]
+    static const double U = PI*PI*PI*PI*PI*2.0/3.0;
+    static const double V = d*d*d*d*d*d;
+
+    const double W     = (n*n - 1.0) / (n*n + 2.0);
+    const double sigma = 1.0e-14 * U*V*W*W / (lambda*lambda*lambda*lambda);   // [cm^2]
+
+    return sigma;
+  }
+
+
+  /**
+   * Rayleigh scattering length.
+   *
+   * \param  n          index of refraction
+   * \param  lambda     wavelength of light [nm]
+   * \return            scattering length   [m]
+   */
+  inline const double getRayleighScatteringLength(const double n,
+						  const double lambda)
+  {
+    static const double amu = 18.01528; // H20 mass in Atomic mass units
+
+    const double sigma = getRayleighCrossSection(n, lambda);
+    const double ls    = 1.0e-2 / (DENSITY_SEA_WATER * AVOGADRO * sigma / amu);     // [m]
+
+    return ls;
+  }
+}
+
+#endif

src/jpp/JPhysics/JPDFTransformer.hh

+#ifndef __JPHYSICS__JPDFTRANSFORMER__
+#define __JPHYSICS__JPDFTRANSFORMER__
+
+#include <cmath>
+
+#include "JLang/JCC.hh"
+#include "JIO/JSerialisable.hh"
+#include "JPhysics/JConstants.hh"
+#include "JTools/JMultiMapTransformer.hh"
+#include "JTools/JFunction1D_t.hh"
+#include "JTools/JGrid.hh"
+#include "JPhysics/JGeant_t.hh"
+
+
+/**
+ * \author mdejong
+ */
+
+namespace JPHYSICS {}
+namespace JPP { using namespace JPHYSICS; }
+
+namespace JPHYSICS {
+
+  using JIO::JReader;
+  using JIO::JWriter;
+  using JTOOLS::JMultiMapTransformer;
+
+
+  /**
+   * Transformer for the 1D probability density function (PDF) of the time response of a PMT to a muon.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# distance of closest approach of the muon to the PMT [m]
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# effective attenuation length
+   */
+  template<class JArgument_t>
+  class JPDFTransformer_t :
+    public JMultiMapTransformer<1, JArgument_t>
+  {
+  public:
+
+    typedef JMultiMapTransformer<1, JArgument_t>                     JMultiMapTransformer_t;
+
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+
+    using JMultiMapTransformer_t::getWeight;
+
+    /**
+     * Get shortest distance of approach.
+     *
+     * \return              distance [m]
+     */
+    static double getRmin()
+    {
+      return 0.01;
+    }
+
+
+    /**
+     * Default constructor.
+     */
+    JPDFTransformer_t() :
+      ln   (0.0),
+      alpha(0),
+      kmin (0.0),
+      kmax (0.0)
+    {}
+
+
+    /**
+     * Constructor.
+     *  
+     * \param  ln           Effective attenuation length [m]
+     * \param  alpha        Distance dependence (power term)
+     * \param  kmin         Minimal kappa
+     * \param  kmax         Maximal kappa
+     */
+    JPDFTransformer_t(const double ln,
+		      const int    alpha,
+		      const double kmin,
+		      const double kmax) :
+      ln   (ln),
+      alpha(alpha),
+      kmin (kmin),
+      kmax (kmax)
+    {}
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPDFTransformer_t(*this);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {R_m}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      using namespace JTOOLS;
+
+      const double R = buffer[0];
+
+      double x = xn;
+
+      const double t0 =  R * getTanThetaC() * getInverseSpeedOfLight();
+      const double t1 =  R * kmin * getInverseSpeedOfLight();
+
+      x -= t1 - t0;
+      
+      if (kmax > kmin) {
+	x /= R * (kmax - kmin) * getInverseSpeedOfLight();
+      }
+
+      return x;
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {R_m}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      using namespace JTOOLS;
+
+      const double R = buffer[0];
+
+      double x = xn;
+
+      if (kmax > kmin) {
+	x *= R * (kmax - kmin) * getInverseSpeedOfLight();
+      }
+
+      const double t0 =  R * getTanThetaC() * getInverseSpeedOfLight();
+      const double t1 =  R * kmin * getInverseSpeedOfLight();
+      
+      x += t1 - t0;
+
+      return x;
+    }
+
+
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {R_m}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    {
+      using namespace JTOOLS;
+
+      const double R = buffer[0];
+
+      const double n   = getIndexOfRefraction();
+      const double ct0 = 1.0 / n;
+      const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
+    
+      const double d  = sqrt(getRmin()*getRmin() + R*R) / st0;
+    
+      return exp(-d/ln) / pow(d,alpha);
+    }
+
+
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> ln;
+      in >> alpha;
+      in >> kmin;
+      in >> kmax;
+
+      return in;
+    }
+
+
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << ln;
+      out << alpha;
+      out << kmin;
+      out << kmax;
+
+      return out;
+    }
+
+
+    /**
+     * Print PDF transformer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      using namespace std;
+
+      out << "Effective attenuation length [m]   " << ln    << endl;
+      out << "Distance dependence (power term)   " << alpha << endl;
+      out << "Minimal kappa                      " << kmin  << endl;
+      out << "Maximal kappa                      " << kmax  << endl;
+
+      return out;
+    }
+
+
+    double ln;                 //!< Effective attenuation length [m]
+    int    alpha;              //!< Distance dependence (power term)
+    double kmin;               //!< minimal kappa
+    double kmax;               //!< maximal kappa
+  };
+
+
+  /**
+   * Transformer for the 1D probability density function (PDF) of the time response of a PMT due to a point source.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# distance between point source and PMT [m]
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# effective attenuation length
+   */
+  template<class JArgument_t>
+  class JPD0Transformer_t : 
+    public JMultiMapTransformer<1, JArgument_t>
+  {
+  public:
+
+    typedef JMultiMapTransformer<1, JArgument_t>                     JMultiMapTransformer_t;
+
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+
+    using JMultiMapTransformer_t::getWeight;
+
+    /**
+     * Get shortest distance.
+     *
+     * \return              distance [m]
+     */
+    static double getDmin()
+    {
+      return 0.01;
+    }
+
+
+    /**
+     * Default constructor.
+     */
+    JPD0Transformer_t() :
+      ln   (0.0),
+      alpha(0),
+      kmin (0.0),
+      kmax (0.0)
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  ln         Effective attenuation length [m]
+     * \param  alpha      Distance dependence (power term)
+     * \param  kmin       Minimal kappa
+     * \param  kmax       Maximal kappa
+     */
+    JPD0Transformer_t(const double  ln,
+		      const int     alpha,
+		      const double  kmin,
+		      const double  kmax) :
+      ln   (ln),
+      alpha(alpha),
+      kmin (kmin),
+      kmax (kmax)
+    {}
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPD0Transformer_t(*this);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      using namespace JTOOLS;
+
+      const double D  = buffer[0];
+
+      double x = xn;
+
+      const double t0 =  D * getIndexOfRefraction() * getInverseSpeedOfLight();
+      const double t1 =  D * kmin * getInverseSpeedOfLight();
+
+      x -= t1 - t0;
+
+      if (kmax > kmin) {
+        x /= D * (kmax - kmin) * getInverseSpeedOfLight();
+      }
+
+      return x;
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      using namespace JTOOLS;
+
+      const double D  = buffer[0];
+
+      double x = xn;
+
+      if (kmax > kmin) {
+        x *= D * (kmax - kmin) * getInverseSpeedOfLight();
+      }
+
+      const double t0 =  D * getIndexOfRefraction() * getInverseSpeedOfLight();
+      const double t1 =  D * kmin * getInverseSpeedOfLight();
+
+      x += t1 - t0;
+
+      return x;
+    }
+
+
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {D_m}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    {
+      using namespace JTOOLS;
+
+      const double D  = buffer[0];
+
+      const double d  = sqrt(getDmin()*getDmin() + D*D);
+
+      return exp(-d/ln) / pow(d,alpha);
+    }
+
+
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> ln;
+      in >> alpha;
+      in >> kmin;
+      in >> kmax;
+
+      return in;
+    }
+
+
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << ln;
+      out << alpha;
+      out << kmin;
+      out << kmax;
+
+      return out;
+    }
+
+
+    /**
+     * Print PDF transformer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      using namespace std;
+
+      out << "Effective attenuation length [m]   " << ln    << endl;
+      out << "Distance dependence (power term)   " << alpha << endl;
+      out << "Minimal kappa                      " << kmin  << endl;
+      out << "Maximal kappa                      " << kmax  << endl;
+
+      return out;
+    }
+
+
+    double ln;                 //!< Effective attenuation length [m]
+    int    alpha;              //!< Distance dependence (power term)
+    double kmin;               //!< minimal kappa
+    double kmax;               //!< maximal kappa
+  };
+
+
+  /**
+   * Transformer for the 2D probability density function (PDF) of the time response of a PMT due to an EM shower.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# distance between EM shower and PMT [m]
+   *   -# cosine angle EM shower direction and EM shower - PMT position
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# effective attenuation length
+   *   -# emission profile of the photons
+   */
+  template<class JArgument_t>
+  class JPDGTransformer_t : 
+    public JMultiMapTransformer<2, JArgument_t>
+  {
+  public:
+
+    typedef JPD0Transformer_t      <JArgument_t>                     JFunction1DTransformer_t; 
+    typedef JMultiMapTransformer<2, JArgument_t>                     JMultiMapTransformer_t;
+
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+
+    using JMultiMapTransformer_t::getWeight;
+
+
+    /**
+     * Default constructor.
+     */
+    JPDGTransformer_t() :
+      transformer(),
+      getShowerProbability()
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  ln         Effective attenuation length [m]
+     * \param  alpha      Distance dependence (power term)
+     * \param  kmin       Minimal kappa
+     * \param  kmax       Maximal kappa
+     * \param  geant      Function photon emission from EM-shower
+     * \param  bmin       Baseline photon emission from EM-shower
+     */
+    JPDGTransformer_t(const double    ln,
+		      const int       alpha,
+		      const double    kmin,
+		      const double    kmax,
+		      const JGeant_t& geant,
+		      const double    bmin) :
+      transformer(ln, alpha, kmin, kmax),
+      getShowerProbability(geant)
+    {
+      getShowerProbability.add(bmin);
+      getShowerProbability.compile();
+    }
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPDGTransformer_t(*this);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m, cd}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.putXn(buffer, xn);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m, cd}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.getXn(buffer, xn);
+    }
+
+
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {D_m, cd}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    {
+      using namespace JTOOLS;
+
+      //const double D  = buffer[0];
+      const double cd = buffer[1];
+
+      return transformer.getWeight(buffer) * getShowerProbability(getIndexOfRefractionPhase(), cd);
+    }
+
+
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> transformer;
+      in >> getShowerProbability;
+
+      return in;
+    }
+
+
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << transformer;
+      out << getShowerProbability;
+
+      return out;
+    }
+
+
+    /**
+     * Print PDF transformer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      return transformer.print(out);
+    }
+
+
+    JFunction1DTransformer_t transformer;
+    JGeant_t                 getShowerProbability;
+  };
+
+
+  /**
+   * Template definition of transformer of the probability density function (PDF) of the time response of a PMT.\n
+   * The actual implementation follows from the number of dimensions.
+   */
+  template<unsigned int N, class JArgument_t>
+  class JPDFTransformer;
+
+
+  /**
+   * Template specialisation of transformer of the 2D probability density function (PDF) of the time response of a PMT due to a bright point.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# distance between bright point and PMT [m]
+   *   -# cosine PMT angle
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# effective attenuation length
+   */
+  template<class JArgument_t>
+  class JPDFTransformer<2, JArgument_t> : 
+    public JMultiMapTransformer<2, JArgument_t>
+  {
+  public:
+
+    typedef JPD0Transformer_t      <JArgument_t>                     JFunction1DTransformer_t;
+    typedef JMultiMapTransformer<2, JArgument_t>                     JMultiMapTransformer_t;
+
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+
+
+    /**
+     * Default constructor.
+     */
+    JPDFTransformer() :
+      transformer()
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  ln         Effective attenuation length [m]
+     * \param  alpha      Distance dependence (power term)
+     * \param  kmin       Minimal kappa
+     * \param  kmax       Maximal kappa
+     */
+    JPDFTransformer(const double ln,
+		    const int    alpha,
+		    const double kmin,
+		    const double kmax) :
+      transformer(ln, alpha, kmin, kmax)
+    {}
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPDFTransformer(*this);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m, ct}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.putXn(buffer, xn);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m, ct}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.getXn(buffer, xn);
+    }
+
+
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {D_m, ct}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    {
+      return transformer.getWeight(buffer);
+    }
+
+
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> transformer;
+
+      return in;
+    }
+
+
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << transformer;
+
+      return out;
+    }
+
+
+    /**
+     * Print PDF transfomer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      return transformer.print(out);
+    }
+
+
+    JFunction1DTransformer_t transformer;
+  };
+
+
+  /**
+   * Template specialisation of transformer of the 3D probability density function (PDF) of the time response of a PMT due to a muon.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# distance of closest approach of the muon to the PMT [m]
+   *   -# zenith    angle of the PMT
+   *   -# azimuthal angle of the PMT
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# effective attenuation length
+   *   -# angular acceptance of PMT
+   */
+  template<class JArgument_t>
+  class JPDFTransformer<3, JArgument_t> : 
+    public JMultiMapTransformer<3, JArgument_t>
+  {
+  public:
+
+    typedef JPDFTransformer_t      <JArgument_t>                     JFunction1DTransformer_t; 
+    typedef JMultiMapTransformer<3, JArgument_t>                     JMultiMapTransformer_t;
+ 
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+
+    typedef JTOOLS::JGridPolint1Function1D_t                         JFunction1D_t;
+
+    using JMultiMapTransformer_t::getWeight;
+
+
+    /**
+     * Default constructor.
+     */
+    JPDFTransformer() :
+      transformer(),
+      getAngularAcceptance()
+    {}
+
+
+    /**
+     * Constructor.
+     *  
+     * \param  ln         Effective attenuation length [m]
+     * \param  alpha      Distance dependence (power term)
+     * \param  kmin       Minimal kappa
+     * \param  kmax       Maximal kappa
+     * \param  pmt        Function angular acceptance of PMT
+     * \param  amin       Baseline angular acceptance of PMT
+     */
+    template<class T>
+    JPDFTransformer(const double ln,
+		    const int    alpha,
+		    const double kmin,
+		    const double kmax,
+		    T            pmt,
+		    const double amin) :
+      transformer(ln, alpha, kmin, kmax),
+      getAngularAcceptance()
+    {
+      getAngularAcceptance.configure(JTOOLS::make_grid(1000, -1.0, +1.0), pmt);
+      getAngularAcceptance.add(amin);
+      getAngularAcceptance.compile();
+      getAngularAcceptance.setExceptionHandler(new JFunction1D_t::JDefaultResult(0.0));
+    }
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPDFTransformer(*this);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {R_m, theta, phi}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.putXn(buffer, xn);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {R_m, theta, phi}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.getXn(buffer, xn);
+    }
+
+
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {R_m, theta, phi}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    {
+      using namespace JTOOLS;
+
+      //const double R     = buffer[0];
+      const double theta = buffer[1];
+      const double phi   = buffer[2];
+
+      const double n   = getIndexOfRefraction();
+      const double ct0 = 1.0 / n;
+      const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
+    
+      const double px = sin(theta)*cos(phi);
+      //const double py = sin(theta)*sin(phi);
+      const double pz = cos(theta);
+    
+      const double ct = st0*px + ct0*pz;
+    
+      return transformer.getWeight(buffer) * getAngularAcceptance(ct);
+    }
+
+
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> transformer;
+      in >> getAngularAcceptance;
+
+      getAngularAcceptance.compile();
+
+      return in;
+    }
+
+
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << transformer;
+      out << getAngularAcceptance;
+
+      return out;
+    }
+
+
+    /**
+     * Print PDF transformer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      return transformer.print(out);
+    }
+
+
+    JFunction1DTransformer_t transformer;
+    JFunction1D_t            getAngularAcceptance;
+  };
+
+
+  /**
+   * Template specialisation of transformer of the 4D probability density function (PDF) of the time response of a PMT due to an EM shower.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# distance between EM shower and PMT [m]
+   *   -# cosine angle EM shower direction and EM shower - PMT position
+   *   -# zenith    angle of the PMT
+   *   -# azimuthal angle of the PMT
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# effective attenuation length
+   *   -# emission profile of the photons
+   *   -# angular acceptance of PMT
+   */
+  template<class JArgument_t>
+  class JPDFTransformer<4, JArgument_t> : 
+    public JMultiMapTransformer<4, JArgument_t>
+  {
+  public:
+
+    typedef JPDGTransformer_t      <JArgument_t>                     JFunction2DTransformer_t;
+    typedef JMultiMapTransformer<4, JArgument_t>                     JMultiMapTransformer_t;
+
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+
+    typedef JTOOLS::JGridPolint1Function1D_t                         JFunction1D_t;
+
+    using JMultiMapTransformer_t::getWeight;
+
+
+    /**
+     * Default constructor.
+     */
+    JPDFTransformer() :
+      transformer(),
+      getAngularAcceptance()
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  ln         Effective attenuation length [m]
+     * \param  alpha      Distance dependence (power term)
+     * \param  kmin       Minimal kappa
+     * \param  kmax       Maximal kappa
+     * \param  geant      Function photon emission from EM-shower
+     * \param  bmin       Baseline photon emission from EM-shower
+     * \param  pmt        Function angular acceptance of PMT
+     * \param  amin       Baseline angular acceptance of PMT
+     */
+    template<class T>
+    JPDFTransformer(const double    ln,
+		    const int       alpha,
+		    const double    kmin,
+		    const double    kmax,
+		    const JGeant_t& geant,
+		    const double    bmin,
+		    T               pmt,
+		    const double    amin) :
+      transformer(ln, alpha, kmin, kmax, geant, bmin),
+      getAngularAcceptance()
+    {
+      getAngularAcceptance.configure(JTOOLS::make_grid(1000, -1.0, +1.0), pmt);
+      getAngularAcceptance.add(amin);
+      getAngularAcceptance.compile();
+      getAngularAcceptance.setExceptionHandler(new JFunction1D_t::JDefaultResult(0.0));
+    }
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPDFTransformer(*this);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m, cd, theta, phi}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.putXn(buffer, xn);
+    }
+
+
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {D_m, cd, theta, phi}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.getXn(buffer, xn);
+    }
+
+
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {D_m, cd, theta, phi}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    {
+      const double cd    = buffer[1];
+      const double theta = buffer[2];
+      const double phi   = buffer[3];
+
+      const double ct0 = (cd > -1.0 ? cd < +1.0 ? cd : +1.0 : -1.0);
+      const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
+    
+      const double px = sin(theta)*cos(phi);
+      //const double py = sin(theta)*sin(phi);
+      const double pz = cos(theta);
+    
+      const double ct = st0*px + ct0*pz;
+
+      return transformer.getWeight(buffer) * getAngularAcceptance(ct);
+    }
+
+
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> transformer;
+      in >> getAngularAcceptance;
+
+      getAngularAcceptance.compile();
+
+      return in;
+    }
+
+
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << transformer;
+      out << getAngularAcceptance;
+
+      return out;
+    }
+
+
+    /**
+     * Print PDF transfomer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      return transformer.print(out);
+    }
+
+
+    JFunction2DTransformer_t transformer;
+    JFunction1D_t            getAngularAcceptance;
+  };
+
+
+  /**
+   * Template specialisation of transformer of the 5D probability density function (PDF) of the time response of a PMT due to an EM shower.
+   *
+   * PDFs are evaluated by interpolation for:
+   *   -# energy of the EM shower
+   *   -# distance between EM shower and PMT [m]
+   *   -# cosine angle EM shower direction and EM shower - PMT position
+   *   -# zenith    angle of the PMT
+   *   -# azimuthal angle of the PMT
+   *   -# arrival time [ns]
+   *
+   * The evaluation of the weights is based on:
+   *   -# energy of the EM shower
+   *   -# effective attenuation length
+   *   -# emission profile of the photons
+   *   -# angular acceptance of PMT
+   */
+  template<class JArgument_t>
+  class JPDFTransformer<5, JArgument_t> : 
+    public JMultiMapTransformer<5, JArgument_t>
+  {
+  public:
+    
+    typedef JPDFTransformer     <4, JArgument_t>                     JFunction4DTransformer_t;
+    typedef JMultiMapTransformer<5, JArgument_t>                     JMultiMapTransformer_t;
+    
+    typedef typename JMultiMapTransformer_t::clone_type              clone_type;
+    typedef typename JMultiMapTransformer_t::argument_type           argument_type;
+    typedef typename JMultiMapTransformer_t::const_array_type        const_array_type;
+    
+    
+    /**
+     * Default constructor.
+     */
+    JPDFTransformer() :
+      transformer()
+    {}    
+
+    
+    /**
+     * Constructor.
+     *
+     * \param  transformer    transformer
+     */
+    JPDFTransformer(const JFunction4DTransformer_t& transformer) :
+      transformer(transformer)
+    {}
+
+
+    /**
+     * Constructor.
+     *
+     * \param  ln         Effective attenuation length [m]
+     * \param  alpha      Distance dependence (power term)
+     * \param  kmin       Minimal kappa
+     * \param  kmax       Maximal kappa
+     * \param  geant      Function photon emission from EM-shower
+     * \param  bmin       Baseline photon emission from EM-shower
+     * \param  pmt        Function angular acceptance of PMT
+     * \param  amin       Baseline angular acceptance of PMT
+     */
+    
+    template<class T>
+    JPDFTransformer(const double    ln,
+		    const int       alpha,
+		    const double    kmin,
+		    const double    kmax,
+		    const JGeant_t& geant,
+		    const double    bmin,
+		    T               pmt,
+		    const double    amin) :
+      transformer(ln, alpha, kmin, kmax, geant, bmin, pmt, amin)
+    {}
+
+
+    /**
+     * Clone object.
+     * 
+     * \return            pointer to newly created transformer
+     */
+    virtual clone_type clone() const override 
+    {
+      return new JPDFTransformer(*this);
+    }
+    
+    
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {E_GeV, D_m, cd, theta, phi}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.putXn(buffer.pop_front(), xn);
+    }
+    
+    
+    /**
+     * Evaluate arrival time.
+     *
+     * \param  buffer     {E_GeV, D_m, cd, theta, phi}
+     * \param  xn         old t_ns
+     * \return            new t_ns
+     */
+    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const override 
+    {
+      return transformer.getXn(buffer.pop_front(), xn);
+    }
+    
+    
+    /**
+     * Weight function.
+     *
+     * \param  buffer     {E_GeV, D_m, cd, theta, phi}
+     * \return            weight
+     */
+    virtual double getWeight(const_array_type& buffer) const override 
+    { 
+      const double E     = buffer[0];
+
+      return transformer.getWeight(buffer.pop_front()) / E;
+    }
+    
+    
+    /**
+     * Read PDF transformer from input.
+     *
+     * \param  in         reader
+     * \return            reader
+     */
+    virtual JReader& read(JReader& in) override 
+    {
+      in >> transformer;
+
+      return in;
+    }
+    
+  
+    /**
+     * Write PDF transformer to output.
+     *
+     * \param  out        writer
+     * \return            writer
+     */
+    virtual JWriter& write(JWriter& out) const override 
+    {
+      out << transformer;
+      
+      return out;
+    }
+    
+    
+    /**
+     * Print PDF transfomer to output stream.
+     *
+     * \param  out        output stream
+     * \return            output stream
+     */
+    std::ostream& print(std::ostream& out) const
+    {
+      return transformer.print(out);
+    }
+    
+    JFunction4DTransformer_t transformer;
+  };
+}
+
+#endif