#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstring>
#include <stdlib.h>
using namespace std;
#ifdef _WIN32
#include <float.h> // _finite
#define finite _finite
#endif

#define TOTALTIME
#include "Stopwatch.h"
#include <vector>

#ifdef DEBUG_THREADS
  // const int coeff = 1;
const int MaxNTracks = 20000;//20000
const int NFits = 1000;//10;
const int Ntimes = 1;//int(3000.*coeff);//100;
#else // if DEBUG_THREADS
const int NFits = 10;

#if defined(TBB) || defined(OMP)
const int MaxNTracks = 80000;
const int Ntimes = 100;
#else
const int MaxNTracks = 20000;
const int Ntimes = 100;//int(3000.*coeff);//100;
#endif // TBB
#endif //  ifndef DEBUG_THREADS


#if defined(CONVENTIONAL)
#ifdef DAF
#include "DAFC.h"
typedef  DAFC FitInterface;
#else
#include "FitC.h"
typedef  FitC FitInterface;
#endif
#elif defined(SQUARE_ROOT1S) || defined(SQUARE_ROOT1) || defined(SQUARE_ROOT2)
#ifdef DAF
#include "DAFSR.h"
typedef  DAFSR FitInterface;
#else
#include "FitSR.h"
typedef  FitSR FitInterface;
#endif
#elif defined(UD_FILTERING_V) || defined(UD_FILTERING_S)
#ifdef DAF
#include "DAFUD.h"
typedef  DAFUD FitInterface;
#else
#include "FitUD.h"
typedef  FitUD FitInterface;
#endif
#endif


  /* Use Intel TBB for multithreading */
#ifdef TBB
  //  #define DEBUG_THREADS  // tbb with output 
#include "openlab_mod/multithread.h"
#include "openlab_mod/parallel_for_simpleEqNThr.h"
#include "openlab_mod/parallel_for_simpleHT.h"
#include "openlab_mod/parallel_for_simple.h"

#include "tbb/task_scheduler_init.h"
#include "tbb/parallel_for.h"
#include "tbb/blocked_range.h"
using namespace tbb;
#endif // TBB


  // default
#if !(defined(ANALYTIC_EXTRAPOLATION) || defined(RK4_EXTRAPOLATION))
#define ANALYTIC_EXTRAPOLATION
#endif
#if !(defined(CONVENTIONAL) || defined(SQUARE_ROOT1) || defined(SQUARE_ROOT1S) || defined(SQUARE_ROOT2) || defined(UD_FILTERING_V) || defined(UD_FILTERING_S))
#define CONVENTIONAL
#endif

#define MUTE

#if defined(TBB) || defined(OMP)
#define COPY // copy NTracksPerThread tracks as more times as it's needed
const int NTracksPerThread = 1000;
#endif

#ifdef OMP
#include <omp.h>
#include "pthread.h"
#include "errno.h"
#define handle_error_en(en, msg) do { errno = en; perror(msg); exit(EXIT_FAILURE); } while (0)
#endif

#include <ctime>
#include <string>
string dataDir = "./data/";

int tasks = 16;  /* #threads <= #tasks */

#ifdef TBB
const bool useObserver = 1;  // 1 for use thread-cpu correspondent
const float speedUpHT = 1.04; // = time1logCore / time2logCoreOn1PhysCore. Need for optimize proc. usage.
#endif // TBB



Station* vStations;

Track vTracks[MaxNTracks];
MCTrack vMCTracks[MaxNTracks];
int NStations = 0;
int NTracks = 0;
int NTracksV = 0;



FieldRegion field0;
FitInterface fitter;

