// HitVector.h

#ifndef HitVector_H
#define HitVector_H

// Class HitVector contains the parameters describing a hit.
// Class HitError contains the corresponding error matrix.
// Class HitDerivative is the derivative with respect to a track
// vector, i.e. is an nx5 matrix where n = hit dimension.
//
// Any object constructed without specified values should be intitialized
// with zeroes.

#include "TrfVector.h"
#include <cassert>
#include <iosfwd>

//**********************************************************************

namespace trf {

class HitVector;
class HitError;
class HitDerivative;
class TrackError;

//**********************************************************************

// Friends.

//**********************************************************************

// output stream
std::ostream& operator<<(std::ostream& stream, const HitVector& rhs);

// vector + vector
HitVector operator+(const HitVector& lhs, const HitVector& rhs);

// vector - vector
HitVector operator-(const HitVector& lhs, const HitVector& rhs);

// equality
// must have the same length and the same values
bool operator==(const HitVector& lhs, const HitVector& rhs);

// inequality
bool operator!=(const HitVector& lhs, const HitVector& rhs);

// equality with tolerance
// must have the same length and the same values
bool is_equal(const HitVector& lhs, const HitVector& rhs);

// error + error
HitError operator+(const HitError& lhs, const HitError& rhs);

// error - error
HitError operator-(const HitError& lhs, const HitError& rhs);

// equality
// must have the same dimension and the same values
bool operator==(const HitError& lhs, const HitError& rhs);

// inequality
bool operator!=(const HitError& lhs, const HitError& rhs);

// equality with tolerance
// must have the same dimension and the same values
bool is_equal(const HitError& lhs, const HitError& rhs);

// output stream
std::ostream& operator<<(std::ostream& stream, const HitError& rhs);

//**********************************************************************

class HitVector {

private:

  trf::TrfVector _vec;

  // construct from implementation
  explicit HitVector(const trf::TrfVector& hvec);

public:

  // default constructor
  explicit HitVector(int size =0);

  // constructor from an array
  explicit HitVector(int size, double* arr);

  // constructor from one value
  explicit HitVector(double x1);

  // constructor from two values
  HitVector(double x1, double x2);

  // constructor from three values
  HitVector(double x1, double x2, double x3);

  // copy constructor
  HitVector (const HitVector& hvec);

  // assignment
  // left side must be unassigned or have the same length
  HitVector& operator=( const HitVector& hvec ) {
    assert( size() == 0 || size() == hvec.size() );
    _vec = hvec._vec;
    return *this;
  }

  // Return the underlying vector.
  const trf::TrfVector& get_vector() const { return _vec; };

  // return the dimension of the hit
  int size() const { return _vec.length(); };

#ifndef TRF_USING_LinearAlgebra
  // accessor
  double& operator()(int i) {
    assert( i >= 0 );
    assert( i < size() );
    return _vec(i);
  }
#endif

  // const accessor
  double operator()(int i) const {
    assert( i >= 0 );
    assert( i < size() );
    return _vec(i);
  }

  // minimum
  double min() const { return _vec.min(); };

  // maximum
  double max() const { return _vec.max(); };

  // absolute minimum
  double amin() const { return _vec.amin(); };

  // absolute maximum
  double amax() const { return _vec.amax(); };

  // +=
  HitVector& operator+=(const HitVector& rhs) {
    _vec += rhs._vec;
    return *this;
  }

  // -=
  HitVector& operator-=(const HitVector& rhs) {
    _vec -= rhs._vec;
    return *this;
  }

  // vector + vector
  friend HitVector
  operator+(const HitVector& lhs, const HitVector& rhs);

  // vector - vector
  friend HitVector 
  operator-(const HitVector& lhs, const HitVector& rhs);

  // equality
  // must have the same length and the same values
  friend bool operator==(const HitVector& lhs, const HitVector& rhs);

  // inequality
  friend bool operator!=(const HitVector& lhs, const HitVector& rhs) {
    return ! ( lhs == rhs );
  }

  // equality with tolerance
  // must have the same length and the same values
  friend bool is_equal(const HitVector& lhs, const HitVector& rhs);

  // output stream
  friend std::ostream&
  operator<<(std::ostream& stream, const HitVector& rhs);

};

//**********************************************************************

class HitError {

private:
  trf::TrfSMatrix _err;

  // constructor from implementation
  explicit HitError(const trf::TrfSMatrix& err);

public:

  // constructor
  explicit HitError(int size =0);

  // constructor from an array
  // order is lower triangle (00, 10, 11, 20, ...)
  HitError(int size, double* arr);

  // constructor for 1D from values
  explicit HitError(double e11);

  // constructor for 2D from values
  HitError(double e11, double e12, double e22);

  // constructor for 3D from values
  HitError(double e11, double e21, double e22, 
           double e31, double e32, double e33);

  // constructor from a track error and a track derivative:
  // E_hit = dhit_dtrack * E_track * dhit_dtrack_transpose
  HitError(const HitDerivative& dhit_dtrack, 
           const TrackError& trkerr);

  // copy constructor
  HitError(const HitError& herr);

  // assignment
  // left side must be unassigned or have the same length
  HitError& operator=(const HitError& herr) {
    assert( size() == 0 || size() == herr.size() );
    _err = herr._err;
    return *this;
  }

