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;
  using JTOOLS::JGridPolint0Function1D_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 JGridPolint0Function1D_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 = infty

      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
      
      u = (y1 - y0) / y0 / (n1 - n0);
      z = 1.0 + u*n0;
    }
 

    /**
     * Function compilation.
     */
    virtual void compile()
    {
      JTOOLS::JGridPolint0Function1D_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 = JTOOLS::JGridPolint0Function1D_t::operator()(ct - 1.0/n);

      return y * (z - u*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_max = std::min(xmax - 1.0/n, buffer.rbegin()->getX());
      const double x_min = std::max(xmin - 1.0/n, buffer. begin()->getX());

      double y = 0.0;

      y += buffer(x_max);
      y -= buffer(x_min);

      return y * (z - u*n);
    }


    /**
     * Read geant from input.
     *
     * \param  in       reader
     * \param  geant    geant
     * \return          reader
     */
    friend inline JReader& operator>>(JReader& in, JGeant& geant)
    {
      in >> geant.z;
      in >> geant.u;
      in >> static_cast<JTOOLS::JGridPolint0Function1D_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.z;
      out << geant.u;
      out << static_cast<const JTOOLS::JGridPolint0Function1D_t&>(geant);

      return out;
    }


  private:
    double z;     //!< offset of the normalisation dependence
    double u;     //!< slope  of the normalisation dependence
    JTOOLS::JGridPolint0Function1D_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