#ifndef __JRECONSTRUCTION__JEVT__
#define __JRECONSTRUCTION__JEVT__

#include <istream>
#include <ostream>
#include <iomanip>
#include <limits>
#include <vector>
#include <algorithm>
#include <cmath>

#include <TROOT.h>
#include <TObject.h>

#include "JLang/JType.hh"
#include "JROOT/JRoot.hh"
#include "JROOT/JTreeParameters.hh"

#include "JReconstruction/JFitStatus.hh"
#include "JReconstruction/JHistory.hh"

/**
 * \author mdejong
 */

namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }

namespace JFIT {


  /**
   * Data structure for track fit results with history and optional associated values.
   */
  class JFit : 
    public TObject,
    public JHistory
  {
  public:
    /**
     * Default constructor.
     *
     * Parameters are initialized with non physical values.
     */
    JFit():
      __x(std::numeric_limits<double>::lowest()),
      __y(std::numeric_limits<double>::lowest()),
      __z(std::numeric_limits<double>::lowest()),
      __dx(0.0),
      __dy(0.0),
      __dz(0.0),
      __t(std::numeric_limits<double>::lowest()),
      __Q(std::numeric_limits<double>::lowest()),
      __NDF(-1),
      __E(0.0),
      __status(JRECONSTRUCTION::UNDEFINED)
    {}


    /**
     * Constructor.
     * 
     * \param  history  history
     * \param  x        X-position
     * \param  y        Y-position
     * \param  z        Z-position
     * \param  dx       X-slope
     * \param  dy       Y-slope
     * \param  dz       Z-slope
     * \param  t        time
     * \param  Q        quality
     * \param  NDF      number of degrees of freedom
     * \param  E        energy
     * \param  status   status
     */
    JFit(const JHistory& history,
	 const double x,
	 const double y,
	 const double z,
	 const double dx,
	 const double dy,
	 const double dz,
	 const double t,
	 const double Q,
	 const int    NDF,
	 const double E      = 0.0,
	 const int    status = JRECONSTRUCTION::SINGLE_STAGE) :
      JHistory(history)
    {
      __x      = x;
      __y      = y;
      __z      = z;
      __dx     = dx;
      __dy     = dy;
      __dz     = dz;
      __t      = t;
      __Q      = Q;
      __NDF    = NDF;
      __E      = E;
      __status = status;
    }


    /**
     * Constructor for storing position only.
     *
     * Note that the type of fit can be obtained via the history information.
     *
     * \param  history  history
     * \param  x        X-position
     * \param  y        Y-position
     * \param  z        Z-position 
     * \param  status   status
     */
    JFit(const JHistory& history,
	 const double x,
	 const double y,
	 const double z,
	 const int    status = JRECONSTRUCTION::SINGLE_STAGE):
      JHistory(history)
    {
      __x      = x;
      __y      = y;
      __z      = z;
      __dx     = 0.0;
      __dy     = 0.0;
      __dz     = 0.0;
      __t      = 0.0;
      __Q      = 0.0;
      __NDF    = -1;
      __E      = 0.0;
     __status  = status;
    }
    

    /**
     * Add event to history.
     *
     * \param  type          application type
     * \return               this fit
     */
    JFit& add(const int type)
    {
      getHistory().add(type);

      return *this;
    }


    double getX()       const { return __x; }           //!< Get X-position.
    double getY()       const { return __y; }           //!< Get Y-position.
    double getZ()       const { return __z; }           //!< Get Z-position.
    double getDX()      const { return __dx; }          //!< Get X-slope.
    double getDY()      const { return __dy; }          //!< Get Y-slope.
    double getDZ()      const { return __dz; }          //!< Get Z-slope.
    double getT()       const { return __t; }           //!< Get time.
    double getQ()       const { return __Q; }           //!< Get quality.
    int    getNDF()     const { return __NDF; }         //!< Get number of degrees of freedom.
    double getE()       const { return __E; }           //!< Get energy.
    int    getStatus()  const { return __status; }      //!< Get status of the fit; negative values should refer to a bad fit.


    /**
     * Set status of the fit.
     *
     * \param  status        status
     */
    void setStatus(const int status)
    {
      this->__status = status;
    }


    /**
     * Move vertex along this track with given velocity.
     *
     * \param  step          step
     * \param  velocity      velocity
     */
    void move(const double step, const double velocity)
    {
      __x += step * __dx;
      __y += step * __dy;
      __z += step * __dz;
      __t += step / velocity;
    }


    /**
     * Set energy.
     *
     * \param  E        energy
     */
    void setE(const double E)
    {
      __E = E;
    }


    /**
     * Get dimension of error matrix.
     *
     * \return          dimension
     */
    size_t getDimensionOfErrorMatrix() const
    {
      return ((size_t) sqrt(1.0 + 8.0*V.size()) - 1) / 2;
    }


    /**
     * Get error matrix.
     *
     * Note that only the lower-half of the matrix is returned.
     *
     * \return          error matrix
     */
    const std::vector<double>& getV() const
    {
      return V;
    }


    /**
     * Get element of error matrix.
     *
     * \param  row      row (<tt>0 <= row < dimension</tt>)
     * \param  col      col (<tt>0 <= col < dimension</tt>)
     * \return          value
     */
    double getV(const size_t row, const size_t col) const
    {
      using namespace std;

      const size_t __row = max(row, col) + 1;
      const size_t __col = min(row, col);

      return V.at((__row * (__row + 1)) / 2 - __row + __col);
    }


    /**
     * Set error matrix.
     *
     * The given size corresponds to the dimension of a 2D-array of which
     * the elements should be accessible via the usual array operators.\n
     * Note that only the lower-half of the given matrix is stored.
     *
     * \param  size     size
     * \param  data     matrix
     */
    template<class T>
    void setV(const size_t size, const T& data)
    {
      V.clear();

      for (size_t row = 0; row != size; ++row) {
        for (size_t col = 0; col <= row; ++col) {
          V.push_back(data[row][col]);
        }
      }
    }
 

    /**
     * Get associated values.
     *
     * \return          values
     */
    const std::vector<double>& getW() const
    {
      return this->W;
    }


    /**
     * Set associated values.
     *
     * \param  W        values
     */
    void setW(const std::vector<double>& W)
    {
      this->W = W;
    }


    /**
     * Get number of associated values.
     *
     * \return          number of values
     */
    int getN() const 
    { 
      return W.size();
    }
  

    /**
     * Check availability of value.
     *
     * \param  i        index
     * \return          true if available; else false
     */
    bool hasW(const int i) const
    { 
      return (i >= 0 && i < (int) W.size());
    }
  

    /**
     * Get associated value.
     *
     * \param  i        index
     * \return          value
     */
    double getW(const int i) const
    { 
      return W.at(i);
    }
  

    /**
     * Get associated value.
     *
     * \param  i        index
     * \param  value    default value
     * \return          value
     */
    double getW(const int i, const double value) const
    { 
      if (hasW(i))
	return W.at(i);
      else
	return value;
    }


    /**
     * Set associated value.
     *
     * \param  i        index
     * \param  value    value
     */
    void setW(const int i, const double value)
    { 
      if (i >= (int) W.size()) {
	W.resize(i + 1, 0.0);
      }

      W[i] = value;
    }


    ClassDef(JFit, 7);

  protected:
    double  __x;
    double  __y;
    double  __z;
    double  __dx;
    double  __dy;
    double  __dz;
    double  __t;
    double  __Q;
    int     __NDF;
    std::vector<double> V;
    std::vector<double> W;
    double  __E;
    int __status;
  };


