#ifndef __JCALIBRATE_JHVINTERPOLATOR__
#define __JCALIBRATE_JHVINTERPOLATOR__

#include "JLang/JException.hh"

#include "JTools/JRange.hh"

#include "JCalibrate/JFitToT.hh"

#include "TKey.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TAttMarker.h"
#include "TColor.h"


/**
 * \author bjung
 */

namespace JCALIBRATE {

  using JTOOLS::JRange;
  using JLANG::JNoValue;
  using JLANG::JValueOutOfRange;
  using JLANG::JNullPointerException;
  
  static const std::string HVINTERPOLATOR_SUFFIX = ".HVxG";
  
  static const char*       HVINTERPOLATOR_DATA   = "HVINTERPOLATOR_DATA";
  static const char*       HVINTERPOLATOR_RESULT = "HVINTERPOLATOR_RESULT";
  
  /**
   * Auxiliary data structure to store high-voltage versus gain data and interpolate the nominal high-voltage.
   */
  struct JHVInterpolator
  {
    /**
     * Constructor.
     *
     * \param  object             TMultiGraph object
     */
    JHVInterpolator(TMultiGraph& object) :
      data(&object)
    {
      TGraphErrors* g0 = new TGraphErrors();

      g0->SetName(HVINTERPOLATOR_RESULT);
      g0->SetLineColor  (kRed);
      g0->SetMarkerColor(kRed);
      g0->SetMarkerStyle(kFullDotSmall);
      g0->Set(0);

      data->Add(g0);
      
      TGraphErrors* g1 = new TGraphErrors();
      
      g1->SetName(HVINTERPOLATOR_DATA);
      g1->SetMarkerStyle(kFullDotSmall);
      g1->Set(0);
      
      data->Add(g1);
    }
    
    
    /**
     * Add point to diagram.
     *
     * \param  HV                 high-voltage        [V]
     * \param  gain               gain value
     * \param  gainError          error on gain value
     */
    void AddPoint(Double_t HV, Double_t gain, Double_t gainError)
    {      
      TGraphErrors* g1 = getData();
      
      const Int_t   N  = g1->GetN();
      
      g1->SetPoint     (N, fabs(HV), gain);
      g1->SetPointError(N, 0.0, gainError);
    }
    
    
    /**
     * Set point with index i.
     *
     * \param  i                  index of point
     * \param  HV                 high-voltage        [V]
     * \param  gain               gain value
     * \param  gainError          error on gain value
     */
    void SetPoint(Int_t i, Double_t HV, Double_t gain, Double_t gainError)
    {
      TGraphErrors* g1 = getData();

      g1->SetPoint     (i, fabs(HV), gain);
      g1->SetPointError(i, 0.0, gainError);
    }


