JDeltaRays.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JDELTARAYS__
#define __JPHYSICS__JDELTARAYS__
 
#include "JPhysics/JConstants.hh"
#include "JPhysics/JPhysicsToolkit.hh"
#include "JPhysics/JEnergyRange.hh"
#include "JPhysics/JSeaWater.hh"
 
/**
 * \author mdejong
 */
 
namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }
 
namespace JPHYSICS {
 
  /**
   * Auxiliary data structure for delta-rays.
   */
  struct JDeltaRays {
 
    static constexpr double TMIN_GEV  = 0.000915499;             //!< Minimum allowed delta-ray kinetic energy [GeV]
    static constexpr double TMAX_GEV  = 1.0e10;                  //!< Maximum allowed delta-ray kinetic energy [GeV]
    static constexpr double K         = 0.307075;                //!< internal constant [MeV mol^-1 cm^2]
 
 
    /**
     * Get ratio charge to mass number for given atomic parameters.
     *
     * \param  parameters    atomic parameters
     * \return               ratio charge to mass number
     */
    static constexpr double getZoverA(const JSeaWater::atom_type& parameters)
    {
      return parameters() * parameters.Z / parameters.A;
    }
 
 
    /**
     * Get average ratio charge to mass number for sea water.
     *
     * \return               ratio charge to mass number
     */
    static constexpr double getZoverA()
    {
      return (getZoverA(JSeaWater::H) +
              getZoverA(JSeaWater::O) +
              getZoverA(JSeaWater::Na) +
              getZoverA(JSeaWater::Cl));
    }
 
 
    /**
     * Get minimum delta-ray kinetic energy.
     *
     * \return               minimum delta-ray kinetic energy [GeV]
     */
    static inline double getTmin()
    {
      const double Emin = MASS_ELECTRON / getSinThetaC();        // delta-ray Cherenkov threshold [GeV]
 
      return getKineticEnergy(Emin, MASS_ELECTRON);
    }
 
 
    /**
     * Get maximum delta-ray kinetic energy for given lepton energy and mass.\n
     * This formula is taken from reference
     * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
     * \n\n
     * N.B.: This function should not be used for electrons.\n
     *       For electrons use `0.5 * getKineticEnergy(E, MASS_ELECTRON)` instead.
     *
     * \param  E             particle energy                  [GeV]
     * \param  M             particle mass                    [GeV]
     * \return               maximum delta-ray kinetic energy [GeV]
     */
    static inline double getTmax(const double E,
                                 const double M)
    {
      const double ratio = MASS_ELECTRON / M;
 
      const double gamma = E / M;
      const double beta  = getBeta(gamma);
 
      return ( (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) /
               (1.0 + 2.0*gamma*ratio + ratio*ratio) );
    }                        
 
 
    /**
     * Auxiliary class for 2nd order polynomial form factor.
     */
    struct JFormFactor
    {
      /**
       * Constructor.
       *
       * \param  a           2nd order polynomial coefficient [GeV^-2]
       * \param  b           1st order polynomial coefficient [GeV^-1]
       * \param  c           0th order polynomial coefficient [unit]
       */
      JFormFactor(const double a,
                  const double b,
                  const double c) :
        a(a),
        b(b),
        c(c)
      {}
 
 
      /**
       * Get form factor for given delta-ray kinetic energy.
       *
       * \param  T           delta-ray kinetic energy [GeV]
       * \return             form factor              [unit]
       */
      double operator()(const double T) const
      {
        return c + T*(b + T*a);
      }
 
 
    private:
      double a;          //!< 2nd order polynomial coefficient [GeV^-2]
      double b;          //!< 1st order polynomial coefficient [GeV^-1]
      double c;          //!< 0th order polynomial coefficient [unit]
    };
 
 
    /**
     * Get number of delta-rays per unit track length for an ionising particle with given energy and given mass.
     *
     * The template parameter corresponds to a class which contains an `operator()(const double)`\n
     * to compute the form factor corresponding to a given delta-ray kinetic energy.
     *
     * \param  E      particle energy                                     [GeV]
     * \param  M      particle mass                                       [GeV]
     * \param  Tmin   minimum delta-ray kinetic energy                    [GeV]
     * \param  Tmax   maximum delta-ray kinetic energy                    [GeV]
     * \param  Z      atomic number                                       [unit]
     * \param  A      atomic mass                                         [g/mol]
     * \param  F      form factor functor
     * \param  N      number of points for numeric integration
     * \return        delta-ray count                                     [g^-1 cm^2]
     */
    template<class JFormFactor_t>
    static inline double getCount(const double         E,
                                  const double         M,
                                  const double         Tmin,
                                  const double         Tmax,
                                  const double         Z,
                                  const double         A,
                                  const JFormFactor_t& F,
                                  const int            N = 1000000)
    {
      if (Tmax > Tmin) {
 
        const double gamma = E / M;
        const double beta  = getBeta(gamma);
 
        const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)) * 1.0e-3;        // [GeV g^-1 cm^2]
 
