#ifndef __JPOINT4DREGRESSOR__
#define __JPOINT4DREGRESSOR__ 

#include <cmath>

#include "JPhysics/JConstants.hh"

#include "JGeometry3D/JVector3D.hh"

#include "JLang/JSharedPointer.hh"
#include "JFit/JTimeRange.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JPoint4D.hh"

#include "Jeep/JMessage.hh"

/**
 * \author adomi
 */

namespace JFIT {

  /**
   * Regressor function object for JPoint4D fit using JSimplex minimiser.
   */
  template<>
  struct JRegressor<JPoint4D, JSimplex> :
    public JAbstractRegressor<JPoint4D, JSimplex>
  {
    using JAbstractRegressor<JPoint4D, JSimplex>::operator();

    /**
     * Constructor.
     */
    JRegressor(double sigma)
    {
      this->sigma = sigma;
    }

    /*  Fit Function
     *  This method is used to determine the chi2 of given hit with respect
     *  to Shower's vertex (brightest point) 
     * 
     *  \param vx           Shower's vertex      
     *  \param hit          hit
     *  \return             chi2
     */
    template<class JHit_t>
    double operator()(const JPoint4D& vx, const JHit_t& hit) const
    {
      using namespace JPP;
      
      const double dt = hit.getT() - vx.getT(hit.getPosition());
      
      const double u  = dt / sigma;
      
      return estimator->getRho(u) * hit.getW();
    }
    
    JLANG::JSharedPointer<JMEstimator> estimator;  //!< M-Estimator function    
    double                             sigma;      //!< Time resolution [ns]
  };
  template<>
  struct JRegressor<JPoint4D, JGandalf> :
    public JAbstractRegressor<JPoint4D, JGandalf>
  {
    using JAbstractRegressor<JPoint4D, JGandalf>::operator();

    /**
     * Constructor.
     */
    JRegressor(double sigma)
    {
      this->sigma = sigma;
    }

    /*  Fit Function
     *  This method is used to determine the chi2 of given hit with respect
     *  to Shower's vertex (brightest point)
     *
     *  \param vx           Shower's vertex
     *  \param hit          hit
     *  \return             chi2
     */
    template<class JHit_t>
    result_type operator()(const JPoint4D& vx, const JHit_t& hit) const
    {
      using namespace JPP;

      const double dt = hit.getT() - vx.getT(hit.getPosition());

      const double u  = dt / sigma;
      result_type result;
      result.chi2 = estimator->getRho(u) * hit.getW();

      JVector3D d = hit.getPosition() - vx.getPosition();
      
      double weight = getIndexOfRefraction() * getInverseSpeedOfLight()/(d.getLength() * sigma);
      
      result.gradient = hit.getW() * JPoint4D( d * weight,
					       -1/sigma);
      result.gradient.mul(0.5 * estimator->getPsi(u));
      
      return result;
    }

    JLANG::JSharedPointer<JMEstimator> estimator;  //!< M-Estimator function
    double                             sigma;      //!< Time resolution [ns]
  };
}

#endif