    /**
     * Checks whether high-voltage is within range
     *
     * \param  HV                 high-voltage [V]
     * \return                    true if high-voltage is within range; else false
     */
    const bool checkHV(const double HV) const
    {
      return (HV > hvRange.getLowerLimit() &&
	      HV < hvRange.getUpperLimit());
    }
    
    
    /**
     * Checks whether two high-voltage values are different.
     *
     * \param  HV1                first  high-voltage [V]
     * \param  HV2                second high-voltage [V]
     * \return                    true if absolute difference is greater than minimal required difference; else false
     */
    const bool checkHV(const double HV1, const double HV2) const
    {
      return (fabs(HV1 - HV2) > dHVmin);
    }
    
    
    /**
     * Checks if gain is within range.
     *
     * \param   gain              gain value
     * \return                    true if gain within allowed range; else false
     */
    const bool checkGain(const double gain) const
    {
      return (gain > gainRange.getLowerLimit() &&
	      gain < gainRange.getUpperLimit());
    }
    
    
    /**
     * Checks whether the gains of two points are strictly increasing as function of their absolute high-voltage
     *
     * \param  i                  index of first  point
     * \param  j                  index of second point
     * \return                    true if gains of given points are strictly increasing\n 
     *                            as function of the absolute high-voltage; else false
     */
    const bool areIncreasing(const Int_t i, const Int_t j) const
    {
      TGraphErrors* g1 = getData();
      
      if ((i >= 0 && i < g1->GetN()) &&
	  (j >= 0 && i < g1->GetN())) {
	  
	return ((g1->GetPointX(i) > g1->GetPointX(j) && g1->GetPointY(i) > g1->GetPointY(j)) ||
		(g1->GetPointX(i) < g1->GetPointX(j) && g1->GetPointY(i) < g1->GetPointY(j)));
	
      } else {

	return false;
      }
    }
    
    
    /**
     * Checks whether two points are valid for inter-/extrapolation.
     *
     * \param  i                  index of first  point
     * \param  j                  index of second point
     * \return                    true if valid for inter-/extrapolation; else false
     */
    const bool areValid(const Int_t i, const Int_t j) const
    {
      TGraphErrors* g1 = getData();
      
      if ((i >= 0 && i < g1->GetN()) &&
	  (j >= 0 && j < g1->GetN())) {

	const double HV_i   = -fabs(g1->GetPointX(i));
	const double HV_j   = -fabs(g1->GetPointX(j));
	const double gain_i =       g1->GetPointY(i);
	const double gain_j =       g1->GetPointY(j);
	
	return (areIncreasing(   i,    j) && checkHV(HV_i) && checkGain(gain_i) &&
		checkHV      (HV_i, HV_j) && checkHV(HV_j) && checkGain(gain_j));
	
      } else {
	
	return false;
      }
    }
    
    
    /**
     * Inter-/Extrapolate the high-voltage value corresponding to the target gain value.
     *
     * \param  gainTarget         target gain value for inter-/extrapolation
     * \return                    true if inter-/extrapolation successful; else false
     */
    bool interpolateHV(const double gainTarget)
    {
      TGraphErrors* g1 = getData();
      
      if (g1->GetN() < 2 || !checkGain(gainTarget)) {
	return false;
      }
      
      // Search for valid inter-/extrapolation points
      
      Int_t i = 0;
      Int_t j = 1;
      
      for (Int_t k = 2; k < g1->GetN(); ++k) {

	const double dGain_i = fabs(g1->GetPointY(i) - gainTarget);	
	const double dGain_j = fabs(g1->GetPointY(j) - gainTarget);
	const double dGain_k = fabs(g1->GetPointY(k) - gainTarget);

	if         (dGain_k < dGain_j) {
	  i = (areValid(j,k) ? j : i);
	  j = k;
	} else if ((dGain_k < dGain_i || !areValid(i,j)) && areValid(j,k)) {
	  i = j;
	  j = k;
	}
      }

      // Inter-/Extrapolate high-voltage corresponding to given gain
      
      if (areValid(i, j)) {
	
	const double logHV0 = log(fabs(g1->GetPointX(i)));
	const double logHV1 = log(fabs(g1->GetPointX(j)));
	
	const double logG0  = log(g1->GetPointY(i));
	const double logG1  = log(g1->GetPointY(j));
	const double elogG0 = g1->GetErrorY(i) / g1->GetPointY(i);
	const double elogG1 = g1->GetErrorY(j) / g1->GetPointY(j);
	
	const double dlogG0 = log(gainTarget) - logG0;
	const double dlogG1 = log(gainTarget) - logG1;
	
	const double slope  = (logG1 - logG0) / (logHV1 - logHV0);
	
	const double HVnom  = exp(dlogG0 / slope + logHV0);
	const double eHVnom = HVnom * sqrt(dlogG1 * dlogG1 * elogG0 * elogG0 +
					   dlogG0 * dlogG0 * elogG1 * elogG1) / fabs(slope * (logG1 - logG0));
	
	const        double distance = fabs((log(HVnom) - logHV0) / (logHV0 - logHV1));
	static const double maxDist  = 2.0;
	
	if (checkHV(-HVnom) && distance < maxDist) {
	  
	  TGraphErrors* g0 = getResult();
	  
	  g0->SetPoint     (0,  HVnom, gainTarget);
	  g0->SetPointError(0, eHVnom,        0.0);
	  
	  return true;
	}
      }
      
      return false;
    }
    
    
    /**
     * Get interpolated high-voltage.
     *
     * \return                    inter-/extrapolated high-voltage\n
     *                            corresponding to target gain value [V]
     */
    double getHV() const
    {
      TGraphErrors* g0 = getResult();
      
      if (g0->GetN() > 0) {
	
	return g0->GetPointX(0);
	
      } else {

	THROW(JNoValue, "JHVInterpolator::getHV(): Missing HV inter-/extrapolation point. Please call JHVInterpolator::interpolateHV() first.");
      }
    }
    
    
    /**
     * Get error estimate on interpolated high-voltage.
     *
     * \return                    error on inter-/extrapolated high-voltage\n
     *                            corresponding to the target gain value [V]
     */
    double getHVError() const
    {
      TGraphErrors* g0 = getResult();
      
      if (g0->GetN() > 0) {

	return g0->GetErrorX(0);
	
      } else {

	THROW(JNoValue, "JHVInterpolator::getHVError(): Missing HV inter-/extrpolation point. Please call JHVInterpolator::interpolateHV() first.");	
      }
    }


