JOscProbInterpolator.hh

Go to the documentation of this file.

#ifndef __JOSCPROB__JOSCPROBINTERPOLATOR__
#define __JOSCPROB__JOSCPROBINTERPOLATOR__
 
#include "Jeep/JMessage.hh"
 
#include "JIO/JSerialisable.hh"
#include "JIO/JObjectBinaryIO.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JObjectIO.hh"
#include "JLang/JException.hh"
#include "JLang/JLoadProperty.hh"
 
#include "JTools/JPolint.hh"
#include "JTools/JMapList.hh"
#include "JTools/JCollection.hh"
#include "JTools/JFunctional.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JFunctionalMap_t.hh"
 
#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbToolkit.hh"
#include "JOscProb/JOscProbInterface.hh"
#include "JOscProb/JBaselineComputer.hh"
 
 
/**
 * \author bjung, mdejong
 */
 
namespace JOSCPROB {}
namespace JPP { using namespace JOSCPROB; }
 
namespace JOSCPROB {
 
  using JEEP::JMessage;
 
  using JLANG::JClonable;
  using JLANG::JLoadProperty;
 
  using JIO::JSerialisable;
  using JIO::JObjectBinaryIO;
 
  using JTOOLS::JFunctional;
  using JTOOLS::JMultiFunction;
  
 
  /** 
   * Template definition of a multi-dimensional oscillation probability interpolation table.
   */
  template<template<class, class> class JCollection_t = JTOOLS::JCollection,
           class JFunction1D_t                        = JTOOLS::JPolintFunction1D <1,
                                                                                   JTOOLS::JElement2D<double, JTOOLS::JArray<NUMBER_OF_OSCCHANNELS, double> >,
                                                                                   JCollection_t,
                                                                                   JTOOLS::JArray<NUMBER_OF_OSCCHANNELS, double> >,
           class JFunctionalMaplist_t                 = JTOOLS::JMAPLIST          <JTOOLS::JPolint1FunctionalMap,
                                                                                   JTOOLS::JPolint1FunctionalMap,
                                                                                   JTOOLS::JPolint1FunctionalMap,
                                                                                   JTOOLS::JPolint1FunctionalMap,
                                                                                   JTOOLS::JPolint1FunctionalMap,
                                                                                   JTOOLS::JPolint1FunctionalMap,
                                                                                   JTOOLS::JPolint2FunctionalMap>::maplist >
  class JOscProbInterpolator final :
    public JMultiFunction                                                   <JFunction1D_t, JFunctionalMaplist_t>,
    public JClonable<JOscProbInterface, JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >,
    public JMessage                    <JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >,
    public JObjectBinaryIO             <JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >,
    public JSerialisable
  {
  public:
    
    typedef JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t>                interpolator_type;
    typedef JMultiFunction<JFunction1D_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 typename JOscProbInterface::JOscParametersHelper_t                                 JOscParametersHelper_t;    
    typedef typename JOscProbInterface::JOscParameters_t                                             JOscParameters_t;
    typedef typename JOscProbInterface::JOscParameterReferences_t                           JOscParameterReferences_t;
    typedef typename JOscProbInterface::JOscParameter_t                                               JOscParameter_t;
    typedef typename JOscProbInterface::JParameter_t                                                     JParameter_t;
 
    
    /**
     * Default constructor.
     */
    JOscProbInterpolator() :
      table(*this)
    {
      JOscParameters<typename JOscParameterReferences_t::argument_type> parameters(this->buffer);
 
      JOscParametersHelper_t::configure(parameters);
    }
 
 
    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     */
    JOscProbInterpolator(const char* fileName) :
      table(*this, fileName)
    {
      JOscParameters<typename JOscParameterReferences_t::argument_type> parameters(this->buffer);
      
      JOscParametersHelper_t::configure(parameters);
      
      load();
    }
 
 
    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     * \param  parameters           oscillation parameters
     */
    JOscProbInterpolator(const char*             fileName,
                         const JOscParameters_t& parameters) :
      JOscProbInterpolator(fileName)
    {
      JOscParametersHelper_t::set(parameters);
    }
 
 
    /**
     * Copy constructor.
     *
     * \param  interpolator         oscillation probability interpolator
     */
    JOscProbInterpolator(const interpolator_type& interpolator) :
      JOscProbInterpolator(interpolator.table.c_str())
    {
      this->baselineComputer   = interpolator.baselineComputer;
      this->getMultiFunction() = interpolator.getMultiFunction();
 
      this->buffer = interpolator.buffer;
 
      this->compile();
    }
 
 
    /**
     * Get cosine zenith angle for a given baseline.
     *
     * \param  L                    baseline [km]
     * \return                      cosine zenith angle
     */
    double getCosth(const double L) const override final
    {
      return baselineComputer.getCosth(L);
    }
 
 
    /**
     * Get baseline for a given cosine zenith angle.
     *
     * \param  costh                cosine zenith angle
     * \return                      baseline [km]
     */
    double getBaseline(const double costh) const override final
    {
      return baselineComputer.getBaseline(costh);
    }
    
 
    /**
     * Get oscillation probability for a given oscillation channel.
     *
     * \param  channel              oscillation channel
     * \param  E                    neutrino energy [GeV]
     * \param  costh                cosine zenith angle
     * \return                      oscillation probability
     */
    double getP(const JOscChannel&    channel,
                const double          E,
                const double          costh) const override final
    {
      using namespace std;
      using namespace JPP;
 
      const JOscChannel* p = find(getOscChannel, getOscChannel + NUMBER_OF_OSCCHANNELS, channel);
 
      if (p != end(getOscChannel)) {
        
        const double L = this->getBaseline(costh);
 
        this->buffer[NUMBER_OF_DIMENSIONS-2] = L/E;
        this->buffer[NUMBER_OF_DIMENSIONS-1] = costh;
 
        const double* arguments = this->buffer.data();
        
        const size_t       index         = distance(getOscChannel, p);
        const result_type& probabilities = this->evaluate(arguments);
 
        const double& P = probabilities[index];
 
        return min(max(P, 0.0), 1.0);
        
      } else {
        
        THROW(JValueOutOfRange, "JOscProbInterpolator<...>::getP(): Invalid oscillation channel " << channel << endl);
      }
    }
 
 
    /**
     * Load oscillation probability table.
     */
    void load()
    {
      using namespace std;
      using namespace JPP;
        
      if (!table.empty()) {
 
        NOTICE("Loading oscillation probability table from file " << table << "... " << flush);
 
        JObjectBinaryIO<interpolator_type>::load(table.c_str());
        
        NOTICE("OK" << endl);
      }
    }
 
 
    /**
     * Load oscillation probability table.
     *
     * \param  table                oscillation probability table filename
     */
    void load(const char* table)
    {
      this->table = table;
 
      load();
    }
    
    
    /**
     * Read from input.
     *
     * \param  in                   reader
     * \return                      reader
     */
    JReader& read(JReader& in) override final 
    {
      multifunction_type& multifunction = this->getMultiFunction();
 
      in >> baselineComputer;
      in >> multifunction;
 
      multifunction.compile();
 
      return in;
    }
    
    
    /**
     * Write from input.
     *
     * \param  out                  writer
     * \return                      writer
     */
    JWriter& write(JWriter& out) const override final 
    {
      out << baselineComputer;
      out << static_cast<const multifunction_type&>(*this);
 
      return out;
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars               equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JOscProbInterface::getEquationParameters()) override final
    {
      return JOscProbInterpolatorHelper(*this, eqpars);
    }
    
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars               equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JOscProbInterface::getEquationParameters()) const override final
    {
      return JOscProbInterpolatorHelper(*this, eqpars);
    }
 
    
  private:
    
