// $Id: AliHLTTPCCAGlobalPerformance.cxx,v 1.11 2010/08/18 20:46:09 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 "AliHLTTPCCounters.h"

#include "AliHLTTPCCATrackPerformanceBase.h"
#include "AliHLTTPCCAGlobalPerformance.h"


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

#include "AliHLTTPCCADisplay.h"

#ifdef DRAW
#define DRAW_GLOBALPERF
#endif

#ifdef DRAW_GLOBALPERF
#include "AliHLTTPCCADisplay.h"
#include "AliHLTTPCCAPerformance.h"
#endif

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

void AliHLTTPCCAGlobalPerformance::SetNewEvent(const AliHLTTPCCAGBTracker * const tracker,
                            AliHLTResizableArray<AliHLTTPCCAHitLabel> *hitLabels,
                            AliHLTResizableArray<AliHLTTPCCAMCTrack> *mcTracks,
                            AliHLTResizableArray<AliHLTTPCCALocalMCPoint> *localMCPoints)
{
  AliHLTTPCCATrackPerformanceBase::SetNewEvent(tracker, hitLabels, mcTracks, localMCPoints);


  nRecoTracks = fTracker->NTracks();


} // void AliHLTTPCCAGlobalPerformance::SetNewEvent

void AliHLTTPCCAGlobalPerformance::CheckMCTracks()
{

  for ( int imc = 0; imc < nMCTracks; imc++ ) (*fMCTracks)[imc].SetNHits( 0 );

  for ( int ih = 0; ih < (*fHitLabels).Size(); ih++ ) { // TODO: do we need to calculate consequtive hits??
    const AliHLTTPCCAHitLabel &l = (*fHitLabels)[ih];
    if ( l.fLab[0] >= 0 ) (*fMCTracks)[l.fLab[0]].SetNHits( (*fMCTracks)[l.fLab[0]].NHits() + 1 );
    if ( l.fLab[1] >= 0 ) (*fMCTracks)[l.fLab[1]].SetNHits( (*fMCTracks)[l.fLab[1]].NHits() + 1 );
    if ( l.fLab[2] >= 0 ) (*fMCTracks)[l.fLab[2]].SetNHits( (*fMCTracks)[l.fLab[2]].NHits() + 1 );
  }

  mcData.resize(nMCTracks);
  for ( int imc = 0; imc < nMCTracks; imc++ ) {
    AliHLTTPCCAMCTrack &mc = (*fMCTracks)[imc];
    AliHLTTPCCAPerformanceMCTrackData &mcTrackData = mcData[imc];
    mc.SetSet( 0 );
    mc.SetNReconstructed( 0 );
    mc.SetNTurns( 1 );
    if ( mc.NHitRows() >= PParameters::MinimumHitsForMCTrack && mc.NHits() == mc.NHitRows() )
//    if ( mc.NMCPoints() >= PParameters::MinimumMCPointsForMCTrack )
    if ( mc.NMCRows() >= PParameters::MinimumMCPointsForMCTrack )
      mcTrackData.SetAsReconstructable();

     // sets of tracks 0-OutSet, 1-ExtraSet, 2-RefSet, 3-ExtraSecSet, 4-ExtraPrimSet, 5-RefSecSet, 6-RefPrimSet, 7-LongRefPrimSet
    if ( mc.P() >= PParameters::ExtraThreshold ) {
      if ( mc.P() >= PParameters::RefThreshold ) {
        mc.SetSet( 2 );
        mcTrackData.SetSet( 2 );
        if ( mc.MotherId() != -1 ) {
          mc.SetSet( 5 );
          mcTrackData.SetSet( 5 );
          
        }
        else {
          mc.SetSet( 6 );
          mcTrackData.SetSet( 6 );
          if ( mc.NHits() == fTracker->NStations() ) {
            mc.SetSet( 7 );
            mcTrackData.SetSet( 7 );
          }
        }
      }
      else{
        mc.SetSet( 1 );
        mcTrackData.SetSet( 1 );
        if ( mc.MotherId() != -1 ) {
          mc.SetSet( 3 );
          mcTrackData.SetSet( 3 );
          
        }
        else {
          mc.SetSet( 4 );
          mcTrackData.SetSet( 4 );
        }
      }
    }
  } // for iMC

} // void AliHLTTPCCAGlobalPerformance::CheckMCTracks()


