#ifndef __PREFIT__
#define __PREFIT__

#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <sstream>
#include <stdexcept>

#include "trigger/util.hh"
#include "nr/prob.hh"

#include "geometry.hh"
#include "tools.hh"


/**
 * \file
 * 1 dimensional prefit.
 */

namespace PREFIT {

  using namespace UTIL;
  using namespace GEOMETRY;
  using namespace TOOLS;
  using namespace NR;

  static const int NUMBER_OF_PARAMETERS = 3;    //!< Number of free parameters in pre-fit

  namespace $$$ {
    static const float RMIN = 1.0; //!< minimal distance for weighed centre-of-gravity [m]
    static const float TMIN = RMIN * C_INVERSE;
  }


  /**
   * exception handling
   */
  class not_enough_points : public std::exception {
  public: virtual const char *what() { return "not enough data points"; }
  };


  /**
   * 1 dimensional prefit in Z-direction
   *
   * The template class <hitIterator> refers to the definition of a hit iterator.
   * It should contain (at least) the following methods:
   *
   *         x()                 x position of hit [m]
   *         y()                 y position of hit [m]
   *         z()                 z position of hit [m]
   *         t()                 time of hit [ns]
   *         a()                 amplitude of hit [p.e.]
   *         sigma()             time resolution of hit [ns]
   */
  class prefit1Z : public track1Z { 
  public:
    /**
     * Constructor
     * \param begin_     begin of hit data
     * \param end_       end of hit data
     */
    template<class hitIterator>
    prefit1Z(const hitIterator& begin_,
	     const hitIterator& end_) :
      track1Z(),
      value(*this),
      error()
    {
      if (std::distance(begin_,end_) < NUMBER_OF_PARAMETERS)
	throw not_enough_points();

      hitIterator i;
      t_pos  dx, dy, dz, r, minz, maxz;
      t_time dt, sumt, mint, maxt;
      double weight, w, sum, beta_inverse, speed_inverse;
      double sumz, sumz2, sumzt;

      // centre-of-gravity

      sum = 0;
      x_  = 0;
      y_  = 0;
      z_  = 0;
      minz =  FLT_MAX;
      maxz = -FLT_MAX;
      mint =  DBL_MAX;
      maxt = -DBL_MAX;

      for (i = begin_; i != end_; ++i) {

	//w = 1.0 + 0.5*i->a()*i->a();
	w = 1.0;
		  
	sum += w;
	x_  += i->x() * w;
	y_  += i->y() * w;
	z_  += i->z() * w;
	dz   = i->z();
	dt   = i->t();
	if( dz < minz ) minz = dz;
	if( dz > maxz ) maxz = dz;
	if( dt < mint ) mint = dt;
	if( dt > maxt ) maxt = dt;
      }

      weight = 1.0/sum;

      x_ *= weight;
      y_ *= weight;
      z_ *= weight;

      speed_inverse = (maxt - mint) / (maxz - minz);
      //cout << "first estimate beta" << C_INVERSE / speed_inverse << endl;

      // first estimate of t
      
      sum = 0;
      t_  = 0;

      for (i = begin_; i != end_; ++i) {

	dx = i->x() - x();
	dy = i->y() - y();
	dz = i->z() - z();
	
	r  = sqrt(dx*dx + dy*dy);
	//if(r > 50)
	//continue;
	//w  = 1.0/(i->sigma()*i->sigma());
	w = 1.0;
	
	sum += w;
	t_  += ( i->t() - dz*speed_inverse ) * w;
      }

      t_ /= sum;
      sumt = sum;

      // weighed centre-of-gravity

      sum = 0;
      x_  = 0;
      y_  = 0;
      
      for (i = begin_; i != end_; ++i) {
	
	dz = i->z() - z();
	dt = i->t() - dz * speed_inverse - t();
	
	if (fabs(dt) < $$$::TMIN)  dt = $$$::TMIN;
	
	w = 1.0/(dt*dt);
	  
	sum += w;
	x_  += i->x() * w;
	y_  += i->y() * w;
      }
      
      x_ /= sum;
      y_ /= sum;

      // (x,y) error estimate **??**

      weight = C * getTanThetaC() / sqrt(sum);

      // error estimate

      error = point(weight, weight, 0, 1.0/sqrt(sumt));


      // determine sums

      dx = dy = dz = r = 0;
      dt = sumt = 0;
      sum = sumz = sumz2 = sumzt = 0;

      for (i = begin_; i != end_; ++i)
	{
	  dx = i->x() - x();
	  dy = i->y() - y();
	  dz = i->z() - z();
	  r  = sqrt( dx*dx + dy*dy );
	  //if(r > 50)
	  //continue;
	  dt = i->t();
	  
	  sumt  += dt;
	  sumz  += dz;
	  sumz2 += dz * dz;
	  sumzt += dz * dt;

	  sum++;
	}

      // determine parameters

      t_    = ( sumt * sumz2 - sumz * sumzt ) / ( sum * sumz2 - sumz * sumz );
      beta_ = C_INVERSE * ( sum * sumz2 - sumz * sumz ) / ( sum * sumzt - sumz * sumt );
      beta_inverse = 1/beta_;


      // evaluation of chi2

      dx = dy = dz = r = 0;
      dt = 0;
      chi2_ = 0;

      for(i = begin_; i != end_; ++i)
	{
	  //dx = i->x() - x();
	  //dy = i->y() - y();
	  dz = i->z() - z();
	  //r  = sqrt( dx*dx + dy*dy );
	  dt = i->t() - dz * beta_inverse * C_INVERSE - t();
	  w  = 1.0 / ( i->sigma() * i->sigma() );
	  
	  chi2_ += dt * dt * w;
	}
      
      chi2ndf_ = chi2_ / (sum - 2);

    }