    /**
     * Auxiliary class for I/O of oscillation probability interpolator.
     */
    struct JOscProbInterpolatorHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  interpolator       oscillation probability interpolator
       * \param  eqpars             equation parameters 
       */
      template<class JOscProbInterpolator_t>
      JOscProbInterpolatorHelper(JOscProbInterpolator_t&    interpolator,
                                 const JEquationParameters& eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JOscProbInterface::getTypeKey()] = "interpolator";
        
        this->insert(gmake_property(interpolator.table));
        
        this->join(interpolator.getParameters().getProperties());
      }
    };
 
    
    JBaselineComputer baselineComputer;                         //!< Baseline computer
    
    JLoadProperty<interpolator_type, std::string> table;        //!< Oscillation probability table filename
 
    mutable JArray<NUMBER_OF_DIMENSIONS, argument_type> buffer; //!< Buffer for interpolation arguments
 
    using JMessage<interpolator_type>::debug;                   //!< Debug level
  };
}
 
 
namespace JEEP {
 
  /**
   * JMessage template specialization for oscillation probability interpolators.
   */
  template<template<class, class> class JCollection_t, class JFunction1D_t, class JFunctionalMaplist_t>
  struct JMessage<JOSCPROB::JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >
  {
    static int debug;
  };
 
 
  /**
   * Default verbosity for oscillation probability interpolators.
   */
  template<template<class, class> class JCollection_t, class JFunction1D_t, class JFunctionalMaplist_t>
  int JMessage<JOSCPROB::JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >::debug = (int) notice_t;
  
}
 
#endif