#ifndef FitClasses_H
#define FitClasses_H 1

  /// data classes, which are used for fit

const int MaxNStations = 10;



#if defined( X87 )

#include "openlab_mod/x87.h"
typedef double Fvec_t;
typedef bool Fvec_m_t;
const int vecN = 1;
typedef double  Single_t;
# define ALIGNMENT_NONE

#elif defined( VALARRAY )

#include <valarray>
#include "openlab_mod/ValVec.h"
const int vecN = ValVec::VEC_SIZE;
typedef ValVec Fvec_t;
typedef float Single_t;
# define ALIGNMENT_NONE

#elif defined( LRB )

#include "openlab_mod/F32LrbVec.h"
typedef F32LrbVec Fvec_t;
const int vecN = F32LrbVec::VECSIZE;
typedef float Single_t;
# define ALIGNMENT_64

#elif defined( SIMPLSIMD ) // header

#ifdef DOUBLE_PRECISION

typedef float Fvar_t;
#include "openlab_mod/P4_F64vec2.h"
typedef F64vec2 Fvec_t;
const int vecN = 2;
typedef double  Single_t;
# define ALIGNMENT_16

#else // DOUBLE_PRECISION

#include "openlab_mod/P4_F32vec4.h"
typedef F32vec4 Fvec_t;
const int vecN = 4;
typedef float  Single_t;
# define ALIGNMENT_16

#endif // DOUBLE_PRECISION

typedef Fvec_t Fvec_m_t;

#elif defined( NOTSIMD ) // pseudo header

#ifdef DOUBLE_PRECISION

typedef float Fvar_t;
#include "openlab_mod/PSEUDO_F64vec2.h"
typedef F64vec2 Fvec_t;
const int vecN = 2;
typedef double  Single_t;

#else

#include "openlab_mod/PSEUDO_F32vec4.h"
typedef F32vec4 Fvec_t;
const int vecN = 4;
typedef float  Single_t;

#endif // DOUBLE_PRECISION

typedef Fvec_t Fvec_m_t;

#else // if Vc

#define VC
#include "Vctypedef.h"
typedef float Single_t;
typedef float_m Fvec_m_t;

#endif

typedef int    Int_t;

struct FieldVector{
  Fvec_t X, Y, Z;
  void Combine( FieldVector &H, const Fvec_t &w ){
    X+= w*(H.X-X);
    Y+= w*(H.Y-Y);
    Z+= w*(H.Z-Z);
  }
} __attribute__ ((aligned(sizeof(Fvec_t))));


struct FieldSlice{
  Fvec_t X[21], Y[21], Z[21]; // polinom coeff.

  FieldSlice(){ for( int i=0; i<21; i++ ) X[i]=Y[i]=Z[i]=0; }

  void GetField( const Fvec_t &x, const Fvec_t &y, Fvec_t &Hx, Fvec_t &Hy, Fvec_t &Hz ){

    Fvec_t x2 = x*x;
    Fvec_t y2 = y*y;
    Fvec_t xy = x*y;
    Fvec_t x3 = x2*x;
    Fvec_t y3 = y2*y;
    Fvec_t xy2 = x*y2;
    Fvec_t x2y = x2*y;

    Fvec_t x4 = x3*x;
    Fvec_t y4 = y3*y;
    Fvec_t xy3 = x*y3;
    Fvec_t x2y2 = x2*y2;
    Fvec_t x3y = x3*y;

    Fvec_t x5 = x4*x;
    Fvec_t y5 = y4*y;
    Fvec_t xy4 = x*y4;
    Fvec_t x2y3 = x2*y3;
    Fvec_t x3y2 = x3*y2;
    Fvec_t x4y = x4*y;

    Hx = X[0] +X[1]*x +X[2]*y +X[3]*x2 +X[4]*xy +X[5]*y2 +X[6]*x3 +X[7]*x2y +X[8]*xy2 +X[9]*y3
       + X[10]*x4 + X[11]*x3y +X[12]*x2y2 +X[13]*xy3 + X[14]*y4
       + X[15]*x5 + X[16]*x4y +X[17]*x3y2 +X[18]*x2y3 +X[19]*xy4 +X[20]*y5;

    Hy = Y[0] +Y[1]*x +Y[2]*y +Y[3]*x2 +Y[4]*xy +Y[5]*y2 +Y[6]*x3 +Y[7]*x2y +Y[8]*xy2 +Y[9]*y3
       + Y[10]*x4 + Y[11]*x3y +Y[12]*x2y2 +Y[13]*xy3 + Y[14]*y4
       + Y[15]*x5 +Y[16]*x4y +Y[17]*x3y2 +Y[18]*x2y3 +Y[19]*xy4 +Y[20]*y5;

    Hz = Z[0] +Z[1]*x +Z[2]*y +Z[3]*x2 +Z[4]*xy +Z[5]*y2 +Z[6]*x3 +Z[7]*x2y +Z[8]*xy2 +Z[9]*y3
       + Z[10]*x4 + Z[11]*x3y +Z[12]*x2y2 +Z[13]*xy3 + Z[14]*y4
       + Z[15]*x5 +Z[16]*x4y +Z[17]*x3y2 +Z[18]*x2y3 +Z[19]*xy4 +Z[20]*y5;
  }

