// **************************************************************************
//                                                                          *
// 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 "AliHLTTPCTopoPerformance.h"

#include "AliHLTTPCCounters.h"

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



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

#include "AliHLTTPCCATracker.h"

#include "AliHLTTPCCADisplay.h"


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

AliHLTTPCTopoPerformance::AliHLTTPCTopoPerformance():fTopoReconstructor(0),fPrimVertices()
{
}

void AliHLTTPCTopoPerformance::SetNewEvent(
                            const AliHLTTPCCAGBTracker * const tracker,
                            AliHLTResizableArray<AliHLTTPCCAHitLabel> *hitLabels,
                            AliHLTResizableArray<AliHLTTPCCAMCTrack> *mcTracks,
                            AliHLTResizableArray<AliHLTTPCCALocalMCPoint> *localMCPoints)
{
  assert( AliHLTTPCCAPerformance::Instance().GetSubPerformance("Global Performance") != 0 );
  
  AliHLTTPCParticlePerformanceBase::SetNewEvent(tracker, hitLabels, mcTracks, localMCPoints);

  nRecoTracks = fTracker->NTracks();
} // void AliHLTTPCTopoPerformance::SetNewEvent

void AliHLTTPCTopoPerformance::SetNewEvent2( const AliHLTTPCTopoReconstructor * const TopoReconstructor)
{
  assert( fTracker != 0 ); // should be called after SetNewEvent
  
  fTopoReconstructor = TopoReconstructor;
} // void AliHLTTPCTopoPerformance::SetNewEvent2

void AliHLTTPCTopoPerformance::CheckMCTracks()
{
  mcData = AliHLTTPCCAPerformance::Instance().GetSubPerformance("Global Performance")->GetMCData();

    // find prim vertex
  AliHLTTPCCAMCTrack *iPrimMC = (*fMCTracks).Data();
  for( ;(iPrimMC->MotherId() != -1) && (iPrimMC < (*fMCTracks).Data() + nMCTracks); iPrimMC++);
  if (iPrimMC < (*fMCTracks).Data() + nMCTracks) {
    AliHLTTPCCAMCVertex primVertex;
    primVertex.SetX( iPrimMC->X() );
    primVertex.SetY( iPrimMC->Y() );
    primVertex.SetZ( iPrimMC->Z() );
    fPrimVertices.push_back(primVertex);
  }
} // void AliHLTTPCTopoPerformance::CheckMCTracks()


void AliHLTTPCTopoPerformance::MatchTracks()
{
  recoData = AliHLTTPCCAPerformance::Instance().GetSubPerformance("Global Performance")->GetRecoData();
} // void AliHLTTPCTopoPerformance::MatchTracks()

