JOscillogram.hh

Go to the documentation of this file.

#ifndef __JOSCPROB__JOSCILLOGRAM__
#define __JOSCPROB__JOSCILLOGRAM__

#include <iostream>
#include <iomanip>
#include <string>

#include "JLang/JException.hh"

#include "JTools/JGrid.hh"

#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscProbHelper.hh"
#include "JOscProb/JOscProbToolkit.hh"
#include "JOscProb/JOscProbInterface.hh"
#include "JOscProb/JOscProbInterpolator.hh"


/**
 * \author bjung, mdejong
 */

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

namespace JOSCPROB {

  using JLANG::JValueOutOfRange;
  using JLANG::JNullPointerException;

  using JTOOLS::JGrid;
  

  /**
   * Auxiliary class for defining an oscillogram axis.
   */
  struct JOscillogramAxis :
    public JGrid<double>
  {
    typedef JGrid<double>     JGrid_t;


    /**
     * Constructor.
     *
     * \param  type               oscillation variable type
     * \param  binning            oscillation variable binning
     */
    JOscillogramAxis(const int     type,
                     const JGrid_t binning) :
      JGrid_t(binning),
      type   (type)
    {}

    
    /**
     * Constructor.
     *
     * \param  type               oscillogram variable type
     * \param  N                  number of bins
     * \param  min                minimum axis value
     * \param  max                maximum axis value
     */
    JOscillogramAxis(const int    type,
                     const int    N,
                     const double min,
                     const double max) :
      JGrid_t(N, min, max),
      type   (type)
    {}

    
    int type; //!< axis type
  };
  
  
  /**
   * Auxiliary class for creating oscillograms.
   */
  struct JOscillogram :
    public JOscProbHelper
  {
    /**
     * Constructor.
     *
     * \param  abscissa           oscillogram abscissa axis definition
     * \param  ordinate           oscillogram ordinate axis definition
     * \param  channel            oscillation channel
     * \param  oscprob            oscillation probability calculator
     */
    JOscillogram(const JOscillogramAxis&  abscissa,
                 const JOscillogramAxis&  ordinate,
                 const JOscChannel&       channel,
                 const JOscProbInterface& oscprob) :
      JOscProbHelper(oscprob),
      abscissa      (abscissa),
      ordinate      (ordinate),
      channel       (channel)
    {}

    
    /**
     * Constructor.
     *
     * \param  abscissaName       oscillogram abscissa variable name
     * \param  abscissaBinning    oscillogram abscissa binning
     * \param  ordinateName       oscillogram ordinate variable name
     * \param  ordinateBinning    oscillogram ordinate binning
     * \param  channel            oscillation channel
     * \param  oscprob            oscillation probability calculator
     */
    JOscillogram(const std::string&       abscissaName,
                 const JGrid<double>      abscissaBinning,
                 const std::string&       ordinateName, 
                 const JGrid<double>      ordinateBinning,       
                 const JOscChannel&       channel,
                 const JOscProbInterface& oscprob) :
      JOscProbHelper(oscprob),
      abscissa      (JOscVars::getType(abscissaName), abscissaBinning),
      ordinate      (JOscVars::getType(ordinateName), ordinateBinning),
      channel       (channel)
    {}

    
    /**
     * Get energy corrresponding to the given bin indices.
     *
     * \param  i                  abscissa bin index
     * \param  j                  ordinate bin index
     * \return                    energy [GeV]
     */
    double getEnergy(const int i, const int j) const
    {
      const std::pair<double, double>& result = getTransformation(i, j);

      return result.first;
    }


    /**
     * Get cosine zenith angle corrresponding to the given bin indices.
     *
     * \param  i                  abscissa bin index
     * \param  j                  ordinate bin index
     * \return                    energy [GeV]
     */
    double getCosth(const int i, const int j) const
    {
      const std::pair<double, double>& result = getTransformation(i, j);

      return result.second;
    }


    /**
     * Get oscillation probability for given bin indices.
     *
     * \param  i                  abscissa bin index
     * \param  j                  ordinate bin index
     * \return                    oscillation probability
     */
    double getP(const int i, const int j) const
    {
      const double& energy = getEnergy(i, j);
      const double& costh  = getCosth (i, j);
      
      return JOscProbHelper::getP(channel, energy, costh);
    }

    
  private:
    

    /**
     * Get energy and cosine zenith angle corresponding to the given bin indices.
     *
     * \param  i                  abscissa bin index
     * \param  j                  ordinate bin index
     * \return                    transformed coordinate (E [GeV], costh)
     */
    std::pair<double, double> getTransformation(const int i, const int j) const
    {
      using namespace std;
      
      static int ix = -1;
      static int iy = -1;

      static pair<double, double> result = {0.0, 0.0}; // Cache result

      if (i != ix || j != iy) {

        ix = i;
        iy = j;
        
        const double x = abscissa.getX(i);
        const double y = ordinate.getX(j);

        switch (ordinate.type) {
        case (int) JOscVars::COSTH:
          result.second = y;
          break;
        case (int) JOscVars::SINTH:
          result.second = sqrt((1 + y) * (1 - y));
          break;
        case (int) JOscVars::BASELINE:
          result.second = JOscProbHelper::getCosth(y);
          break;
        default:
          THROW(JValueOutOfRange, "JOscillogram::getCosth(const double): Invalid ordinate type " << ordinate.type);
        }

        switch (abscissa.type) {
        case (int) JOscVars::ENERGY:
          result.first = x;
          break;
        case (int) JOscVars::LOG10E:
          result.first = pow(10.0, x);
          break;
        case (int) JOscVars::LOE:
          result.first = (x > 0.0 ? this->getBaseline(result.second) / x : 0.0);
          break;
        default:
          THROW(JValueOutOfRange, "JOscillogram::getEnergy(const double, const double): Invalid abscissa type " << abscissa.type);      
        }
      }

      return result;
    }
    
    
    JOscillogramAxis abscissa; //!< Abscissa axis
    JOscillogramAxis ordinate; //!< Ordinate axis

    JOscChannel      channel;  //!< Oscillation channel
  };
}

#endif