void AliHLTTPCCAGlobalPerformance::MatchTracks()
{
  recoData.resize(nRecoTracks);
  for ( int itr = 0; itr < nRecoTracks; itr++ ) {
    int traLabels = -1;
    double traPurity = 0;
    const AliHLTTPCCAGBTrack &tCA = fTracker->Track( itr );
    const int nhits = tCA.NHits();
    int *lb = new int[nhits*3];
    int nla = 0;
    for ( int ihit = 0; ihit < nhits; ihit++ ) {
      const int index = fTracker->TrackHit( tCA.FirstHitRef() + ihit );
      const AliHLTTPCCAHitLabel &l = (*fHitLabels)[fTracker->Hit( index ).ID()];
      if ( l.fLab[0] >= 0 ) lb[nla++] = l.fLab[0];
      if ( l.fLab[1] >= 0 ) lb[nla++] = l.fLab[1];
      if ( l.fLab[2] >= 0 ) lb[nla++] = l.fLab[2];
    }
    sort( lb, lb + nla );
    int labmax = -1, labcur = -1, lmax = 0, lcurr = 0;
    for ( int i = 0; i < nla; i++ ) {
      if ( lb[i] != labcur ) {
        if ( labcur >= 0 && lmax < lcurr ) {
          lmax = lcurr;
          labmax = labcur;
        }
        labcur = lb[i];
        lcurr = 0;
      }
      lcurr++;
    }
    if ( labcur >= 0 && lmax < lcurr ) {
      lmax = lcurr;
      labmax = labcur;
    }
    lmax = 0;
    for ( int ihit = 0; ihit < nhits; ihit++ ) {
      const int index = fTracker->TrackHit( tCA.FirstHitRef() + ihit );
      const AliHLTTPCCAHitLabel &l = (*fHitLabels)[fTracker->Hit( index ).ID()];
      if ( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax
         ) lmax++;
    }
    traLabels = labmax;
    traPurity = ( ( nhits > 0 ) ? double( lmax ) / double( nhits ) : 0 );
    if ( lb ) delete[] lb;

    recoData[itr].SetMCTrack(traLabels, traPurity, nhits);
    if ( recoData[itr].IsReco(PParameters::MinTrackPurity, PParameters::MinimumHitsForRecoTrack) ) mcData[traLabels].AddReconstructed();
  } // for iReco
} // void AliHLTTPCCAGlobalPerformance::MatchTracks()


void AliHLTTPCCAGlobalPerformance::EfficiencyPerformance( )
{
  for ( int iRTr = 0; iRTr < nRecoTracks; iRTr++ ) {
    if (  recoData[iRTr].IsGhost(PParameters::MinTrackPurity) )
      fEff.ghosts++;
  }

  for ( int iMCTr = 0; iMCTr < nMCTracks; iMCTr++ ) {
    AliHLTTPCCAPerformanceMCTrackData &mc = mcData[iMCTr];
    if ( !mc.IsReconstructable() ) continue;
    const bool reco = mc.IsReconstructed();
    const int clones = mc.GetNClones();

    if ( mc.GetSet() == 0){ // rest, out track
      fEff.Inc(reco,clones,"rest");
    }
    else{ // good
      fEff.Inc(reco,clones,"total");
      switch( mc.GetSet() ){
        case 1:
          fEff.Inc(reco,clones,"extra");
          break;
        case 2:
          fEff.Inc(reco,clones,"ref");
          break;
        case 3:
          fEff.Inc(reco,clones,"extra");
          fEff.Inc(reco,clones,"extra_sec");
          break;
        case 4:
          fEff.Inc(reco,clones,"extra");
          fEff.Inc(reco,clones,"extra_prim");
          break;
        case 5:
          fEff.Inc(reco,clones,"ref");
          fEff.Inc(reco,clones,"ref_sec");
          break;
        case 6:
          fEff.Inc(reco,clones,"ref");
          fEff.Inc(reco,clones,"ref_prim");
          break;
            break;
        case 7:
          fEff.Inc(reco,clones,"ref");
          fEff.Inc(reco,clones,"ref_prim");
          fEff.Inc(reco,clones,"ref_prim_long");
          break;
      }
    }
  } // for iMCTr

  AliHLTTPCCATrackPerformanceBase::EfficiencyPerformance();
} // void AliHLTTPCCAGlobalPerformance::EfficiencyPerformance( )