        const double xmin = 1.0/Tmax;
        const double xmax = 1.0/Tmin;
        const double dx   = (xmax - xmin) / ((double) N);
 
        double count = 0.0;
 
        for (double x = xmin; x <= xmax; x += dx) {
 
          const double T = 1.0/x;
          const double y = W * F(T) * dx;                                     // [g^-1 cm^2]
 
          count += y;
        }
 
        return count;                                                         // [g^-1 cm^2]
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get equivalent EM-shower energy loss due to delta-rays per unit track length\n
     * for an ionising particle with given energy and given mass and for a given form factor.
     *
     * The template parameter corresponds to a class which contains an `operator()(const double)`\n
     * to compute the form factor corresponding to a given delta-ray kinetic energy.
     *
     * \param  E      particle energy                                     [GeV]
     * \param  M      particle mass                                       [GeV]
     * \param  Tmin   minimum delta-ray kinetic energy                    [GeV]
     * \param  Tmax   maximum delta-ray kinetic energy                    [GeV]
     * \param  Z      atomic number                                       [unit]
     * \param  A      atomic mass                                         [g/mol]
     * \param  F      form factor functor
     * \param  N      number of points for numeric integration
     * \return        equivalent energy loss                              [GeV g^-1 cm^2]
     */
    template<class JFormFactor_t>
    static inline double getEnergyLoss(const double         E,
                                       const double         M,
                                       const double         Tmin,
                                       const double         Tmax,
                                       const double         Z,
                                       const double         A,
                                       const JFormFactor_t& F,
                                       const int            N = 1000000)
    {
      if (Tmax > Tmin) {
 
        const double gamma = E / M;
        const double beta  = getBeta(gamma);
 
        const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta));                 // [MeV g^-1 cm^2]
 
        const double xmin = log(Tmin);
        const double xmax = log(Tmax);
        const double dx   = (xmax - xmin) / ((double) N);
        
        double weight = 0.0;
 
        for (double x = xmin; x <= xmax; x += dx) {
 
          const double T = exp(x);
          const double y = W * F(T) * dx;                                     // [MeV g^-1 cm^2]
          
          weight += y;
        }
 