  void GetField( const Fvec_t &x, const Fvec_t &y, FieldVector &H ){
    GetField( x, y, H.X, H.Y, H.Z );
  }
} __attribute__ ((aligned(sizeof(Fvec_t))));


struct FieldRegion{
  Fvec_t x0, x1, x2 ; // Hx(Z) = x0 + x1*(Z-z) + x2*(Z-z)^2
  Fvec_t y0, y1, y2 ; // Hy(Z) = y0 + y1*(Z-z) + y2*(Z-z)^2
  Fvec_t z0, z1, z2 ; // Hz(Z) = z0 + z1*(Z-z) + z2*(Z-z)^2
  Fvec_t z;

  friend ostream& operator<<( ostream &os, const FieldRegion &a ) {
    os << a.x0 << endl
       << a.x1 << endl
       << a.x2 << endl
       << a.y0 << endl
       << a.y1 << endl
       << a.y2 << endl
       << a.z0 << endl
       << a.z1 << endl
       << a.z2 << endl
       << a.z;
    return os;
  }
  
  FieldRegion(){
    x0 = x1 = x2 = y0 = y1 = y2 = z0 = z1 = z2 = z = 0.;
  }

  void Get(const Fvec_t z_, Fvec_t* B) const{
    Fvec_t dz = (z_ - z);
    Fvec_t dz2 = dz*dz;
    B[0] = x0 + x1*dz + x2*dz2;
    B[1] = y0 + y1*dz + y2*dz2;
    B[2] = z0 + z1*dz + z2*dz2;
  }
  
  void Set( const FieldVector &H0, const Fvec_t &H0z, 
	    const FieldVector &H1, const Fvec_t &H1z, 
	    const FieldVector &H2, const Fvec_t &H2z){
    z = H0z;
    Fvec_t dz1 = H1z-H0z, dz2 = H2z-H0z;
    Fvec_t det = rcp(dz1*dz2*(dz2-dz1));
    Fvec_t w21 = -dz2*det;
    Fvec_t w22 = dz1*det;
    Fvec_t w11 = -dz2*w21;
    Fvec_t w12 = -dz1*w22;
    
    Fvec_t dH1 = H1.X - H0.X;
    Fvec_t dH2 = H2.X - H0.X;
    x0 = H0.X;
    x1 = dH1*w11 + dH2*w12 ;
    x2 = dH1*w21 + dH2*w22 ;
      
    dH1 = H1.Y - H0.Y;
    dH2 = H2.Y - H0.Y;
    y0 = H0.Y;
    y1 = dH1*w11 + dH2*w12 ;
    y2 = dH1*w21 + dH2*w22  ;

    dH1 = H1.Z - H0.Z;
    dH2 = H2.Z - H0.Z;
    z0 = H0.Z;
    z1 = dH1*w11 + dH2*w12 ;
    z2 = dH1*w21 + dH2*w22 ;         
  }
  