void AliHLTTPCCAGlobalPerformance::FillHistos()
{
  AliHLTTPCCATrackPerformanceBase::FillHistos();


    //

  const int NMCTracks = (*fMCTracks).Size();
  vector<int> mcTrackNRecoHits;
  vector<int> nHitsVsRow;
  vector<int> nMCPointsVsRow;
  const int Multiplicity = (*fMCTracks).Size();
  
  mcTrackNRecoHits.resize(NMCTracks, 0);
  nHitsVsRow.resize(fTracker->NStations());
  nMCPointsVsRow.resize(fTracker->NStations());
  for(int iH=0; iH < fTracker->NHits(); iH++){
    const AliHLTTPCCAGBHit &hit = fTracker->Hit( iH );
    
    nHitsVsRow[hit.IRow()]++;
    
    const AliHLTTPCCAHitLabel &l = (*fHitLabels)[hit.ID()];
    if ( l.fLab[0] >= 0 ) mcTrackNRecoHits[l.fLab[0]]++;
    if ( l.fLab[1] >= 0 ) mcTrackNRecoHits[l.fLab[1]]++;
    if ( l.fLab[2] >= 0 ) mcTrackNRecoHits[l.fLab[2]]++;
  }
  for(int i=0; i < NMCTracks; i++){
    AliHLTTPCCAMCTrack &mcT = (*fMCTracks)[i];
        
    GetHisto("mcTrackNRecoHits")->Fill( mcTrackNRecoHits[i] );
    GetHisto("nMCPointsVsMCMom")->Fill( mcT.P(), mcT.NMCPoints() );
    
    if ( mcT.NMCPoints() > 0 ) {
      double mcEx = mcT.Px();
      double mcEy = mcT.Py();
      double mcEz = mcT.Pz();
      double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );

      const double dZdS = mcEz/mcEt;

      GetHisto("nHitsOverNMCPointsVsMCMom")->Fill( mcT.P(), float(mcTrackNRecoHits[i])/float(mcT.NMCPoints()) );
      GetHisto("nHitsOverNMCPointsVsMCDzDs")->Fill( dZdS, float(mcTrackNRecoHits[i])/float(mcT.NMCPoints()) );
    }
  }

  for(int i=0; i < (*fLocalMCPoints).Size(); i++){
    nMCPointsVsRow[(*fLocalMCPoints)[i].IRow()]++;
  }
  for(int i=0; i < fTracker->NStations(); i++){
    if ( nMCPointsVsRow[i] > 0 ) {
      GetHisto("nHitsOverNMCPointsVsRow")->Fill( i, float(nHitsVsRow[i])/float(nMCPointsVsRow[i]) );
      GetHisto("nHitsOverNMCPointsVsNMCTracks")->Fill( Multiplicity, float(nHitsVsRow[i])/float(nMCPointsVsRow[i]) );
    }
  }
  
  for(int iRTr=0; iRTr < nRecoTracks; iRTr++){  // TODO: make common
    AliHLTTPCCAPerformanceRecoTrackData &recoD = recoData[iRTr];

    const AliHLTTPCCAGBTrack &recoTr = fTracker->Track( iRTr );  // TODO: make common
    AliHLTTPCCAMCTrack &mcTr = (*fMCTracks)[ recoD.GetMCTrackId() ];

    //AliHLTTPCCATrackParam param = t.EndPoint();
    double RecoPt  = 1. / fabs(recoTr.InnerParam().QPt());
    double RecoMom = RecoPt * sqrt(1. + recoTr.InnerParam().DzDs()*recoTr.InnerParam().DzDs());
      //     fNVsMom->Fill( param.GetY());
      //     fLengthVsMom->Fill( param.GetY(), t.NHits());

    GetHisto("purity")->Fill( recoData[iRTr].GetPurity() );
    if (  recoD.IsGhost(PParameters::MinTrackPurity) ) {
      GetHisto("ghostsLength")->Fill( recoTr.NHits() );
      GetHisto("ghostsRMom")->Fill( RecoMom );
      GetHisto("ghostsMCMom")->Fill( mcTr.P() );
      GetHisto("ghostsRPt")->Fill( RecoPt );
      GetHisto("ghostsMCPt")->Fill( mcTr.Pt() );
      GetHisto("ghostsLengthAndMCMom")->Fill( recoTr.NHits(), mcTr.P() );
      GetHisto("ghostsLengthAndRMom")->Fill( recoTr.NHits(), RecoMom );
      GetHisto("ghostsChi2")->Fill( recoTr.InnerParam().GetChi2() );
      GetHisto("ghostsProb")->Fill( TMath::Prob(recoTr.InnerParam().GetChi2(),recoTr.InnerParam().GetNDF()));
    }
    else {
      GetHisto("recosLength")->Fill( recoTr.NHits() );
      GetHisto("recosRMom")->Fill( RecoMom );
      GetHisto("recosMCMom")->Fill( mcTr.P() );
      GetHisto("recosRPt")->Fill( RecoPt );
      GetHisto("recosMCPt")->Fill( mcTr.Pt() );
      GetHisto("recosLengthAndMCMom")->Fill( recoTr.NHits() , mcTr.P() );
      GetHisto("recosLengthAndRMom")->Fill( recoTr.NHits() , RecoMom );
      GetHisto("recosChi2")->Fill( recoTr.InnerParam().GetChi2() );
      GetHisto("recosProb")->Fill( TMath::Prob(recoTr.InnerParam().GetChi2(),recoTr.InnerParam().GetNDF()));
      // GetHisto("recosChi2")->Fill( recoTr.OuterParam().GetChi2() );
      // GetHisto("recosProb")->Fill( TMath::Prob(recoTr.OuterParam().GetChi2(),recoTr.OuterParam().GetNDF()));
      if( mcTr.P() > 0.5 ) GetHisto("nHitsRecoTOverNHitsMCT")->Fill( float(recoTr.NHits()) / mcTrackNRecoHits[recoD.GetMCTrackId()] );
    }
  }
  
  // global tracker performance
  {
    for ( int itr = 0; itr < nRecoTracks; itr++ ) {
      const int iMC = recoData[itr].GetMCTrackId();
      if ( recoData[itr].IsGhost(PParameters::MinTrackPurity) ) continue;
      AliHLTTPCCAMCTrack &mc = (*fMCTracks)[iMC];

      // if ( mc.MotherId() > -1 ) continue;
      if ( mc.NHits() != 6 ) continue;
      if ( mc.Pt() < 1 ) continue;
      // if ( mc.DzDs() < 0.5 ) continue;
        
      int nFirstMC = mc.FirstMCPointID();
      int nMCPoints = mc.NMCPoints();

      AliHLTTPCCALocalMCPoint *points = &((*fLocalMCPoints).Data()[nFirstMC]);

      const AliHLTTPCCAGBTrack &t = fTracker->Track( itr );
//      AliHLTTPCCATrackParam p = t.Param();
       AliHLTTPCCATrackParam p = t.OuterParam();
      if( !p.IsValid() ) continue;
      
      int nHits = t.NHits();
      if(nHits < 3) continue;
// int nHits = t.NHits();
// std::cout << "hits:" << std::endl;
// for(int iH=0; iH<nHits; iH++)
//   std::cout << fTracker->Hit( fTracker->TrackHit(t.FirstHitRef() + iH) ).X() << " " <<
//                fTracker->Hit( fTracker->TrackHit(t.FirstHitRef() + iH) ).Y() << " " << 
//                fTracker->Hit( fTracker->TrackHit(t.FirstHitRef() + iH) ).Z() << " " << std::endl; 
// std::cout << "MCPoints:" << std::endl;
// for(int iMCPoint=0; iMCPoint<nMCPoints; iMCPoint++)
//   std::cout << points[iMCPoint].X() << " " <<
//                points[iMCPoint].Y() << " " << 
//                points[iMCPoint].Z() << " " << std::endl; 
// 
// std::cout << "track:" << std::endl;
// std::cout << p.X() << " " << p.Y() << " " << p.Z() << std::endl;
      int MCindex=-1;
      for(int iMCPoint=0; iMCPoint<nMCPoints; iMCPoint++)
      {
        double xPoint = points[iMCPoint].X();
        double yPoint = points[iMCPoint].Y();
        double xLoc =  xPoint*cos(p.Angle()) + yPoint*sin(p.Angle());
        // double yLoc = -xPoint*sin(p.Angle()) + yPoint*cos(p.Angle());

        // cout << xLoc << " " << p.X() << " " << p.Y()<< " " << yLoc << " " << p.Z()<< " " << points[iMCPoint].Z() << endl;
        if(fabs(xLoc - p.X()) > 2.f) continue;
        
          // find hit
        double y = -1e10, z = -1e10;
        for(int iH=0; iH<nHits; iH++) {
          int ith = t.FirstHitRef() + iH;
          int ih = fTracker->TrackHits()[ith];
          const AliHLTTPCCAGBHit& hit = fTracker->Hit( ih );
          float xtmp, ytmp;
//           hit.GetLocalXY( xtmp, ytmp );
          xtmp =  hit.X()*cos(p.Angle()) + hit.Y()*sin(p.Angle());
          ytmp = -hit.X()*sin(p.Angle()) + hit.Y()*cos(p.Angle());
          if(fabs(xtmp - p.X()) > 2.f) continue;
          y = hit.Y();
          z = hit.Z();
        }
          
        // if(fabs(p.Y() - yLoc)<2 && fabs(p.Z() - points[iMCPoint].Z())<2)
        if(fabs(y - yPoint)<2 && fabs(z - points[iMCPoint].Z())<2)
          MCindex = iMCPoint;
        
      }
      if(MCindex == -1)
      {
        continue;
      }

      // track resolutions
      while ( 1/*mc.Set() == 2 && TMath::Abs( mc.TPCPar()[0] ) + TMath::Abs( mc.TPCPar()[1] ) > 1*/ ) {
/*        double cosA = TMath::Cos( t.Angle() );
        double sinA = TMath::Sin( t.Angle() );
        double mcX =  mc.TPCPar()[0] * cosA + mc.TPCPar()[1] * sinA;
        double mcY = -mc.TPCPar()[0] * sinA + mc.TPCPar()[1] * cosA;
        double mcZ =  mc.TPCPar()[2];
        double mcEx =  mc.TPCPar()[3] * cosA + mc.TPCPar()[4] * sinA;
        double mcEy = -mc.TPCPar()[3] * sinA + mc.TPCPar()[4] * cosA;
        double mcEz =  mc.TPCPar()[5];
        double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );
        if ( TMath::Abs( mcEt ) < 1.e-4 ) break;
        double mcSinPhi = mcEy / mcEt;
        double mcDzDs   = mcEz / mcEt;
        double mcQPt = mc.TPCPar()[6] / mcEt;
        if ( TMath::Abs( mcQPt ) < 1.e-4 ) break;
        double mcPt = 1. / TMath::Abs( mcQPt );
        if ( mcPt < Parameters::RefThreshold ) break;
        if ( t.NHits() <  PParameters::MinimumHitsForMCTrack ) break;
        double bz = fTracker->Slice( 0 ).Param().Bz();*/

        double angle = p.Angle();
//         double mcX =  points[MCindex].X();
//         double mcY =  points[MCindex].Y();

        double mcX0 =  points[MCindex].X();
        double mcY0 =  points[MCindex].Y();
        double mcX =  mcX0*cos(angle) +  mcY0*sin(angle);
        double mcY = -mcX0*sin(angle) +  mcY0*cos(angle);
        double mcZ =  points[MCindex].Z();
        double mcPX0 =  points[MCindex].Px();
        double mcPY0 =  points[MCindex].Py();
        double mcPX =  mcPX0*cos(angle) +  mcPY0*sin(angle);
        double mcPY = -mcPX0*sin(angle) +  mcPY0*cos(angle);

//        double mcZ =  points[MCindex].Z();
        double mcQP = points[MCindex].QP();
        double mcEx = mcPX*(mcQP);
        double mcEy = mcPY*(mcQP);
        double mcEz = points[MCindex].Pz()*(mcQP);
        double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );

        if ( TMath::Abs( mcEt ) < 1.e-4 ) break;
        double mcSinPhi = mcEy / mcEt;
        double mcDzDs   = mcEz / mcEt;
        double mcQPt =  mcQP / mcEt;

        if ( TMath::Abs( mcQPt ) < 1.e-6 ) break;
        double mcPt = 1. / TMath::Abs( mcQPt );
//        if ( mcPt < Parameters::RefThreshold ) break;
//	if ( t.NHits() <  PParameters::MinimumHitsForMCTrack ) break;
        double bz = fTracker->GetParameters().cBz();

        if ( !p.TransportToXWithMaterial( mcX, bz ) ) break;

        bool dirOk = 1;
        if( (mcX > 0) && (mcEx < 0) ) dirOk =0;
        if( (mcX < 0) && (mcEx > 0) ) dirOk =0;
       if ( p.GetCosPhi()*mcEx < 0 ) { // change direction
//         if(!dirOk) {
          mcSinPhi = -mcSinPhi;
          mcDzDs = -mcDzDs;
//           mcQPt = -mcQPt;
        }
//         const double kCLight = 0.000299792458;
//         double k2QPt = 100;
//         if ( TMath::Abs( bz ) > 1.e-4 ) k2QPt = 1. / ( bz * kCLight );
/*        double qPt = p.GetKappa( bz ) * k2QPt;
        double pt = 100;*/
	double qPt = p.GetQPt(); 
        // double pt = 1. / TMath::Abs( qPt );
        if ( mcPt < 0.010 ) break;
//        if ( TMath::Abs( qPt ) > 1.e-4 ) pt = 1. / TMath::Abs( qPt );
        GetHisto("resY")->Fill( p.GetY() - mcY );
        GetHisto("resZ")->Fill( p.GetZ() - mcZ );
        GetHisto("resSinPhi")->Fill( p.GetSinPhi() - mcSinPhi );
        GetHisto("resDzDs")->Fill( p.GetDzDs() - mcDzDs );
        if(CAMath::Abs(qPt) > 1.e-8){
          GetHisto("resPt")->Fill( (qPt - mcQPt)/mcQPt );
        }
        if ( p.GetErr2Y() > 0 ) GetHisto("pullY")->Fill( ( p.GetY() - mcY ) / TMath::Sqrt( p.GetErr2Y() ) );
        if ( p.GetErr2Z() > 0 ) GetHisto("pullZ")->Fill( ( p.GetZ() - mcZ ) / TMath::Sqrt( p.GetErr2Z() ) );

        if ( p.GetErr2SinPhi() > 0 ) GetHisto("pullSinPhi")->Fill( ( p.GetSinPhi() - mcSinPhi ) / TMath::Sqrt( p.GetErr2SinPhi() ) );
        if ( p.GetErr2DzDs() > 0 ) GetHisto("pullDzDs")->Fill( ( p.DzDs() - mcDzDs ) / TMath::Sqrt( p.GetErr2DzDs() ) );
        if(CAMath::Abs(qPt) > 1.e-7 && p.GetErr2QPt()>0 ) GetHisto("pullQPt")->Fill( (qPt - mcQPt)/TMath::Sqrt(p.GetErr2QPt()) );

        break;
      }
    }

  }

    // hits pulls
  { // TODO
    const int nHits = fTracker->NHits();
    for ( int ih = 0; ih < nHits; ih++ ) {
      const AliHLTTPCCAGBHit &hit = fTracker->Hit( ih );
      const AliHLTTPCCAHitLabel &l = (*fHitLabels)[hit.ID()];
//       fhHitErrY->Fill( hit.ErrY() );
//       fhHitErrZ->Fill( hit.ErrZ() );

      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++)
      {
        double xPoint = points[iMCPoint].X();
        double yPoint = points[iMCPoint].Y();
        double xLoc =  xPoint*cos(hit.Angle()) + yPoint*sin(hit.Angle());
        
        if(fabs(xLoc - hit.X()) > 2.f) continue;
 
        double y = hit.Y();
        double z = hit.Z();
  
          // if(fabs(hit.Y() - yLoc)<2 && fabs(hit.Z() - points[iMCPoint].Z())<2)
        if(fabs(y - yPoint)<2 && fabs(z - points[iMCPoint].Z())<2)
          MCindex = iMCPoint;
      }
      
      if(MCindex == -1)
      {
        continue;
      }


      const AliHLTTPCCALocalMCPoint& point = points[MCindex];
      
      GetHisto("resXHit")->Fill( hit.X() - point.X() );
      GetHisto("resYHit")->Fill( hit.Y() - point.Y() );
      GetHisto("resZHit")->Fill( hit.Z() - point.Z() );
      GetHisto("resXHitVsZ")->Fill( hit.Z(), hit.X() - point.X() );
      GetHisto("resYHitVsZ")->Fill( hit.Z(), hit.Y() - point.Y() );
      GetHisto("resZHitVsZ")->Fill( hit.Z(), hit.Z() - point.Z() );
      GetHisto("resXHitVsX")->Fill( hit.X(), hit.X() - point.X() );
      GetHisto("resYHitVsX")->Fill( hit.X(), hit.Y() - point.Y() );
      GetHisto("resZHitVsX")->Fill( hit.X(), hit.Z() - point.Z() );

      GetHisto("xMCPoint")->Fill( point.X() );
      GetHisto("rMCPoint")->Fill( sqrt(point.Y()*point.Y() + point.X()*point.X()) );

        // get errors
      double mcEx = point.Px();
      double mcEy = point.Py();
      double mcEz = point.Pz();
      double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );
      
      sfloat_v Err2Y = 0, Err2Z = 0;
      const AliHLTTPCCAParam par;
      TrackParamVector t;
      t.SetZ( point.Z() );
      t.SetSinPhi(mcEy / mcEt);
      t.SetDzDs(mcEz / mcEt);

      // GetHisto("pullXHit")->Fill( hit.X() - point.X() );
      GetHisto("pullYHit")->Fill( (hit.Y() - point.Y())/sqrt(hit.Err2Y()) );
      GetHisto("pullZHit")->Fill( (hit.Z() - point.Z())/sqrt(hit.Err2Z()) );
    }
  }


  
} // void AliHLTTPCCAGlobalPerformance::FillHistos()