    /**
     * Constructor
     * \param point_     new starting point
     */
    prefit1Z(const point& point_) :
      track1Z(point_),
      value(*this),
      error()
    {}

    /**
     * assignment operator
     * \param point_     new starting point
     */
    prefit1Z& operator=(const point& point_) 
    {
      value = point_;
      error = point();
      return *this;
    }

    /**
     * assignment operator
     * \param prefit_    new fit values
     */
    prefit1Z& operator=(const prefit1Z& prefit_) 
    {
      value = prefit_.value;
      error = prefit_.error;
      return *this;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    /*
    std::ostream& write(std::ostream& out) const
    {
      track1Z::write(out);
      out << ' ' << beta_ << ' ' << chi2_ << ' ' << chi2ndf_;
      return out;
    }
    */

  public:
    point& value; //!< reference to actual values
    point  error; //!< error esitmate
    float  beta_;  //!< velocity as fraction of c
    float  chi2_;  //!< chi square of prefit
    float  chi2ndf_;
  };


  /**
   * 1 dimensional prefit in Z-direction for muons
   *
   * The template class <hitIterator> refers to the definition of a hit iterator.
   * It should contain (at least) the following methods:
   *
   *         x()                 x position of hit [m]
   *         y()                 y position of hit [m]
   *         z()                 z position of hit [m]
   *         t()                 time of hit [ns]
   *         a()                 amplitude of hit [p.e.]
   *         sigma()             time resolution of hit [ns]
   */
  class prefit1Zmu : public track1Z { 
  public:
    /**
     * Constructor
     * \param begin_     begin of hit data
     * \param end_       end of hit data
     */
    template<class hitIterator>
    prefit1Zmu(const hitIterator& begin_,
	       const hitIterator& end_) :
      track1Z(),
      value(*this),
      error()
    {
      if (std::distance(begin_,end_) < NUMBER_OF_PARAMETERS)
	throw not_enough_points();
      
      hitIterator i;
      double weight, w, sum;
      t_pos  dx, dy, dz, r;
      t_time tt;

      // centre-of-gravity

      sum = 0;
      x_  = 0;
      y_  = 0;
      z_  = 0;

      for (i = begin_; i != end_; ++i) {

	//w = 1.0 + 0.5*i->a()*i->a();
	w = 1.0;
		  
	sum += w;
	x_  += i->x() * w;
	y_  += i->y() * w;
	z_  += i->z() * w;
      }

      weight = 1.0/sum;

      x_ *= weight;
      y_ *= weight;
      z_ *= weight;


      // first estimate of t
      
      sum = 0;
      t_  = 0;

      for (i = begin_; i != end_; ++i) {

	dx = i->x() - x();
	dy = i->y() - y();
	dz = i->z() - z();
	
	r  = sqrt(dx*dx + dy*dy);
	w  = 1.0/(i->sigma()*i->sigma());
	//w = 1.0;
	
	sum += w;
	t_  += (i->t() - (dz + r*getTanThetaC())*C_INVERSE) * w;
      }

      t_ /= sum;


      // weighed centre-of-gravity

      sum = 0;
      x_  = 0;
      y_  = 0;
      
      for (i = begin_; i != end_; ++i) {
	
	dz = i->z() - z();
	tt = i->t() - dz * C_INVERSE - t();
	
	if (fabs(tt) < $$$::TMIN)  tt = $$$::TMIN;
	
	w = 1.0/(tt*tt);
	  
	sum += w;
	x_  += i->x() * w;
	y_  += i->y() * w;
      }
      
      x_ /= sum;
      y_ /= sum;

      // (x,y) error estimate

      weight = C * getTanThetaC() / sqrt(sum);


      // final estimate of t

      sum = 0;
      t_  = 0;

      for (i = begin_; i != end_; ++i)
	{
	  dx = i->x() - x();
	  dy = i->y() - y();
	  dz = i->z() - z();

	  r  = sqrt( dx*dx + dy*dy );
	  w  = 1.0/(i->sigma()*i->sigma());
	  
	  sum += w;
	  t_  += (i->t() - (dz + r*getTanThetaC())*C_INVERSE) * w;
	}

      t_ /= sum;

      // error estimate

      error = point(weight, weight, 0, 1.0/sqrt(sum));


      // evaluation of chi2

      sum   = 0;
      chi2_ = 0;

      for(i = begin_; i != end_; ++i)
	{
	  dx = i->x() - x();
	  dy = i->y() - y();
	  dz = i->z() - z();
	  r  = sqrt( dx*dx + dy*dy );
	  
	  tt = i->t() - (dz + r*getTanThetaC()) * C_INVERSE - t();
	  w  = 1.0 / ( i->sigma() * i->sigma() );
	  
	  chi2_ += tt * tt * w;
	  sum++;
	}
      
      chi2ndf_ = chi2_ / (sum - 1);

    }

    /**
     * Constructor
     * \param point_     new starting point
     */
    prefit1Zmu(const point& point_) :
      track1Z(point_),
      value(*this),
      error()
    {}

    /**
     * assignment operator
     * \param point_     new starting point
     */
    prefit1Zmu& operator=(const point& point_) 
    {
      value = point_;
      error = point();
      return *this;
    }

    /**
     * assignment operator
     * \param prefit_    new fit values
     */
    prefit1Zmu& operator=(const prefit1Zmu& prefit_) 
    {
      value = prefit_.value;
      error = prefit_.error;
      return *this;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    /*
    std::ostream& write(std::ostream& out) const
    {
      track1Z::write(out);
      out << ' ' << chi2_ << ' ' << chi2ndf_;
      return out;
    }
    */

  public:
    point& value; //!< reference to actual values
    point  error; //!< error esitmate
    float  chi2_;  //!< chi square of prefit
    float  chi2ndf_;
  };

}

#endif