  void Shift( const Fvec_t &Z0){
    Fvec_t dz = Z0-z;
    Fvec_t x2dz = x2*dz;
    Fvec_t y2dz = y2*dz;
    Fvec_t z2dz = z2*dz;
    z = Z0;
    x0+= (x1 + x2dz)*dz;
    x1+= x2dz+x2dz; 
    y0+= (y1 + y2dz)*dz;
    y1+= y2dz+y2dz; 
    z0+= (z1 + z2dz)*dz;
    z1+= z2dz+z2dz; 
  }

} __attribute__ ((aligned(sizeof(Fvec_t))));


struct HitInfo{ // strip info
  Fvec_t cos_phi, sin_phi, sigma2, sigma216;
};

struct HitXYInfo{
  Fvec_t C00, C10, C11;
};

struct Station{
  Fvec_t z, thick, zhit, RL,  RadThick, logRadThick, 
    SyF, SyL; //  field intergals with respect to First(last) station 

  HitInfo UInfo, VInfo; // front and back
  HitXYInfo XYInfo;
  
  FieldSlice Map;
} __attribute__ ((aligned(sizeof(Fvec_t))));


struct Hit{
  Single_t x, y;
  Int_t ista;
  Single_t tmp1;
} __attribute__ ((aligned(sizeof(Fvec_t))));

struct MCPoint{
  Single_t x, y, z;
  Single_t px, py, pz;
  Int_t ista;
} __attribute__ ((aligned(sizeof(Fvec_t))));

struct MCTrack{
  Single_t MC_x, MC_y, MC_z, MC_px, MC_py, MC_pz, MC_q;
  MCPoint vPoints[MaxNStations*2];
  int NMCPoints;
} __attribute__ ((aligned(sizeof(Fvec_t))));

struct Track{
  Int_t NHits;
  Hit vHits[MaxNStations];

#ifdef DAF
  Single_t Ts[MaxNStations][6];
  Single_t Cs[MaxNStations][15];

  Single_t Chi2s[MaxNStations];
#endif
  
  Single_t T[6]; // x, y, tx, ty, qp, z  
  Single_t C[15]; // cov matr.
  Single_t Chi2;
  Int_t NDF;
} __attribute__ ((aligned(sizeof(Fvec_t))));

struct HitV{
  Fvec_t x, y, w;
  FieldVector H;
} __attribute__ ((aligned(sizeof(Fvec_t))));



#if defined(CONVENTIONAL)

struct CovV{
  Fvec_t C00, 
    C10, C11, 
    C20, C21, C22, 
    C30, C31, C32, C33, 
    C40, C41, C42, C43, C44;

#ifdef DAF
  CovV & operator += (const CovV &A)
  {
    Fvec_t *c = &C00;
    const Fvec_t *a = &A.C00;
    for(int i=0; i<15; i++) c[i] += a[i];
    return *this;
  }
  CovV operator + (const CovV &A) const
  {
    CovV R = *this;
    R += A;
    return R;
  }
  CovV & operator -= (const CovV &A)
  {
    Fvec_t *c = &C00;
    const Fvec_t *a = &A.C00;
    for(int i=0; i<15; i++) c[i] -= a[i];
    return *this;
  }
  CovV operator - (const CovV &A) const
  {
    CovV R = *this;
    R -= A;
    return R;
  }

  CovV & operator *= (const Single_t &w)
  {
    Fvec_t *c = &C00;
    for(int i=0; i<15; i++) c[i] *= w;
    return *this;
  }

  CovV & operator *= (const Fvec_t &w)
  {
    Fvec_t *c = &C00;
    for(int i=0; i<15; i++) c[i] *= w;
    return *this;
  }

  const Fvec_t &GetC00() const { return C00; };
  const Fvec_t &GetC10() const { return C10; };
  const Fvec_t &GetC11() const { return C11; };
#endif // DAF

  friend ostream& operator<<( ostream &os, const CovV &a ) {
    os << a.C00 << endl
       << a.C10 << endl
       << a.C11 << endl
       << a.C20 << endl
       << a.C21 << endl
       << a.C22 << endl
       << a.C30 << endl
       << a.C31 << endl
       << a.C32 << endl
       << a.C33 << endl
       << a.C40 << endl
       << a.C41 << endl
       << a.C42 << endl
       << a.C43 << endl
       << a.C44;
    return os;
  }
} __attribute__ ((aligned(sizeof(Fvec_t))));