void ReadInput(){

  fstream FileGeo, FileTracks, FileMCTracks;

  FileGeo.open( (dataDir+"geo.dat").c_str(), ios::in );
  FileTracks.open( (dataDir+"tracks.dat").c_str(), ios::in );
  FileMCTracks.open( (dataDir+"mctracksin.dat").c_str(), ios::in );
  {
    FieldVector H[3];
    Fvec_t Hz[3];
    for( int i=0; i<3; i++){
      Single_t Bx, By, Bz, z;
      FileGeo >> z >> Bx >> By >> Bz;
      Hz[i] = z; H[i].X = Bx;   H[i].Y = By; H[i].Z = Bz; 
#ifndef MUTE
      cout<<"Input Magnetic field:"<<z<<" "<<Bx<<" "<<By<<" "<<Bz<<endl;
#endif
    }                   
    field0.Set(H[0],Hz[0], H[1],Hz[1], H[2],Hz[2] );
  }
  FileGeo >> NStations;
#ifndef MUTE
  cout<<"Input "<<NStations<<" Stations:"<<endl;
#endif
  
  for( int i=0; i<NStations; i++ ){
    int ist;    
    FileGeo >> ist;
    if( ist!=i ) break;
    Station &st = vStations[i];

    FileGeo >> st.z >> st.thick >> st.RL >> st.UInfo.sigma2 >> st.VInfo.sigma2;
    st.UInfo.sigma2 *= st.UInfo.sigma2;
    st.VInfo.sigma2 *= st.VInfo.sigma2;
    st.UInfo.sigma216 = st.UInfo.sigma2*16.f;
    st.VInfo.sigma216 = st.VInfo.sigma2*16.f;
    
    if( i < 2 ){ // mvd //TODO From Geo File!!!
      st.UInfo.cos_phi = 1.f;
      st.UInfo.sin_phi = 0.f;
      st.VInfo.cos_phi = 0.f;
      st.VInfo.sin_phi = 1.f;
    }
    else{
      st.UInfo.cos_phi = 1.f;           // 0 degree
      st.UInfo.sin_phi = 0.f;
      st.VInfo.cos_phi = 0.9659258244f; // 15 degree
      st.VInfo.sin_phi = 0.2588190521f;
    }

    Fvec_t idet = st.UInfo.cos_phi*st.VInfo.sin_phi - st.UInfo.sin_phi*st.VInfo.cos_phi;
    idet = 1.f/(idet*idet);
    st.XYInfo.C00 = ( st.VInfo.sin_phi*st.VInfo.sin_phi*st.UInfo.sigma2 +
                      st.UInfo.sin_phi*st.UInfo.sin_phi*st.VInfo.sigma2 )*idet;
    st.XYInfo.C10 =-( st.VInfo.sin_phi*st.VInfo.cos_phi*st.UInfo.sigma2 +
                      st.UInfo.sin_phi*st.UInfo.cos_phi*st.VInfo.sigma2 )*idet;
    st.XYInfo.C11 = ( st.VInfo.cos_phi*st.VInfo.cos_phi*st.UInfo.sigma2 +
                      st.UInfo.cos_phi*st.UInfo.cos_phi*st.VInfo.sigma2 )*idet;
    
#ifndef MUTE
#ifdef X87
    cout<<"    "<<st.z <<" "<<st.thick<<" "<<st.RL<<", ";
#else
    cout<<"    "<<st.z[0] <<" "<<st.thick[0]<<" "<<st.RL[0]<<", ";
#endif
#endif
    st.zhit = st.z;
    st.RadThick = st.thick/st.RL;
    st.logRadThick = log(st.RadThick);

    int N = 0;
    FileGeo >> N;
#ifndef MUTE
    cout<<N<<" field coeff."<<endl;
#endif
    for( int i=0; i<N; i++ ) FileGeo >> st.Map.X[i];
    for( int i=0; i<N; i++ ) FileGeo >> st.Map.Y[i];
    for( int i=0; i<N; i++ ) FileGeo >> st.Map.Z[i];
  }
  {
      // field intergals with respect to Last station
    Fvec_t z0  = vStations[NStations-1].z;
    Fvec_t sy = 0.f, Sy = 0.f;
    for( int i=NStations-1; i>=0; i-- ){
      Station &st = vStations[i];
      Fvec_t dz = st.z-z0;
      Fvec_t Hy = vStations[i].Map.Y[0];
      Sy += dz*sy + dz*dz*Hy/2.f;
      sy += dz*Hy;
      st.SyL = Sy;
      z0 = st.z;
    }
      // field intergals with respect to First station
    z0 = vStations[0].z;
    sy = 0.f, Sy = 0.f;
    for( int i=0; i<NStations; i++ ){
      Station &st = vStations[i];
      Fvec_t dz = st.z-z0;
      Fvec_t Hy = vStations[i].Map.Y[0];
      Sy += dz*sy + dz*dz*Hy/2.f;
      sy += dz*Hy;
      st.SyF = Sy;
      z0 = st.z;
    }
  }

  FileGeo.close();

  NTracks = 0;
  int TrackIndex[MaxNTracks];
  while( !FileTracks.eof() ){
    
    int itr;
    FileTracks>>itr;
      //if( itr!=NTracks ) break;
    if( NTracks>=MaxNTracks ) break;

    Track &t = vTracks[NTracks];
    MCTrack &mc = vMCTracks[NTracks];
    FileTracks >> mc.MC_x   >> mc.MC_y  >> mc.MC_z 
               >> mc.MC_px >> mc.MC_py >> mc.MC_pz >> mc.MC_q
               >> t.NHits;
    for( int i=0; i<t.NHits; i++ ){
      int ist;
      FileTracks >> ist;
      t.vHits[i].ista = ist;
      FileTracks >> t.vHits[i].x >> t.vHits[i].y;
    }
    TrackIndex[NTracks] = itr;
    if( t.NHits==NStations )   NTracks++;
  }
  int NMCTracks=0;
  int iPoint = 0;
  while( !FileMCTracks.eof() ){
    
    int itr;
    FileMCTracks>>itr;
      //if( itr!=NTracks ) break;
    if( NMCTracks>=MaxNTracks ) break;
    MCTrack &mc = vMCTracks[NMCTracks];
    Single_t temp;
    int NMCPoints;
    FileMCTracks >> temp   >> temp  >> temp
                 >> temp >> temp >> temp >> temp
                 >> NMCPoints;
    mc.NMCPoints = NMCPoints;
    for( int i=0; i<NMCPoints; i++ ){
      int ist;
      FileMCTracks >> ist;
      mc.vPoints[i].ista = ist;
      FileMCTracks >> mc.vPoints[i].x >> mc.vPoints[i].y>>mc.vPoints[i].z>>mc.vPoints[i].px>>mc.vPoints[i].py>>mc.vPoints[i].pz;

#ifdef DAF
        // if(ist == 5) { // shift  for DAF
        //   const float coeff = 10.f;
        //   // vTracks[NMCTracks].vHits[ist].x = mc.vPoints[i].x + sqrt(vStations[ist].UInfo.sigma2[0]) * coeff;
        //   // cout << sqrt(vStations[ist].UInfo.sigma2[0]) << " U " << sqrt(vStations[ist].VInfo.sigma2[0]) << endl;
        //   // cout << sqrt(vStations[ist].XYInfo.C00[0]) << " X " << sqrt(vStations[ist].XYInfo.C11[0]) << endl;
        //   vTracks[NMCTracks].vHits[ist].x = mc.vPoints[i].x + sqrt(vStations[ist].XYInfo.C00[0]) * coeff;
        //   vTracks[NMCTracks].vHits[ist].y = mc.vPoints[i].y + sqrt(vStations[ist].XYInfo.C11[0]) * coeff;
        // }
#endif
    }
#ifdef DAF
      // Track &t = vTracks[NMCTracks];
      // for( int i=0; i<NMCPoints; i++ ){
      //   int ist = mc.vPoints[i].ista;
      //   for( int j=0; j<t.NHits; j++ ){
      //     if(t.vHits[j].ista == ist)
      //     {
      //       t.vHits[j].x = mc.vPoints[i].x;
      //       t.vHits[j].y = mc.vPoints[i].y;
      //       if(ist == 5) t.vHits[j].x += sqrt(vStations[ist].UInfo.sigma2[0]) * 20.f;
      //     }
      //   }
      // }
#endif
    
    iPoint = 0; // compare paraments at the first station
      // iPoint = NMCPoints-1;
    mc.MC_x = mc.vPoints[iPoint].x;
    mc.MC_y = mc.vPoints[iPoint].y;
    mc.MC_z = mc.vPoints[iPoint].z;
    mc.MC_px = mc.vPoints[iPoint].px;
    mc.MC_py = mc.vPoints[iPoint].py;
    mc.MC_pz = mc.vPoints[iPoint].pz;

    if( itr ==  TrackIndex[NMCTracks]) NMCTracks++;
  }
    // 	cout << NTracks <<" "<<NMCTracks<< " reco and Mc tracks have been read" << endl;
  FileTracks.close();
  FileMCTracks.close();
  
#ifdef OMP
  int threadNumberToCpuMap[80];
    /*    for (int i=0; i<8; i++){
          threadNumberToCpuMap[2*i+0] = i;
          threadNumberToCpuMap[2*i+1] = i+8;
          }*/
  for ( int iProc = 0; iProc < 4; iProc++ ) {
    for ( int i = 0; i < 8; i++ ) {
      threadNumberToCpuMap[2 * i + 0 + iProc * 20] = 4 * i + iProc;
      threadNumberToCpuMap[2 * i + 1 + iProc * 20] = 4 * i + 32 + iProc;
    }
    for ( int i = 0; i < 2; i++ ) {
      threadNumberToCpuMap[2 * i + 0 + 16 + iProc * 20] = 4 * i + iProc + 64;
      threadNumberToCpuMap[2 * i + 1 + 16 + iProc * 20] = 4 * i + 8 + iProc + 64;
    }
  }
#endif

#ifdef COPY
  NTracks = 0;
#pragma omp parallel reduction(+:NTracks) num_threads(tasks)
  {
#ifdef OMP
    int s;
    cpu_set_t cpuset;
    int cpuId = threadNumberToCpuMap[omp_get_thread_num()];
    pthread_t thread  = pthread_self();
    CPU_ZERO( &cpuset );
    CPU_SET( cpuId, &cpuset );
    s = pthread_setaffinity_np( thread, sizeof( cpu_set_t ), &cpuset );
    if ( s != 0 ) {
      cout << " pthread_setaffinity_np " << endl;
      handle_error_en( s, "pthread_setaffinity_np" );
    }
#endif

#pragma omp for
    for ( int j = 0; j < tasks; ++j ) {
      for ( int i = 0; i < NTracksPerThread; ++i ) {
        if ( j * NTracksPerThread + i > MaxNTracks ) {
          continue;
        }
        vTracks[j * NTracksPerThread + i] = vTracks[i];
        vMCTracks[j * NTracksPerThread + i] = vMCTracks[i];
#pragma omp atomic
        NTracks++;
      }
    }
  }
#endif // COPY

  NTracksV = NTracks/vecN;
  NTracks =  NTracksV*vecN;
}

