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/JFileStreamIO.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JObjectIO.hh"
#include "JLang/JException.hh"
 
#include "JTools/JPolint.hh"
#include "JTools/JMapList.hh"
#include "JTools/JCollection.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 JIO::JSerialisable;
 
  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 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()
    {
      JOscParameters<typename JOscParameterReferences_t::argument_type> parameters(multifunction_type::buffer);
 
      JOscParametersHelper_t::configure(parameters);
    }
 
 
    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     */
    JOscProbInterpolator(const char* fileName) :
      JOscProbInterpolator()
    {
      this->load(fileName);
    }
 
 
    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     * \param  parameters           oscillation parameters
     */
    JOscProbInterpolator(const char*             fileName,
                         const JOscParameters_t& parameters) :
      JOscProbInterpolator()
    {
      JOscParametersHelper_t::set(parameters);
      
      this->load(fileName);
    }
 
 
    /**
     * Copy constructor.
     *
     * \param  interpolator         oscillation probability interpolator
     */
    JOscProbInterpolator(const interpolator_type& interpolator) :
      JOscProbInterpolator()
    {
      static_cast<multifunction_type&>(*this) = static_cast<const multifunction_type&>(interpolator);
    }
    
 
    /**
     * Load oscillation probability table from file.
     *
     * \param  fileName             oscillation probability table filename
     */
    void load(const char* const fileName)
    {
      using namespace std;
      using namespace JPP;
      
      try {
 
        NOTICE("Loading oscillation probability table from file " << fileName << "... " << flush);
        
        JLANG::load<JIO::JFileStreamReader>(fileName, static_cast<interpolator_type&>(*this));
 
        NOTICE("OK" << endl);
      }
      catch(const JException& error) {
        THROW(JFileReadException, "JOscProbInterpolator::load(): Error reading file " << fileName);
      }
    }
 
 
    /**
     * 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 (P > 1.0 ? 1.0 : (P < 0.0 ? 0.0 : P));
        
      } else {
        
        THROW(JValueOutOfRange, "JOscProbInterpolator<...>::getP(): Invalid oscillation channel " << channel << endl);
      }
    }
    
    
    /**
     * Read from input.
     *
     * \param  in                   reader
     * \return                      reader
     */
    JReader& read(JReader& in) override final 
    {
      in >> this->Lmin >> this->Lmax;
      in >> static_cast<multifunction_type&>(*this);
      
      this->compile();
      
      return in;
    }
    
    
    /**
     * Write from input.
     *
     * \param  out                  writer
     * \return                      writer
     */
    JWriter& write(JWriter& out) const override final 
    {
      out << this->Lmin << this->Lmax;
      out << static_cast<const multifunction_type&>(*this);
 
      return out;
    }
 
    
  private:
 
    using JMessage<interpolator_type>::debug;
  };
}
 
 
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