    /**
     * Get graph with the input data for the interpolation.
     *
     * \return                   pointer to graph with input data
     */
    TGraphErrors* getData() const
    {
      TList* list = data->GetListOfGraphs();

      if (list != NULL && list->FindObject(HVINTERPOLATOR_DATA) != NULL) {
	
	return (TGraphErrors*)list->FindObject(HVINTERPOLATOR_DATA);
	
      } else if (list == NULL) {

	THROW(JNullPointerException, "JHVInterpolator()::getData(): Emtpy list of graphs.");

      } else {

	THROW(JNullPointerException, "JHVInterpolator()::getData(): No " << HVINTERPOLATOR_DATA << " graph.");
      }
    }


    /**
     * Get graph with the interpolation result.
     *
     * \return                   pointer to graph with the interpolation result
     */
    TGraphErrors* getResult() const
    {
      TList* list = data->GetListOfGraphs();

      if (list != NULL && list->FindObject(HVINTERPOLATOR_RESULT) != NULL) {
	
	return (TGraphErrors*)list->FindObject(HVINTERPOLATOR_RESULT);
	
      } else if (list == NULL) {

	THROW(JNullPointerException, "JHVInterpolator::getResult(): Empty list of graphs.");

      } else {
	
	THROW(JNullPointerException, "JHVInterpolator::getResult(): No " << HVINTERPOLATOR_RESULT << " graph.");
      }
    }

    
    /**
     * Set minimal separation distance for high-voltage.
     *
     * \param  minDist            minimal separation distance for high-voltage
     */
    static void setMinHVDistance(const double minDist)
    {
      dHVmin = minDist;
    }
    

    /**
     * Set valid gain range.
     *
     * \param  range              valid high-voltage range [V]
     */
    static void setHVRange(const JRange<double> range)
    {
      hvRange = range;
    }

    
    /**
     * Set valid gain range.
     *
     * \param  range              valid gain range
     */
    static void setGainRange(const JRange<double> range)
    {
      gainRange = range;
    }
    
    
  private:

    TMultiGraph*          data;           //!< HV-versus-gain data

    static double         dHVmin;         //!< Minimal high-voltage difference between two points [V]
    static JRange<double> hvRange;        //!< Allowed high-voltage range [V]
    static JRange<double> gainRange;      //!< Allowed gain range
  };


  /**
   * Default values.
   */
  double         JHVInterpolator::dHVmin    = 2 * 3.14;
  JRange<double> JHVInterpolator::hvRange   = JRange<double>(-1500, -80);
  JRange<double> JHVInterpolator::gainRange = JRange<double>(FITTOT_GAIN_MIN,
							     FITTOT_GAIN_MAX);
}

#endif