JOscProbInterpolator.hh

Go to the documentation of this file.

#ifndef __JOSCPROB__JOSCPROBINTERPOLATOR__
#define __JOSCPROB__JOSCPROBINTERPOLATOR__

#include "Jeep/JMessage.hh"
#include "Jeep/JProperties.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 :
    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::JOscParameter_t                                               JOscParameter_t;
    typedef typename JOscProbInterface::JParameter_t                                                     JParameter_t;

    
    /**
     * Default constructor.
     */
    JOscProbInterpolator()
    {
      JOscParameters<double&> parameters(multifunction_type::buffer);

      JOscParametersHelper_t::configure(parameters);
    }


    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     */
    JOscProbInterpolator(const char* fileName)
    {
      JOscParameters<double&> parameters(multifunction_type::buffer);

      JOscParametersHelper_t::configure(parameters);
      
      this->load(fileName);
    }


    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     * \param  parameters           oscillation parameters
     */
    JOscProbInterpolator(const char*                            fileName,
                         const JOscParametersInterface<double>& parameters)
    {
      JOscParameters<double&> pars(multifunction_type::buffer);

      JOscParametersHelper_t::configure(pars);

      JOscParametersHelper_t::set(parameters);
      
      this->load(fileName);
    }


    /**
     * Load oscillation probability table from file.
     *
     * \param  fileName             oscillation probability table filename
     */
    void load(const char* fileName)
    {
      using namespace std;
      using namespace JPP;
      
      try {

        NOTICE("loading oscillation probability table from file " << fileName << "... " << flush);
        
        JLANG::load<JIO::JFileStreamReader>(fileName, static_cast<JOscProbInterpolator&>(*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
    {
      using namespace std;
      using namespace JPP;

      const JOscChannel* p = find(getOscChannel, getOscChannel + NUMBER_OF_OSCCHANNELS, channel);
        
      if (p != end(getOscChannel)) {
        
        const double L = baselineComputer.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);
      }
    }


    /**
     * Get cosine zenith angle for a given baseline.
     *
     * \param  L              baseline [km]
     * \return                cosine zenith angle
     */
    double getCosth(const double L) const override
    {
      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
    {
      return baselineComputer.getBaseline(costh);
    }
    
    
    /**
     * Read from input.
     *
     * \param  in                   reader
     * \return                      reader
     */
    JReader& read(JReader& in) override 
    {
      using namespace JPP;

      in >> baselineComputer;
      in >> static_cast<multifunction_type&>(*this);
      
      this->compile();
      
      return in;
    }
    
    
    /**
     * Write from input.
     *
     * \param  out                  writer
     * \return                      writer
     */
    JWriter& write(JWriter& out) const override 
    {
      out << baselineComputer;
      out << static_cast<const multifunction_type&>(*this);

      return out;
    }

    
  private:

    JBaselineComputer baselineComputer;
    
    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