// $Id: AliHLTTPCCAPerformance.cxx,v 1.13 2010/09/01 10:38:27 ikulakov Exp $
// **************************************************************************
// This file is property of and copyright by the ALICE HLT Project          *
// ALICE Experiment at CERN, All rights reserved.                           *
//                                                                          *
// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
//                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
//                  for The ALICE HLT Project.                              *
//                                                                          *
// Developed by:   Igor Kulakov <I.Kulakov@gsi.de>                          *
//                 Maksym Zyzak <M.Zyzak@gsi.de>                            *
//                                                                          *
// Permission to use, copy, modify and distribute this software and its     *
// documentation strictly for non-commercial purposes is hereby granted     *
// without fee, provided that the above copyright notice appears in all     *
// copies and that both the copyright notice and this permission notice     *
// appear in the supporting documentation. The authors make no claims       *
// about the suitability of this software for any purpose. It is            *
// provided "as is" without express or implied warranty.                    *
//                                                                          *
//***************************************************************************
#ifdef DO_TPCCATRACKER_EFF_PERFORMANCE

#include "AliHLTTPCCAPerformance.h"

#include "AliHLTTPCCounters.h"
#include "AliHLTTPCPerformanceBase.h"
#include "AliHLTTPCCAGlobalPerformance.h"
#include "AliHLTTPCCAParam.h"

#include "AliHLTTPCCAGBHit.h"
#include "AliHLTTPCCAMCTrack.h"
#ifndef HLTCA_STANDALONE
#include "AliHLTTPCCAMCPoint.h"
#endif
#include "AliHLTTPCCAGBTrack.h"
#include "AliHLTTPCCAGBTracker.h"

#include "AliHLTTPCCADisplay.h"

#include "TRandom3.h" //dbg
#include "TMath.h"
#include "TROOT.h"
#include "Riostream.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TStyle.h"

AliHLTTPCCAPerformance::AliHLTTPCCAPerformance()
{
  static bool first_call = true;

  if (first_call){
    typedef TSubPerformance TSP;
    
      /// Just define here all sub-performances
      /// TSP(new __ClassName__               , __Name__      ),
#ifndef STAR_STANDALONE
    const int NSPerfo = 1;
#else // HLTCA_STANDALONE
    const int NSPerfo = 1;
#endif // HLTCA_STANDALONE
    const TSP perfos[NSPerfo] = {
      TSP(new AliHLTTPCCAGlobalPerformance, "Global Performance"),
    };

    subPerformances.resize(NSPerfo);
    for (int iP = 0; iP < NSPerfo; iP++){
      subPerformances[iP] = perfos[iP];
    }
  }
  first_call = false;
  //* constructor
}

AliHLTTPCCAPerformance::~AliHLTTPCCAPerformance()
{
  //* destructor
  fHitLabels.Resize( 0 );
  fMCTracks.Resize( 0 );
  for (unsigned int iPerf = 0; iPerf < subPerformances.size(); iPerf++){ // TODO: why we can't do this in subPerformances.~destructor() ??
    if (subPerformances[iPerf].perf) delete subPerformances[iPerf].perf;
  }
}

AliHLTTPCCAPerformance &AliHLTTPCCAPerformance::Instance()
{
  // reference to static object
  static AliHLTTPCCAPerformance gAliHLTTPCCAPerformance;
  return gAliHLTTPCCAPerformance;
}

bool AliHLTTPCCAPerformance::SetNewEvent(AliHLTTPCCAGBTracker* const tracker, string mcTracksFile, string mcPointsFile)
{
  fTracker = tracker;
    // Read MC info from file
  FILE *file = std::fopen( mcTracksFile.data(), "rb" );
  if ( !file ) {
    return false;
  }
  ReadMCEvent( file );
  fclose( file );
  
  file = std::fopen( mcPointsFile.data(), "rb" );
  if ( !file ) {
    return false;
  }
  ReadLocalMCPoints( file );
  fclose( file );

  return true;
} // void AliHLTTPCCAPerformance::SetNewEvent