        return weight * 1.0e-3;                                               // [GeV g^-1 cm^2]
        
      } else {
 
        return 0.0;
      }
    }
 
    
    /**
     * Get number of delta-rays per unit track length for an ionising particle with given energy and given mass.
     *
     * N.B: For this function a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$ is assumed, where:
     *      - \f$T\f$ corresponds to the delta-ray kinetic energy,
     *      - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,
     *      - \f$E\f$ to the ionising particle's energy and
     *      - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
     *
     * \param  E      particle energy                  [GeV]
     * \param  M      particle mass                    [GeV]
     * \param  Tmin   minimum delta-ray kinetic energy [GeV]
     * \param  Tmax   maximum delta-ray kinetic energy [GeV]
     * \param  Z      atomic number                    [unit]
     * \param  A      atomic mass                      [g/mol]
     * \return        delta-ray count                  [g^-1 cm^2]
     */
    static inline double getCount(const double E,
                                  const double M,
                                  const double Tmin,
                                  const double Tmax,
                                  const double Z,
                                  const double A)
    {
      if (Tmax > Tmin) {
 
        const double gamma = E / M;
        const double beta  = getBeta(gamma);
 
        const double a = 0.5/(E*E);
        const double b = beta*beta/Tmax;
        const double c = 1.0;
 
        const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta));                 // [MeV g^-1 cm^2]
 
        const double dT = Tmax - Tmin;
        const double rT = Tmax / Tmin;
        const double mT = Tmax * Tmin;
 
        return W * (a*dT - b*log(rT) + c*dT/mT) * 1.0e-3;                     // [g^-1 cm^2]
 
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Get equivalent EM-shower energy loss due to delta-rays per unit track length\n
     * for an ionising particle with given energy and given mass.\n\n
     *
     * N.B: For this function a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$ is assumed, where:
     *      - \f$T\f$ corresponds to the delta-ray kinetic energy,
     *      - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,
     *      - \f$E\f$ to the ionising particle's energy and
     *      - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
     *
     * \param  E      particle energy                  [GeV]
     * \param  M      particle mass                    [GeV]
     * \param  Tmin   minimum delta-ray kinetic energy [GeV]
     * \param  Tmax   maximum delta-ray kinetic energy [GeV]
     * \param  Z      atomic number                    [unit]
     * \param  A      atomic mass                      [g/mol]
     * \return        equivalent energy loss           [GeV g^-1 cm^2]
     */
    static inline double getEnergyLoss(const double E,
                                       const double M,
                                       const double Tmin,
                                       const double Tmax,
                                       const double Z,
                                       const double A)
    {      
      if (Tmax > Tmin) {
 
        const double gamma = E / M;
        const double beta  = getBeta(gamma);
      
        const double a = 0.25/(E*E);
        const double b = beta*beta/Tmax;
        const double c = 1.0;
 
        const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta));                 // [MeV g^-1 cm^2]
 
        const double sT = Tmax + Tmin;
        const double dT = Tmax - Tmin;
        const double rT = Tmax / Tmin;
        
        return W * (a*sT*dT - b*dT + c*log(rT)) * 1.0e-3;                     // [GeV g^-1 cm^2]
        
      } else {
 
        return 0.0;
      }
    }
 
    
    /**
     * Get number of delta-rays per unit track length for an ionising particle with given energy and given mass in sea water.
     *
     * \param  E      particle energy                  [GeV]
     * \param  M      particle mass                    [GeV]
     * \param  Tmin   minimum delta-ray kinetic energy [GeV]
     * \param  Tmax   maximum delta-ray kinetic energy [GeV]
     * \return        delta-ray count                  [g^-1 cm^2]
     */
    static inline double getCount(const double E,
                                  const double M,
                                  const double Tmin,
                                  const double Tmax)
    {
      return getCount(E, M, Tmin, Tmax, 1.0, 1.0) * getZoverA();
    }
 
 
    /**
     * Get equivalent EM-shower energy loss due to delta-rays per unit track length\n
     * for an ionising particle with given energy and given mass in sea water.
     *
     * \param  E      particle energy                  [GeV]
     * \param  M      particle mass                    [GeV]
     * \param  Tmin   minimum delta-ray kinetic energy [GeV]
     * \param  Tmax   maximum delta-ray kinetic energy [GeV]
     * \return        equivalent energy loss           [GeV g^-1 cm^2]
     */
    static inline double getEnergyLoss(const double E,
                                       const double M,
                                       const double Tmin,
                                       const double Tmax)
    {
      return getEnergyLoss(E, M, Tmin, Tmax, 1.0, 1.0) * getZoverA();
    }
 
 
    /**
     * Equivalent EM-shower energy loss due to delta-rays per unit electron track length in sea water.
     *
     * \param  E       electron energy                            [GeV]
     * \param  T_GeV   allowed kinetic energy range of delta-rays [GeV]
     * \return         equivalent energy loss                     [GeV/m]
     */  
    static inline double getEnergylossFromElectron(const double       E,
                                                   const JEnergyRange T_GeV = JEnergyRange(TMIN_GEV, TMAX_GEV))
    {
      const double Tmin = T_GeV.constrain(getTmin());
      const double Tmax = T_GeV.constrain(0.5 * getKineticEnergy(E, MASS_ELECTRON));
    
      return getEnergyLoss(E, MASS_ELECTRON, Tmin, Tmax) * DENSITY_SEA_WATER * 1.0e2;
    }
  
  
    /**
     * Equivalent EM-shower energy loss due to delta-rays per unit muon track length in sea water.
     *
     * \param  E       muon energy                                [GeV]
     * \param  T_GeV   allowed kinetic energy range of delta-rays [GeV]
     * \return         equivalent energy loss                     [GeV/m]
     */  
    static inline double getEnergyLossFromMuon(const double       E,
                                               const JEnergyRange T_GeV = JEnergyRange(TMIN_GEV, TMAX_GEV))
    {
      const double Tmin = T_GeV.constrain(getTmin());
      const double Tmax = T_GeV.constrain(getTmax(E, MASS_MUON));
    
      return getEnergyLoss(E, MASS_MUON, Tmin, Tmax) * DENSITY_SEA_WATER * 1.0e2;
    }
 
  
    /**
     * Equivalent EM-shower energy loss due to delta-rays per unit tau track length in sea water.
     *
     * \param  E       tau energy                                 [GeV]
     * \param  T_GeV   allowed kinetic energy range of delta-rays [GeV]
     * \return         equivalent energy loss                     [GeV/m]
     */  
    static inline double getEnergyLossFromTau(const double       E,
                                              const JEnergyRange T_GeV = JEnergyRange(TMIN_GEV, TMAX_GEV))
    {
      const double Tmin = T_GeV.constrain(getTmin());
      const double Tmax = T_GeV.constrain(getTmax(E, MASS_TAU));
    
      return getEnergyLoss(E, MASS_TAU, Tmin, Tmax) * DENSITY_SEA_WATER * 1.0e2;
    }
 
 
    /**
     * Auxiliary data structure to encapsulate energy loss methods based on fit.
     */
    struct FIT {
 
      /**
       * Equivalent EM-shower energy loss due to delta-rays per unit muon track length.
       *
       * \param  E       muon energy                                [GeV]
       * \return         equivalent energy loss                     [GeV/m]
       */  
      static inline double getEnergyLossFromMuon(const double E)
      {
        static const double a =  3.195e-01;
        static const double b =  3.373e-01;
        static const double c = -2.731e-02;
        static const double d =  1.610e-03;
        static const double Emin = 0.13078; // [GeV]
     
        if (E > Emin) {
      
          const double x = log10(E);
          const double y = a + x*(b + x*(c + x*(d)));                           // [MeV g^-1 cm^2]
      
          if (y > 0.0) {
            return y * DENSITY_SEA_WATER * 1.0e-1;                              // [GeV/m]
          }
        }
    
        return 0.0;
      }
 
 
      /**
       * Equivalent EM-shower energy loss due to delta-rays per unit tau track length.
       *
       * \param  E       tau energy                                 [GeV]
       * \return         equivalent energy loss                     [GeV/m]
       */  
      static inline double getEnergyLossFromTau(const double E)
      {
        static const double a = -2.178e-01;
        static const double b =  4.995e-01;
        static const double c = -3.906e-02;
        static const double d =  1.615e-03;
        static const double Emin = 2.19500; // [GeV]
       
        if (E > Emin) {
   
          const double x = log10(E);
          const double y = a + x*(b + x*(c + x*(d)));                           // [MeV g^-1 cm^2]
   
          if (y > 0) {
            return y * DENSITY_SEA_WATER * 1.0e-1;                              // [GeV/m]
          }
        }
   
        return 0.0;
      }
    };
 
 
    /**
     * Emission profile of photons from delta-rays.
     *
     * Profile is taken from reference ANTARES-SOFT-2002-015, J.\ Brunner (fig.\ 3).
     *
     * \param  x      cosine emission angle
     * \return        probability
     */
    static inline double getProbability(const double x)
    {
      //return 1.0 / (4.0 * PI);
      return 0.188 * exp(-1.25 * pow(fabs(x - 0.90), 1.30));
    }
  };
}
 
#endif