#ifndef __JACOUSTICS__JKATOOMBA_T__
#define __JACOUSTICS__JKATOOMBA_T__

#include "JFit/JEstimator.hh"

#include "JAcoustics/JKatoomba.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JGeometry.hh"
#include "JAcoustics/JModel.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JTimeRange.hh"


/**
 * \file
 *
 * Fit function of acoustic model.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }

namespace JACOUSTICS {

  /**
   * Worker class for complete fit procedure of acoustic model.
   */
  struct JKatoomba_t {
    /**
     * Constructor.
     *
     * \param  geometry           detector geometry
     * \param  V                  sound velocity
     * \param  parameters         parameters
     */
    JKatoomba_t(const JGeometry& geometry, const JSoundVelocity& V, const JFitParameters& parameters) :
      parameters(parameters),
      estimator(geometry, V, parameters.option),
      evaluator(geometry, V, parameters.option),
      gandalf  (geometry, V, parameters.option)
    {
      using namespace JPP;

      evaluator.estimator.reset(getMEstimator(parameters.mestimator));
      gandalf  .estimator.reset(getMEstimator(parameters.mestimator));
    }


    /**
     * Auxiliary data structure for return value of fit.
     */
    template<class T>
    struct result_type {
      JModel& value;
      double  chi2;
      double  ndf;
      T       begin;
      T       end;

      /**
       * Get number of hits used in fit (after outlier removal).
       *
       * \return                size
       */
      int size() const
      {
	return std::distance(begin, end);
      }
      
      /**
       * Get time range of data.
       *
       * \return                time range
       */
      JTimeRange getTimeRange() const
      {
	JTimeRange range(JTimeRange::DEFAULT_RANGE());

	for (T i = begin; i != end; ++i) {
	  range.include(i->getValue());
	}

	return range;
      }
    };


    /**
     * Fit.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     * \return                fit result
     */
    template<class T>
    result_type<T> operator()(T __begin, T __end)
    {
      using namespace std;

      JModel  model;
      JModel& result = gandalf.value;                                 // start value

      result = estimator(__begin, __end);                             // pre-fit

      if (parameters.stdev > 0.0) {                                   // remove outliers

	double chi2[] = { numeric_limits<double>::max(),
			  numeric_limits<double>::max() };

	for ( ; ; ) {

	  T out = __end;

	  double xmax = 0.0;

	  for (T hit = __begin; hit != __end; ++hit) {

	    const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->getSigma();

	    if (x > xmax) {
	      xmax = x;
	      out  = hit;
	    }
	  }

	  if (xmax > parameters.stdev) {

	    if (chi2[0] == numeric_limits<double>::max()) {
	      chi2[0] = evaluator(result, __begin, __end);
	    }

	    model = result;

	    iter_swap(out, --__end);

	    result  = estimator(__begin, __end);
	    chi2[1] = evaluator(result, __begin, __end);

	    if (chi2[0] - chi2[1]  >  parameters.stdev * parameters.stdev) {

	      chi2[0] = chi2[1];

	      continue;

	    } else {

	      iter_swap(out, __end++);

	      model = result;
	    }
	  }

	  break;
	}
      }

      const double chi2  =  gandalf  (__begin, __end)  /  gandalf.estimator->getRho(1.0);
      const double ndf   =  getWeight(__begin, __end)  -  gandalf.value.getN();

      return { gandalf.value, chi2, ndf, __begin, __end };            // return values
    }


    JFitParameters                parameters;
    JKatoomba<JEstimator>         estimator;
    JKatoomba<JAbstractMinimiser> evaluator;
    JKatoomba<JGandalf>           gandalf;
  };
}

#endif