void AliHLTTPCCAPerformance::InitSubPerformances()
{
    // Init subperformances
  static bool first_call = true;

  for (unsigned int iPerf = 0; iPerf < subPerformances.size(); iPerf++){
    subPerformances[iPerf]->SetNewEvent(fTracker, &fHitLabels, &fMCTracks, &fLocalMCPoints);
  }

  if (first_call) CreateHistos();
  
  first_call = false;
} // void AliHLTTPCCAPerformance::InitSubPerformances

void AliHLTTPCCAPerformance::CreateHistos()
{
  for (unsigned int iPerf = 0; iPerf < subPerformances.size(); iPerf++){
    if(!(subPerformances[iPerf]->IsHistoCreated())) 
      subPerformances[iPerf]->CreateHistos(subPerformances[iPerf].name, fOutputFile);
  }
}

bool AliHLTTPCCAPerformance::CreateHistos(string name)
{
  unsigned i = 0;
  for( ; i < (subPerformances.size()) && (subPerformances[i].name != name); i++);
  if(!(subPerformances[i]->IsHistoCreated()) && i != subPerformances.size()) 
    subPerformances[i]->CreateHistos(subPerformances[i].name, fOutputFile);
  if ( i == subPerformances.size() ) return 0;
  return 1;
}


void AliHLTTPCCAPerformance::WriteDir2Current( TObject *obj )
{
  //* recursive function to write down all created in the file histos
  if ( !obj->IsFolder() ) obj->Write();
  else {
    TDirectory *cur = gDirectory;
    ((TDirectory*)obj)->cd();
    TList *listSub = ( ( TDirectory* )obj )->GetList();
    TIter it( listSub );
    while ( TObject *obj1 = it() ) WriteDir2Current( obj1 );
    cur->cd();
  }
}


void AliHLTTPCCAPerformance::WriteHistos()
{
  if(fOutputFile) WriteDir2Current(fOutputFile);
}

void AliHLTTPCCAPerformance::ExecPerformance()
{
  fStatNEvents++;
  
  for (unsigned int iPerf = 0; iPerf < subPerformances.size(); iPerf++){
    if(subPerformances[iPerf].IsGlobalPerf) subPerformances[iPerf]->Exec(0);
  }
#if 1  // current event efficiencies
  for (unsigned int iPerf = 0; iPerf < subPerformances.size(); iPerf++){    
    cout << endl
        << " ---- " << subPerformances[iPerf].name << " event " << fStatNEvents << " ---- "<< endl;
    subPerformances[iPerf]->PrintEfficiency();
  }
  cout << endl << " ============================== " << endl;
#endif //0
  for (unsigned int iPerf = 0; iPerf < subPerformances.size(); iPerf++){  
    cout << endl
        << " ---- " << subPerformances[iPerf].name << " " << fStatNEvents << " events Statistic ---- " << endl;
    if(subPerformances[iPerf].IsGlobalPerf) subPerformances[iPerf]->PrintEfficiencyStatistic();
  };

}

AliHLTTPCPerformanceBase* AliHLTTPCCAPerformance::GetSubPerformance(string name)
{
  unsigned i = 0;
  for( ; (i < subPerformances.size()) && (subPerformances[i].name != name); i++);
  //ASSERT ( i != subPerformances.size() , " Incorrect name of subPerformance used.");
  if ( i == subPerformances.size() ) return 0;
  return subPerformances[i].perf;
} // AliHLTTPCPerformanceBase* AliHLTTPCCAPerformance::GetSubPerformance(string name)


  /// -------------------- Read\write MC information ---------------------
void AliHLTTPCCAPerformance::WriteMCEvent( FILE *file ) const
{
  // write MC information to the file
  int n = fMCTracks.Size();
  std::fwrite( &n, sizeof( int ), 1, file );
  n = fHitLabels.Size();
  std::fwrite( &n, sizeof( int ), 1, file );
  std::fwrite( fMCTracks.Data(), sizeof( AliHLTTPCCAMCTrack ), fMCTracks.Size(), file );
  std::fwrite( fHitLabels.Data(), sizeof( AliHLTTPCCAHitLabel ), fHitLabels.Size(), file );
}