void AliHLTTPCTopoPerformance::FillHistos()
{
  AliHLTTPCParticlePerformanceBase::FillHistos();

  const int nParticles = fTopoReconstructor->fRTrackIds.size();

  for(int iP = 0; iP < nParticles; iP++){

    const int itr = fTopoReconstructor->fRTrackIds[iP];
    const int iMC = recoData[itr].GetMCTrackId();

    if ( recoData[itr].IsGhost(PParameters::MinTrackPurity) ) continue;

    AliHLTTPCCAMCTrack &mc = (*fMCTracks)[iMC];


      //      const AliHLTTPCCAGBTrack &t = fTracker->Track( itr );
      //      AliHLTTPCCATrackParam p = t.Param();
    AliKFParticle &p = fTopoReconstructor->fKFPTopoReconstructor->GetParticle(iP);

    for(;;){ // technical DO // track pulls and residuals at track end
        // find point which corresponds to our track end - mcPoint
      int nFirstMC = mc.FirstMCPointID();
      int nMCPoints = mc.NMCPoints();
      AliHLTTPCCALocalMCPoint *points = &((*fLocalMCPoints).Data()[nFirstMC]);
        
      int MCindex=-1;
      AliHLTTPCCALocalMCPoint mcPoint;
      for(int iMCPoint=0; iMCPoint<nMCPoints; iMCPoint++)
      {
          // get global CS
        mcPoint = points[iMCPoint];
        const float alpha = - fTracker->Slice(mcPoint.ISlice()).SliceAlpha();
        mcPoint.RotateXY( alpha );
          
          // check if the point correspond to our track end
        if (fabs(p.X() - mcPoint.X()) < 2.f) {
          if (fabs(p.Y() - mcPoint.Y()) < 2.f) {
            if (fabs(p.Z() - mcPoint.Z()) < 2.f) {
              MCindex = iMCPoint;
              break;
            }
          }
        }
      }
      if(MCindex == -1) break;

        // get mc info
      const float mcX =  mcPoint.X();
      const float mcY =  mcPoint.Y();
      const float mcZ =  mcPoint.Z();
      const float mcPx = mcPoint.Px();
      const float mcPy = mcPoint.Py();
      const float mcPz = mcPoint.Pz();

      if ( CAMath::Sqrt( mcPx*mcPx + mcPy*mcPy /*+ mcPz*mcPz*/ ) < 0.010 ) break;
      
      AliKFParticle &pTmp = p;
      
      GetHisto("resX")->Fill ( pTmp.GetX()  - mcX  );
      GetHisto("resY")->Fill ( pTmp.GetY()  - mcY  );
      GetHisto("resZ")->Fill ( pTmp.GetZ()  - mcZ  );
      GetHisto("resPx")->Fill( pTmp.GetPx() - mcPx );
      GetHisto("resPy")->Fill( pTmp.GetPy() - mcPy );
      GetHisto("resPz")->Fill( pTmp.GetPz() - mcPz );
        
      if ( pTmp.GetErrX()  > 0 ) GetHisto("pullX")->Fill ( ( pTmp.GetX()  - mcX  ) / pTmp.GetErrX()  );
      if ( pTmp.GetErrY()  > 0 ) GetHisto("pullY")->Fill ( ( pTmp.GetY()  - mcY  ) / pTmp.GetErrY()  );
      if ( pTmp.GetErrZ()  > 0 ) GetHisto("pullZ")->Fill ( ( pTmp.GetZ()  - mcZ  ) / pTmp.GetErrZ()  );
      if ( pTmp.GetErrPx() > 0 ) GetHisto("pullPx")->Fill( ( pTmp.GetPx() - mcPx ) / pTmp.GetErrPx() );
      if ( pTmp.GetErrPy() > 0 ) GetHisto("pullPy")->Fill( ( pTmp.GetPy() - mcPy ) / pTmp.GetErrPy() );
      if ( pTmp.GetErrPz() > 0 ) GetHisto("pullPz")->Fill( ( pTmp.GetPz() - mcPz ) / pTmp.GetErrPz() );

      break;
    }
    
    for(;;){ // technical DO // track pulls and residuals at vertex
    
        // get mc info
      const float mcX =  mc.X();
      const float mcY =  mc.Y();
      const float mcZ =  mc.Z();
      const float mcPx = mc.Px();
      const float mcPy = mc.Py();
      const float mcPz = mc.Pz();

        // transport to vertex
      AliKFParticle pTmp = p;
      const float vertex[3] = { mcX, mcY, mcZ };
      pTmp.TransportToDS(  pTmp.GetDStoPoint(vertex) );
      
      if ( mcX*mcX+mcY*mcY < 50*50 ) break; // dbg // TPC-vertices only
//      if ( pTmp.GetPt() < 1.f ) break; // dbg
//      if ( mc.MotherId() != -1 ) break; // dbg // primary only !!
      
        // fill histos

      if ( CAMath::Sqrt( mcPx*mcPx + mcPy*mcPy /*+ mcPz*mcPz*/ ) < 0.010 ) break;
        //        if ( TMath::Abs( qPt ) > 1.e-4 ) pt = 1. / TMath::Abs( qPt );

      // std::cout << " MC " << mcX << " " << mcY << std::endl;
      // std::cout << " diff " << pTmp.GetX()  - mcX << " " << pTmp.GetY()  - mcY << std::endl;
      GetHisto("resAtProdPointX")->Fill ( pTmp.GetX()  - mcX  );
      GetHisto("resAtProdPointY")->Fill ( pTmp.GetY()  - mcY  );
      GetHisto("resAtProdPointZ")->Fill ( pTmp.GetZ()  - mcZ  );
      GetHisto("resAtProdPointPx")->Fill( pTmp.GetPx() - mcPx );
      GetHisto("resAtProdPointPy")->Fill( pTmp.GetPy() - mcPy );
      GetHisto("resAtProdPointPz")->Fill( pTmp.GetPz() - mcPz );
        
      if ( pTmp.GetErrX()  > 0 ) GetHisto("pullAtProdPointX")->Fill ( ( pTmp.GetX()  - mcX  ) / pTmp.GetErrX()  );
      if ( pTmp.GetErrY()  > 0 ) GetHisto("pullAtProdPointY")->Fill ( ( pTmp.GetY()  - mcY  ) / pTmp.GetErrY()  );
      if ( pTmp.GetErrZ()  > 0 ) GetHisto("pullAtProdPointZ")->Fill ( ( pTmp.GetZ()  - mcZ  ) / pTmp.GetErrZ()  );
      if ( pTmp.GetErrPx() > 0 ) GetHisto("pullAtProdPointPx")->Fill( ( pTmp.GetPx() - mcPx ) / pTmp.GetErrPx() );
      if ( pTmp.GetErrPy() > 0 ) GetHisto("pullAtProdPointPy")->Fill( ( pTmp.GetPy() - mcPy ) / pTmp.GetErrPy() );
      if ( pTmp.GetErrPz() > 0 ) GetHisto("pullAtProdPointPz")->Fill( ( pTmp.GetPz() - mcPz ) / pTmp.GetErrPz() );

      break;
    }

    for(;;){ // technical DO // vertex pulls and residuals
      if (fPrimVertices.size() == 0) continue; // no prim tracks
      
        // get mc info
      const float mcX =  fPrimVertices[0].X();
      const float mcY =  fPrimVertices[0].Y();
      const float mcZ =  fPrimVertices[0].Z();

        // prim vertex
      AliKFParticle &pTmp = fTopoReconstructor->fKFPTopoReconstructor->GetPrimVertex();

        // fill histos

      GetHisto("resVertexX")->Fill ( pTmp.GetX()  - mcX  );
      GetHisto("resVertexY")->Fill ( pTmp.GetY()  - mcY  );
      GetHisto("resVertexZ")->Fill ( pTmp.GetZ()  - mcZ  );
        
      if ( pTmp.GetErrX()  > 0 ) GetHisto("pullVertexX")->Fill ( ( pTmp.GetX()  - mcX  ) / pTmp.GetErrX()  );
      if ( pTmp.GetErrY()  > 0 ) GetHisto("pullVertexY")->Fill ( ( pTmp.GetY()  - mcY  ) / pTmp.GetErrY()  );
      if ( pTmp.GetErrZ()  > 0 ) GetHisto("pullVertexZ")->Fill ( ( pTmp.GetZ()  - mcZ  ) / pTmp.GetErrZ()  );

      break;
    }
    
  }
  
} // void AliHLTTPCTopoPerformance::FillHistos()

#endif //DO_TPCCATRACKER_EFF_PERFORMANCE