JGeanz.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JGEANZ__
#define __JPHYSICS__JGEANZ__
 
#include <istream>
#include <ostream>
#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;
    }
 
 
    /**
     * Derivative of PDF.
     *
     * \param  E       EM-shower energy [GeV]
     * \param  z       z position [m] (> 0)
     * \return         dP/dz
     */
    double getDerivative(const double E, 
                         const double z) const
    {
      if (E > Emin) {
 
        const double a = a0 + a1 * log(E);
 
        return ((a - 1.0)/z - 1.0/b) * getProbability(E, z);
      }
 
      return 0.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-5) const
    {
      if (E > Emin) {
 
        double z0 = 0.0;
        double z1 = getMaximum(E);
        
        for (int i = 1000; i != 0; --i) {
 
          const double p = getIntegral(E, z1);
 
          if (fabs(p-P) < P*eps) {
            return z1;
          }
          
          if (p > P) {
            z1  = 0.5 * (z0 + z1);
          } else {
            z0  = z1;
            z1 *= 1.5;
          }
        }
        
        return z1;
 
      } 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;
    }
 
 
    /**
     * Read longitudinal profile from input.
     *
     * \param  in              input stream
     * \param  object          longitudinal profile
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JGeanz& object)
    {
      in >> object.a0 >> object.a1 >> object.b;
 
      object.Emin = exp(-object.a0/object.a1);
 
      return in;
    }
 
 
    /**
     * Write longitidunal profile to output.
     *
     * \param  out             output stream
     * \param  object          longitudinal profile
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JGeanz& object)
    {
      return out << object.a0 << ' ' << object.a1 << ' ' << object.b;
    }
 
  protected:
    double a0;
    double a1;
    double b;
    double Emin;
  };
 
 
  /**
   * Function object for longitudinal EM-shower profile
   */
  static const JGeanz geanz(1.85, 0.62, 0.54);
}
 
#endif