void WriteOutput(){
  
  fstream Out, Diff;

  Out.open( (dataDir + "fit.dat").c_str(), ios::out );

#ifdef DAF
  Out << "Smoother" << endl;
#else
  Out << "Fitter" << endl;
#endif // DAF

  Out << MaxNTracks << endl;
  
  for( int it=0, itt=0; itt<NTracks; itt++ ){
    Track &t = vTracks[itt];
    MCTrack &mc = vMCTracks[itt];

    bool ok = 1;
    for( int i=0; i<6; i++ ){
      ok = ok && finite(t.T[i]);
    }
    for( int i=0; i<15; i++ ) ok = ok && finite(t.C[i]);
    
    if(!ok){ cout<<" infinite "<<endl; }

#ifdef DAF
    Out << it << endl;
    for( int iPoint=0; iPoint<NStations; iPoint++ ) {
      Out << iPoint << endl << "   "
          << " " << mc.vPoints[iPoint].x  << " " << mc.vPoints[iPoint].y  << " " << mc.vPoints[iPoint].z 
          << " " << mc.vPoints[iPoint].px << " " << mc.vPoints[iPoint].py << " " << mc.vPoints[iPoint].pz 
          << " " << mc.MC_q<<endl;
      Out << "   ";
      for( int i=0; i<6; i++ )
        Out<< " " <<t.Ts[iPoint][i];
      Out << endl;
      Out<<"   ";
      for( int i=0; i<15; i++ )
        Out<< " " <<t.Cs[iPoint][i];
      Out<<endl;
    }
    
#else // DAF
    const int iPoint = 0;
    Out <<it<<endl<<"   "
        << " " << mc.vPoints[iPoint].x  << " " << mc.vPoints[iPoint].y  << " " << mc.vPoints[iPoint].z 
        << " " << mc.vPoints[iPoint].px << " " << mc.vPoints[iPoint].py << " " << mc.vPoints[iPoint].pz 
        << " " << mc.MC_q<<endl;
    
    Out<<"   ";
    for( int i=0; i<6; i++ ) Out<< " " <<t.T[i];
    Out<<endl<<"   ";
    for( int i=0; i<15; i++ ) Out<< " " <<t.C[i];
    Out<<endl;
#endif // else DAF

    it++;
  }
  Out.close();
}