#ifdef DAF
struct MatV{ // matrix used for multiplication
  Fvec_t C[5][5];
} __attribute__ ((aligned(sizeof(Fvec_t))));
#endif // DAF

typedef CovV CovVConventional;

#else

#if defined(SQUARE_ROOT1) || defined(SQUARE_ROOT2)

struct CovV{
  Fvec_t       S00,
    S10, S11,
    S20, S21, S22,
    S30, S31, S32, S33,
    S40, S41, S42, S43, S44;
  
#ifdef DAF
  Fvec_t GetC00() const { return S00*S00; };
  Fvec_t GetC10() const { return S10*S00; };
  Fvec_t GetC11() const { return S10*S10 + S11*S11; };
#endif // DAF

  friend ostream& operator<<( ostream &os, const CovV &a ) {
    os << a.S00 << endl
       << a.S10 << endl
       << a.S11 << endl
       << a.S20 << endl
       << a.S21 << endl
       << a.S22 << endl
       << a.S30 << endl
       << a.S31 << endl
       << a.S32 << endl
       << a.S33 << endl
       << a.S40 << endl
       << a.S41 << endl
       << a.S42 << endl
       << a.S43 << endl
       << a.S44;
    return os;
  }

  CovV &operator*=( const CovV& a ){
    S40 = S40*a.S00 + S41*a.S10 + S42*a.S20 + S43*a.S30 + S44*a.S40;
    S30 = S30*a.S00 + S31*a.S10 + S32*a.S20 + S33*a.S30;
    S20 = S20*a.S00 + S21*a.S10 + S22*a.S20;
    S10 = S10*a.S00 + S11*a.S10;
    S00 = S00*a.S00;
    S41 =             S41*a.S11 + S42*a.S21 + S43*a.S31 + S44*a.S41;
    S31 =             S31*a.S11 + S32*a.S21 + S33*a.S31;
    S21 =             S21*a.S11 + S22*a.S21;
    S11 =             S11*a.S11;
    S42 =                         S42*a.S22 + S43*a.S32 + S44*a.S42;
    S32 =                         S32*a.S22 + S33*a.S32;
    S22 =                         S22*a.S22;
    S43 =                                     S43*a.S33 + S44*a.S43;
    S33 =                                     S33*a.S33;
    S44 =                                                 S44*a.S44;
    
    return *this;
  }
  
} __attribute__ ((aligned(sizeof(Fvec_t))));

#elif defined(SQUARE_ROOT1S)

struct CovV{
  Fvec_t
    S00, S01, S02, S03, S04,
    S10, S11, S12, S13, S14,
    S20, S21, S22, S23, S24,
    S30, S31, S32, S33, S34,
    S40, S41, S42, S43, S44;
  
#ifdef DAF
  Fvec_t GetC00() const { return S00*S00 + S01*S01 + S02*S02 + S03*S03 + S04*S04; };
  Fvec_t GetC10() const { return S10*S00 + S11*S01 + S12*S02 + S13*S03 + S14*S04; };
  Fvec_t GetC11() const { return S10*S00 + S11*S11 + S12*S12 + S13*S13 + S14*S14; };
#endif // DAF

  friend ostream& operator<<( ostream &os, const CovV &a ) {
    os << a.S00 << endl
       << a.S01 << endl
       << a.S02 << endl
       << a.S03 << endl
       << a.S04 << endl
       << a.S10 << endl
       << a.S11 << endl
       << a.S12 << endl
       << a.S13 << endl
       << a.S14 << endl
       << a.S20 << endl
       << a.S21 << endl
       << a.S22 << endl
       << a.S23 << endl
       << a.S24 << endl
       << a.S30 << endl
       << a.S31 << endl
       << a.S32 << endl
       << a.S33 << endl
       << a.S34 << endl
       << a.S40 << endl
       << a.S41 << endl
       << a.S42 << endl
       << a.S43 << endl
       << a.S44;
     
    return os;
  }

