#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/JOscProbToolkit.hh"
#include "JOscProb/JOscProbInterpolatorInterface.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
  {
    /**
     * Constructor.
     *
     * \param  abscissa           oscillogram abscissa axis definition
     * \param  ordinate           oscillogram ordinate axis definition
     * \param  channel            oscillation channel
     * \param  interpolator       pointer to oscillation probability interpolator
     */
    JOscillogram(const JOscillogramAxis&              abscissa,
		 const JOscillogramAxis&              ordinate,
		 const JOscChannel&                   channel,
		 const JOscProbInterpolatorInterface* pInterpolator) :
      abscissa     (abscissa),
      ordinate     (ordinate),
      channel      (channel),
      pInterpolator(pInterpolator)
    {
      if (pInterpolator == NULL) {
	THROW(JNullPointerException, "JOscillogram::JOscillogram(...): Oscillation probability interpolator is not set");
      }
    }

    
    /**
     * 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  interpolator       pointer to oscillation probability interpolator
     */
    JOscillogram(const std::string&                   abscissaName,
		 const JGrid<double>                  abscissaBinning,
		 const std::string&                   ordinateName,	
		 const JGrid<double>                  ordinateBinning,	 
		 const JOscChannel&                   channel,
		 const JOscProbInterpolatorInterface* pInterpolator) :
      abscissa     (JOscVars::getType(abscissaName), abscissaBinning),
      ordinate     (JOscVars::getType(ordinateName), ordinateBinning),
      channel      (channel),
      pInterpolator(pInterpolator)
    {
      if (pInterpolator == NULL) {
	THROW(JNullPointerException, "JOscillogram::JOscillogram(...): Oscillation probability interpolator is not set");
      }
    }

    
    /**
     * 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 (*pInterpolator)(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 const JBaselineCalculator& baselineCalculator = pInterpolator->getBaselineCalculator();
      
      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 = baselineCalculator.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 ? baselineCalculator.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

    const JOscProbInterpolatorInterface* pInterpolator; //!< Pointer to oscillation probability interpolator
  };
}

#endif