JHondaFluxInterpolator.hh

Go to the documentation of this file.

#ifndef __JAANET__JHONDAFLUXINTERPOLATOR__
#define __JAANET__JHONDAFLUXINTERPOLATOR__
 
#include <iostream>
#include <string>
#include <regex>
 
#include "Jeep/JPrint.hh"
#include "Jeep/JMessage.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JException.hh"
#include "JLang/JObjectStreamIO.hh"
#include "JLang/JLoadProperty.hh"
 
#include "JTools/JRange.hh"
#include "JTools/JPolint.hh"
#include "JTools/JMapList.hh"
#include "JTools/JCollection.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JFunctionalMap_t.hh"
 
#include "JAAnet/JDiffuseFlux.hh"
 
 
/**
 * \author bjung
 */
 
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JEEP::JMessage;
 
  using JLANG::JClonable;
  using JLANG::JObjectStreamIO;
  using JLANG::JLoadProperty;
 
  using JTOOLS::JRange;
  using JTOOLS::JArray;
  
  using JTOOLS::JMAPLIST;
 
  using JTOOLS::JMultiFunction;
  using JTOOLS::JConstantFunction1D;
  using JTOOLS::JPolint1FunctionalMap;
 
  
  static const char* const HONDA_ENERGY_KEYWORD                         = "Enu(GeV)";
  static const char* const HONDA_COSINE_ZENITH_ANGLE_KEYWORD            = "cosZ";
  static const char* const HONDA_AZIMUTH_ANGLE_KEYWORD                  = "phi_Az";
  
  static const char* const HONDA_MUON_NEUTRINO_FLUX_KEYWORD             = "NuMu";
  static const char* const HONDA_MUON_ANTINEUTRINO_FLUX_KEYWORD         = "NuMubar";
  static const char* const HONDA_ELECTRON_NEUTRINO_FLUX_KEYWORD         = "NuE";
  static const char* const HONDA_ELECTRON_ANTINEUTRINO_FLUX_KEYWORD     = "NuEbar";
 
  static const char* const HONDA_BINSPECS_SEPARATOR                     = "--";
  static const char        HONDA_BINSPECS_BEGIN                         = '[';
  static const char        HONDA_BINSPECS_END                           = ']';
  
 
  /**
   * Auxiliary data structure for reading Honda bin ranges.
   */
  struct JHondaBinRange :
    JRange<double>
  {
    typedef JRange<double>        JRange_t;
 
        
    /**
     * Default constructor.
     */
    JHondaBinRange()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  minimum            minimum value
     * \param  maximum            maximum value
     */
    JHondaBinRange(const double minimum,
                   const double maximum) :
      JRange_t(minimum, maximum)
    {}
 
 
    /**
     * Get central value of this range.
     *
     * \return                    central value
     */
    double getCentralValue() const
    {
      return (getLowerLimit() + getUpperLimit()) / 2.0;
    }
 
        
    /**
     * Read bin range from input.
     *
     * \param  in                 input stream
     * \param  object             bin range
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JHondaBinRange& object)
    {
      using namespace std;
        
      double min = JRange_t::getMinimum();
      double max = JRange_t::getMaximum();
 
      string buffer;
          
      if (in >> min >> buffer >> max &&
          min    <= max              &&
          buffer == HONDA_BINSPECS_SEPARATOR) {
        static_cast<JRange_t&>(object) = JRange_t(min, max);
      }
 
      return in;
    }
  };
      
 
  /**
   * Auxiliary data structure for reading angular binning specifications of Honda flux table.
   */
  struct JHondaAngularBinSpecs
  {
    /**
     * Default constructor.
     */
    JHondaAngularBinSpecs() :
      coszRange(-1.0,   1.0),
      phiRange ( 0.0, 360.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  minCosz            minimum cosine zenith angle
     * \param  maxCosz            maximum cosine zenith angle
     * \param  minPhi             minimum azimuth angle
     * \param  maxPhi             maximum azimuth angle
     */
    JHondaAngularBinSpecs(const double minCosz,
                          const double maxCosz,
                          const double minPhi,
                          const double maxPhi) :
      coszRange(minCosz, maxCosz),
      phiRange (minPhi,  maxPhi)
    {}
    
 
    /**
     * Read bin specifications from input.
     *
     * \param  in                 input stream
     * \param  binspecs           bin specifications
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JHondaAngularBinSpecs& binspecs)
    {
      static const JEquationParameters eqpars("=", ",", "./", "#");
          
      JProperties properties(eqpars, 1);
          
      properties[HONDA_COSINE_ZENITH_ANGLE_KEYWORD] = binspecs.coszRange;
      properties[HONDA_AZIMUTH_ANGLE_KEYWORD]       = binspecs.phiRange;
 
      in >> properties;
 
      return in >> properties;
    }
 
 
    JHondaBinRange coszRange; //!< Cosine zenith-angle range
    JHondaBinRange phiRange;  //!< Azimuthal angle range [deg]
  };
    
  
  /** 
   * Template definition for Honda flux table interpolator.
   */
  template<class                               JPhiFunction_t         = JConstantFunction1D<double, JArray<4, double> >,
           template<class, class, class> class JCoszFunctionalMap_t   = JPolint1FunctionalMap,
           template<class, class, class> class JEnergyFunctionalMap_t = JPolint1FunctionalMap>
  class JHondaFluxInterpolator final :
    public JMultiFunction                          <JPhiFunction_t, typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist>,
    public JClonable<JFlux, JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
    public JObjectStreamIO <JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
    public JMessage        <JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
  {
  public:
 
    typedef typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist                                      JFunctionalMaplist_t;
    
    typedef JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t>                             interpolator_type;
    typedef JMultiFunction<JPhiFunction_t, JFunctionalMaplist_t>                                                            multifunction_type;
    
    enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS };
 
    typedef typename multifunction_type::abscissa_type                                                                           abscissa_type;
    typedef typename multifunction_type::argument_type                                                                           argument_type;
    typedef typename multifunction_type::result_type                                                                               result_type;
    typedef typename multifunction_type::value_type                                                                                 value_type;
 
    typedef typename multifunction_type::multimap_type                                                                           multimap_type;
    
    typedef typename multifunction_type::super_const_iterator                                                             super_const_iterator;
    typedef typename multifunction_type::super_iterator                                                                         super_iterator;
    typedef typename multifunction_type::function_type                                                                           function_type;
 
    typedef JRange<abscissa_type>                                                                                                     JRange_t;
    
    
    /**
     * Default constructor.
     */
    JHondaFluxInterpolator() :
      table    (*this),
      coszRange(0.0, -1.0),
      phiRange (0.0, -1.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  filename             Honda flux table filename
     */
    JHondaFluxInterpolator(const char* filename) :
      table    (*this, filename),
      coszRange(0.0, -1.0),
      phiRange (0.0, -1.0)
    {
      load();
    }
 
 
    /**
     * Copy constructor.
     *
     * \param  interpolator    Honda flux interpolator
     */
    JHondaFluxInterpolator(const interpolator_type& interpolator) :
      table    (*this, interpolator.table.c_str()),
      coszRange(interpolator.coszRange)
    {
      this->getMultiFunction() = interpolator.getMultiFunction();
 
      this->compile();
    }
 
 
    /**
     * Check whether this interpolator is valid.
     *
     * \return                   true if valid; else false
     */
    bool is_valid() const override final
    {
      return (!table.empty() ?
              this->size() > 0 && coszRange.getLength() > 0.0 && phiRange.getLength() > 0.0 :
              true);
    }
 
    
    /**
     * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \param  phi               azimuth angle (defined clockwise w.r.t. the North) [deg]
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const int    type,
                     const double log10E,
                     const double costh,
                     const double phi) const
    {
      using namespace std;
      using namespace JPP;
 
      static const double epsilon_costh = 1.0e-3; //!< [--]
      static const double epsilon_phi   = 1.0;    //!< [deg]
 
      static const JRange_t mod(0.0, 360.0);
      
      const multifunction_type& function = static_cast<const multifunction_type&>(*this);
 
      const double _costh = min(max(costh, -1.0), 1.0);
      const double _phi   = mod.mod(180.0 - phi); //!< Honda azimuth angle is defined counterclockwise w.r.t. the South\n
                                                  //!< (c.f. arXiv:1102.2688), so a conversion is necessary.
      
      result_type lnFluxes;
      
      if (coszRange(_costh) && phiRange(_phi)) {
        
        lnFluxes = function(log10E, _costh, _phi);
        
      } else { // Perform linear extrapolation for edges of angular ranges
        
        if        (_costh > coszRange.getUpperLimit() && phiRange(_phi)) {
 
          const abscissa_type x1 = coszRange.getUpperLimit();
          const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh;
            
          const result_type&  y1 = function(log10E, x1, _phi);  
          const result_type&  y2 = function(log10E, x2, _phi);
 
          lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
          
        } else if (_costh < coszRange.getLowerLimit() && phiRange(_phi)) {
        
          const abscissa_type x1 = coszRange.getLowerLimit();
          const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh;
            
          const result_type&  y1 = function(log10E, x1, _phi);  
          const result_type&  y2 = function(log10E, x2, _phi);
 
          lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
        
        } else if (coszRange(_costh) && _phi > phiRange.getUpperLimit()) {
 
          const abscissa_type x1 = phiRange.getUpperLimit();
          const abscissa_type x2 = phiRange.getUpperLimit() - epsilon_costh;
            
          const result_type&  y1 = function(log10E, _costh, x1);        
          const result_type&  y2 = function(log10E, _costh, x2);
 
          lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1);
        
        } else if (coszRange(_costh) && _phi < phiRange.getLowerLimit()) {
 
          const abscissa_type x1 = phiRange.getLowerLimit();
          const abscissa_type x2 = phiRange.getLowerLimit() + epsilon_costh;
            
          const result_type&  y1 = function(log10E, _costh, x1);        
          const result_type&  y2 = function(log10E, _costh, x2);
 
          lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1);
        
        } else if (_costh > coszRange.getUpperLimit() && _phi > phiRange.getUpperLimit()) {
 
          const abscissa_type c1 = coszRange.getUpperLimit();
          const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh;
          
          const abscissa_type p1 = phiRange .getUpperLimit();
          const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi;
 
          const result_type&  y1 = function(log10E, c1, p1);
          const result_type&  y2 = function(log10E, c2, p2);
 
          lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
          
        } else if (_costh > coszRange.getUpperLimit() && _phi < phiRange.getUpperLimit()) {
 
          const abscissa_type c1 = coszRange.getUpperLimit();
          const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh;
          
          const abscissa_type p1 = phiRange .getLowerLimit();
          const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi;
 
          const result_type&  y1 = function(log10E, c1, p1);
          const result_type&  y2 = function(log10E, c2, p2);
 
          lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
          
        } else if (_costh < coszRange.getLowerLimit() && phi > phiRange.getUpperLimit()) {
 
          const abscissa_type c1 = coszRange.getLowerLimit();
          const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh;
          
          const abscissa_type p1 = phiRange .getUpperLimit();
          const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi;
 
          const result_type&  y1 = function(log10E, c1, p1);
          const result_type&  y2 = function(log10E, c2, p2);
 
          lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
          
        } else {
 
          const abscissa_type c1 = coszRange.getLowerLimit();
          const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh;
          
          const abscissa_type p1 = phiRange .getLowerLimit();
          const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi;
 
          const result_type&  y1 = function(log10E, c1, p1);
          const result_type&  y2 = function(log10E, c2, p2);
 
          lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
        }
      }
 
      switch(type) {
      case TRACK_TYPE_NUMU:
        return exp(lnFluxes[0]);
      case TRACK_TYPE_ANTINUMU:
        return exp(lnFluxes[1]);
      case TRACK_TYPE_NUE:
        return exp(lnFluxes[2]);
      case TRACK_TYPE_ANTINUE:
        return exp(lnFluxes[3]);
      default:
        return 0.0;
      }
    }
 
 
    /**
     * Get flux of given event.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \param  phi               azimuth angle (defined clockwise w.r.t. the North) [deg]
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFlux(const int    type,
                   const double log10E,
                   const double costh,
                   const double phi) const
    {
      return getFactor(type, log10E, costh, phi);
    }
 
 
    /**
     * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \param  phi               azimuth angle (defined clockwise w.r.t. the North) [deg]
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double operator()(const int    type,
                      const double log10E,
                      const double costh,
                      const double phi) const
    {
      return getFactor(type, log10E, costh, phi);
    }
 
 
    /**
     * Get flux of given event.
     *
     * \param  evt               event
     * \return                   flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const Evt& evt) const override final
    {
      const Trk& neutrino = get_neutrino(evt);
      
      const double log10E = log10(neutrino.E);
      const double costh  = neutrino.dir.z;
      
      //!< The gSeaGen azimuth angle is defined clockwise w.r.t. the North (c.f. arXiv:2003.14040)
      const double phi = atan(neutrino.dir.y / neutrino.dir.x) * 180.0 / M_PI;
 
      return getFactor(neutrino.type, log10E, costh, phi);
    }
 
 
    /**
     * Get flux of given event.
     *
     * \param  event             event
     * \return                   flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFlux(const Evt& event) const
    {
      return getFactor(event);
    }
 
 
    /**
     * Load Honda flux table.
     */
    void load()
    {
      using namespace std;
      using namespace JPP;
      
      if (!table.empty()) {
 
        NOTICE("Loading Honda flux table from file " << table << "... " << flush);
        
        JObjectStreamIO<interpolator_type>::load(table.c_str());
        
        this->check_validity();
 
        NOTICE("OK" << endl);
      }
    }
 
 
    /**
     * Load Honda flux table.
     *
     * \param  table           Honda flux table filename
     */
    void load(const char* table)
    {
      this->table = table;
 
      load();
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }
    
 
    /**
     * Read Honda atmospheric neutrino flux interpolator from file.
     *
     * \param  ifs             input file stream
     * \param  interpolator    Honda flux interpolator
     * \return                 input file stream 
     */
    friend std::ifstream& operator>>(std::ifstream& ifs, interpolator_type& interpolator)
    {
      using namespace std;
      
      multimap_type& zmap = static_cast<multimap_type&>(interpolator);
      
      JHondaAngularBinSpecs binspecs;
        
      for (string line; getline(ifs, line); ) {
 
        double energy = 0.0;
        
        result_type values;
 
        istringstream iss1(line);
 
        const double cosz = binspecs.coszRange.getCentralValue();
        const double phi  = binspecs.phiRange .getCentralValue();
 
        interpolator.coszRange.include(cosz);
        interpolator.phiRange .include(phi);
          
        if (iss1 >> energy >> values) {
            
          const double log10E = log10(energy);
 
          for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) {
            *i = log(*i);
          }
 
          zmap[log10E][cosz][phi] = values;
            
        } else {
 
          const size_t p1 = line.find (HONDA_BINSPECS_BEGIN);
          const size_t p2 = line.rfind(HONDA_BINSPECS_END);
        
          if (p1 != string::npos && p2 != string::npos) {
 
            const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1);
            
            istringstream iss2(binspecstr);
            
            iss2 >> binspecs;
          }
        }
      }
 
      interpolator.compile();
 
      return ifs;
    }
 
 
  private:
    
    /**
     * Auxiliary class for I/O of Honda flux interpolator.
     */
    struct JHondaFluxInterpolatorHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  interpolator  Honda flux interpolator
       * \param  eqpars        equation parameters 
       */
      template<class JHondaFluxInterpolator_t>
      JHondaFluxInterpolatorHelper(JHondaFluxInterpolator_t&  interpolator,
                                   const JEquationParameters& eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator";
 
        this->insert(gmake_property(interpolator.table));
      }
    };
 
    
    JLoadProperty<interpolator_type, std::string> table;  //!< Table filename
    
    JRange_t coszRange; //!< Zenith  angle range
    JRange_t phiRange;  //!< Azimuth angle range [deg]
    
    using JMessage<interpolator_type>::debug; //!< Debug level
  };
 
 
  /** 
   * Template specialisation for interpolator of azimuth-averaged Honda flux table.
   */
  template<template<class, class, class> class JCoszFunctionalMap_t,
           template<class, class, class> class JEnergyFunctionalMap_t>
  class JHondaFluxInterpolator<JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> final :
    public JMultiFunction                                 <JConstantFunction1D<double, JArray<4, double> >, typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist>,
    public JClonable<JDiffuseFlux, JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
    public JObjectStreamIO        <JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
    public JMessage               <JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
  {
  public:
    
    typedef typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist                                      JFunctionalMaplist_t;
    
    typedef JConstantFunction1D<double, JArray<4, double> >                                                                     JPhiFunction_t;
    
    typedef JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t>                             interpolator_type;
    typedef JMultiFunction<JPhiFunction_t, JFunctionalMaplist_t>                                                            multifunction_type;
    
    enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS };
 
    typedef typename multifunction_type::abscissa_type                                                                           abscissa_type;
    typedef typename multifunction_type::argument_type                                                                           argument_type;
    typedef typename multifunction_type::result_type                                                                               result_type;
    typedef typename multifunction_type::value_type                                                                                 value_type;
 
    typedef typename multifunction_type::multimap_type                                                                           multimap_type;
 
    typedef typename multifunction_type::super_const_iterator                                                             super_const_iterator;
    typedef typename multifunction_type::super_iterator                                                                         super_iterator;
    typedef typename multifunction_type::function_type                                                                           function_type;
 
    typedef JRange<abscissa_type>                                                                                                     JRange_t;
    
    
    /**
     * Default constructor.
     */
    JHondaFluxInterpolator() :
      table    (*this),
      coszRange(0.0, -1.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  filename        Honda flux table filename
     */
    JHondaFluxInterpolator(const char* filename) :
      table    (*this, filename),
      coszRange(0.0, -1.0)
    {
      load();
    }
 
 
    /**
     * Copy constructor.
     *
     * \param  interpolator    Honda flux interpolator
     */
    JHondaFluxInterpolator(const interpolator_type& interpolator) :
      table    (*this, interpolator.table.c_str()),
      coszRange(interpolator.coszRange)
    {
      this->getMultiFunction() = interpolator.getMultiFunction();
 
      this->compile();
    }
 
 
    /**
     * Check whether this interpolator is valid.
     *
     * \return                 true if valid; else false
     */
    bool is_valid() const override final
    {
      return (!table.empty() ? this->size() > 0 && coszRange.getLength() : true);
    }
    
    
    /**
     * Get diffuse flux for given particle PDG-identifier, energy and zenith-angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double dNdEdOmega(const int    type,
                      const double log10E,
                      const double costh) const override final
    {
      using namespace std;
      using namespace JPP;
 
      static const double epsilon_costh = 1.0e-3;
      
      const multifunction_type& function = static_cast<const multifunction_type&>(*this);
 
      const double _costh = min(max(costh, -1.0), 1.0);
 
      result_type lnFluxes;
      
      if (coszRange(_costh)) {
        
        lnFluxes = function(log10E, _costh);
        
      } else { // Perform linear extrapolation for edges of cosine zenith angle range
        
        if        (_costh > coszRange.getUpperLimit()) {
 
          const abscissa_type x1 = coszRange.getUpperLimit();
          const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh;
            
          const result_type&  y1 = function(log10E, x1);        
          const result_type&  y2 = function(log10E, x2);
 
          lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
          
        } else {
        
          const abscissa_type x1 = coszRange.getLowerLimit();
          const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh;
            
          const result_type&  y1 = function(log10E, x1);        
          const result_type&  y2 = function(log10E, x2);
 
          lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
        }
      }
 
      switch(type) {
      case TRACK_TYPE_NUMU:
        return exp(lnFluxes[0]);
      case TRACK_TYPE_ANTINUMU:
        return exp(lnFluxes[1]);
      case TRACK_TYPE_NUE:
        return exp(lnFluxes[2]);
      case TRACK_TYPE_ANTINUE:
        return exp(lnFluxes[3]);
      default:
        return 0.0;
      }
    }
 
 
    /**
     * Load Honda flux table.
     */
    void load()
    {
      using namespace std;
      using namespace JPP;
      
      if (!table.empty()) {
        
        NOTICE("Loading Honda flux table from file " << table << "... " << flush);
        
        JObjectStreamIO<interpolator_type>::load(table.c_str());
        
        this->check_validity();
 
        NOTICE("OK" << endl);   
      }
    }
 
 
    /**
     * Load Honda flux table.
     *
     * \param  table           Honda flux table filename
     */
    void load(const char* table)
    {
      this->table = table;
 
      load();
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }
 
 
    /**
     * Read Honda atmospheric neutrino flux interpolator from file.
     *
     * \param  ifs             input file stream
     * \param  interpolator    Honda flux interpolator
     * \return                 input file stream 
     */
    friend std::ifstream& operator>>(std::ifstream& ifs, interpolator_type& interpolator)
    {
      using namespace std;
      
      multimap_type& zmap = static_cast<multimap_type&>(interpolator);
        
      JHondaAngularBinSpecs binspecs;
      
      for (string line; getline(ifs, line); ) {
 
        double energy = 0.0;
        
        result_type values;
 
        istringstream iss1(line);
 
        const double cosz = binspecs.coszRange.getCentralValue();
        
        interpolator.coszRange.include(cosz);
 
        if (iss1 >> energy >> values) {
 
          const double log10E = log10(energy);
          
          for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) {
            *i = log(*i);
          }
 
          zmap[log10E][cosz] = values;
            
        } else {
 
          const size_t p1 = line.find (HONDA_BINSPECS_BEGIN);
          const size_t p2 = line.rfind(HONDA_BINSPECS_END);
        
          if (p1 != string::npos && p2 != string::npos) {
 
            const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1);
            
            istringstream iss2(binspecstr);
            
            iss2 >> binspecs;
          }
        }
      }
 
      interpolator.compile();
 
      return ifs;
    }
    
    
  private:
    
    /**
     * Auxiliary class for I/O of Honda flux interpolator.
     */
    struct JHondaFluxInterpolatorHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  interpolator  Honda flux interpolator
       * \param  eqpars        equation parameters 
       */
      template<class JHondaFluxInterpolator_t>
      JHondaFluxInterpolatorHelper(JHondaFluxInterpolator_t&  interpolator,
                                   const JEquationParameters& eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator";
 
        this->insert(gmake_property(interpolator.table));
      }
    };
 
    
    JLoadProperty<interpolator_type, std::string> table; //!< Honda flux table filename
    
    JRange_t coszRange; //!< Cosine zenith angle range
    
    using JMessage<interpolator_type>::debug; //!< Debug level
  };
 
 
  /**
   * Alias template definition for 2-dimensional Honda flux interpolator.
   */
  template<template<class, class, class> class JCoszFunctionalMap_t   = JPolint1FunctionalMap,
           template<class, class, class> class JEnergyFunctionalMap_t = JPolint1FunctionalMap>
  using JHondaFluxInterpolator2D = JHondaFluxInterpolator<JConstantFunction1D<double, JArray<4, double> >,
                                                          JCoszFunctionalMap_t,
                                                          JEnergyFunctionalMap_t>;
}
 
 
namespace JEEP {
 
  /**
   * JMessage template specialization for Honda atmospheric neutrino flux interpolators.
   */
  template<class                               JPhiFunction_t,
           template<class, class, class> class JCoszFunctionalMap_t,
           template<class, class, class> class JEnergyFunctionalMap_t>  
  struct JMessage<JAANET::JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
  {
    static int debug;
  };
 
 
  /**
   * Default verbosity for Honda atmospheric neutrino flux interpolators.
   */
  template<class                               JPhiFunction_t,
           template<class, class, class> class JCoszFunctionalMap_t,
           template<class, class, class> class JEnergyFunctionalMap_t>  
  int JMessage<JAANET::JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >::debug = (int) notice_t;
  
}
 
#endif