void FitTracksV(){

  double TimeTable[Ntimes];

  TrackV *TracksV = new TrackV[MaxNTracks/vecN+1];
  Fvec_t *Z0      = new Fvec_t[MaxNTracks/vecN+1]; // mc-z, used for result comparison
  Fvec_t *Z0s[MaxNStations];
  for( int is = 0; is < NStations; ++is )
    Z0s[is] = new Fvec_t[MaxNTracks/vecN+1]; 
  
#ifdef VC
  Fvec_t::Memory Z0mem;
  Fvec_t::Memory Z0smem[MaxNStations];
#endif
#ifndef MUTE
  cout<<"Prepare data..."<<endl;
#endif
  Stopwatch timer1;

  for( int iV=0; iV<NTracksV; iV++ ){ // loop on set of 4 tracks
#ifndef MUTE
    if( iV*vecN%100==0 ) cout<<iV*vecN<<endl;
#endif
    TrackV &t = TracksV[iV];
    for( int ist=0; ist<NStations; ist++ ){
      HitV &h = t.vHits[ist];
    
      h.x = 0.;
      h.y = 0.;
      h.w = 0.;
      h.H.X = 0.;
      h.H.Y = 0.;
      h.H.Z = 0.;
    }

#ifdef VC
    Fvec_t::Memory hxmem[NStations],hymem[NStations],hwmem[NStations];
    for( int it=0; it<vecN; it++ ){
      Track &ts = vTracks[iV*vecN+it];

      Z0mem[it] = vMCTracks[iV*vecN+it].MC_z;
      for( int is = 0; is < NStations; ++is )
        Z0smem[is][it] = vMCTracks[iV*vecN+it].vPoints[is].z;

      for( int ista=0, ih=0; ista<NStations; ista++ ){
        Hit &hs = ts.vHits[ih];
        if (hs.ista != ista) continue;
        ih++;

        hxmem[ista][it] = hs.x;
        hymem[ista][it] = hs.y;
        hwmem[ista][it] = 1.;
      }

    }
    for( int ista=0; ista<NStations; ista++ ){
      Fvec_t hxtemp( hxmem[ista] );
      Fvec_t hytemp( hymem[ista] );
      Fvec_t hwtemp( hwmem[ista] );
      t.vHits[ista].x = hxtemp;
      t.vHits[ista].y = hytemp;
      t.vHits[ista].w = hwtemp;
    }


    Fvec_t Z0temp(Z0mem);
    Z0[iV] = Z0temp;

    for( int is = 0; is < NStations; ++is ) {
      Fvec_t Z0stemp(Z0smem[is]);
      Z0s[is][iV] = Z0stemp;
    }
#else // VC
    
    for( int it=0; it<vecN; it++ ){
      Track &ts = vTracks[iV*vecN+it];
#ifdef X87
      Z0[iV] = vMCTracks[iV*vecN+it].MC_z;
      for( int is = 0; is < NStations; ++is )
        Z0s[is][iV] = vMCTracks[iV*vecN+it].vPoints[is].z;
#else
      Z0[iV][it] = vMCTracks[iV*vecN+it].MC_z;
      for( int is = 0; is < NStations; ++is )
        Z0s[is][iV][it] = vMCTracks[iV*vecN+it].vPoints[is].z;
#endif
      for( int ih=0; ih<ts.NHits; ih++ ){
        Hit &hs = ts.vHits[ih];
        HitV &h = t.vHits[hs.ista];
#ifdef X87
        h.x = hs.x;
        h.y = hs.y;
        h.w = 1.;
#else
        h.x[it] = hs.x;
        h.y[it] = hs.y;
        h.w[it] = 1.;
#endif
      }
    }
#endif // else VC


    if (0){    // output for check
      cout << "track " << iV << "  ";
      for( int ista=0; ista<NStations; ista++ )
        cout << t.vHits[ista].x << " ";
      cout << endl;
    }


    for( int ist=0; ist<NStations; ist++ ){
      HitV &h = t.vHits[ist];
      vStations[ist].Map.GetField(h.x, h.y, h.H);
    }
  }
  timer1.Stop();
#ifndef MUTE
  cout<<"Start fit..."<<endl;
#endif
  std::vector<double> fitSpeedCPU(tasks,0);
  std::vector<double> fitSpeedReal(tasks,0);
  std::vector<int> nFittedTracks(tasks,0);

  Stopwatch timerLocal;

  Stopwatch timer;
  Stopwatch timer2;
    //   Stopwatch timer_test;
  timer.Start();
  for( int times=0; times<Ntimes; times++){
    timer2.Start();
#ifdef TBB
#ifdef DEBUG_THREADS
    cout << NTracksV << " tracks by " << NTracksV / tasks << " tracks in group are run." << endl;
#endif // DEBUG_THREADS

      // FIXME: using affinity_partitioner only gives performance gain with #CPUs > 2
      // use affinity_partitioner
      //   static simple_partitioner ap;
      //   static affinity_partitioner ap;  // standart
      //   parallel_for(blocked_range<int>(0, NTracksV, NTracksV / tasks), ApplyFit(TracksV, vStations, NStations, NFits), ap);
      // use auto partitioning. Good then tasks != 2^N
      //   static auto_partitioner ap;
      //   parallel_for(blocked_range<int>(0, NTracksV), ApplyFit(TracksV, vStations, NStations, NFits), ap);
    {      // use special partitioning.  Good used with fixed thread-cpu correspondent by observer
        //     if (speedUpHT >= 0){
        //       parallel_for_simpleHT(NTracksV, tasks, ApplyFit(TracksV, vStations, NStations, NFits),speedUpHT);
        //     }
        //     else{
        //       parallel_for_simple(NTracksV, 10, tasks, ApplyFit(TracksV, vStations, NStations, NFits));
        //     };

      {   // fix nTask on thread
        const int nTracksVOnThread = NTracksPerThread/vecN;
        parallel_for_simpleEqNThr(nTracksVOnThread, tasks, ApplyFit(TracksV, vStations, NStations, NFits, &fitSpeedCPU, &nFittedTracks));
        NTracks = nTracksVOnThread*vecN*tasks; // for time counters
      }
    }
#else 
    int ifit;
    int iV;

    int nFittedTracks = 0;
#ifndef OMP
    tasks = 1;
#endif
  
#pragma omp parallel num_threads(tasks) private(iV,ifit) firstprivate(NTracksV,NStations, nFittedTracks,timerLocal)
    {
      timerLocal.Start();
#pragma omp for nowait //schedule(static) 
      for( iV=0; iV<NTracksV; iV++ ){ // loop on set of 4 tracks
        for( ifit=0; ifit<NFits; ifit++){
#ifdef DAF
#if 0
          fitter.Fit( TracksV[iV], vStations, NStations );
          for( int ip=0; ip<6; ip++ )
            TracksV[iV].Ts[0][ip] = TracksV[iV].T[ip];
          TracksV[iV].Cs[0] = TracksV[iV].C;
#elif 0
          fitter.Smooth( TracksV[iV], vStations, NStations, 1 );
#elif 0
          for( int ih = 0; ih < NStations; ++ih )
            fitter.SmoothHit( ih, TracksV[iV], vStations, NStations, 0 );
#else
          fitter.RunDAF( TracksV[iV], vStations, NStations );
#endif // 0
#else
          fitter.Fit( TracksV[iV], vStations, NStations );
#endif // DAF
        }
        nFittedTracks += vecN;
      }
      timerLocal.Stop();
      double cpuTimeLocal = timerLocal.CpuTime();
      double realTimeLocal = timerLocal.RealTime();
    
      if(fabs(cpuTimeLocal)<1.e-8) cpuTimeLocal = realTimeLocal;
    
#ifdef OMP
      int curThredNum = omp_get_thread_num();
#else
      int curThredNum = 0;
#endif
      if(cpuTimeLocal > 0)
        fitSpeedCPU[curThredNum] = nFittedTracks/cpuTimeLocal;
      else
        fitSpeedCPU[curThredNum] = 0;
    
      if(realTimeLocal > 0)
        fitSpeedReal[curThredNum] = nFittedTracks/realTimeLocal;
      else
        fitSpeedReal[curThredNum] = 0;
    }
#endif
    timer2.Stop();
    TimeTable[times]=timer2.RealTime();
  }
  timer.Stop();
        
#ifdef DAF
    // print results
  int NP[MaxNStations] = {0,0,0,0,0,0,0,0,0,0}; // TODO move out
  for(int iSt=0; iSt < NStations; iSt++) {
    for( int iV=0; iV<NTracksV; iV++ ){ // loop on set of 4 tracks
      for(int i=0; i < vecN; i++)
        if(TracksV[iV].w[iSt][i] < 0.5) NP[iSt]++;
    }
    cout << " ist  " << iSt << "  " << (double)NP[iSt]/(double)(NTracksV*vecN)*100. << endl;
  }
  
#endif
  
  for( int iV=0; iV<NTracksV; iV++ ){ // loop on set of 4 tracks
    TrackV &t = TracksV[iV];
    fitter.ExtrapolateALight( t.T, t.C, Z0[iV], TracksV[iV].T[4], t.f );
#ifdef DAF
    for( int is = 0; is < NStations; ++is )
      fitter.ExtrapolateALight( t.Ts[is], t.Cs[is], Z0s[is][iV], TracksV[iV].Ts[is][4], t.fs[is] );
#endif
  }
  
  double realtime=0;
  fstream TimeFile;
  TimeFile.open( (dataDir + "time.dat").c_str(), ios::out );
  for( int times=0; times<Ntimes; times++ ){
    TimeFile << TimeTable[times]*1.e6/(NTracks*NFits)<<endl;
    realtime += TimeTable[times]*1.e6/(NTracks*NFits);
  }  
  TimeFile.close();
  realtime /= Ntimes;

#ifndef MUTE
  cout<<"Preparation time/track = "<<timer1.CpuTime()*1.e6/NTracks/NFits<<" [us]"<<endl;
  cout<<"CPU  fit time/track = "<<timer.CpuTime()*1.e6/(NTracks*NFits)/Ntimes<<" [us]"<<endl;
  cout<<"Real fit time/track = "<<realtime <<" [us]"<<endl;
  cout<<"Total fit time = "<<timer.CpuTime()<<" [sec]"<<endl;
  cout<<"Total fit real time = "<<timer.RealTime()<<" [sec]"<<endl;
#else
#ifdef TOTALTIME
  
#ifdef TBB
  for(int iTask=0; iTask<tasks; iTask++)
  {
    fitSpeedCPU[iTask] = nFittedTracks[iTask]*vecN/fitSpeedCPU[iTask];
    fitSpeedReal[iTask] = fitSpeedCPU[iTask];
  }
#endif

  double realSpeedTotal = 0.;
  double cpuSpeedTotal = 0.;
  
  for(int iThread=0; iThread<tasks; iThread++)
  {
    realSpeedTotal += fitSpeedReal[iThread];
    cpuSpeedTotal += fitSpeedCPU[iThread];
  }
  
  cout<<"Prep[us], CPU fit/tr[us], Real fit/tr[us], CPU[sec], Real[sec] = "<<timer1.CpuTime()*1.e6/NTracks/NFits<<"\t";
  cout<< (1.e6/cpuSpeedTotal/NFits) <<"\t";
  cout<< (1.e6/realSpeedTotal/NFits) <<"\t";
    //   cout<<timer.CpuTime()*1.e6/(NTracks*NFits)/Ntimes<<"\t";
    //   cout<<realtime <<"\t";
  cout<<timer.CpuTime()<<"\t";
  cout<<timer.RealTime()<<endl;
#else //TOTALTIME
  cout<<"Prep[us], CPU fit/tr[us], Real fit/tr[us], CPU[sec], Real[sec] = "<<timer1.CpuTime()*1.e6/NTracks/NFits<<"\t";
  cout<<timer.CpuTime()*1.e6/(NTracks*NFits)/Ntimes<<"\t";
  cout<<realtime <<"\t";
  cout<<timer.CpuTime()<<"\t";
  cout<<timer.RealTime()<<endl;
#endif //OMP
#endif
  
  for( int iV=0; iV<NTracksV; iV++ ){ // loop on set of 4 tracks
    TrackV &t = TracksV[iV];
    
    for( int it=0; it < vecN; it++ ){
      Track &ts = vTracks[iV*vecN+it];
#if defined(CONVENTIONAL)
      Fvec_t *C = (Fvec_t*)&t.C;
#ifdef DAF
      Fvec_t *Cs[MaxNStations];
      for( int i = 0; i < NStations; ++i ) {
        Cs[i] = (Fvec_t*)&(t.Cs[i]);
      }
#endif // DAF
#else // CONVENTIONAL
      CovVConventional C_tmp = t.C;
      Fvec_t *C = (Fvec_t*)&C_tmp;

#ifdef DAF
      Fvec_t *Cs[MaxNStations];
      CovVConventional Cs_tmp[MaxNStations];
      for( int i = 0; i < NStations; ++i ) {
        Cs_tmp[i] = t.Cs[i];
        Cs[i] = (Fvec_t*)&(Cs_tmp[i]);
      }
#endif // DAF
#endif // else CONVENTIONAL

#ifdef X87
      for( int i=0; i<6; i++ )
        ts.T[i] = t.T[i];
      for( int i=0; i<15; i++ ) 
        ts.C[i] = C[i];
#ifdef DAF
      for( int is = 0; is < NStations; ++is ) {
        for( int i=0; i<6; i++ )
          ts.Ts[is][i] = t.Ts[is][i];
        for( int i=0; i<15; i++ ) 
          ts.Cs[is][i] = Cs[is][i];
      }
#endif // DAF
#else // else X87
      for( int i=0; i<6; i++ )
        ts.T[i] = t.T[i][it];
      for( int i=0; i<15; i++ )
        ts.C[i] = C[i][it];
      
#ifdef DAF
      for( int is = 0; is < NStations; ++is ) {
        for( int i=0; i<6; i++ )
          ts.Ts[is][i] = t.Ts[is][i][it];
        for( int i=0; i<15; i++ )
          ts.Cs[is][i] = Cs[is][i][it];
      }
#endif // DAF
#endif // X87
    }    
  }

  delete [] Z0;
  for( int is = 0; is < NStations; ++is )
    delete [] Z0s[is];
  delete [] TracksV;
}
 