  // CovV &operator*=( const CovV& a ){
  //   S40 = S40*a.S00 + S41*a.S10 + S42*a.S20 + S43*a.S30 + S44*a.S40;
  //   S30 = S30*a.S00 + S31*a.S10 + S32*a.S20 + S33*a.S30;
  //   S20 = S20*a.S00 + S21*a.S10 + S22*a.S20;
  //   S10 = S10*a.S00 + S11*a.S10;
  //   S00 = S00*a.S00;
  //   S41 =             S41*a.S11 + S42*a.S21 + S43*a.S31 + S44*a.S41;
  //   S31 =             S31*a.S11 + S32*a.S21 + S33*a.S31;
  //   S21 =             S21*a.S11 + S22*a.S21;
  //   S11 =             S11*a.S11;
  //   S42 =                         S42*a.S22 + S43*a.S32 + S44*a.S42;
  //   S32 =                         S32*a.S22 + S33*a.S32;
  //   S22 =                         S22*a.S22;
  //   S43 =                                     S43*a.S33 + S44*a.S43;
  //   S33 =                                     S33*a.S33;
  //   S44 =                                                 S44*a.S44;
    
  //   return *this;
  // }
  
} __attribute__ ((aligned(sizeof(Fvec_t))));

#elif defined(UD_FILTERING_V) || defined(UD_FILTERING_S)

struct CovV{
  Fvec_t U01, U02, U03, U04,
    U12, U13, U14,
    U23, U24,
    U34,
    D00, D11, D22, D33, D44;

#ifdef DAF
  Fvec_t GetC00() const { return D00 + U01*D11*U01 + U02*D22*U02 + U03*D33*U03 + U04*D44*U04; };
  Fvec_t GetC10() const { return           D11*U01 + U12*D22*U02 + U13*D33*U03 + U14*D44*U04; };
  Fvec_t GetC11() const { return           D11     + U12*D22*U12 + U13*D33*U13 + U14*D44*U14; };
#endif // DAF

  friend ostream& operator<<( ostream &os, const CovV &a ) {
    os << a.D00 << endl
       << a.D11 << endl
       << a.D22 << endl
       << a.D33 << endl
       << a.D44 << endl
       << a.U01 << endl
       << a.U02 << endl
       << a.U03 << endl
       << a.U04 << endl
       << a.U12 << endl
       << a.U13 << endl
       << a.U14 << endl
       << a.U23 << endl
       << a.U24 << endl
       << a.U34;
    return os;
  }
  
} __attribute__ ((aligned(sizeof(Fvec_t))));

#endif

struct CovVConventional{ // used for convertion
  Fvec_t  C00, 
    C10, C11, 
    C20, C21, C22, 
    C30, C31, C32, C33, 
    C40, C41, C42, C43, C44;
  CovVConventional() {};

  friend ostream& operator<<( ostream &os, const CovVConventional &a ) {
    os << a.C00 << endl
       << a.C10 << endl
       << a.C11 << endl
       << a.C20 << endl
       << a.C21 << endl
       << a.C22 << endl
       << a.C30 << endl
       << a.C31 << endl
       << a.C32 << endl
       << a.C33 << endl
       << a.C40 << endl
       << a.C41 << endl
       << a.C42 << endl
       << a.C43 << endl
       << a.C44;
    return os;
  }

  
#if defined(SQUARE_ROOT1) || defined(SQUARE_ROOT2)
  CovVConventional( const CovV& T ) {
      // C = S*S'
    C00 = T.S00*T.S00;
    C10 = T.S10*T.S00;
    C11 = T.S10*T.S10 + T.S11*T.S11;
    C20 = T.S20*T.S00;
    C21 = T.S20*T.S10 + T.S21*T.S11;
    C22 = T.S20*T.S20 + T.S21*T.S21 + T.S22*T.S22;
    C30 = T.S30*T.S00;
    C31 = T.S30*T.S10 + T.S31*T.S11;
    C32 = T.S30*T.S20 + T.S31*T.S21 + T.S32*T.S22;
    C33 = T.S30*T.S30 + T.S31*T.S31 + T.S32*T.S32 + T.S33*T.S33;
    C40 = T.S40*T.S00;
    C41 = T.S40*T.S10 + T.S41*T.S11;
    C42 = T.S40*T.S20 + T.S41*T.S21 + T.S42*T.S22;
    C43 = T.S40*T.S30 + T.S41*T.S31 + T.S42*T.S32 + T.S43*T.S33;
    C44 = T.S40*T.S40 + T.S41*T.S41 + T.S42*T.S42 + T.S43*T.S43 + T.S44*T.S44;
  }

