// $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 * // Ivan Kisel * // for The ALICE HLT Project. * // * // Developed by: Igor Kulakov * // Maksym Zyzak * // * // 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 *hitLabels, AliHLTResizableArray *mcTracks, AliHLTResizableArray *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.NHits() >= PParameters::MinimumHitsForMCTrack ) // 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 mcTrackNRecoHits; vector nHitsVsRow; vector 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 RecoP = 1. / fabs(recoTr.InnerParam().QP()); double RecoMom = RecoP; // 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("ghostsRP")->Fill( RecoP ); GetHisto("ghostsMCP")->Fill( mcTr.P() ); GetHisto("ghostsLengthAndMCMom")->Fill( recoTr.NHits(), mcTr.P() ); GetHisto("ghostsLengthAndRMom")->Fill( recoTr.NHits(), RecoMom ); GetHisto("ghostsChi2")->Fill( recoTr.InnerParam().Chi2() ); GetHisto("ghostsProb")->Fill( TMath::Prob(recoTr.InnerParam().Chi2(),recoTr.InnerParam().NDF())); } else { GetHisto("recosLength")->Fill( recoTr.NHits() ); GetHisto("recosRMom")->Fill( RecoMom ); GetHisto("recosMCMom")->Fill( mcTr.P() ); GetHisto("recosRP")->Fill( RecoP ); GetHisto("recosMCP")->Fill( mcTr.P() ); GetHisto("recosLengthAndMCMom")->Fill( recoTr.NHits() , mcTr.P() ); GetHisto("recosLengthAndRMom")->Fill( recoTr.NHits() , RecoMom ); GetHisto("recosChi2")->Fill( recoTr.InnerParam().Chi2() ); GetHisto("recosProb")->Fill( TMath::Prob(recoTr.InnerParam().Chi2(),recoTr.InnerParam().NDF())); // 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.P() < 1 ) continue; int nFirstMC = mc.FirstMCPointID(); int nMCPoints = mc.NMCPoints(); AliHLTTPCCALocalMCPoint *points = &((*fLocalMCPoints).Data()[nFirstMC]); const AliHLTTPCCAGBTrack &t = fTracker->Track( itr ); #if 0 AliHLTTPCCATrackParam p = t.Param(); #else AliHLTTPCCATrackParam p = t.OuterParam(); #endif 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; iHHit( 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 2.f) continue; // find hit double y = -1e10, z = -1e10; for(int iH=0; iHTrackHits()[ith]; const AliHLTTPCCAGBHit& hit = fTracker->Hit( ih ); // hit.GetLocalXY( xtmp, ytmp ); float xtmp = hit.X()*cos(p.Angle()) + hit.Y()*sin(p.Angle()); // float 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 mcPZ = points[MCindex].Pz(); // double mcZ = points[MCindex].Z(); double mcQP = points[MCindex].QP(); // double mcP = fabs(1/mcQP); double mcP = sqrt( mcPX*mcPX + mcPY*mcPY + mcPZ*mcPZ ); double mcTx = mcPX/mcPZ; double mcTy = mcPY/mcPZ; // if ( mcPt < Parameters::RefThreshold ) break; // if ( t.NHits() < PParameters::MinimumHitsForMCTrack ) break; if ( !p.TransportToX0Line( mcZ ) ) break; double qp = p.QP(); if ( mcP < 0.01 ) break; // if ( TMath::Abs( qPt ) > 1.e-4 ) pt = 1. / TMath::Abs( qPt ); GetHisto("resX")->Fill( (p.X() - mcX)*10000 ); // mum GetHisto("resY")->Fill( (p.Y() - mcY)*10000 ); GetHisto("resTx")->Fill( (p.Tx() - mcTx)*1000 ); // 10^-3 GetHisto("resTy")->Fill( (p.Ty() - mcTy)*1000 ); if(CAMath::Abs(qp) > 1.e-8){ GetHisto("resP")->Fill( (fabs(1.f/qp) - mcP)/mcP*100 ); // % } if ( p.Err2X() > 0 ) GetHisto("pullX")->Fill( ( p.X() - mcX ) / TMath::Sqrt( p.Err2X() ) ); if ( p.Err2Y() > 0 ) GetHisto("pullY")->Fill( ( p.Y() - mcY ) / TMath::Sqrt( p.Err2Y() ) ); if ( p.Err2Tx() > 0 ) GetHisto("pullTx")->Fill( ( p.Tx() - mcTx ) / TMath::Sqrt( p.Err2Tx() ) ); if ( p.Err2Ty() > 0 ) GetHisto("pullTy")->Fill( ( p.Ty() - mcTy ) / TMath::Sqrt( p.Err2Ty() ) ); if(CAMath::Abs(qp) > 1.e-7 && p.Err2QP()>0 ) GetHisto("pullQP")->Fill( (qp - mcQP)/TMath::Sqrt(p.Err2QP()) ); break; } } } // distribution of cluster errors { 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() ); int nmc = 0; for ( int il = 0; il < 3; il++ ) if ( l.fLab[il] >= 0 ) nmc++; // if ( nmc == 1 ) fhHitShared->Fill( hit.IRow(), 0 ); // else if ( nmc > 1 ) fhHitShared->Fill( hit.IRow(), 1 ); } } } // 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