JRadiationFunction.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JRADIATIONFUNCTION__
#define __JPHYSICS__JRADIATIONFUNCTION__

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JGrid.hh"

#include "JPhysics/JRadiation.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  /**
   * Fast implementation of class JRadiation.
   *
   * In this, the methods
   *   - JRadiation::TotalCrossSectionEErad,
   *   - JRadiation::TotalCrossSectionGNrad,
   *   - JRadiation::CalculateACoeff and
   *   - JRadiation::IntegralofG
   *
   * are reimplemented using lookup tables.
   */
  class JRadiationFunction : 
    public JRadiation
  {
    typedef JTOOLS::JGridPolint1Function1D_t                                                       JFunction1D_t;
    typedef JTOOLS::JMultiFunction<JFunction1D_t, 
                                   JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap> >           JFunction2D_t;

  public:
    /**
     * Constructor.
     *
     * \param  radiation       JRadiation object
     * \param  number_of_bins  number of bins
     * \param  Emin            minimal muon energy [GeV]
     * \param  Emax            maximal muon energy [GeV]
     */
    JRadiationFunction(const JRadiation&  radiation,
                       const unsigned int number_of_bins,
                       const double       Emin,
                       const double       Emax) :
      JRadiation(radiation)
    {
      using namespace JTOOLS;

      const double xmin = log(Emin);
      const double xmax = log(Emax);

    
      integral.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction2D_t::iterator i = integral.begin(); i != integral.end(); ++i) {

        const double x    = i->getX();
        const double E    = exp(x);

        const double ymin = log(EminEErad/E);
        const double ymax = 0.0;

        if (ymax > ymin) {

          i->getY().configure(make_grid(number_of_bins, ymin, ymax));
          
          for (JFunction1D_t::iterator j = i->getY().begin(); j != i->getY().end(); ++j) {
            
            const double y   = j->getX();
            const double eps = exp(y) * E;
            const double z   = JRadiation::IntegralofG(E, eps);

            j->getY() = z;
          }
        }
      }
    
      integral.compile();
      integral.setExceptionHandler(new JFunction1D_t::JDefaultResult(0.0));


      sigmaEE.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction1D_t::iterator i = sigmaEE.begin(); i != sigmaEE.end(); ++i) {

        const double E = exp(i->getX());
        const double y = JRadiation::TotalCrossSectionEErad(E);

        i->getY() = y;
      }
                    
      sigmaEE.compile();


      sigmaGN.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction1D_t::iterator i = sigmaGN.begin(); i != sigmaGN.end(); ++i) {

        const double E = exp(i->getX());
        const double y = JRadiation::TotalCrossSectionGNrad(E);

        i->getY() = y;
      }
                    
      sigmaGN.compile();


      Acoeff.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction1D_t::iterator i = Acoeff.begin(); i != Acoeff.end(); ++i) {
     
        const double E = exp(i->getX());
        const double y = JRadiation::CalculateACoeff(E);

        i->getY() = y;
      }
    
      Acoeff.compile();
    }
       

    /**
     * Pair production cross section.
     *
     * \param  E     muon    energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double TotalCrossSectionEErad(const double E) const override 
    {
      const double x = log(E);

      if (x >= sigmaEE. begin()->getX() &&
          x <= sigmaEE.rbegin()->getX()) {

        try {
          return sigmaEE(x);
        }
        catch(std::exception& error) {}
      }

      return JRadiation::TotalCrossSectionEErad(E);
    }
  

    /**
     * Photo-nuclear cross section.
     *
     * \param  E     muon    energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double TotalCrossSectionGNrad(const double E) const override 
    {
      const double x = log(E);

      if (x >= sigmaGN. begin()->getX() &&
          x <= sigmaGN.rbegin()->getX()) {

        try {
          return sigmaGN(x);
        }
        catch(std::exception& error) {}
      }

      return JRadiation::TotalCrossSectionGNrad(E);
    }


    /**
     * Ionization a parameter.
     *
     * \param  E       muon energy [GeV]
     * \return         ionization coefficient [GeV*m^2/g]
     */
    virtual double CalculateACoeff(const double E) const override 
    {
      const double x = log(E);

      if (x >= Acoeff. begin()->getX() &&
          x <= Acoeff.rbegin()->getX()) {

        try {
          return Acoeff(x);
        }
        catch(std::exception& error) {}
      }

      return JRadiation::CalculateACoeff(E);
    }

  protected:
    virtual double IntegralofG(const double E,
                               const double eps) const override
    {
      try {
        
        const double x = log(E); 
        const double y = log(eps/E); 
        const double z = integral(x, y);

        return z;
      }
      catch(std::exception& error) {}

      return JRadiation::IntegralofG(E, eps);
    }

    JFunction1D_t sigmaEE;
    JFunction1D_t sigmaGN;
    JFunction1D_t Acoeff;
    JFunction2D_t integral;
  };
}

#endif