  // Return the underlying matrix.
  const trf::TrfSMatrix& get_matrix() const { return _err; };

  // return the dimension of the matrix
  int size() const { return _err.nrow(); };

#ifndef TRF_USING_LinearAlgebra
  // accessor
  double& operator()( int i, int j ) {
    assert( i >= 0 );
    assert( i < size() );
    assert( j >= 0 );
    assert( j < size() );
    return _err(i,j);
  }
#endif

  // const accessor
  double operator()(int i, int j) const {
    assert( i >= 0 );
    assert( i < size() );
    assert( j >= 0 );
    assert( j < size() );
    return _err(i,j);
  }

  // minimum
  double min() const { return _err.min(); };

  // maximum
  double max() const { return _err.max(); };

  // absolute minimum
  double amin() const { return _err.amin(); };

  // absolute maximum
  double amax() const { return _err.amax(); };

  // invert -- return 0 for success
  int invert();

  // +=
  HitError& operator+=(const HitError& rhs) {
    _err += rhs._err;
    return *this;
  }

  // -=
  HitError& operator-=(const HitError& rhs) {
    _err -= rhs._err;
    return *this;
  }

  // error + error
  friend HitError
  operator+(const HitError& lhs, const HitError& rhs) {
    return HitError(lhs._err + rhs._err );
  }

  // error - error
  friend HitError
  operator-(const HitError& lhs, const HitError& rhs) {
    return HitError( lhs._err - rhs._err );
  }

  // equality
  // must have the same dimension and the same values
  friend bool operator==(const HitError& lhs, const HitError& rhs) {
    if ( lhs.size() != rhs.size() ) return false;
    return lhs._err == rhs._err;
  }

  // inequality
  friend bool operator!=(const HitError& lhs, const HitError& rhs) {
    return ! ( lhs == rhs );
  }

  // equality with tolerance
  // must have the same dimension and the same values
  friend bool is_equal(const HitError& lhs, const HitError& rhs) {
    if ( lhs.size() != rhs.size() ) return false;
    return is_equal(lhs._err, rhs._err);
  }

  // output stream
  friend std::ostream&
  operator<<(std::ostream& stream, const HitError& rhs);

};

//**********************************************************************

class HitDerivative {

friend class HitError;

private:
  trf::TrfMatrix _mtx;

  // constructor from implementation
  explicit HitDerivative(const trf::TrfMatrix& mtx);

public:

  // constructor
  explicit HitDerivative(int size =0);

  // constructor from an array
  // order is by completed rows (00, 01, 02, ..., 10, 12, ...)
  HitDerivative(int size, double* arr);

  // copy constructor
  HitDerivative( const HitDerivative& hder );

  // assignment
  // left side must be unassigned or have the same length
  HitDerivative& operator=(const HitDerivative& hder) {
    assert( size() == 0 || size() == hder.size() );
    _mtx = hder._mtx;
    return *this;
  }

  // Return the underlying matrix.
  const trf::TrfMatrix& get_matrix() const { return _mtx; };

  // return the # rows in the matrix
  int size() const { return _mtx.nrow(); };

#ifndef TRF_USING_LinearAlgebra
  // accessor
  double& operator()( int i, int j ) {
    assert( i >= 0 );
    assert( i < size() );
    assert( j >= 0 );
    assert( j < 5 );
    return _mtx(i,j);
  }
#endif

  // const accessor
  double operator()(int i, int j) const {
    assert( i >= 0 );
    assert( i < size() );
    assert( j >= 0 );
    assert( j < 5 );
    return _mtx(i,j);
  }

  // minimum
  double min() const { return _mtx.min(); };

  // maximum
  double max() const { return _mtx.max(); };

  // absolute minimum
  double amin() const { return _mtx.amin(); };

  // absolute maximum
  double amax() const { return _mtx.amax(); };

  // +=
  HitDerivative& operator+=(const HitDerivative& rhs) {
    _mtx += rhs._mtx;
    return *this;
  }

  // -=
  HitDerivative& operator-=(const HitDerivative& rhs) {
    _mtx -= rhs._mtx;
    return *this;
  }

  // derivative + derivative
  friend HitDerivative
  operator+(const HitDerivative& lhs, const HitDerivative& rhs) {
    return HitDerivative( lhs._mtx + rhs._mtx );
  }

  // derivative - derivative
  friend HitDerivative
  operator-(const HitDerivative& lhs, const HitDerivative& rhs) {
    return HitDerivative( lhs._mtx - rhs._mtx );
  }

  // equality
  // must have the same dimension and the same values
  friend bool
  operator==(const HitDerivative& lhs, const HitDerivative& rhs) {
    if ( lhs.size() != rhs.size() ) return false;
    return lhs._mtx == rhs._mtx;
  }

  // inequality
  friend bool
  operator!=(const HitDerivative& lhs, const HitDerivative& rhs) {
    return ! ( lhs == rhs );
  }

  // equality with tolerance
  friend bool
  is_equal(const HitDerivative& lhs, const HitDerivative& rhs) {
    if ( lhs.size() != rhs.size() ) return false;
    return is_equal(lhs._mtx, rhs._mtx);
  }

  // output stream
  friend std::ostream&
  operator<<(std::ostream& stream, const HitDerivative& rhs);

};

} //  end namespace trf

#endif