void AliHLTTPCCAPerformance::ReadMCEvent( FILE *file )
{
  // read mc info from the file
  int read;
  int n;

  read = std::fread( &n, sizeof( int ), 1, file );
  assert( read == 1 );
  fMCTracks.Resize( n );
  read = std::fread( &n, sizeof( int ), 1, file );
  assert( read == 1 );
  fHitLabels.Resize( n );

  read = std::fread( fMCTracks.Data(), sizeof( AliHLTTPCCAMCTrack ), fMCTracks.Size(), file );
  assert( read == fMCTracks.Size() );

  read = std::fread( fHitLabels.Data(), sizeof( AliHLTTPCCAHitLabel ), fHitLabels.Size(), file );
  assert( read == fHitLabels.Size() );
  UNUSED_PARAM1(read);
//   for (int i = 0; i < fMCTracks.Size(); i++){
//     AliHLTTPCCAMCTrack& mc = fMCTracks[i];
//         std::cout << mc.PDG() << " ";
//         std::cout << std::endl;
//         for (int iPar = 0; iPar < 7; iPar++)
//           std::cout << mc.Par(iPar) << " ";
//         std::cout << std::endl;
//         for (int iPar = 0; iPar < 7; iPar++)
//           std::cout << mc.TPCPar(iPar) << " ";
//         std::cout << std::endl;
// //         std::cout << mc.P() << " ";
//   };
//   for (int i = 0; i < fHitLabels.Size(); i++){
//     AliHLTTPCCAHitLabel& l = fHitLabels[i];
//     for (int iPar = 0; iPar < 3; iPar++)
//      std::cout << l.fLab[iPar] << " ";
//     std::cout << std::endl;
//   };
}

void AliHLTTPCCAPerformance::ReadLocalMCPoints( FILE *file )
{
  
  int read;
  int n;

  read = std::fread( &n, sizeof( int ), 1, file );
  assert( read == 1 );
  fLocalMCPoints.Resize( n );

  read = std::fread( fLocalMCPoints.Data(), sizeof( AliHLTTPCCALocalMCPoint ), fLocalMCPoints.Size(), file );
  assert( read == fLocalMCPoints.Size() );
  UNUSED_PARAM1(read);
} // void AliHLTTPCCAPerformance::ReadLocalMCPoints( FILE *file )


void AliHLTTPCCAPerformance::SetMCTracks(vector<AliHLTTPCCAMCTrack>& mcTracks)
{
  const int N = mcTracks.size();
  fMCTracks.Resize(N);
  for(int i = 0; i < N; i++){
    fMCTracks[i] = mcTracks[i];
  }
}

void AliHLTTPCCAPerformance::SetMCPoints(vector<AliHLTTPCCALocalMCPoint>& mcPoints)
{
  const int N = mcPoints.size();
  fLocalMCPoints.Resize(N);
  for(int i = 0; i < N; i++){
    fLocalMCPoints[i] = mcPoints[i];
  }
}

void AliHLTTPCCAPerformance::SetHitLabels(vector<AliHLTTPCCAHitLabel>& hitLabels)
{
  const int N = hitLabels.size();
  fHitLabels.Resize(N);
  for(int i = 0; i < N; i++){
    fHitLabels[i] = hitLabels[i];
  }					   
}

