JOscProbInterpolator2D.hh

Go to the documentation of this file.

#ifndef __JOSCPROB__JOSCPROBINTERPOLATOR2D__
#define __JOSCPROB__JOSCPROBINTERPOLATOR2D__

#include <string>

#include "km3net-dataformat/offline/Evt.hh"

#include "Jeep/JMessage.hh"

#include "JLang/JException.hh"

#include "JTools/JDistance.hh"

#include "JOscProb/JOscProbTable_t.hh"
#include "JOscProb/JOscProbTable.hh"
#include "JOscProb/JOscProbToolkit.hh"
#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscParameters.hh"


namespace JOSCPROB {};
namespace JPP { using namespace JOSCPROB; }

namespace JOSCPROB {

  using JLANG::JValueOutOfRange;  
  using JLANG::JNullPointerException;
  
  
  /**
   * Auxiliary class for interpolating an oscillation probability table in 8 dimensions in terms of:
   *
   * 1. \f$ E                     \left[\mathrm{GeV}\right]  \f$
   * 2. \f$ \sin(\theta)\f$
   */
  template<class JDistance_t = JTOOLS::JDistance<typename JOscProbTable2D_t::JFunction2D_t::argument_type> >
  struct JOscProbInterpolator2D
  {
    typedef JOscProbTable2D<JDistance_t>             table_type;
    typedef typename table_type::data_type            data_type;
    
    
    /**
     * Default constructor.
     */
    JOscProbInterpolator2D() :
      pTable()
    {}


    /**
     * Constructor.
     *
     * \param  pTable             pointer to 2D oscillation probability table
     */
    JOscProbInterpolator2D(const table_type* pTable) :
      pTable(pTable)
    {}    

    
    /**
     * Get oscillation probability for given oscillation channel.
     *
     * \param  channel            oscillation channel
     * \param  E                  neutrino energy     [GeV]
     * \param  costh              cosine zenith angle
     * \return                    oscillation probability
     */
    double operator()(const JOscChannel& channel,
                      const double       E,
                      const double       costh) const
    {
      using namespace std;
      using namespace JPP;

      if (pTable != NULL) {
      
        const JOscChannel* p = find(getOscChannel, getOscChannel + NUMBER_OF_OSCCHANNELS, channel);
        
        if (p != end(getOscChannel)) {

          const double L     = getL(costh);
          const double sinth = (fabs(costh) < 1.0 ? sqrt((1 - costh) * (1 + costh)) : 0.0);

          const size_t      index         = distance(getOscChannel, p);
          const data_type&  probabilities = (*pTable)(L/E, sinth);

          return probabilities[index];
        
        } else {
        
          THROW(JValueOutOfRange, "JOscProbInterpolator2D::operator(): Invalid oscillation channel " << channel << endl);
        }
        
      } else {
        
        THROW(JNullPointerException, "JOscProbInterpolator2D::operator(): Unspecified oscillation probability table." << endl);
      }
    }

    
  private:
    
    const table_type* pTable; //!< Pointer to 2D oscillation probability table
  };
}


#endif