  operator CovV(){
    CovV r;
    
    r.S00 = sqrt( C00 );
    const Fvec_t iS00 = 1.f/r.S00;
    r.S10 = iS00*C10;

    r.S11 = sqrt( C11 - r.S10*r.S10 );
    const Fvec_t iS11 = 1.f/r.S11;
    r.S20 = iS00*C20;
    r.S21 = iS11*( C21 - r.S20*r.S10 );

    r.S22 = sqrt( C22 - r.S20*r.S20 - r.S21*r.S21 );
    const Fvec_t iS22 = 1.f/r.S22;
    r.S30 = iS00*C30;
    r.S31 = iS11*( C31 - r.S30*r.S10 );
    r.S32 = iS22*( C32 - r.S30*r.S20 - r.S31*r.S21 );

    r.S33 = sqrt( C33 - r.S30*r.S30 - r.S31*r.S31 - r.S32*r.S32 );
    const Fvec_t iS33 = 1.f/r.S33;
    r.S40 = iS00*C40;
    r.S41 = iS11*( C41 - r.S40*r.S10 );
    r.S42 = iS22*( C42 - r.S40*r.S20 - r.S41*r.S21 );
    r.S43 = iS33*( C43 - r.S40*r.S20 - r.S41*r.S21 - r.S42*r.S32 );

    r.S44 = sqrt( C44 - r.S40*r.S40 - r.S41*r.S41 - r.S42*r.S42 - r.S43*r.S43 );

    return r;
  }
#elif defined(SQUARE_ROOT1S)
  CovVConventional( const CovV& T ) {
      // C = S*S'
#define M(i,j) C##i##j = T.S##i##0*T.S##j##0 + T.S##i##1*T.S##j##1 + T.S##i##2*T.S##j##2 + T.S##i##3*T.S##j##3 + T.S##i##4*T.S##j##4

    M(0,0);
    M(1,0);
    M(1,1);
    M(2,0);
    M(2,1);
    M(2,2);
    M(3,0);
    M(3,1);
    M(3,2);
    M(3,3);
    M(4,0);
    M(4,1);
    M(4,2);
    M(4,3);
    M(4,4);
#undef M
  }

#elif defined(UD_FILTERING_V) || defined(UD_FILTERING_S)
  CovVConventional( const CovV& T ) {
      // DU = D*U'
      // C = U*D*U'= U*DU
    C44 = T.D44;
    C43 = T.D44*T.U34;
    C42 = T.D44*T.U24;
    C41 = T.D44*T.U14;
    C40 = T.D44*T.U04;

    const Fvec_t &DU43 = C43;
    const Fvec_t &DU42 = C42;
    const Fvec_t &DU41 = C41;
    const Fvec_t &DU40 = C40;
#define DU_DEF(i,j) const Fvec_t DU##i##j = T.D##i##i*T.U##j##i
    DU_DEF(3,2);
    DU_DEF(3,1);
    DU_DEF(3,0);
    DU_DEF(2,1);
    DU_DEF(2,0);
    DU_DEF(1,0);
#undef DU_DEF
  
    C33 =                                         T.D33 + T.U34*DU43;
    C32 =                                          DU32 + T.U34*DU42;
    C31 =                                          DU31 + T.U34*DU41;
    C30 =                                          DU30 + T.U34*DU40;
    C22 =                            T.D22 + T.U23*DU32 + T.U24*DU42;
    C21 =                             DU21 + T.U23*DU31 + T.U24*DU41;
    C20 =                             DU20 + T.U23*DU30 + T.U24*DU40;
    C11 =               T.D11 + T.U12*DU21 + T.U13*DU31 + T.U14*DU41;
    C10 =                DU10 + T.U12*DU20 + T.U13*DU30 + T.U14*DU40;
    C00 =  T.D00 + T.U01*DU10 + T.U02*DU20 + T.U03*DU30 + T.U04*DU40;
  }