void AliHLTTPCCAPerformance::SaveDataInFiles(string prefix) const
{

  {
    ofstream ofile((prefix+"hitLabels.data").data(),ios::out|ios::app);
    const int Size = fHitLabels.Size();
    ofile << Size << std::endl;
    for (int i = 0; i < fHitLabels.Size(); i++){
      const AliHLTTPCCAHitLabel &l = fHitLabels[i];
      ofile << l;
    }
    ofile.close();
  }

  {
    ofstream ofile((prefix+"MCTracks.data").data(),ios::out|ios::app);
    const int Size = fMCTracks.Size();
    ofile << Size << std::endl;
    for (int i = 0; i < fMCTracks.Size(); i++){
      const AliHLTTPCCAMCTrack &l = fMCTracks[i];
      ofile << l;
    }
    ofile.close();
  }

  {
    ofstream ofile((prefix+"MCPoints.data").data(),ios::out|ios::app);
    const int Size = fLocalMCPoints.Size();
    ofile << Size << std::endl;
    for (int i = 0; i < fLocalMCPoints.Size(); i++){
      const AliHLTTPCCALocalMCPoint &l = fLocalMCPoints[i];
      ofile << l;
    }
    ofile.close();
  }
}

struct CbmL1MCPoint
{
  double x, y, z, px, py, pz;
  double xIn,  yIn,  zIn,  pxIn,  pyIn,  pzIn;
  double xOut, yOut, zOut, pxOut, pyOut, pzOut;
  double p, q, mass;
  int pdg, ID, mother_ID;
  int iStation;
};

class CbmL1MCTrack
{
 public:
  double mass, q, p, x, y, z, px, py, pz;
  int ID, mother_ID, pdg;

  int nHits, nMCPoints;
  
  int nMCContStations;   // number of consecutive stations with mcPoints
  int nHitContStations;  // number of consecutive stations with hits
  int maxNStaMC;         // max number of mcPoints on station
  int maxNSensorMC;         // max number of mcPoints with same z
  int maxNStaHits;       // max number of hits on station
   
  int nStations;         // number of stations with hits

  bool isReconstructable;
};

struct TMcPointsSortHelper{
  int oldI, newI, trackId, z;
  static bool compare1( const TMcPointsSortHelper& a, const TMcPointsSortHelper& b ) {
    return (a.trackId < b.trackId) || ( (a.trackId == b.trackId) && (a.z < b.z) ); }
  static bool compare2( const TMcPointsSortHelper& a, const TMcPointsSortHelper& b ) { return a.oldI < b.oldI; }
};

