JGeant.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JGEANT__
#define __JPHYSICS__JGEANT__
 
#include <cmath>
#include <algorithm>
 
#include "JTools/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JIO/JSerialisable.hh"
#include "JPhysics/JGeanx.hh"
 
/**
 * \file
 * Photon emission profile EM-shower.
 * \author mdejong
 */
 
namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }
 
namespace JPHYSICS {
 
  using JIO::JReader;
  using JIO::JWriter;
  using JTOOLS::getSpeedOfLight;
  using JTOOLS::getIndexOfRefraction;
 
  typedef JTOOLS::JGridPolint1Function1D_t  JGeantFunction1D_t;
 
 
  /**
   * Function object for the probability density function of photon emission from EM-shower 
   * as a function of 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 : 
    public JGeantFunction1D_t
  {
  public:
    /**
     * Default constructor.
     */
    JGeant() 
    {}
 
 
    /**
     * Constructor.
     * 
     * \param  geanx   PDF of EM-shower
     * \param  dx      step size for interpolation of function values
     */
    JGeant(const JGeanx& geanx,
           const double  dx = 0.0001) 
    {
      using JTOOLS::PI;
 
      const double n0 = geanx.n - 0.02;
      const double n1 = geanx.n + 0.02;
 
      const JGeanx g0(geanx.a, geanx.b, n0);
      const JGeanx g1(geanx.a, geanx.b, n1);
 
      const double y0 = g0.evaluate(-1.0, +1.0);
      const double y1 = g1.evaluate(-1.0, +1.0);
 
      const double c  = 1.0 / (y0 * 2*PI);
 
      // following range corresponds to n = 1 to n = infinity
 
      for (double x = -2.0; x < 1.0; x += dx) {
        (*this)[x] = c * g0.evaluate(x + 1.0/n0);
      }
 
      compile();
 
      // linear approximation of dependence normalisation constant on index of refraction
      
      a1 = (y1 - y0) / y0 / (n1 - n0);
      a0 = 1.0 + a1*n0;
    }
 
 
    /**
     * Function compilation.
     */
    virtual void compile()
    {
      JGeantFunction1D_t::compile();
 
      buffer.clear();
 
      JTOOLS::integrate(*this, buffer);
 
      buffer.compile();
    }
 
 
    /**
     * 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& 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& geant)
    {
      out << geant.a0;
      out << geant.a1;
      out << static_cast<const JGeantFunction1D_t&>(geant);
 
      return out;
    }
 
 
  private:
    double a0;     //!< offset of the normalisation dependence
    double a1;     //!< slope  of the normalisation dependence
    JGeantFunction1D_t buffer;
  };
 
 
  /**
   * Function object for the number of photons from EM-shower as a function of emission angle.
   */
  static const JGeant geant(geanx, 0.0001);
}
 
#endif