void AliHLTTPCCAGlobalPerformance::Draw()
{
#ifdef DRAW_GLOBALPERF
  // AliHLTTPCCAPerformance::Instance().Init();
  AliHLTTPCCAPerformance& gbPerfo = AliHLTTPCCAPerformance::Instance();
  AliHLTTPCCADisplay &disp = AliHLTTPCCADisplay::Instance();
  // disp.SetGB( gbPerfo.GetTracker() );
  disp.SetTPCView();
  disp.DrawTPC();
  disp.DrawGBHits( *gbPerfo.GetTracker(), -1, 0.1 );

  for ( int imc = 0; imc < nMCTracks; imc++ ) {
    AliHLTTPCCAPerformanceMCTrackData &mc = mcData[imc];
    bool doDraw = true;
    doDraw &= (mc.GetSet() == 7); // long ref prim
    doDraw &= (*fMCTracks)[imc].P() > 1;
    doDraw &= mc.IsReconstructable();
    // doDraw &= mc.IsReconstructed();
    if ( doDraw ) {
      disp.DrawMCTrack( imc, kRed+2, 0 );
    }
  }

  for( int iT = 0; iT < nRecoTracks; ++iT ) {
    AliHLTTPCCAPerformanceRecoTrackData &r = recoData[iT];
    const bool IsGhost = r.IsGhost(PParameters::MinTrackPurity);
    if ( IsGhost ) {
      disp.DrawRecoTrack( iT, kGreen+1, 0 );
    }
    else {
      const int imc = r.GetMCTrackId();
      AliHLTTPCCAPerformanceMCTrackData &mc = mcData[imc];
      bool doDraw = true;
      doDraw &= (mc.GetSet() == 7); // long ref prim
      doDraw &= (*fMCTracks)[imc].P() > 1;
      doDraw &= mc.IsReconstructable();
        // doDraw &= mc.IsReconstructed();
      if ( doDraw ) disp.DrawRecoTrack( iT, kBlue+1, 0 );
    }
  }

  disp.SaveCanvasToFile( "DrawGlobalPerformance.pdf" );
  disp.Ask();
  
#endif // DRAW_GLOBALPERF
} // void AliHLTTPCCAGlobalPerformance::Draw()

#endif //DO_TPCCATRACKER_EFF_PERFORMANCE