bool AliHLTTPCCAPerformance::ReadDataFromFiles(ifstream& ifile)
{
  const int iVerbose = 0;
  static int nEvent_ = 1;
  
  {
    vector<CbmL1MCPoint> vMCPoints;
    vector<CbmL1MCTrack> vMCTracks;
    vector<int> vHitMCRef;
       
    vMCPoints.clear();
    vMCTracks.clear();
    vHitMCRef.clear();


      // check if it is right position in file
    char s[] = "EVENT:     ";  // buffer
    int nEv=0;                // event number
    ifile >> s;
    ifile >> nEv;

    if (nEv != nEvent_)  cout << "-E- CbmL1: Performance: can't read event number " << nEvent_ << " from file " << "data_perfo.txt" << endl;
      // vMCPoints
    int n;  // number of elements
    ifile >> n;
    for (int i = 0; i < n; i++){
      CbmL1MCPoint element;

      ifile >> element.xIn;
      ifile >> element.yIn;
      ifile >> element.zIn;
      ifile >> element.pxIn;
      ifile >> element.pyIn;
      ifile >> element.pzIn;
            
      ifile >> element.xOut;
      ifile >> element.yOut;
      ifile >> element.zOut;
      ifile >> element.pxOut;
      ifile >> element.pyOut;
      ifile >> element.pzOut;

      element.x = (element.xIn + element.xOut)/2;
      element.y = (element.yIn + element.yOut)/2;
      element.z = (element.zIn + element.zOut)/2;
      element.px = (element.pxIn + element.pxOut)/2;
      element.py = (element.pyIn + element.pyOut)/2;
      element.pz = (element.pzIn + element.pzOut)/2;
      
      ifile >> element.p;
      ifile >> element.q;
      ifile >> element.mass;
      
      ifile >> element.pdg;
      ifile >> element.ID;
      ifile >> element.mother_ID;
      ifile >> element.iStation;

      int nhits;
      ifile >> nhits;
      for (int k = 0; k < nhits; k++){
        int helement;
        ifile >> helement;
        // element.hitIds.push_back(helement);
      };
      
      vMCPoints.push_back(element);
    };
    if (iVerbose >= 4) cout << "vMCPoints[" << n << "]" << " have been read." << endl;

      // vMCTracks . without Points
    ifile >> n;
    for (int i = 0; i < n; i++){
      CbmL1MCTrack element;

      ifile >> element.x;
      ifile >> element.y;
      ifile >> element.z;
      ifile >> element.px;
      ifile >> element.py;
      ifile >> element.pz;
      ifile >> element.p;
      ifile >> element.q;
      ifile >> element.mass;

      ifile >> element.pdg;
      ifile >> element.ID;
      ifile >> element.mother_ID;

      ifile >> element.nHits;
      for (int k = 0; k < element.nHits; k++){
        int helement;
        ifile >> helement;
        // element.StsHits.push_back(helement);
      };
      ifile >> element.nMCPoints;
      for (int k = 0; k < element.nMCPoints; k++){
        int helement;
        ifile >> helement;
        // element.Points.push_back(helement);
      };

      ifile >> element.nMCContStations;
      ifile >> element.nHitContStations;
      ifile >> element.maxNStaMC;
      ifile >> element.maxNSensorMC;
      ifile >> element.maxNStaHits;
      ifile >> element.nStations;

      // element.CalculateIsReconstructable();
      vMCTracks.push_back(element);
    };
    if (iVerbose >= 4) cout << "vMCTracks[" << n << "]" << " have been read." << endl;
    
        // vHitMCRef
    ifile >> n;
    for (int i = 0; i < n; i++){
      int element;
      ifile >> element;
      vHitMCRef.push_back(element);
    };
    if (iVerbose >= 4) cout << "vHitMCRef[" << n << "]" << " have been read." << endl;

      // vHitStore
    ifile >> n;
    int tmpI;
    double tmpD;
    for (int i = 0; i < n; i++){
      // CbmL1HitStore element;
      ifile >> tmpI;//element.ExtIndex;
      ifile >> tmpI;//element.iStation;

      ifile >> tmpD;//element.x;
      ifile >> tmpD;//element.y;
      // vHitStore.push_back(element);
//  cout <<   element.ExtIndex << " "<<   element.iStation << " "<<   element.x << " "<<  element.y << " " <<  endl;
    };
    if (iVerbose >= 4) cout << "vHitStore[" << n << "]" << " have been read." << endl;
    
          // vStsHits
    ifile >> n;

    for (int i = 0; i < n; i++){
      // CbmL1StsHit element;
      ifile >> tmpI;//element.hitId;
      ifile >> tmpI;//element.extIndex;

      int nPoints;
      ifile >> nPoints;
      // element.mcPointIds.clear();
      for (int k = 0; k < nPoints; k++){
        int id;
        ifile >> id;
        // element.mcPointIds.push_back(id);
      };
      // vStsHits.push_back(element);
    };
    if (iVerbose >= 4) cout << "vStsHits[" << n << "]" << " have been read." << endl;


//    if (nEvent == maxNEvent) { // file open on begin of all work class and close at end
//       ifile.close();
//       cout << " -I- Performance: data read from file " << "data_perfo.txt" << " successfully"<< endl;
//     }
    if (iVerbose >= 2) cout << "-I- CbmL1: L1Performance data for event " << nEvent_ << " has been read  successfully." << endl;



    std::map< int, int > mapTrackIdToIndex;
      
    fMCTracks.Resize(vMCTracks.size());
    for (unsigned int i = 0; i < vMCTracks.size(); i++){
      AliHLTTPCCAMCTrack &t = fMCTracks[i];
      CbmL1MCTrack &ti = vMCTracks[i];
      t.SetId(ti.ID);
      mapTrackIdToIndex[ti.ID] = i;
      t.SetMotherId(ti.mother_ID);
      t.SetPDG(ti.pdg);
      t.SetP(ti.p);
      t.SetPt(ti.p); // TODO tmp
      t.SetPar( 0, ti.x );
      t.SetPar( 1, ti.y );
      t.SetPar( 2, ti.z );
      t.SetPar( 6, ti.q/ti.p );
      t.SetPar( 3, ti.px*t.Par(6) );
      t.SetPar( 4, ti.py*t.Par(6));
      t.SetPar( 5, ti.pz*t.Par(6) );

      t.SetNHits( ti.nHits );
      t.SetNMCPoints( ti.nMCPoints );

      t.SetNMCRows( ti.nStations );
    }

      // save mc points in TrackId-order
    fLocalMCPoints.Resize(vMCPoints.size());
    
    vector<TMcPointsSortHelper> sortHelper(vMCPoints.size());
    for (unsigned int i = 0; i < vMCPoints.size(); i++){
      CbmL1MCPoint &pi = vMCPoints[i];

      sortHelper[i].oldI = i;
      sortHelper[i].trackId = pi.ID;
      sortHelper[i].z = pi.z;
    }

    sort(sortHelper.begin(), sortHelper.end(), TMcPointsSortHelper::compare1); // sort by trackId
    
    for( unsigned int i = 0; i < sortHelper.size(); ++i ) {
      sortHelper[i].newI = i; // remember order
    }
    sort(sortHelper.begin(), sortHelper.end(), TMcPointsSortHelper::compare2); // sort back

    for (unsigned int i = 0; i < vMCPoints.size(); i++){
      AliHLTTPCCALocalMCPoint &p = fLocalMCPoints[sortHelper[i].newI];
      CbmL1MCPoint &pi = vMCPoints[i]; // take points according to trackId-order
      
      p.SetX( pi.x );
      p.SetY( pi.y );
      p.SetZ( pi.z );
      p.SetPx( pi.px );
      p.SetPy( pi.py);
      p.SetPz( pi.pz );
      p.SetQP( pi.q/pi.p );
      p.SetIRow( pi.iStation );
      int trackId = pi.ID;
      int trackI = mapTrackIdToIndex[trackId];
      p.SetTrackI( trackI );
      p.SetTrackID( trackId );
    }

      // find FirstMCPointID
    int oldTI = -1;
    for (unsigned int i = 0; i < vMCPoints.size(); i++){
      AliHLTTPCCALocalMCPoint &p = fLocalMCPoints[i];
      const int ti = p.TrackI();
      if ( ti == oldTI ) continue;
      fMCTracks[ti].SetFirstMCPointID( i );
      oldTI = ti;
    }
    
    fHitLabels.Resize(vHitMCRef.size());

    for (unsigned int i = 0; i < vHitMCRef.size(); i++){
      AliHLTTPCCAHitLabel &l = fHitLabels[i];
      l.fLab[0] = -1;
      l.fLab[1] = -1;
      l.fLab[2] = -1;
      int li =  vHitMCRef[i];
      if (li < 0) continue;
      l.fLab[0] = mapTrackIdToIndex[vMCPoints[li].ID];
    }
  
  }
  nEvent_++;
  

  //   // calculate needed additional info
  // {
  //   for (int i = 0; i < fMCTracks.Size(); ++i ){
  //     AliHLTTPCCAMCTrack& t = fMCTracks[i];

  //     vector<int> rows;
  //     rows.resize( AliHLTTPCCAParameters::NStations, 0 );
  //     for (int j = 0, iP = t.FirstMCPointID(); j < t.NMCPoints(); iP++, j++){
  //       rows[fLocalMCPoints[iP].IRow()]++;
  //     }

  //     int nRows = 0;
  //     for (unsigned int j = 0 ; j < rows.size(); j++)
  //       if( rows[j] > 0 ) nRows++;

  //     t.SetNMCRows( nRows );
  //   }
  // }


  return 1;
}

