///
  /// @primary authors: S.Gorbunov; I.Kisel
  /// @        authors: H.Pabst et al. (Intel); M.Zyzak; I.Kulakov
  ///

#ifndef FitBase_H
#define FitBase_H

  /// common for all KF-approaches functions

#include "fit_interface_functional.h"


template<typename U, typename W>
class FitBase: public virtual FitFunctional<U,W> { // base class for all approaches
 public:

    /// Fit tracks
  void Fit( TracksCt<U> &ts, StationsCt<U> &ss, FieldRegionCt<U> &f, FitResults<U>& results ) const;

 private:
  void Combine( const dense<U> H[3], const U &w, dense<U> des[3] ) const;
};



template<typename U, typename W>
void FitBase<U,W>::Combine( const dense<U> H[3], const U &w, dense<U> des[3] ) const
{
    des[0] += w * (H[0] - des[0]);
    des[1] += w * (H[1] - des[1]);
    des[2] += w * (H[2] - des[2]);
}




//! main procedure of Ct track fitting
template<typename U, typename W>
void FitBase<U,W>::Fit( TracksCt<U> &ts, StationsCt<U> &ss, FieldRegionCt<U> &f, FitResults<U>& results ) const
{
    // upstream 
    GuessVec( ts, ss, results, 0 );

    // downstream
    //FieldRegionCt<U> f;
    dense<U>* T = results.VecT;
    dense<U> qp0 = T[4];

    isize i = vStations.nStations - 1;

    FilterFirst( ts, ss, i, results );
   
    AddMaterial( ss, i, qp0, results, false );

    U z1 = ss.z[i];
    dense<U> H0[3], H1[3], H2[3], HH[3];
    ss.field(i, T[0], T[1], H1);

    U w = 1.0f;   //h.w, since it's always 1.0.
    ss.field(i, ts.hitsX2.row( i ), ts.hitsY2.row( i ), HH);
    Combine( HH, w, H1 );
    
    U z2 = ss.z[i - 2];
    U dz = z2-z1;
    ss.field( i - 2, T[0] + T[2] * dz, T[1] + T[3] * dz, H2);
    ss.field( i - 2, ts.hitsX2.row( i - 2 ), ts.hitsY2.row( i - 2), HH);
    Combine( HH, w, H2 );

    _for( i -= 1, i >= 0, i-- ){
        U z0 = ss.z[i];
        dz = z1 - z0;
        ss.field( i, T[0] - T[2] * dz, T[1] - T[3] * dz, H0 );
        ss.field( i, ts.hitsX2.row( i ), ts.hitsY2.row( i ), HH );
        Combine( HH, w, H0 );
        //! note: FieldRegionCt f sets values here, needn't pass parameters
        f.set( H0, z0, H1, z1, H2, z2);

        // Extrapolate( results, ss.zhit[i], qp0, f, fill<FTYPE>(0,qp0.length()) );
        // AddMaterial( ts, ss, i, qp0, results );
        _if( i == 1 )
        {// 2nd MVD
          ExtrapolateWithMaterial( results, ss.zhit[i], qp0, f, ss, i, true, fill<FTYPE>(0,qp0.length()) ); // add pipe
        }
        _else {
          ExtrapolateWithMaterial( results, ss.zhit[i], qp0, f, ss, i, false, fill<FTYPE>(0,qp0.length()) );
        }_end_if;

        dense<U> u = ts.hitsX2.row( i )*ss.UInfo.cos_phi[i] + ts.hitsY2.row( i )*ss.UInfo.sin_phi[i]; 
 	dense<U> v = ts.hitsX2.row( i )*ss.VInfo.cos_phi[i] + ts.hitsY2.row( i )*ss.VInfo.sin_phi[i]; 
        Filter( ts, ss, ss.UInfo, i, u, fill<FTYPE>(1,qp0.length()), results );
        Filter( ts, ss, ss.VInfo, i, v, fill<FTYPE>(1,qp0.length()), results );
        for( int j = 0; j < 3; j ++ ){
            H2[j] = H1[j];
            H1[j] = H0[j];
        }
        z2 = z1;
        z1 = z0;
    }_end_for;
}


#endif