JGeane.hh

Go to the documentation of this file.

#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.
   *
   * \return        equivalent muon track length [m/Gev]
   */
  inline double geanc()
  {
    return 4.7;                   // 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 Bremstrahlung [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 Bremstrahlung [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 Bremstrahlung [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 Bremstrahlung [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) {

        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 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) {

        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