#include "AliHLTTPCCADisplay.h"
  // Use smeared MCposition instead of hits
void AliHLTTPCCAPerformance::ShiftHitsToMC()
{ // TODO
  AliHLTResizableArray<AliHLTTPCCAGBHit>& hits = const_cast<AliHLTTPCCAGBTracker *>(fTracker)->fHits;

  std::sort( hits.Data(), hits.Data() + const_cast<AliHLTTPCCAGBTracker *>(fTracker)->NHits(), AliHLTTPCCAGBHit::Compare ); // has an influence on the procedure. CHECKME why?
  
  TRandom3 rand;
  {
    int nHits = fTracker->NHits();

    for ( int ih = 0; ih < nHits; ih++ ) {
      AliHLTTPCCAGBHit &hit = hits[ih];
      const AliHLTTPCCAHitLabel &l = (fHitLabels)[hit.ID()];

      const int iMC = l.fLab[0];
      if( iMC<0 ) continue;
      AliHLTTPCCAMCTrack &mc = (fMCTracks)[iMC];


      int MCindex = -1;
      int nFirstMC = mc.FirstMCPointID();
      int nMCPoints = mc.NMCPoints();
      AliHLTTPCCALocalMCPoint *points = &((fLocalMCPoints).Data()[nFirstMC]);

      for(int iMCPoint=0; iMCPoint<nMCPoints; iMCPoint++)
      {
        if(points[iMCPoint].IRow() == hit.IRow())
        {
          if(fabs(hit.Y() - points[iMCPoint].Y())<2 && fabs(hit.Z() - points[iMCPoint].Z())<2) // dbg +
            MCindex = iMCPoint;
        }
      }
      if(MCindex == -1)
      {
        continue;
      }

      const AliHLTTPCCALocalMCPoint& point = points[MCindex];

// #ifdef DRAW
//       AliHLTTPCCADisplay::Instance().SetTPCView();
//       AliHLTTPCCADisplay::Instance().DrawTPC();
//       for(int iMCPoint=0; iMCPoint<nMCPoints; iMCPoint++)
//       {
//         AliHLTTPCCADisplay::Instance().DrawGBPoint((float)points[iMCPoint].X(),
//                                                    (float)points[iMCPoint].Y(),
//                                                    (float)points[iMCPoint].Z(), 1, (Size_t)1);
//       }
//       AliHLTTPCCADisplay::Instance().DrawGBPoint((float)hit.X(), 
//                                                  (float)hit.Y(),
//                                                  (float)hit.Z(),
//                                                  -1, (Size_t)0.5);
//       AliHLTTPCCADisplay::Instance().DrawGBPoint((float)point.X(), 
//                                                  (float)point.Y(),
//                                                  (float)point.Z(),
//                                                   2, (Size_t)0.5);
//       std::cout << hit.IRow() << "    " << hit.X() << " " << hit.Y() << " " << hit.Z() << std::endl;
//       std::cout << point.IRow() << "    " << point.X() << " " << point.Y() << " " << point.Z() << std::endl;
// 
//       AliHLTTPCCADisplay::Instance().Ask();
// #endif

      sfloat_v Err2Y = 0.01, Err2Z = 0.01;

      double mcX =  point.X();
      double mcY =  point.Y();
      double mcZ =  point.Z();

      double mcXTmp  =  mcX*cos( hit.Angle() ) +  mcY*sin( hit.Angle() );
      double mcYTmp  = -mcX*sin( hit.Angle() ) +  mcY*cos( hit.Angle() );
      mcYTmp -= sqrt(Err2Y[0])*rand.Gaus();
      mcZ -= sqrt(Err2Z[0])*rand.Gaus();
      mcX = mcXTmp*cos( hit.Angle() ) -  mcYTmp*sin( hit.Angle() );
      mcY = mcXTmp*sin( hit.Angle() ) +  mcYTmp*cos( hit.Angle() );
      // fix hits
      hit.SetX( mcX );
      hit.SetY( mcY );
      hit.SetZ( mcZ );

      hit.SetErr2Y( Err2Y[0] );
      hit.SetErr2Y( Err2Z[0] );
    }
  }
}

  // Create new hits spreading MCPositions
