/******************************************************************************
* Header file for workload Track Fitting based on Kalman Filter.
* C classes and prodecure defined here.
* 
* $Header$
*
******************************************************************************/

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

#ifndef __CT_EXAMPLE_TRACK_FITTING_FIT_C_H_
#define __CT_EXAMPLE_TRACK_FITTING_FIT_C_H_

#ifdef SINGLE_PRECISION
typedef float ftype;
#else
typedef double ftype;
#endif //SINGLE_PRECISION
#include <cstring>

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

    FieldRegion(){
        x0 = x1 = x2 = y0 = y1 = y2 = z0 = z1 = z2 = z = 0.;
    }

    void set( const ftype h0[3], const ftype &h0z, 
              const ftype h1[3], const ftype &h1z, 
              const ftype h2[3], const ftype &h2z){
            z = h0z;
            ftype dz1 = h1z-h0z, dz2 = h2z-h0z;
            ftype det = 1.0f/(dz1*dz2*(dz2-dz1));
            ftype w21 = -dz2*det;
            ftype w22 = dz1*det;
            ftype w11 = -dz2*w21;
            ftype w12 = -dz1*w22;

            ftype dh1 = h1[0] - h0[0];
            ftype dh2 = h2[0] - h0[0];
            x0 = h0[0];
            x1 = dh1*w11 + dh2*w12 ;
            x2 = dh1*w21 + dh2*w22 ;

            dh1 = h1[1] - h0[1];
            dh2 = h2[1] - h0[1];
            y0 = h0[1];
            y1 = dh1*w11 + dh2*w12 ;
            y2 = dh1*w21 + dh2*w22  ;

            dh1 = h1[2] - h0[2];
            dh2 = h2[2] - h0[2];
            z0 = h0[2];
            z1 = dh1*w11 + dh2*w12 ;
            z2 = dh1*w21 + dh2*w22 ;  
    }

    void Shift( const ftype &newZ){
        ftype dz = newZ-z;
        ftype x2dz = x2*dz;
        ftype y2dz = y2*dz;
        ftype z2dz = z2*dz;
        z = newZ;
        x0+= (x1 + x2dz)*dz;
        x1+= x2dz+x2dz; 
        y0+= (y1 + y2dz)*dz;
        y1+= y2dz+y2dz; 
        z0+= (z1 + z2dz)*dz;
        z1+= z2dz+z2dz; 
    }

};

struct HitInfo{ // strip info
  void init( int ns );
  ~HitInfo();
  
  ftype *cos_phi, *sin_phi, *sigma2, *sigma216;
};

struct HitXYInfo{
  void init( int ns );
  ~HitXYInfo();
  
  ftype *C00, *C10, *C11;
};

//! Define stations (tracking detector) (SOA)
struct Stations{
    int nStations;

    ftype *z, *thick, *zhit, *RL,  *RadThick, *logRadThick, *SyF, *SyL;
    ftype (**mapX), (**mapY), (**mapZ); // polinom coeff.

   HitInfo UInfo, VInfo; // front and back
   HitXYInfo XYInfo;
  
    void initMap( int i, ftype *x, ftype *y, ftype *z ){
        memcpy( mapX[i], x, sizeof( ftype ) * MaxNFieldCoeff );
        memcpy( mapY[i], y, sizeof( ftype ) * MaxNFieldCoeff );
        memcpy( mapZ[i], z, sizeof( ftype ) * MaxNFieldCoeff );
    }

    Stations();
    Stations( int ns );

    void init( int ns );
    
    ~Stations();

    void field( int i, const ftype x, const ftype y, ftype H[3] ){
        ftype x2 = x*x;
        ftype y2 = y*y;
        ftype xy = x*y;
        ftype x3 = x2*x;
        ftype y3 = y2*y;
        ftype xy2 = x*y2;
        ftype x2y = x2*y;

        H[0] = mapX[i][0] +mapX[i][1]*x +mapX[i][2]*y +mapX[i][3]*x2 +mapX[i][4]*xy +mapX[i][5]*y2 +mapX[i][6]*x3 +mapX[i][7]*x2y +mapX[i][8]*xy2 +mapX[i][9]*y3;
        H[1] = mapY[i][0] +mapY[i][1]*x +mapY[i][2]*y +mapY[i][3]*x2 +mapY[i][4]*xy +mapY[i][5]*y2 +mapY[i][6]*x3 +mapY[i][7]*x2y +mapY[i][8]*xy2 +mapY[i][9]*y3;
        H[2] = mapZ[i][0] +mapZ[i][1]*x +mapZ[i][2]*y +mapZ[i][3]*x2 +mapZ[i][4]*xy +mapZ[i][5]*y2 +mapZ[i][6]*x3 +mapZ[i][7]*x2y +mapZ[i][8]*xy2 +mapZ[i][9]*y3;
    }

private:
    ftype *mapXPool, *mapYPool, *mapZPool;  //! to make the mem address contineous (for Ct copyin)
};


struct MCPoint{
  ftype x, y, z;
  ftype px, py, pz;
  int ista;
};

struct MCTracks
{
  MCPoint vPoints[maxNStations*2]; // won't be binded to ct-class, so isn't SOA
  int NMCPoints;
};

//! Define tracks (SOA)
class Tracks
{
public:
    Tracks();
    void init( int ns, int nt );
    ~Tracks();

    void setMC( int i, ftype * MC_para );
    void setHits( int i, int *iSta, ftype *hitsX, ftype *hitsY );

//private:
    unsigned int nStations;
    unsigned int nTracks;

    ftype *MC_x, *MC_y, *MC_z, *MC_px, *MC_py, *MC_pz, *MC_q;
//    ftype (*T)[6];    //for each track, T[6] = x, y, tx, ty, qp, z
//    ftype (*C)[15];   //for each track, C[15] is the covariance matrix

  
    int *nHits, *NDF;
    ftype **hitsX, **hitsY;
    ftype *hitsX2, *hitsY2;
    int **hitsIsta;

private:
    //! define three pool for hitsX, hitsY and hitsIsta to make the memory location is continueous.
    ftype *memXPool, *memYPool;
    int *memIstaPool;
};

void fitTracks( ftype (*T)[6], ftype (*C)[15] );

#endif //__CT_EXAMPLE_TRACK_FITTING_FIT_CT_H_