  /**
   * Data structure for set of track fit results.
   */
  class JEvt : 
    public TObject,
    public std::vector<JFit>
  {
  public:
    /**
     * Default constructor.
     */
    JEvt()
    {}


    /**
     * Select fits.
     *
     * \param  selector         fit selector
     */
    template<class JSelector_t>
    void select(const JSelector_t& selector)
    {
      using namespace std;

      if (!empty()) {

	iterator __end = partition(this->begin(), this->end(), selector);

	this->erase(__end, this->end());
      }
    }


    /**
     * Select fits.
     *
     * \param  number_of_fits   maximal number of best fits to select
     * \param  comparator       quality comparator
     */
    template<class JComparator_t>
    void select(const size_t         number_of_fits,
		const JComparator_t& comparator)
    {
      using namespace std;

      if (!empty()) {

	iterator __end = this->end();

	if (number_of_fits > 0 && number_of_fits < this->size()) { 

	  advance(__end = this->begin(), number_of_fits);

	  partial_sort(this->begin(), __end, this->end(), comparator);

	} else {

	  sort(this->begin(), __end, comparator);
	}

	this->erase(__end, this->end());
      }
    }
      
    
    /**
     * Write event to output.
     *
     * \param  out      output stream
     * \param  event    event
     * \return          output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JEvt& event)
    {
      using namespace std;
    
      out << "Event: " << endl;

      for (JEvt::const_iterator fit = event.begin(); fit != event.end(); ++fit) {
	out << *fit;
      }
    
      return out;
    }


    ClassDef(JEvt, 5);
  };
}


namespace JRECONSTRUCTION {
  typedef JFIT::JFit   JFit;
  typedef JFIT::JEvt   JEvt;
}


/**
 * Write fit results to output.
 *
 * \param  out      output stream
 * \param  fit      fit results
 * \return          output stream
 */
std::ostream& operator<<(std::ostream& out, const JRECONSTRUCTION::JFit& fit);


/**
 * Get TTree parameters for given data type.
 *
 * \return          TTree parameters
 */
inline JROOT::JTreeParameters getTreeParameters(JLANG::JType<JRECONSTRUCTION::JEvt>)
{
  return JROOT::JTreeParameters("EVT", "evt", "", 0);
}

#endif