void AliHLTTPCCAPerformance::ResimulateHits()
{ // TODO
  // AliHLTResizableArray<AliHLTTPCCAGBHit>& hits = const_cast<AliHLTTPCCAGBTracker *>(fTracker)->fHits;

  // const int NHits = fLocalMCPoints.Size();
  // const_cast<AliHLTTPCCAGBTracker *>(fTracker)->fNHits = NHits;
  // // hits.Clear();
  // hits.Resize(NHits);
  // fHitLabels.Resize(NHits);

  
  // TRandom3 rand;
  // for ( int ih = 0; ih < NHits; ih++ ) {
  //   AliHLTTPCCALocalMCPoint &mcP = fLocalMCPoints[ih];

  //   AliHLTTPCCAGBHit &hit = hits[ih];
  //   hit.SetID(ih);
  //   AliHLTTPCCAHitLabel &l = fHitLabels[ih];
  //   l.fLab[0] = mcP.TrackI();
  //   l.fLab[1] = -1;
  //   l.fLab[2] = -1;


  //     // get errors
  //   double mcEx = mcP.Px();
  //   double mcEy = mcP.Py();
  //   double mcEz = mcP.Pz();
  //   double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );
      
  //   sfloat_v Err2Y = 0.12*0.12, Err2Z = 0.16*0.16;
  //     //float Err2Y = 0, Err2Z = 0;
  //   const AliHLTTPCCAParam par;
  //   TrackParamVector t;
  //   t.SetZ( mcP.Z() );
  //   t.SetSinPhi(mcEy / mcEt);
  //   t.SetDzDs(mcEz / mcEt);
      
  //     // fix hits
  //   hit.SetZ( mcP.Z() - sqrt(Err2Z[0])*rand.Gaus() );
  //   hit.SetY( mcP.Y() - sqrt(Err2Y[0])*rand.Gaus() );
  //   hit.SetX( mcP.X() );

  //   hit.SetIRow( mcP.IRow() );
  //   hit.SetErr2Y( 0.12*0.12 );
  //   hit.SetErr2Z( 0.16*0.16 );
  // }
}

  // Match hits with closest MCPoint
