JRadiationSource.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JRADIATIONSOURCE__
#define __JPHYSICS__JRADIATIONSOURCE__

#include <utility>

#include "JLang/JSharedPointer.hh"
#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(...) and
   *   JRadiation::IntegralofG(...)
   *
   * are replaced by lookup tables.
   */
  class JRadiationFunction : 
    public JRadiation
  {
    //typedef JTOOLS::JGridSplineFunction1D_t                                                        JFunction1D_t;
    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]
     */
    explicit 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();
    }
  

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

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

        try {
          return sigmaEE(x);
        }
        catch(std::exception& error) {
          std::cerr << "JRadiation::TotalCrossSectionEErad() " << E << ' ' << error.what() << std::endl;
        }
      }

      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
    {
      const double x = log(E);

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

        try {
          return sigmaGN(x);
        }
        catch(std::exception& error) {
          std::cerr << "JRadiation::TotalCrossSectionGNrad() " << E << ' ' << error.what() << std::endl;
        }
      }

      return JRadiation::TotalCrossSectionGNrad(E);
    }

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

        return z;
      }
      catch(std::exception& error) {
        std::cerr << "JRadiation::IntegralofG() " << E << ' ' << eps << ' ' << error.what() << std::endl;
      }

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

    JFunction1D_t sigmaEE;
    JFunction1D_t sigmaGN;
    JFunction2D_t integral;
  };
  

  /**
   * Interface for calculation of inverse interaction length and shower energy.
   */
  class JRadiationInterface {
  public:

    /**
     * Virtual destructor.
     */
    virtual ~JRadiationInterface()
    {}


    /**
     * Get inverse interaction length.
     *
     * \param  E      muon energy [GeV]
     * \return        inverse interaction length [m^-1]
     */
    virtual double getInverseInteractionLength(const double E) const = 0;


    /**
     * Get energy of shower.
     *
     * \param  E      muon   energy [GeV]
     * \return        shower energy [GeV]
     */
    virtual double getEnergyOfShower(const double E) const = 0;
  };
  

  /**
   * Source of radiation.
   * This class implements the JRadiationInterface interface.
   * N.B: This class owns the object pointed to using JSharedPointer<>.
   */
  class JRadiationSource : 
    public JRadiationInterface
  {
  public:

    typedef double (JRadiation::*sigma_type)(const double) const;   // total cross section method
    typedef double (JRadiation::*eloss_type)(const double) const;   // energy loss method

    typedef std::pair<sigma_type, eloss_type> source_type;

 
    /**
     * Constructor.
     *
     * \param  radiation       pointer to valid JRadition object
     * \param  density         Mass density of radiation material [gr/cm³]
     * \param  source          pointers to total cross section and energy loss methods
     */
    JRadiationSource(const JLANG::JSharedPointer<JRadiation>& radiation,
                     const double                             density,
                     const source_type                        source) :
      rad  (radiation),
      rho  (density),
      sigma(source.first),
      eloss(source.second)
    {}
    
    
    /**
     * Get inverse interaction length.
     *
     * \param  E      muon energy [GeV]
     * \return        inverse interaction length [m^-1]
     */
    virtual double getInverseInteractionLength(const double E) const 
    { 
      return (*rad.*sigma)(E) * rho * 1.0e6;
    }
    

    /**
     * Get energy of shower.
     *
     * \param  E      muon   energy [GeV]
     * \return        shower energy [GeV]
     */
    virtual double getEnergyOfShower(const double E) const 
    { 
      return (*rad.*eloss)(E);
    }

  protected:
    JLANG::JSharedPointer<JRadiation> rad;
    const double  rho;
    sigma_type    sigma;
    eloss_type    eloss;
  };
}    

#endif