int main(int argc, char *argv[]){

    //   time_t seconds;
    //   seconds = time (NULL);
    //   cout << "Begin NowIs(h:m:s) " << seconds/3600 << ":"  << (seconds/60)%60 << ":" << seconds%60 << endl;

    // RRMOD
  vStations = new Station[MaxNStations];

    // parse cmdline and set #tasks if TBB or OpenMP
  if(string(argv[0]).find("omp") != string::npos || string(argv[0]).find("tbb") != string::npos){
    if(argc != 2){
      cout << "Please specify #parallel tasks" << endl;
      return 1;
    }
    tasks = atoi( argv[1] );
  }

#ifdef TBB
  task_scheduler_init init(tasks);
  TMyObserver obs;
  obs.FInit();    // set cpu-threads correspondence
  obs.observe(useObserver); // Turn on observer. It's turn off by default.
#endif

  ReadInput();
  FitTracksV();
  WriteOutput();
#ifndef MUTE
  cout<<"track size = "<<sizeof(TrackV)<<endl;
  cout<<"char size = "<<sizeof(char)<<endl;
#endif

    // RRMOD
  delete[] vStations;

    //   seconds = time (NULL);
    //   cout << "End NowIs(h:m:s) " << seconds/3600 << ":"  << (seconds/60)%60 << ":" << seconds%60 << endl;
  return 0;
}