void AliHLTTPCCAPerformance::RematchHits(){
  AliHLTResizableArray<AliHLTTPCCAGBHit>& hits = const_cast<AliHLTTPCCAGBTracker *>(fTracker)->fHits;

  const int NHits = fTracker->NHits();

  for ( int ih = 0; ih < NHits; ih++ ) {
    AliHLTTPCCAGBHit &hit = hits[ih];
    AliHLTTPCCAHitLabel &l = fHitLabels[ih];

    int MCindex = -1;
    float R2 = 10e30;
    for(int iMCPoint=0; iMCPoint < fLocalMCPoints.Size(); iMCPoint++)
    {
      if(fLocalMCPoints[iMCPoint].IRow()   != hit.IRow()) continue;

      const float dY = hit.Y() - fLocalMCPoints[iMCPoint].Y();
      const float dZ = hit.Z() - fLocalMCPoints[iMCPoint].Z();
      const float r2 = dY*dY + dZ*dZ;
      
      if( r2 < R2 ) {
        MCindex = iMCPoint;
        R2 = r2;
      }
    }
    if(MCindex == -1)
    {
      l.fLab[0] = -1;
      l.fLab[1] = -1;
      l.fLab[2] = -1;
      continue;
    }

    
    l.fLab[0] = fLocalMCPoints[MCindex].TrackI();
    l.fLab[1] = -1;
    l.fLab[2] = -1;
  }
}


#endif //DO_TPCCATRACKER_EFF_PERFORMANCE