JDIS.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JDIS__
#define __JPHYSICS__JDIS__
 
#include <cmath>
 
#include <TRandom3.h>
 
#include "JPhysics/JConstants.hh"
 
/**
 * \file
 * Deep-inelastic muon-nucleon scattering.
 *
 * \author mdejong
 */
 
namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }
 
namespace JPHYSICS {
 
  /**
   * Deep-inelastic muon-nucleon scattering.
   *
   * For cross section, see reference:\n
   * D.A. Timashkov and A.A. Petrukhin, 29th International Cosmic Ray Conference Pune (2005) 9, 89-92.
   *
   * For kinematics, see reference:\n
   * A. van Ginneken, Nucl.\ Instr.\ and Meth.\ A251 (1986) 21.
   */
  class JDIS {
  public:
    /**
     * Default constructor.
     */
    JDIS()
    {}
 
 
    /**
     * Get cross section.
     *
     * \param  E         muon energy   [GeV]
     * \return           cross section [cm^2]
     */
    inline double getCrossSection(const double E) const
    {
      const double x = log10(E*1.0e+1);
 
      if (E > JDIS_t::E0)
        return 0.35e-30 * pow(1.0 - JDIS_t::E0/E, 2.0) * exp(-2.10/x) * pow(10.0, 0.125*x);
      else
        return 0.0;
    }
 
 
    /**
     * Get probability of given energy fraction.
     *
     * \param  E         muon energy [GeV]
     * \param  v         energy fraction
     * \return           probability
     */
    double getP(const double E, const double v) const
    {
      const JDIS_t dis(E);
 
      return dis.getP(v);
    }
  
 
    /**
     * Get shower energy.
     *
     * \param  E         muon energy [GeV]
     * \return           shower energy [GeV]
     */
    double getE(const double E) const
    {
      const JDIS_t dis(E);
 
      return dis.getE();
    }
 
 
    /**
     * Get RMS of scattering angle.
     *
     * \param  E         muon energy [GeV]
     * \param  v         energy loss fraction
     * \return           angle [rad]
     */
    double getThetaRMS(const double E, const double v) const
    {
      const double u = 1.0 - v;
 
      if (E*v > 0.135)
        return (0.39 / (E*u)) * pow(sqrt(E)*u*v, 0.17) * (1.0 - 0.135/(E*v));
      else
        return  0.0;
    }
 
 
    /**
     * Get breakpoint.
     *
     * \param  E         muon energy [GeV]
     * \return           energy fraction
     */
    double getV(const double E) const
    {
      return JDIS_t::E1/E;
    }
 
 
    /**
     * Auxiliary helper class for kinematics of deep-inelastic muon-nucleon scattering at fixed muon energy.
     */
    class JDIS_t {
    public:
      /**
       * Constructor.
       *
       * \param  E         muon energy [GeV]
       */
      JDIS_t(const double E) :
        E (E),
        E2(E - MASS_MUON)
      {
        k2 = pow(E1, p0 - p1) / pow(1.0 - E1/E, 2);
 
        const double A = getA(E1/E, false) - getA(E0/E, false);
        const double B = getB(E2/E, false) - getB(E1/E, false);
 
        k1 = 1.0 / (A + k2 * B);
      }
 
 
      /**
       * Get probability of given energy fraction.
       *
       * \param  v         energy fraction
       * \return           probability
       */
      double getP(const double v) const
      {
        const double x = E*v;
        const double u = 1.0 - v; 
 
        if      (x > E2)
          return 0.0;
        else if (x > E1)
          return k1 * k2 * pow(x, p1) * u*u;
        else if (x > E0)
          return k1 *      pow(x, p0);
        else
          return 0.0;
      }
  
 
      /**
       * Get shower energy.
       *
       * \return           shower energy [GeV]
       */
      double getE() const
      {
        const double A = getA(E1/E) - getA(E0/E);
        const double B = getB(E2/E) - getB(E1/E);
 
        for ( ; ; ) {
 
          double y = gRandom->Rndm() * (A + B);
 
          if (y <= A) {
 
            y += getA(E0/E);
            y /= k1 * pow(E,p0) / (p0 + 1);
 
            const double v = pow(y, 1.0 / (p0 + 1));
 
            return E * v;
 
          } else {
 
            y -= A;
            y *= (getb(E2/E) - getb(E1/E)) / B;
            y +=  getb(E1/E);
            y /= k1 * k2 * pow(E,p1) / (p1 + 1);
          
            const double v = pow(y, 1.0 / (p1 + 1));
            const double u = 1.0 - v; 
 
            if (gRandom->Rndm() < u*u) {
              return E * v;
            }
          }
        }
      }
 
      const double E;                        //!< actual  energy [GeV]
      const double E2;                       //!< maximal energy [GeV]
 
      static constexpr double E0 = 0.144;    //!< minimal energy [GeV]
      static constexpr double E1 = 0.35;     //!< breakpoint     [GeV]
 
      static constexpr double p0 =  2.0;     //!< spectral index upto breakpoint
      static constexpr double p1 = -1.11;    //!< spectral index from breakpoint
 
    private:
      /**
       * Integral upto breakpoint.
       *
       * \param  v         energy fraction
       * \param  option    option (if false, suppress normalisation constant)
       * \return           integral
       */    
      double getA(const double v, const bool option = true) const
      {
        return (option ? k1 : 1.0)      * pow(E, p0) * pow(v, p0 + 1) * (1.0 / (p0 + 1));
      }
 
 
      /**
       * Integral from breakpoint.
       *
       * \param  v         energy fraction
       * \param  option    option (if false, suppress normalisation constant)
       * \return           integral
       */    
      double getB(const double v, const bool option = true) const
      {
        return (option ? k1 * k2 : 1.0) * pow(E, p1) * pow(v, p1 + 1) * (1.0 / (p1 + 1)  -
                                                                         2*v / (p1 + 2)  +
                                                                         v*v / (p1 + 3));
      }
 
 
      /**
       * Integral from breakpoint without suppression factor.
       *
       * \param  v         energy fraction
       * \return           integral
       */    
      double getb(const double v) const
      {
        return k1 * k2 * pow(E, p1) * pow(v, p1 + 1) / (p1 + 1);
      }
 
      double k1;                             //!< normalisation constant upto breakpoint
      double k2;                             //!< normalisation constant from breakpoint
    };
  };
}
 
#endif