  operator CovV(){
    CovV r;
    
    r.D44 = C44;
    const Fvec_t iD44 = 1.f/r.D44;
    const Fvec_t &DU43 = C43;
    const Fvec_t &DU42 = C42;
    const Fvec_t &DU41 = C41;
    const Fvec_t &DU40 = C40;
    r.U34 = DU43*iD44;
    r.U24 = DU42*iD44;
    r.U14 = DU41*iD44;
    r.U04 = DU40*iD44;

    r.D33 = C33 - r.U34*DU43;
    const Fvec_t iD33 = 1.f/r.D33;
    const Fvec_t DU32 = C32 - r.U34*DU42;
    const Fvec_t DU31 = C31 - r.U34*DU41;
    const Fvec_t DU30 = C30 - r.U34*DU40;
    r.U23 = DU32*iD33;
    r.U13 = DU31*iD33;
    r.U03 = DU30*iD33;

    r.D22 = C22 - r.U23*DU32 - r.U24*DU42;
    const Fvec_t iD22 = 1.f/r.D22;
    const Fvec_t DU21 = C21 - r.U23*DU31 - r.U24*DU41;
    const Fvec_t DU20 = C20 - r.U23*DU30 - r.U24*DU40;
    r.U12 = DU21*iD22;
    r.U02 = DU20*iD22;

    r.D11 = C11 - r.U12*DU21 - r.U13*DU31 - r.U14*DU41;
    const Fvec_t DU10 = C10 - r.U12*DU20 - r.U13*DU30 - r.U14*DU40;
    r.U01 = DU10/r.D11;

    r.D00 = C00 - r.U01*DU10 - r.U02*DU20 - r.U03*DU30 - r.U04*DU40;
    
    return r;
  }
  
#endif

  CovVConventional& operator+=(const CovVConventional &A){
    Fvec_t *c = &C00;
    const Fvec_t *a = &A.C00;
    for(int i=0; i<15; i++) c[i] += a[i];
    return *this;
  }
  CovVConventional operator+(CovVConventional &A){
    CovVConventional R = *this;
    R += A;
    return R;
  }
  CovVConventional& operator-=(const CovVConventional &A){
    Fvec_t *c = &C00;
    const Fvec_t *a = &A.C00;
    for(int i=0; i<15; i++) c[i] -= a[i];
    return *this;
  }
} __attribute__ ((aligned(sizeof(Fvec_t))));


#endif // !defined(CONVENTIONAL)

struct TrackV{
  HitV vHits[MaxNStations];

#ifdef DAF
  Fvec_t Ts[MaxNStations][6];
  CovV   Cs[MaxNStations];
  FieldRegion fs[MaxNStations];

  Fvec_t Chi2s[MaxNStations];

  Fvec_t w[MaxNStations]; // weight = 1 / ( 1 + exp ( ( chi2 - chi2_cut ) / 2T ) )
#endif
  
  Fvec_t T[6]; // x, y, tx, ty, qp, z
  CovV   C;    // cov matr.

  Fvec_t Chi2;
  Fvec_t NDF;
  
  FieldRegion f; // field at first hit (needed for extrapolation to MC and check of results)

#ifdef DAF
  friend ostream& operator<<( ostream &os, const TrackV &a ) {
    for ( int is = 0; is < MaxNStations; ++is ) {
      os << "sta: " << is << endl
         << a.Ts[is][0] << endl
         << a.Ts[is][1] << endl
         << a.Ts[is][2] << endl
         << a.Ts[is][3] << endl
         << a.Ts[is][4] << endl
         << a.Ts[is][5] << endl
         << a.Cs[is] << endl;
    }
    return os;
  }
#endif
  
} __attribute__ ((aligned(sizeof(Fvec_t))));

#endif