HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.cxx

   1 // $Id$
   2 // **************************************************************************
   3 // This file is property of and copyright by the ALICE HLT Project          *
   4 // ALICE Experiment at CERN, All rights reserved.                           *
   5 //                                                                          *
   6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
   7 //                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
   8 //                  for The ALICE HLT Project.                              *
   9 //                                                                          *
  10 // Permission to use, copy, modify and distribute this software and its     *
  11 // documentation strictly for non-commercial purposes is hereby granted     *
  12 // without fee, provided that the above copyright notice appears in all     *
  13 // copies and that both the copyright notice and this permission notice     *
  14 // appear in the supporting documentation. The authors make no claims       *
  15 // about the suitability of this software for any purpose. It is            *
  16 // provided "as is" without express or implied warranty.                    *
  17 //                                                                          *
  18 //***************************************************************************
  19
  20
  21 #include "AliHLTTPCCAPerformance.h"
  22 #include "AliHLTTPCCAMCTrack.h"
  23 #include "AliHLTTPCCAMCPoint.h"
  24 #include "AliHLTTPCCAOutTrack.h"
  25 #include "AliHLTTPCCATracker.h"
  26 #include "AliHLTTPCCATracklet.h"
  27 #include "AliHLTTPCCAStandaloneFramework.h"
  28 #include "AliHLTTPCCASliceTrack.h"
  29 #include "AliHLTTPCCASliceOutput.h"
  30 #include "AliHLTTPCCAMergerOutput.h"
  31 #include "AliHLTTPCCAMergedTrack.h"
  32
  33 #include "TMath.h"
  34 #include "TROOT.h"
  35 #include "Riostream.h"
  36 #include "TFile.h"
  37 #include "TH1.h"
  38 #include "TH2.h"
  39 #include "TProfile.h"
  40 #include "TRandom.h"
  41 #include <cmath>
  42
  43 AliHLTTPCCAPerformance &AliHLTTPCCAPerformance::Instance()
  44 {
  45   // reference to static object
  46   static AliHLTTPCCAPerformance gAliHLTTPCCAPerformance;
  47   return gAliHLTTPCCAPerformance;
  48 }
  49
  50 AliHLTTPCCAPerformance::AliHLTTPCCAPerformance()
  51     :
  52     fHitLabels( 0 ),
  53     fNHits( 0 ),
  54     fMCTracks( 0 ),
  55     fNMCTracks( 0 ),
  56     fMCPoints( 0 ),
  57     fNMCPoints( 0 ),
  58     fDoClusterPulls( 0 ),
  59     fStatNEvents( 0 ),
  60     fStatTime( 0 ),
  61     fStatSeedNRecTot( 0 ),
  62     fStatSeedNRecOut( 0 ),
  63     fStatSeedNGhost( 0 ),
  64     fStatSeedNMCAll( 0 ),
  65     fStatSeedNRecAll( 0 ),
  66     fStatSeedNClonesAll( 0 ),
  67     fStatSeedNMCRef( 0 ),
  68     fStatSeedNRecRef( 0 ),
  69     fStatSeedNClonesRef( 0 ),
  70     fStatCandNRecTot( 0 ),
  71     fStatCandNRecOut( 0 ),
  72     fStatCandNGhost( 0 ),
  73     fStatCandNMCAll( 0 ),
  74     fStatCandNRecAll( 0 ),
  75     fStatCandNClonesAll( 0 ),
  76     fStatCandNMCRef( 0 ),
  77     fStatCandNRecRef( 0 ),
  78     fStatCandNClonesRef( 0 ),
  79     fStatNRecTot( 0 ),
  80     fStatNRecOut( 0 ),
  81     fStatNGhost( 0 ),
  82     fStatNMCAll( 0 ),
  83     fStatNRecAll( 0 ),
  84     fStatNClonesAll( 0 ),
  85     fStatNMCRef( 0 ),
  86     fStatNRecRef( 0 ),
  87     fStatNClonesRef( 0 ),
  88     fStatGBNRecTot( 0 ),
  89     fStatGBNRecOut( 0 ),
  90     fStatGBNGhost( 0 ),
  91     fStatGBNMCAll( 0 ),
  92     fStatGBNRecAll( 0 ),
  93     fStatGBNClonesAll( 0 ),
  94     fStatGBNMCRef( 0 ),
  95     fStatGBNRecRef( 0 ),
  96     fStatGBNClonesRef( 0 ),
  97     fHistoDir( 0 ),
  98     fhResY( 0 ),
  99     fhResZ( 0 ),
 100     fhResSinPhi( 0 ),
 101     fhResDzDs( 0 ),
 102     fhResPt( 0 ),
 103     fhPullY( 0 ),
 104     fhPullZ( 0 ),
 105     fhPullSinPhi( 0 ),
 106     fhPullDzDs( 0 ),
 107     fhPullQPt( 0 ),
 108     fhPullYS( 0 ),
 109     fhPullZT( 0 ),
 110     fhHitErrY( 0 ),
 111     fhHitErrZ( 0 ),
 112     fhHitResY( 0 ),
 113     fhHitResZ( 0 ),
 114     fhHitPullY( 0 ),
 115     fhHitPullZ( 0 ),
 116     fhHitShared( 0 ),
 117     fhHitResY1( 0 ),
 118     fhHitResZ1( 0 ),
 119     fhHitPullY1( 0 ),
 120     fhHitPullZ1( 0 ),
 121     fhCellPurity( 0 ),
 122     fhCellNHits( 0 ),
 123     fhCellPurityVsN( 0 ),
 124     fhCellPurityVsPt( 0 ),
 125     fhEffVsP( 0 ),
 126     fhSeedEffVsP( 0 ),
 127     fhCandEffVsP( 0 ),
 128     fhGBEffVsP( 0 ),
 129     fhGBEffVsPt( 0 ),
 130     fhNeighQuality( 0 ),
 131     fhNeighEff( 0 ),
 132     fhNeighQualityVsPt( 0 ),
 133     fhNeighEffVsPt( 0 ),
 134     fhNeighDy( 0 ),
 135     fhNeighDz( 0 ),
 136     fhNeighChi( 0 ),
 137     fhNeighDyVsPt( 0 ),
 138     fhNeighDzVsPt( 0 ),
 139     fhNeighChiVsPt( 0 ),
 140     fhNeighNCombVsArea( 0 ),
 141     fhNHitsPerSeed ( 0 ),
 142     fhNHitsPerTrackCand( 0 ),
 143     fhTrackLengthRef( 0 ),
 144     fhRefRecoX( 0 ),
 145     fhRefRecoY( 0 ),
 146     fhRefRecoZ( 0 ),
 147     fhRefRecoP( 0 ),
 148     fhRefRecoPt( 0 ),
 149     fhRefRecoAngleY( 0 ),
 150     fhRefRecoAngleZ( 0 ),
 151     fhRefRecoNHits( 0 ),
 152     fhRefNotRecoX( 0 ),
 153     fhRefNotRecoY( 0 ),
 154     fhRefNotRecoZ( 0 ),
 155     fhRefNotRecoP( 0 ),
 156     fhRefNotRecoPt( 0 ),
 157     fhRefNotRecoAngleY( 0 ),
 158     fhRefNotRecoAngleZ( 0 ),
 159     fhRefNotRecoNHits( 0 )
 160 {
 161   //* constructor
 162 }
 163
 164
 165 AliHLTTPCCAPerformance::AliHLTTPCCAPerformance( const AliHLTTPCCAPerformance& )
 166     :
 167     fHitLabels( 0 ),
 168     fNHits( 0 ),
 169     fMCTracks( 0 ),
 170     fNMCTracks( 0 ),
 171     fMCPoints( 0 ),
 172     fNMCPoints( 0 ),
 173     fDoClusterPulls( 0 ),
 174     fStatNEvents( 0 ),
 175     fStatTime( 0 ),
 176     fStatSeedNRecTot( 0 ),
 177     fStatSeedNRecOut( 0 ),
 178     fStatSeedNGhost( 0 ),
 179     fStatSeedNMCAll( 0 ),
 180     fStatSeedNRecAll( 0 ),
 181     fStatSeedNClonesAll( 0 ),
 182     fStatSeedNMCRef( 0 ),
 183     fStatSeedNRecRef( 0 ),
 184     fStatSeedNClonesRef( 0 ),
 185     fStatCandNRecTot( 0 ),
 186     fStatCandNRecOut( 0 ),
 187     fStatCandNGhost( 0 ),
 188     fStatCandNMCAll( 0 ),
 189     fStatCandNRecAll( 0 ),
 190     fStatCandNClonesAll( 0 ),
 191     fStatCandNMCRef( 0 ),
 192     fStatCandNRecRef( 0 ),
 193     fStatCandNClonesRef( 0 ),
 194     fStatNRecTot( 0 ),
 195     fStatNRecOut( 0 ),
 196     fStatNGhost( 0 ),
 197     fStatNMCAll( 0 ),
 198     fStatNRecAll( 0 ),
 199     fStatNClonesAll( 0 ),
 200     fStatNMCRef( 0 ),
 201     fStatNRecRef( 0 ),
 202     fStatNClonesRef( 0 ),
 203     fStatGBNRecTot( 0 ),
 204     fStatGBNRecOut( 0 ),
 205     fStatGBNGhost( 0 ),
 206     fStatGBNMCAll( 0 ),
 207     fStatGBNRecAll( 0 ),
 208     fStatGBNClonesAll( 0 ),
 209     fStatGBNMCRef( 0 ),
 210     fStatGBNRecRef( 0 ),
 211     fStatGBNClonesRef( 0 ),
 212     fHistoDir( 0 ),
 213     fhResY( 0 ),
 214     fhResZ( 0 ),
 215     fhResSinPhi( 0 ),
 216     fhResDzDs( 0 ),
 217     fhResPt( 0 ),
 218     fhPullY( 0 ),
 219     fhPullZ( 0 ),
 220     fhPullSinPhi( 0 ),
 221     fhPullDzDs( 0 ),
 222     fhPullQPt( 0 ),
 223     fhPullYS( 0 ),
 224     fhPullZT( 0 ),
 225     fhHitErrY( 0 ),
 226     fhHitErrZ( 0 ),
 227     fhHitResY( 0 ),
 228     fhHitResZ( 0 ),
 229     fhHitPullY( 0 ),
 230     fhHitPullZ( 0 ),
 231     fhHitShared( 0 ),
 232     fhHitResY1( 0 ),
 233     fhHitResZ1( 0 ),
 234     fhHitPullY1( 0 ),
 235     fhHitPullZ1( 0 ),
 236     fhCellPurity( 0 ),
 237     fhCellNHits( 0 ),
 238     fhCellPurityVsN( 0 ),
 239     fhCellPurityVsPt( 0 ),
 240     fhEffVsP( 0 ),
 241     fhSeedEffVsP( 0 ),
 242     fhCandEffVsP( 0 ),
 243     fhGBEffVsP( 0 ),
 244     fhGBEffVsPt( 0 ),
 245     fhNeighQuality( 0 ),
 246     fhNeighEff( 0 ),
 247     fhNeighQualityVsPt( 0 ),
 248     fhNeighEffVsPt( 0 ),
 249     fhNeighDy( 0 ),
 250     fhNeighDz( 0 ),
 251     fhNeighChi( 0 ),
 252     fhNeighDyVsPt( 0 ),
 253     fhNeighDzVsPt( 0 ),
 254     fhNeighChiVsPt( 0 ),
 255     fhNeighNCombVsArea( 0 ),
 256     fhNHitsPerSeed ( 0 ),
 257     fhNHitsPerTrackCand( 0 ),
 258     fhTrackLengthRef( 0 ),
 259     fhRefRecoX( 0 ),
 260     fhRefRecoY( 0 ),
 261     fhRefRecoZ( 0 ),
 262     fhRefRecoP( 0 ),
 263     fhRefRecoPt( 0 ),
 264     fhRefRecoAngleY( 0 ),
 265     fhRefRecoAngleZ( 0 ),
 266     fhRefRecoNHits( 0 ),
 267     fhRefNotRecoX( 0 ),
 268     fhRefNotRecoY( 0 ),
 269     fhRefNotRecoZ( 0 ),
 270     fhRefNotRecoP( 0 ),
 271     fhRefNotRecoPt( 0 ),
 272     fhRefNotRecoAngleY( 0 ),
 273     fhRefNotRecoAngleZ( 0 ),
 274     fhRefNotRecoNHits( 0 )
 275 {
 276   //* dummy
 277 }
 278
 279 const AliHLTTPCCAPerformance &AliHLTTPCCAPerformance::operator=( const AliHLTTPCCAPerformance& ) const
 280 {
 281   //* dummy
 282   return *this;
 283 }
 284
 285 AliHLTTPCCAPerformance::~AliHLTTPCCAPerformance()
 286 {
 287   //* destructor
 288   StartEvent();
 289 }
 290
 291 void AliHLTTPCCAPerformance::StartEvent()
 292 {
 293   //* clean up arrays
 294   if ( !fHistoDir )  CreateHistos();
 295   if ( fHitLabels ) delete[] fHitLabels;
 296   fHitLabels = 0;
 297   fNHits = 0;
 298   if ( fMCTracks ) delete[] fMCTracks;
 299   fMCTracks = 0;
 300   fNMCTracks = 0;
 301   if ( fMCPoints ) delete[] fMCPoints;
 302   fMCPoints = 0;
 303   fNMCPoints = 0;
 304 }
 305
 306 void AliHLTTPCCAPerformance::SetNHits( int NHits )
 307 {
 308   //* set number of hits
 309   if ( fHitLabels ) delete[] fHitLabels;
 310   fHitLabels = 0;
 311   fHitLabels = new AliHLTTPCCAHitLabel[ NHits ];
 312   fNHits = NHits;
 313 }
 314
 315 void AliHLTTPCCAPerformance::SetNMCTracks( int NumberOfMCTracks )
 316 {
 317   //* set number of MC tracks
 318   if ( fMCTracks ) delete[] fMCTracks;
 319   fMCTracks = 0;
 320   fMCTracks = new AliHLTTPCCAMCTrack[ NumberOfMCTracks ];
 321   fNMCTracks = NumberOfMCTracks;
 322 }
 323
 324 void AliHLTTPCCAPerformance::SetNMCPoints( int NMCPoints )
 325 {
 326   //* set number of MC points
 327   if ( fMCPoints ) delete[] fMCPoints;
 328   fMCPoints = 0;
 329   fMCPoints = new AliHLTTPCCAMCPoint[ NMCPoints ];
 330   fNMCPoints = 0;
 331 }
 332
 333 void AliHLTTPCCAPerformance::ReadHitLabel( int HitID,
 334     int lab0, int lab1, int lab2 )
 335 {
 336   //* read the hit labels
 337   AliHLTTPCCAHitLabel hit;
 338   hit.fLab[0] = lab0;
 339   hit.fLab[1] = lab1;
 340   hit.fLab[2] = lab2;
 341   fHitLabels[HitID] = hit;
 342 }
 343
 344 void AliHLTTPCCAPerformance::ReadMCTrack( int index, const TParticle *part )
 345 {
 346   //* read mc track to the local array
 347   fMCTracks[index] = AliHLTTPCCAMCTrack( part );
 348 }
 349
 350 void AliHLTTPCCAPerformance::ReadMCTPCTrack( int index, float X, float Y, float Z,
 351     float Px, float Py, float Pz )
 352 {
 353   //* read mc track parameters at TPC
 354   fMCTracks[index].SetTPCPar( X, Y, Z, Px, Py, Pz );
 355 }
 356
 357 void AliHLTTPCCAPerformance::ReadMCPoint( int TrackID, float X, float Y, float Z, float Time, int iSlice )
 358 {
 359   //* read mc point to the local array
 360   AliHLTTPCCAMCPoint &p = fMCPoints[fNMCPoints];
 361   p.SetTrackID( TrackID );
 362   p.SetX( X );
 363   p.SetY( Y );
 364   p.SetZ( Z );
 365   p.SetTime( Time );
 366   p.SetISlice( iSlice );
 367   float sx, sy, sz;
 368   AliHLTTPCCAStandaloneFramework::Instance().Param( iSlice ).Global2Slice( X, Y, Z, &sx, &sy, &sz );
 369   p.SetSx( sx );
 370   p.SetSy( sy );
 371   p.SetSz( sz );
 372   if ( X*X + Y*Y > 10. ) fNMCPoints++;
 373 }
 374
 375 void AliHLTTPCCAPerformance::CreateHistos()
 376 {
 377   //* create performance histogramms
 378   TDirectory *curdir = gDirectory;
 379   fHistoDir = gROOT->mkdir( "HLTTPCCATrackerPerformance" );
 380   fHistoDir->cd();
 381
 382   gDirectory->mkdir( "Links" );
 383   gDirectory->cd( "Links" );
 384
 385   fhLinkEff[0] = new TProfile( "fhLinkEffPrimRef", "fhLinkEffPrimRef vs row", 156, 2., 158. );
 386   fhLinkEff[1] = new TProfile( "fhLinkEffPrimExt", "fhLinkEffPrimExt vs row", 156, 2., 158. );
 387   fhLinkEff[2] = new TProfile( "fhLinkEffSecRef", "fhLinkEffSecRef vs row", 156, 2., 158. );
 388   fhLinkEff[3] = new TProfile( "fhLinkEffSecExt", "fhLinkEffSecExt vs row", 156, 2., 158. );
 389   fhLinkAreaY[0] = new TH1D( "fhLinkAreaYPrimRef", "fhLinkAreaYPrimRef", 100, 0, 10 );
 390   fhLinkAreaZ[0] = new TH1D( "fhLinkAreaZPrimRef", "fhLinkAreaZPrimRef", 100, 0, 10 );
 391   fhLinkAreaY[1] = new TH1D( "fhLinkAreaYPrimExt", "fhLinkAreaYPrimExt", 100, 0, 10 );
 392   fhLinkAreaZ[1] = new TH1D( "fhLinkAreaZPrimExt", "fhLinkAreaZPrimExt", 100, 0, 10 );
 393   fhLinkAreaY[2] = new TH1D( "fhLinkAreaYSecRef", "fhLinkAreaYSecRef", 100, 0, 10 );
 394   fhLinkAreaZ[2] = new TH1D( "fhLinkAreaZSecRef", "fhLinkAreaZSecRef", 100, 0, 10 );
 395   fhLinkAreaY[3] = new TH1D( "fhLinkAreaYSecExt", "fhLinkAreaYSecExt", 100, 0, 10 );
 396   fhLinkAreaZ[3] = new TH1D( "fhLinkAreaZSecExt", "fhLinkAreaZSecExt", 100, 0, 10 );
 397   fhLinkChiRight[0] = new TH1D( "fhLinkChiRightPrimRef", "fhLinkChiRightPrimRef", 100, 0, 10 );
 398   fhLinkChiRight[1] = new TH1D( "fhLinkChiRightPrimExt", "fhLinkChiRightPrimExt", 100, 0, 10 );
 399   fhLinkChiRight[2] = new TH1D( "fhLinkChiRightSecRef", "fhLinkChiRightSecRef", 100, 0, 10 );
 400   fhLinkChiRight[3] = new TH1D( "fhLinkChiRightSecExt", "fhLinkChiRightSecExt", 100, 0, 10 );
 401   fhLinkChiWrong[0] = new TH1D( "fhLinkChiWrongPrimRef", "fhLinkChiWrongPrimRef", 100, 0, 10 );
 402   fhLinkChiWrong[1] = new TH1D( "fhLinkChiWrongPrimExt", "fhLinkChiWrongPrimExt", 100, 0, 10 );
 403   fhLinkChiWrong[2] = new TH1D( "fhLinkChiWrongSecRef", "fhLinkChiWrongSecRef", 100, 0, 10 );
 404   fhLinkChiWrong[3] = new TH1D( "fhLinkChiWrongSecExt", "fhLinkChiWrongSecExt", 100, 0, 10 );
 405
 406   gDirectory->cd( ".." );
 407
 408   gDirectory->mkdir( "Neighbours" );
 409   gDirectory->cd( "Neighbours" );
 410
 411   fhNeighQuality = new TProfile( "NeighQuality", "Neighbours Quality vs row", 160, 0., 160. );
 412   fhNeighEff = new TProfile( "NeighEff", "Neighbours Efficiency vs row", 160, 0., 160. );
 413   fhNeighQualityVsPt = new TProfile( "NeighQualityVsPt", "Neighbours Quality vs Pt", 100, 0., 5. );
 414   fhNeighEffVsPt = new TProfile( "NeighEffVsPt", "Neighbours Efficiency vs Pt", 100, 0., 5. );
 415   fhNeighDy = new TH1D( "NeighDy", "Neighbours dy", 100, -10, 10 );
 416   fhNeighDz =  new TH1D( "NeighDz", "Neighbours dz", 100, -10, 10 );
 417   fhNeighChi = new TH1D( "NeighChi", "Neighbours chi", 100, 0, 20 );
 418
 419   fhNeighDyVsPt = new TH2D( "NeighDyVsPt", "NeighDyVsPt", 100, 0, 5, 100, -20, 20 );
 420   fhNeighDzVsPt = new TH2D( "NeighDzVsPt", "NeighDzVsPt", 100, 0, 5, 100, -20, 20 );
 421   fhNeighChiVsPt = new TH2D( "NeighChiVsPt", "NeighChiVsPt", 100, 0, 5, 100, 0, 40 );
 422   fhNeighNCombVsArea = new TH2D( "NeighNCombVsArea", "NeighNCombVsArea", 15, 0, 3, 40, 0, 40 );
 423
 424   gDirectory->cd( ".." );
 425
 426   gDirectory->mkdir( "Tracklets" );
 427   gDirectory->cd( "Tracklets" );
 428
 429   fhNHitsPerSeed = new TH1D( "NHitsPerSeed", "NHitsPerSeed", 160, 0, 160 );
 430   fhSeedEffVsP = new TProfile( "fhSeedEffVsP", "Track Seed Eff vs P", 100, 0., 5. );
 431
 432   gDirectory->cd( ".." );
 433
 434   gDirectory->mkdir( "TrackCandidates" );
 435   gDirectory->cd( "TrackCandidates" );
 436
 437   fhNHitsPerTrackCand = new TH1D( "NHitsPerTrackCand", "NHitsPerTrackCand", 160, 0, 160 );
 438   fhCandEffVsP = new TProfile( "fhCandEffVsP", "Track Candidate Eff vs P", 100, 0., 5. );
 439
 440   gDirectory->cd( ".." );
 441
 442   gDirectory->mkdir( "Tracks" );
 443   gDirectory->cd( "Tracks" );
 444
 445   fhTrackLengthRef = new TH1D( "TrackLengthRef", "TrackLengthRef", 100, 0, 1 );
 446
 447   fhRefRecoX = new TH1D( "fhRefRecoX", "fhRefRecoX", 100, 0, 200. );
 448   fhRefRecoY = new TH1D( "fhRefRecoY", "fhRefRecoY", 100, -200, 200. );
 449   fhRefRecoZ = new TH1D( "fhRefRecoZ", "fhRefRecoZ", 100, -250, 250. );
 450
 451
 452   fhRefRecoP = new TH1D( "fhRefRecoP", "fhRefRecoP", 100, 0, 10. );
 453   fhRefRecoPt = new TH1D( "fhRefRecoPt", "fhRefRecoPt", 100, 0, 10. );
 454   fhRefRecoAngleY = new TH1D( "fhRefRecoAngleY", "fhRefRecoAngleY", 100, -180., 180. );
 455   fhRefRecoAngleZ = new TH1D( "fhRefRecoAngleZ", "fhRefRecoAngleZ", 100, -180., 180 );
 456   fhRefRecoNHits = new TH1D( "fhRefRecoNHits", "fhRefRecoNHits", 100, 0., 200 );
 457
 458   fhRefNotRecoX = new TH1D( "fhRefNotRecoX", "fhRefNotRecoX", 100, 0, 200. );
 459   fhRefNotRecoY = new TH1D( "fhRefNotRecoY", "fhRefNotRecoY", 100, -200, 200. );
 460   fhRefNotRecoZ = new TH1D( "fhRefNotRecoZ", "fhRefNotRecoZ", 100, -250, 250. );
 461
 462
 463   fhRefNotRecoP = new TH1D( "fhRefNotRecoP", "fhRefNotRecoP", 100, 0, 10. );
 464   fhRefNotRecoPt = new TH1D( "fhRefNotRecoPt", "fhRefNotRecoPt", 100, 0, 10. );
 465   fhRefNotRecoAngleY = new TH1D( "fhRefNotRecoAngleY", "fhRefNotRecoAngleY", 100, -180., 180. );
 466   fhRefNotRecoAngleZ = new TH1D( "fhRefNotRecoAngleZ", "fhRefNotRecoAngleZ", 100, -180., 180 );
 467   fhRefNotRecoNHits = new TH1D( "fhRefNotRecoNHits", "fhRefNotRecoNHits", 100, 0., 200 );
 468
 469   gDirectory->cd( ".." );
 470
 471   gDirectory->mkdir( "TrackFit" );
 472   gDirectory->cd( "TrackFit" );
 473
 474   fhResY = new TH1D( "resY", "track Y resoltion [cm]", 30, -.5, .5 );
 475   fhResZ = new TH1D( "resZ", "track Z resoltion [cm]", 30, -.5, .5 );
 476   fhResSinPhi = new TH1D( "resSinPhi", "track SinPhi resoltion ", 30, -.03, .03 );
 477   fhResDzDs = new TH1D( "resDzDs", "track DzDs resoltion ", 30, -.01, .01 );
 478   fhResPt = new TH1D( "resPt", "track relative Pt resoltion", 30, -.2, .2 );
 479   fhPullY = new TH1D( "pullY", "track Y pull", 30, -10., 10. );
 480   fhPullZ = new TH1D( "pullZ", "track Z pull", 30, -10., 10. );
 481   fhPullSinPhi = new TH1D( "pullSinPhi", "track SinPhi pull", 30, -10., 10. );
 482   fhPullDzDs = new TH1D( "pullDzDs", "track DzDs pull", 30, -10., 10. );
 483   fhPullQPt = new TH1D( "pullQPt", "track Q/Pt pull", 30, -10., 10. );
 484   fhPullYS = new TH1D( "pullYS", "track Y+SinPhi chi deviation", 100, 0., 30. );
 485   fhPullZT = new TH1D( "pullZT", "track Z+DzDs chi deviation ", 100, 0., 30. );
 486
 487   gDirectory->cd( ".." );
 488
 489   fhEffVsP = new TProfile( "EffVsP", "Eff vs P", 100, 0., 5. );
 490   fhGBEffVsP = new TProfile( "GBEffVsP", "Global tracker: Eff vs P", 100, 0., 5. );
 491   fhGBEffVsPt = new TProfile( "GBEffVsPt", "Global tracker: Eff vs Pt", 100, 0.2, 5. );
 492
 493   gDirectory->mkdir( "Clusters" );
 494   gDirectory->cd( "Clusters" );
 495
 496   fhHitShared = new TProfile( "fhHitSharedf", "fhHitShared vs row", 160, 0., 160. );
 497
 498   fhHitResY = new TH1D( "resHitY", "Y cluster resoltion [cm]", 100, -2., 2. );
 499   fhHitResZ = new TH1D( "resHitZ", "Z cluster resoltion [cm]", 100, -2., 2. );
 500   fhHitPullY = new TH1D( "pullHitY", "Y cluster pull", 100, -10., 10. );
 501   fhHitPullZ = new TH1D( "pullHitZ", "Z cluster pull", 100, -10., 10. );
 502
 503   fhHitResY1 = new TH1D( "resHitY1", "Y cluster resoltion [cm]", 100, -2., 2. );
 504   fhHitResZ1 = new TH1D( "resHitZ1", "Z cluster resoltion [cm]", 100, -2., 2. );
 505   fhHitPullY1 = new TH1D( "pullHitY1", "Y cluster pull", 100, -10., 10. );
 506   fhHitPullZ1 = new TH1D( "pullHitZ1", "Z cluster pull", 100, -10., 10. );
 507
 508   fhHitErrY = new TH1D( "HitErrY", "Y cluster error [cm]", 100, 0., 3. );
 509   fhHitErrZ = new TH1D( "HitErrZ", "Z cluster error [cm]", 100, 0., 3. );
 510
 511   gDirectory->cd( ".." );
 512
 513   gDirectory->mkdir( "Cells" );
 514   gDirectory->cd( "Cells" );
 515   fhCellPurity = new TH1D( "CellPurity", "Cell Purity", 100, -0.1, 1.1 );
 516   fhCellNHits = new TH1D( "CellNHits", "Cell NHits", 40, 0., 40. );
 517   fhCellPurityVsN = new TProfile( "CellPurityVsN", "Cell purity Vs N hits", 40, 2., 42. );
 518   fhCellPurityVsPt = new TProfile( "CellPurityVsPt", "Cell purity Vs Pt", 100, 0., 5. );
 519   gDirectory->cd( ".." );
 520
 521   curdir->cd();
 522 }
 523
 524 void AliHLTTPCCAPerformance::WriteDir2Current( TObject *obj )
 525 {
 526   //* recursive function to copy the directory 'obj' to the current one
 527   if ( !obj->IsFolder() ) obj->Write();
 528   else {
 529     TDirectory *cur = gDirectory;
 530     TDirectory *sub = cur->mkdir( obj->GetName() );
 531     sub->cd();
 532     TList *listSub = ( ( TDirectory* )obj )->GetList();
 533     TIter it( listSub );
 534     while ( TObject *obj1 = it() ) WriteDir2Current( obj1 );
 535     cur->cd();
 536   }
 537 }
 538
 539 void AliHLTTPCCAPerformance::WriteHistos()
 540 {
 541   //* write histograms to the file
 542   TDirectory *curr = gDirectory;
 543   // Open output file and write histograms
 544   TFile* outfile = new TFile( "HLTTPCCATrackerPerformance.root", "RECREATE" );
 545   outfile->cd();
 546   WriteDir2Current( fHistoDir );
 547   outfile->Close();
 548   curr->cd();
 549 }
 550
 551
 552
 553
 554 void AliHLTTPCCAPerformance::GetMCLabel( std::vector<int> &ClusterIDs, int &Label, float &Purity )
 555 {
 556   // find MC label for the track
 557
 558   Label = -1;
 559   Purity = 0;
 560   int nClusters = ClusterIDs.size();
 561   vector<int> labels;
 562   for ( int i = 0; i < nClusters; i++ ) {
 563     const AliHLTTPCCAHitLabel &l = fHitLabels[ClusterIDs[i]];
 564     if ( l.fLab[0] >= 0 ) labels.push_back( l.fLab[0] );
 565     if ( l.fLab[1] >= 0 ) labels.push_back( l.fLab[1] );
 566     if ( l.fLab[2] >= 0 ) labels.push_back( l.fLab[2] );
 567   }
 568   sort( labels.begin(), labels.end() );
 569   int nMax = 0, labCur = -1, nCur = 0;
 570
 571   for ( unsigned int i = 0; i < labels.size(); i++ ) {
 572     if ( labels[i] != labCur ) {
 573       if ( nMax < nCur ) {
 574         nMax = nCur;
 575         Label = labCur;
 576       }
 577       labCur = labels[i];
 578       nCur = 0;
 579     }
 580     nCur++;
 581   }
 582   if ( nMax < nCur ) Label = labCur;
 583
 584   nMax = 0;
 585   for ( int i = 0; i < nClusters; i++ ) {
 586     const AliHLTTPCCAHitLabel &l = fHitLabels[ClusterIDs[i]];
 587     if ( l.fLab[0] == Label || l.fLab[1] == Label || l.fLab[2] == Label ) nMax++;
 588   }
 589   Purity = ( nClusters > 0 ) ? ( ( double ) nMax ) / nClusters : 0 ;
 590 }
 591
 592
 593 void AliHLTTPCCAPerformance::LinkPerformance( int /*iSlice*/ )
 594 {
 595   // Efficiency and quality of the found neighbours
 596 #ifdef XXX
 597   std::cout << "Link performance..." << std::endl;
 598   if ( !fTracker ) return;
 599   const AliHLTTPCCATracker &slice = fTracker->Slice( iSlice );
 600
 601   AliHLTResizableArray<int> mcType( fNMCTracks );
 602
 603   for ( int imc = 0; imc < fNMCTracks; imc++ ) {
 604     if ( fMCTracks[imc].P() < .2 ) {  mcType[imc] = -1; continue; }
 605     float x = fMCTracks[imc].Par()[0];
 606     float y = fMCTracks[imc].Par()[1];
 607     //float z = fMCTracks[imc].Par()[2];
 608     if ( x*x + y*y < 100. ) {
 609       if ( fMCTracks[imc].P() >= 1 ) mcType[imc] = 0;
 610       else mcType[imc] = 1;
 611     } else {
 612       if ( fMCTracks[imc].P() >= 1 ) mcType[imc] = 2;
 613       else mcType[imc] = 3;
 614     }
 615   }
 616
 617   struct AliHLTTPCCAMCHits {
 618     int fNHits;
 619     int fID[30];
 620   };
 621   AliHLTTPCCAMCHits *mcGbHitsUp = new AliHLTTPCCAMCHits[fNMCTracks];
 622   AliHLTTPCCAMCHits *mcGbHitsDn = new AliHLTTPCCAMCHits[fNMCTracks];
 623
 624   for ( int iRow = 2; iRow < slice.Param().NRows() - 2; iRow++ ) {
 625
 626     const AliHLTTPCCARow &row = slice.Row( iRow );
 627     const AliHLTTPCCARow &rowUp = slice.Row( iRow + 2 );
 628     const AliHLTTPCCARow &rowDn = slice.Row( iRow - 2 );
 629
 630     AliHLTResizableArray<int> gbHits  ( row.NHits() );
 631     AliHLTResizableArray<int> gbHitsUp( rowUp.NHits() );
 632     AliHLTResizableArray<int> gbHitsDn( rowDn.NHits() );
 633
 634     for ( int ih = 0; ih < row.NHits()  ; ih++ ) gbHits  [ih] = fTracker->FirstSliceHit()[iSlice] + slice.HitInputID( row  , ih );
 635     for ( int ih = 0; ih < rowUp.NHits(); ih++ ) gbHitsUp[ih] = fTracker->FirstSliceHit()[iSlice] + slice.HitInputID( rowUp, ih );
 636     for ( int ih = 0; ih < rowDn.NHits(); ih++ ) gbHitsDn[ih] = fTracker->FirstSliceHit()[iSlice] + slice.HitInputID( rowDn, ih );
 637
 638     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
 639       mcGbHitsUp[imc].fNHits = 0;
 640       mcGbHitsDn[imc].fNHits = 0;
 641     }
 642
 643     for ( int ih = 0; ih < rowUp.NHits(); ih++ ) {
 644       AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[gbHitsUp[ih]].ID()];
 645       for ( int il = 0; il < 3; il++ ) {
 646         int imc = l.fLab[il];
 647         if ( imc < 0 ) break;
 648         int &nmc = mcGbHitsUp[imc].fNHits;
 649         if ( nmc >= 30 ) continue;
 650         mcGbHitsUp[imc].fID[nmc] = gbHitsUp[ih];
 651         nmc++;
 652       }
 653     }
 654
 655     for ( int ih = 0; ih < rowDn.NHits(); ih++ ) {
 656       AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[gbHitsDn[ih]].ID()];
 657       for ( int il = 0; il < 3; il++ ) {
 658         int imc = l.fLab[il];
 659         if ( imc < 0 ) break;
 660         int &nmc = mcGbHitsDn[imc].fNHits;
 661         if ( nmc >= 30 ) continue;
 662         mcGbHitsDn[imc].fID[nmc] = gbHitsDn[ih];
 663         nmc++;
 664       }
 665     }
 666
 667     //float dxUp = rowUp.X() - row.X();
 668     //float dxDn = row.X() - rowDn.X();
 669     float tUp = rowUp.X() / row.X();
 670     float tDn = rowDn.X() / row.X();
 671
 672     for ( int ih = 0; ih < row.NHits(); ih++ ) {
 673
 674       int up = slice.HitLinkUpData( row, ih );
 675       int dn = slice.HitLinkDownData( row, ih );
 676
 677       const AliHLTTPCCAGBHit &h = fTracker->Hits()[gbHits[ih]];
 678       AliHLTTPCCAHitLabel &l = fHitLabels[h.ID()];
 679
 680       int isMC = -1;
 681       int mcFound = -1;
 682
 683       float yUp = h.Y() * tUp, zUp = h.Z() * tUp;
 684       float yDn = h.Y() * tDn, zDn = h.Z() * tDn;
 685
 686       for ( int il = 0; il < 3; il++ ) {
 687         int imc = l.fLab[il];
 688         if ( imc < 0 ) break;
 689
 690         bool isMcUp = 0, isMcDn = 0;
 691
 692         float dyMin = 1.e8, dzMin = 1.e8;
 693         for ( int i = 0; i < mcGbHitsUp[imc].fNHits; i++ ) {
 694           const AliHLTTPCCAGBHit &h1 = fTracker->Hits()[mcGbHitsUp[imc].fID[i]];
 695           float dy = TMath::Abs( h1.Y() - yUp );
 696           float dz = TMath::Abs( h1.Z() - zUp );
 697           if ( dy*dy + dz*dz < dyMin*dyMin + dzMin*dzMin ) {
 698             dyMin = dy;
 699             dzMin = dz;
 700           }
 701         }
 702
 703         if ( mcType[imc] >= 0 && mcGbHitsUp[imc].fNHits >= 0 ) {
 704           fhLinkAreaY[mcType[imc]]->Fill( dyMin );
 705           fhLinkAreaZ[mcType[imc]]->Fill( dzMin );
 706         }
 707         if ( dyMin*dyMin + dzMin*dzMin < 100. ) isMcUp = 1;
 708
 709         dyMin = 1.e8;
 710         dzMin = 1.e8;
 711         for ( int i = 0; i < mcGbHitsDn[imc].fNHits; i++ ) {
 712           const AliHLTTPCCAGBHit &h1 = fTracker->Hits()[mcGbHitsDn[imc].fID[i]];
 713           float dy = TMath::Abs( h1.Y() - yDn );
 714           float dz = TMath::Abs( h1.Z() - zDn );
 715           if ( dy*dy + dz*dz < dyMin*dyMin + dzMin*dzMin ) {
 716             dyMin = dy;
 717             dzMin = dz;
 718           }
 719         }
 720
 721         if ( mcType[imc] >= 0 && mcGbHitsDn[imc].fNHits >= 0 ) {
 722           fhLinkAreaY[mcType[imc]]->Fill( dyMin );
 723           fhLinkAreaZ[mcType[imc]]->Fill( dzMin );
 724         }
 725         if ( dyMin*dyMin + dzMin*dzMin < 100. ) isMcDn = 1;
 726
 727         if ( !isMcUp || !isMcDn ) continue;
 728         isMC = imc;
 729
 730         bool found = 0;
 731         if ( up >= 0 && dn >= 0 ) {
 732           //std::cout<<"row, ih, mc, up, dn = "<<iRow<<" "<<ih<<" "<<imc<<" "<<up<<" "<<dn<<std::endl;
 733           const AliHLTTPCCAGBHit &hUp = fTracker->Hits()[gbHitsUp[up]];
 734           const AliHLTTPCCAGBHit &hDn = fTracker->Hits()[gbHitsDn[dn]];
 735           AliHLTTPCCAHitLabel &lUp = fHitLabels[hUp.ID()];
 736           AliHLTTPCCAHitLabel &lDn = fHitLabels[hDn.ID()];
 737           bool foundUp = 0, foundDn = 0;
 738           for ( int jl = 0; jl < 3; jl++ ) {
 739             if ( lUp.fLab[jl] == imc ) foundUp = 1;
 740             if ( lDn.fLab[jl] == imc ) foundDn = 1;
 741             //std::cout<<"mc up, dn = "<<lUp.fLab[jl]<<" "<<lDn.fLab[jl]<<std::endl;
 742           }
 743           if ( foundUp && foundDn ) found = 1;
 744         }
 745         if ( found ) { mcFound = imc; break;}
 746       }
 747
 748       if ( mcFound >= 0 ) {
 749         //std::cout<<" mc "<<mcFound<<" found"<<std::endl;
 750         if ( mcType[mcFound] >= 0 ) fhLinkEff[mcType[mcFound]]->Fill( iRow, 1 );
 751       } else if ( isMC >= 0 ) {
 752         //std::cout<<" mc "<<isMC<<" not found"<<std::endl;
 753         if ( mcType[isMC] >= 0 ) fhLinkEff[mcType[isMC]]->Fill( iRow, 0 );
 754       }
 755
 756     } // ih
 757   } // iRow
 758   delete[] mcGbHitsUp;
 759   delete[] mcGbHitsDn;
 760 #endif
 761 }
 762
 763
 764 void AliHLTTPCCAPerformance::SliceTrackletPerformance( int /*iSlice*/, bool /*PrintFlag*/ )
 765 {
 766   //* calculate slice tracker performance
 767 #ifdef XXX
 768   if ( !fTracker ) return;
 769
 770   int nRecTot = 0, nGhost = 0, nRecOut = 0;
 771   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
 772   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
 773   const AliHLTTPCCATracker &slice = fTracker->Slice( iSlice );
 774
 775   int firstSliceHit = fTracker->FirstSliceHit()[iSlice];
 776   int endSliceHit = fTracker->NHits();
 777   if ( iSlice < fTracker->NSlices() - 1 ) endSliceHit = fTracker->FirstSliceHit()[iSlice+1];
 778
 779   // Select reconstructable MC tracks
 780
 781   {
 782     for ( int imc = 0; imc < fNMCTracks; imc++ ) fMCTracks[imc].SetNHits( 0 );
 783
 784     for ( int ih = firstSliceHit; ih < endSliceHit; ih++ ) {
 785       int id = fTracker->Hits()[ih].ID();
 786       if ( id < 0 || id >= fNHits ) break;
 787       AliHLTTPCCAHitLabel &l = fHitLabels[id];
 788       if ( l.fLab[0] >= 0 ) fMCTracks[l.fLab[0]].SetNHits( fMCTracks[l.fLab[0]].NHits() + 1 );
 789       if ( l.fLab[1] >= 0 ) fMCTracks[l.fLab[1]].SetNHits( fMCTracks[l.fLab[1]].NHits() + 1 );
 790       if ( l.fLab[2] >= 0 ) fMCTracks[l.fLab[2]].SetNHits( fMCTracks[l.fLab[2]].NHits() + 1 );
 791     }
 792
 793     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
 794       AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
 795       mc.SetSet( 0 );
 796       mc.SetNReconstructed( 0 );
 797       mc.SetNTurns( 1 );
 798       if ( mc.NHits() >=  30 && mc.P() >= .05 ) {
 799         mc.SetSet( 1 );
 800         nMCAll++;
 801         if ( mc.NHits() >=  30 && mc.P() >= 1. ) {
 802           mc.SetSet( 2 );
 803           nMCRef++;
 804         }
 805       }
 806     }
 807   }
 808
 809
 810   int traN = slice.NTracklets();
 811   int *traLabels = 0;
 812   double *traPurity = 0;
 813
 814   traLabels = new int[traN];
 815   traPurity = new double[traN];
 816   {
 817     for ( int itr = 0; itr < traN; itr++ ) {
 818       traLabels[itr] = -1;
 819       traPurity[itr] = 0;
 820
 821       int hits[1600];
 822       int nHits = 0;
 823
 824       {
 825         const AliHLTTPCCAHitId &id = slice.TrackletStartHit( itr );
 826         int iRow = id.RowIndex();
 827         int ih =  id.HitIndex();
 828
 829         while ( ih >= 0 ) {
 830           const AliHLTTPCCARow &row = slice.Row( iRow );
 831           hits[nHits] = firstSliceHit + slice.HitInputID( row, ih );
 832           nHits++;
 833           ih = slice.HitLinkUpData( row, ih );
 834           iRow++;
 835         }
 836       }
 837
 838       if ( nHits < 5 ) continue;
 839
 840       int lb[1600*3];
 841       int nla = 0;
 842
 843       for ( int ih = 0; ih < nHits; ih++ ) {
 844         AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[hits[ih]].ID()];
 845         if ( l.fLab[0] >= 0 ) lb[nla++] = l.fLab[0];
 846         if ( l.fLab[1] >= 0 ) lb[nla++] = l.fLab[1];
 847         if ( l.fLab[2] >= 0 ) lb[nla++] = l.fLab[2];
 848       }
 849
 850       sort( lb, lb + nla );
 851       int labmax = -1, labcur = -1, lmax = 0, lcurr = 0;
 852       for ( int i = 0; i < nla; i++ ) {
 853         if ( lb[i] != labcur ) {
 854           if ( labcur >= 0 && lmax < lcurr ) {
 855             lmax = lcurr;
 856             labmax = labcur;
 857           }
 858           labcur = lb[i];
 859           lcurr = 0;
 860         }
 861         lcurr++;
 862       }
 863       if ( labcur >= 0 && lmax < lcurr ) {
 864         lmax = lcurr;
 865         labmax = labcur;
 866       }
 867       lmax = 0;
 868       for ( int ih = 0; ih < nHits; ih++ ) {
 869         AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[hits[ih]].ID()];
 870         if ( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax
 871            ) lmax++;
 872       }
 873       traLabels[itr] = labmax;
 874       traPurity[itr] = ( ( nHits > 0 ) ? double( lmax ) / double( nHits ) : 0 );
 875     }
 876   }
 877
 878   nRecTot += traN;
 879
 880   for ( int itr = 0; itr < traN; itr++ ) {
 881     if ( traPurity[itr] < .9 || traLabels[itr] < 0 || traLabels[itr] >= fNMCTracks ) {
 882       nGhost++;
 883       continue;
 884     }
 885
 886     AliHLTTPCCAMCTrack &mc = fMCTracks[traLabels[itr]];
 887     mc.SetNReconstructed( mc.NReconstructed() + 1 );
 888     if ( mc.Set() == 0 ) nRecOut++;
 889     else {
 890       if ( mc.NReconstructed() == 1 ) nRecAll++;
 891       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
 892       if ( mc.Set() == 2 ) {
 893         if ( mc.NReconstructed() == 1 ) nRecRef++;
 894         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
 895       }
 896     }
 897   }
 898
 899   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
 900     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
 901     if ( mc.Set() > 0 ) fhSeedEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
 902   }
 903
 904   if ( traLabels ) delete[] traLabels;
 905   if ( traPurity ) delete[] traPurity;
 906
 907   fStatSeedNRecTot += nRecTot;
 908   fStatSeedNRecOut += nRecOut;
 909   fStatSeedNGhost  += nGhost;
 910   fStatSeedNMCAll  += nMCAll;
 911   fStatSeedNRecAll  += nRecAll;
 912   fStatSeedNClonesAll  += nClonesAll;
 913   fStatSeedNMCRef  += nMCRef;
 914   fStatSeedNRecRef  += nRecRef;
 915   fStatSeedNClonesRef  += nClonesRef;
 916
 917   if ( nMCAll == 0 ) return;
 918
 919   if ( PrintFlag ) {
 920     cout << "Track seed performance for slice " << iSlice << " : " << endl;
 921     cout << " N tracks : "
 922          << nMCAll << " mc all, "
 923          << nMCRef << " mc ref, "
 924          << nRecTot << " rec total, "
 925          << nRecAll << " rec all, "
 926          << nClonesAll << " clones all, "
 927          << nRecRef << " rec ref, "
 928          << nClonesRef << " clones ref, "
 929          << nRecOut << " out, "
 930          << nGhost << " ghost" << endl;
 931
 932     int nRecExtr = nRecAll - nRecRef;
 933     int nMCExtr = nMCAll - nMCRef;
 934     int nClonesExtr = nClonesAll - nClonesRef;
 935
 936     double dRecTot = ( nRecTot > 0 ) ? nRecTot : 1;
 937     double dMCAll = ( nMCAll > 0 ) ? nMCAll : 1;
 938     double dMCRef = ( nMCRef > 0 ) ? nMCRef : 1;
 939     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
 940     double dRecAll = ( nRecAll + nClonesAll > 0 ) ? nRecAll + nClonesAll : 1;
 941     double dRecRef = ( nRecRef + nClonesRef > 0 ) ? nRecRef + nClonesRef : 1;
 942     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
 943
 944     cout << " EffRef = ";
 945     if ( nMCRef > 0 ) cout << nRecRef / dMCRef; else cout << "_";
 946     cout << ", CloneRef = ";
 947     if ( nRecRef > 0 ) cout << nClonesRef / dRecRef; else cout << "_";
 948     cout << endl;
 949     cout << " EffExtra = ";
 950     if ( nMCExtr > 0 ) cout << nRecExtr / dMCExtr; else cout << "_";
 951     cout << ", CloneExtra = ";
 952     if ( nRecExtr > 0 ) cout << nClonesExtr / dRecExtr; else cout << "_";
 953     cout << endl;
 954     cout << " EffAll = ";
 955     if ( nMCAll > 0 ) cout << nRecAll / dMCAll; else cout << "_";
 956     cout << ", CloneAll = ";
 957     if ( nRecAll > 0 ) cout << nClonesAll / dRecAll; else cout << "_";
 958     cout << endl;
 959     cout << " Out = ";
 960     if ( nRecTot > 0 ) cout << nRecOut / dRecTot; else cout << "_";
 961     cout << ", Ghost = ";
 962     if ( nRecTot > 0 ) cout << nGhost / dRecTot; else cout << "_";
 963     cout << endl;
 964   }
 965 #endif
 966 }
 967
 968
 969
 970
 971 void AliHLTTPCCAPerformance::SliceTrackCandPerformance( int /*iSlice*/, bool /*PrintFlag*/ )
 972 {
 973   //* calculate slice tracker performance
 974 #ifdef XXX
 975   if ( !fTracker ) return;
 976
 977   int nRecTot = 0, nGhost = 0, nRecOut = 0;
 978   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
 979   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
 980   const AliHLTTPCCATracker &slice = fTracker->Slice( iSlice );
 981
 982   int firstSliceHit = fTracker->FirstSliceHit()[iSlice];
 983   int endSliceHit = fTracker->NHits();
 984   if ( iSlice < fTracker->NSlices() - 1 ) endSliceHit = fTracker->FirstSliceHit()[iSlice+1];
 985
 986   // Select reconstructable MC tracks
 987
 988   {
 989     for ( int imc = 0; imc < fNMCTracks; imc++ ) fMCTracks[imc].SetNHits( 0 );
 990
 991     for ( int ih = firstSliceHit; ih < endSliceHit; ih++ ) {
 992       int id = fTracker->Hits()[ih].ID();
 993       if ( id < 0 || id >= fNHits ) break;
 994       AliHLTTPCCAHitLabel &l = fHitLabels[id];
 995       if ( l.fLab[0] >= 0 ) fMCTracks[l.fLab[0]].SetNHits( fMCTracks[l.fLab[0]].NHits() + 1 );
 996       if ( l.fLab[1] >= 0 ) fMCTracks[l.fLab[1]].SetNHits( fMCTracks[l.fLab[1]].NHits() + 1 );
 997       if ( l.fLab[2] >= 0 ) fMCTracks[l.fLab[2]].SetNHits( fMCTracks[l.fLab[2]].NHits() + 1 );
 998     }
 999
1000     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
1001       AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
1002       mc.SetSet( 0 );
1003       mc.SetNReconstructed( 0 );
1004       mc.SetNTurns( 1 );
1005       if ( mc.NHits() >=  30 && mc.P() >= .05 ) {
1006         mc.SetSet( 1 );
1007         nMCAll++;
1008         if ( mc.NHits() >=  30 && mc.P() >= 1. ) {
1009           mc.SetSet( 2 );
1010           nMCRef++;
1011         }
1012       }
1013     }
1014   }
1015
1016   int traN = slice.NTracklets();
1017   int *traLabels = 0;
1018   double *traPurity = 0;
1019   traLabels = new int[traN];
1020   traPurity = new double[traN];
1021   {
1022     for ( int itr = 0; itr < traN; itr++ ) {
1023       traLabels[itr] = -1;
1024       traPurity[itr] = 0;
1025
1026       const AliHLTTPCCATracklet &t = slice.Tracklet( itr );
1027
1028       int nHits = t.NHits();
1029       if ( nHits < 10 ) continue;
1030       int firstRow = t.FirstRow();
1031       int lastRow = t.LastRow();
1032       nHits = 0;
1033
1034       int lb[1600*3];
1035       int nla = 0;
1036
1037       for ( int irow = firstRow; irow <= lastRow; irow++ ) {
1038 #ifdef EXTERN_ROW_HITS
1039         int ih = slice.TrackletRowHits[iRow * *slice.NTracklets() + itr];
1040 #else
1041                 int ih = t.RowHit( irow );
1042 #endif
1043         if ( ih < 0 ) continue;
1044         int index = firstSliceHit + slice.HitInputID( slice.Row( irow ), ih );
1045         AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
1046         if ( l.fLab[0] >= 0 ) lb[nla++] = l.fLab[0];
1047         if ( l.fLab[1] >= 0 ) lb[nla++] = l.fLab[1];
1048         if ( l.fLab[2] >= 0 ) lb[nla++] = l.fLab[2];
1049         nHits++;
1050       }
1051       if ( nHits < 10 ) continue;
1052
1053       sort( lb, lb + nla );
1054       int labmax = -1, labcur = -1, lmax = 0, lcurr = 0;
1055       for ( int i = 0; i < nla; i++ ) {
1056         if ( lb[i] != labcur ) {
1057           if ( labcur >= 0 && lmax < lcurr ) {
1058             lmax = lcurr;
1059             labmax = labcur;
1060           }
1061           labcur = lb[i];
1062           lcurr = 0;
1063         }
1064         lcurr++;
1065       }
1066       if ( labcur >= 0 && lmax < lcurr ) {
1067         lmax = lcurr;
1068         labmax = labcur;
1069       }
1070       lmax = 0;
1071       for ( int irow = firstRow; irow <= lastRow; irow++ ) {
1072 #ifdef EXTERN_ROW_HITS
1073         int ih = slice.TrackletRowHits[iRow * *slice.NTracklets() + itr];
1074 #else
1075                 int ih = t.RowHit( irow );
1076 #endif
1077         if ( ih < 0 ) continue;
1078         int index = firstSliceHit + slice.HitInputID( slice.Row( irow ), ih );
1079         AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
1080         if ( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax
1081            ) lmax++;
1082       }
1083       traLabels[itr] = labmax;
1084       traPurity[itr] = ( ( nHits > 0 ) ? double( lmax ) / double( nHits ) : 0 );
1085     }
1086   }
1087
1088   nRecTot += traN;
1089
1090   for ( int itr = 0; itr < traN; itr++ ) {
1091     if ( traPurity[itr] < .9 || traLabels[itr] < 0 || traLabels[itr] >= fNMCTracks ) {
1092       nGhost++;
1093       continue;
1094     }
1095
1096     AliHLTTPCCAMCTrack &mc = fMCTracks[traLabels[itr]];
1097     mc.SetNReconstructed( mc.NReconstructed() + 1 );
1098     if ( mc.Set() == 0 ) nRecOut++;
1099     else {
1100       if ( mc.NReconstructed() == 1 ) nRecAll++;
1101       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
1102       if ( mc.Set() == 2 ) {
1103         if ( mc.NReconstructed() == 1 ) nRecRef++;
1104         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
1105       }
1106     }
1107   }
1108
1109   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1110     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1111     if ( mc.Set() > 0 ) fhCandEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1112   }
1113
1114   if ( traLabels ) delete[] traLabels;
1115   if ( traPurity ) delete[] traPurity;
1116
1117   fStatCandNRecTot += nRecTot;
1118   fStatCandNRecOut += nRecOut;
1119   fStatCandNGhost  += nGhost;
1120   fStatCandNMCAll  += nMCAll;
1121   fStatCandNRecAll  += nRecAll;
1122   fStatCandNClonesAll  += nClonesAll;
1123   fStatCandNMCRef  += nMCRef;
1124   fStatCandNRecRef  += nRecRef;
1125   fStatCandNClonesRef  += nClonesRef;
1126
1127   if ( nMCAll == 0 ) return;
1128
1129   if ( PrintFlag ) {
1130     cout << "Track candidate performance for slice " << iSlice << " : " << endl;
1131     cout << " N tracks : "
1132          << nMCAll << " mc all, "
1133          << nMCRef << " mc ref, "
1134          << nRecTot << " rec total, "
1135          << nRecAll << " rec all, "
1136          << nClonesAll << " clones all, "
1137          << nRecRef << " rec ref, "
1138          << nClonesRef << " clones ref, "
1139          << nRecOut << " out, "
1140          << nGhost << " ghost" << endl;
1141
1142     int nRecExtr = nRecAll - nRecRef;
1143     int nMCExtr = nMCAll - nMCRef;
1144     int nClonesExtr = nClonesAll - nClonesRef;
1145
1146     double dRecTot = ( nRecTot > 0 ) ? nRecTot : 1;
1147     double dMCAll = ( nMCAll > 0 ) ? nMCAll : 1;
1148     double dMCRef = ( nMCRef > 0 ) ? nMCRef : 1;
1149     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1150     double dRecAll = ( nRecAll + nClonesAll > 0 ) ? nRecAll + nClonesAll : 1;
1151     double dRecRef = ( nRecRef + nClonesRef > 0 ) ? nRecRef + nClonesRef : 1;
1152     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1153
1154     cout << " EffRef = ";
1155     if ( nMCRef > 0 ) cout << nRecRef / dMCRef; else cout << "_";
1156     cout << ", CloneRef = ";
1157     if ( nRecRef > 0 ) cout << nClonesRef / dRecRef; else cout << "_";
1158     cout << endl;
1159     cout << " EffExtra = ";
1160     if ( nMCExtr > 0 ) cout << nRecExtr / dMCExtr; else cout << "_";
1161     cout << ", CloneExtra = ";
1162     if ( nRecExtr > 0 ) cout << nClonesExtr / dRecExtr; else cout << "_";
1163     cout << endl;
1164     cout << " EffAll = ";
1165     if ( nMCAll > 0 ) cout << nRecAll / dMCAll; else cout << "_";
1166     cout << ", CloneAll = ";
1167     if ( nRecAll > 0 ) cout << nClonesAll / dRecAll; else cout << "_";
1168     cout << endl;
1169     cout << " Out = ";
1170     if ( nRecTot > 0 ) cout << nRecOut / dRecTot; else cout << "_";
1171     cout << ", Ghost = ";
1172     if ( nRecTot > 0 ) cout << nGhost / dRecTot; else cout << "_";
1173     cout << endl;
1174   }
1175 #endif
1176 }
1177
1178
1179
1180 void AliHLTTPCCAPerformance::SlicePerformance( int iSlice, bool PrintFlag )
1181 {
1182   //* calculate slice tracker performance
1183
1184   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1185
1186   int nRecTot = 0, nGhost = 0, nRecOut = 0;
1187   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
1188   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
1189   //const AliHLTTPCCATracker &tracker = hlt.SliceTracker( iSlice );
1190   const AliHLTTPCCAClusterData &clusterdata = hlt.ClusterData(iSlice);
1191
1192   // Select reconstructable MC tracks
1193
1194   {
1195     for ( int imc = 0; imc < fNMCTracks; imc++ ) fMCTracks[imc].SetNHits( 0 );
1196
1197     for ( int ih = 0; ih < clusterdata.NumberOfClusters(); ih++ ) {
1198       int id = clusterdata.Id( ih );
1199       if ( id < 0 || id > fNHits ) break;
1200       AliHLTTPCCAHitLabel &l = fHitLabels[id];
1201       if ( l.fLab[0] >= 0 ) fMCTracks[l.fLab[0]].SetNHits( fMCTracks[l.fLab[0]].NHits() + 1 );
1202       if ( l.fLab[1] >= 0 ) fMCTracks[l.fLab[1]].SetNHits( fMCTracks[l.fLab[1]].NHits() + 1 );
1203       if ( l.fLab[2] >= 0 ) fMCTracks[l.fLab[2]].SetNHits( fMCTracks[l.fLab[2]].NHits() + 1 );
1204     }
1205
1206     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
1207       AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
1208       mc.SetSet( 0 );
1209       mc.SetNReconstructed( 0 );
1210       mc.SetNTurns( 1 );
1211       if ( mc.NHits() >=  30 && mc.P() >= .05 ) {
1212         mc.SetSet( 1 );
1213         nMCAll++;
1214         if ( mc.NHits() >=  30 && mc.P() >= 1. ) {
1215           mc.SetSet( 2 );
1216           nMCRef++;
1217         }
1218       }
1219     }
1220   }
1221
1222   //if ( !tracker.Output() ) return;
1223
1224   const AliHLTTPCCASliceOutput &output = hlt.Output(iSlice);
1225
1226   int traN = output.NTracks();
1227
1228   nRecTot += traN;
1229
1230   for ( int itr = 0; itr < traN; itr++ ) {
1231
1232     const AliHLTTPCCASliceTrack &tCA = output.Track( itr );
1233     std::vector<int> clusterIDs;
1234     for ( int i = 0; i < tCA.NClusters(); i++ ) {
1235       clusterIDs.push_back( output.ClusterId( tCA.FirstClusterRef() + i ) );
1236     }
1237     int label;
1238     float purity;
1239     GetMCLabel( clusterIDs, label, purity );
1240
1241     if ( purity < .9 || label < 0 || label >= fNMCTracks ) {
1242       nGhost++;
1243       continue;
1244     }
1245
1246     AliHLTTPCCAMCTrack &mc = fMCTracks[label];
1247     mc.SetNReconstructed( mc.NReconstructed() + 1 );
1248     if ( mc.Set() == 0 ) nRecOut++;
1249     else {
1250       if ( mc.NReconstructed() == 1 ) nRecAll++;
1251       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
1252       if ( mc.Set() == 2 ) {
1253         if ( mc.NReconstructed() == 1 ) nRecRef++;
1254         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
1255       }
1256     }
1257
1258   }
1259
1260
1261   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1262     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1263     if ( mc.Set() > 0 ) fhEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1264   }
1265
1266
1267   fStatNRecTot += nRecTot;
1268   fStatNRecOut += nRecOut;
1269   fStatNGhost  += nGhost;
1270   fStatNMCAll  += nMCAll;
1271   fStatNRecAll  += nRecAll;
1272   fStatNClonesAll  += nClonesAll;
1273   fStatNMCRef  += nMCRef;
1274   fStatNRecRef  += nRecRef;
1275   fStatNClonesRef  += nClonesRef;
1276
1277   if ( nMCAll == 0 ) return;
1278
1279   if ( PrintFlag ) {
1280     cout << "Performance for slice " << iSlice << " : " << endl;
1281     cout << " N tracks : "
1282          << nMCAll << " mc all, "
1283          << nMCRef << " mc ref, "
1284          << nRecTot << " rec total, "
1285          << nRecAll << " rec all, "
1286          << nClonesAll << " clones all, "
1287          << nRecRef << " rec ref, "
1288          << nClonesRef << " clones ref, "
1289          << nRecOut << " out, "
1290          << nGhost << " ghost" << endl;
1291
1292     int nRecExtr = nRecAll - nRecRef;
1293     int nMCExtr = nMCAll - nMCRef;
1294     int nClonesExtr = nClonesAll - nClonesRef;
1295
1296     double dRecTot = ( nRecTot > 0 ) ? nRecTot : 1;
1297     double dMCAll = ( nMCAll > 0 ) ? nMCAll : 1;
1298     double dMCRef = ( nMCRef > 0 ) ? nMCRef : 1;
1299     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1300     double dRecAll = ( nRecAll + nClonesAll > 0 ) ? nRecAll + nClonesAll : 1;
1301     double dRecRef = ( nRecRef + nClonesRef > 0 ) ? nRecRef + nClonesRef : 1;
1302     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1303
1304     cout << " EffRef = ";
1305     if ( nMCRef > 0 ) cout << nRecRef / dMCRef; else cout << "_";
1306     cout << ", CloneRef = ";
1307     if ( nRecRef > 0 ) cout << nClonesRef / dRecRef; else cout << "_";
1308     cout << endl;
1309     cout << " EffExtra = ";
1310     if ( nMCExtr > 0 ) cout << nRecExtr / dMCExtr; else cout << "_";
1311     cout << ", CloneExtra = ";
1312     if ( nRecExtr > 0 ) cout << nClonesExtr / dRecExtr; else cout << "_";
1313     cout << endl;
1314     cout << " EffAll = ";
1315     if ( nMCAll > 0 ) cout << nRecAll / dMCAll; else cout << "_";
1316     cout << ", CloneAll = ";
1317     if ( nRecAll > 0 ) cout << nClonesAll / dRecAll; else cout << "_";
1318     cout << endl;
1319     cout << " Out = ";
1320     if ( nRecTot > 0 ) cout << nRecOut / dRecTot; else cout << "_";
1321     cout << ", Ghost = ";
1322     if ( nRecTot > 0 ) cout << nGhost / dRecTot; else cout << "_";
1323     cout << endl;
1324   }
1325 }
1326
1327
1328
1329 void AliHLTTPCCAPerformance::MergerPerformance()
1330 {
1331   // performance calculation for merged tracks
1332
1333   int nRecTot = 0, nGhost = 0, nRecOut = 0;
1334   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
1335   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
1336
1337   // Select reconstructable MC tracks
1338
1339   {
1340     for ( int imc = 0; imc < fNMCTracks; imc++ ) fMCTracks[imc].SetNHits( 0 );
1341
1342     for ( int ih = 0; ih < fNHits; ih++ ) {
1343       AliHLTTPCCAHitLabel &l = fHitLabels[ih];
1344       if ( l.fLab[0] >= 0 ) fMCTracks[l.fLab[0]].SetNHits( fMCTracks[l.fLab[0]].NHits() + 1 );
1345       if ( l.fLab[1] >= 0 ) fMCTracks[l.fLab[1]].SetNHits( fMCTracks[l.fLab[1]].NHits() + 1 );
1346       if ( l.fLab[2] >= 0 ) fMCTracks[l.fLab[2]].SetNHits( fMCTracks[l.fLab[2]].NHits() + 1 );
1347     }
1348
1349     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
1350       AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
1351       mc.SetSet( 0 );
1352       mc.SetNReconstructed( 0 );
1353       mc.SetNTurns( 1 );
1354       if ( mc.NHits() >=  50 && mc.P() >= .05 ) {
1355         mc.SetSet( 1 );
1356         nMCAll++;
1357         if ( mc.P() >= 1. ) {
1358           mc.SetSet( 2 );
1359           nMCRef++;
1360         }
1361       }
1362     }
1363   }
1364
1365   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1366
1367   if ( !hlt.Merger().Output() ) return;
1368
1369   const AliHLTTPCCAMergerOutput &output = *( hlt.Merger().Output() );
1370
1371   int traN = output.NTracks();
1372
1373   nRecTot += traN;
1374
1375   for ( int itr = 0; itr < traN; itr++ ) {
1376
1377     const AliHLTTPCCAMergedTrack &tCA = output.Track( itr );
1378     std::vector<int> clusterIDs;
1379     for ( int i = 0; i < tCA.NClusters(); i++ ) {
1380       clusterIDs.push_back( output.ClusterId( tCA.FirstClusterRef() + i ) );
1381     }
1382     int label;
1383     float purity;
1384     GetMCLabel( clusterIDs, label, purity );
1385
1386     if ( purity < .9 || label < 0 || label >= fNMCTracks ) {
1387       nGhost++;
1388       continue;
1389     }
1390
1391     AliHLTTPCCAMCTrack &mc = fMCTracks[label];
1392     mc.SetNReconstructed( mc.NReconstructed() + 1 );
1393     if ( mc.Set() == 0 ) nRecOut++;
1394     else {
1395       if ( mc.NReconstructed() == 1 ) nRecAll++;
1396       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
1397       if ( mc.Set() == 2 ) {
1398         if ( mc.NReconstructed() == 1 ) nRecRef++;
1399         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
1400         fhTrackLengthRef->Fill( tCA.NClusters() / ( ( double ) mc.NHits() ) );
1401       }
1402     }
1403
1404     // track resolutions
1405     while ( mc.Set() == 2 && TMath::Abs( mc.TPCPar()[0] ) + TMath::Abs( mc.TPCPar()[1] ) > 1 ) {
1406
1407       if ( purity < .90 ) break;
1408       AliHLTTPCCATrackParam p = tCA.InnerParam();
1409       double cosA = TMath::Cos( tCA.InnerAlpha() );
1410       double sinA = TMath::Sin( tCA.InnerAlpha() );
1411       double mcX =  mc.TPCPar()[0] * cosA + mc.TPCPar()[1] * sinA;
1412       double mcY = -mc.TPCPar()[0] * sinA + mc.TPCPar()[1] * cosA;
1413       double mcZ =  mc.TPCPar()[2];
1414       double mcEx =  mc.TPCPar()[3] * cosA + mc.TPCPar()[4] * sinA;
1415       double mcEy = -mc.TPCPar()[3] * sinA + mc.TPCPar()[4] * cosA;
1416       double mcEz =  mc.TPCPar()[5];
1417       double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );
1418       if ( TMath::Abs( mcEt ) < 1.e-4 ) break;
1419       double mcSinPhi = mcEy / mcEt;
1420       double mcDzDs   = mcEz / mcEt;
1421       double mcQPt = mc.TPCPar()[6] / mcEt;
1422       if ( TMath::Abs( mcQPt ) < 1.e-4 ) break;
1423       double mcPt = 1. / TMath::Abs( mcQPt );
1424
1425       if ( mcPt < 1. ) break;
1426
1427       if ( tCA.NClusters() <  50 ) break;
1428       if ( !p.TransportToXWithMaterial( mcX, hlt.Merger().SliceParam().GetBz( p ) ) ) break;
1429       if ( p.GetCosPhi()*mcEx < 0 ) { // change direction
1430         mcSinPhi = -mcSinPhi;
1431         mcDzDs = -mcDzDs;
1432         mcQPt = -mcQPt;
1433       }
1434
1435       double qPt = p.GetQPt();
1436       double pt = 100;
1437       if ( TMath::Abs( qPt ) > 1.e-4 ) pt = 1. / TMath::Abs( qPt );
1438
1439       fhResY->Fill( p.GetY() - mcY );
1440       fhResZ->Fill( p.GetZ() - mcZ );
1441       fhResSinPhi->Fill( p.GetSinPhi() - mcSinPhi );
1442       fhResDzDs->Fill( p.GetDzDs() - mcDzDs );
1443       fhResPt->Fill( ( pt - mcPt ) / mcPt );
1444
1445       if ( p.GetErr2Y() > 0 ) fhPullY->Fill( ( p.GetY() - mcY ) / TMath::Sqrt( p.GetErr2Y() ) );
1446       if ( p.GetErr2Z() > 0 ) fhPullZ->Fill( ( p.GetZ() - mcZ ) / TMath::Sqrt( p.GetErr2Z() ) );
1447
1448       if ( p.GetErr2SinPhi() > 0 ) fhPullSinPhi->Fill( ( p.GetSinPhi() - mcSinPhi ) / TMath::Sqrt( p.GetErr2SinPhi() ) );
1449       if ( p.GetErr2DzDs() > 0 ) fhPullDzDs->Fill( ( p.DzDs() - mcDzDs ) / TMath::Sqrt( p.GetErr2DzDs() ) );
1450       if ( p.GetErr2QPt() > 0 ) fhPullQPt->Fill( ( qPt - mcQPt ) / TMath::Sqrt( p.GetErr2QPt() ) );
1451       fhPullYS->Fill( TMath::Sqrt( hlt.Merger().GetChi2( p.GetY(), p.GetSinPhi(), p.GetCov()[0], p.GetCov()[3], p.GetCov()[5], mcY, mcSinPhi, 0, 0, 0 ) ) );
1452       fhPullZT->Fill( TMath::Sqrt( hlt.Merger().GetChi2( p.GetZ(), p.GetDzDs(), p.GetCov()[2], p.GetCov()[7], p.GetCov()[9], mcZ, mcDzDs, 0, 0, 0 ) ) );
1453
1454       break;
1455     } // end resolutions
1456
1457   }// end reco tracks
1458
1459
1460   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1461     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1462     if ( mc.Set() > 0 ) fhGBEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1463     if ( mc.Set() > 0 ) fhGBEffVsPt->Fill( mc.Pt(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1464     if ( mc.Set() == 2 ) {
1465       const double *p = mc.TPCPar();
1466       double r = TMath::Sqrt( p[0] * p[0] + p[1] * p[1] );
1467       double cosA = p[0] / r;
1468       double sinA = p[1] / r;
1469
1470
1471       double phipos = TMath::Pi() + TMath::ATan2( -p[1], -p[0] );
1472       double alpha =  TMath::Pi() * ( 20 * ( ( ( ( int )( phipos * 180 / TMath::Pi() ) ) / 20 ) ) + 10 ) / 180.;
1473       cosA = TMath::Cos( alpha );
1474       sinA = TMath::Sin( alpha );
1475
1476       double mcX =  p[0] * cosA + p[1] * sinA;
1477       double mcY = -p[0] * sinA + p[1] * cosA;
1478       double mcZ =  p[2];
1479       double mcEx =  p[3] * cosA + p[4] * sinA;
1480       double mcEy = -p[3] * sinA + p[4] * cosA;
1481       double mcEz =  p[5];
1482       //double mcEt = TMath::Sqrt(mcEx*mcEx + mcEy*mcEy);
1483       double angleY = TMath::ATan2( mcEy, mcEx ) * 180. / TMath::Pi();
1484       double angleZ = TMath::ATan2( mcEz, mcEx ) * 180. / TMath::Pi();
1485
1486       if ( mc.NReconstructed() > 0 ) {
1487         fhRefRecoX->Fill( mcX );
1488         fhRefRecoY->Fill( mcY );
1489         fhRefRecoZ->Fill( mcZ );
1490         fhRefRecoP->Fill( mc.P() );
1491         fhRefRecoPt->Fill( mc.Pt() );
1492         fhRefRecoAngleY->Fill( angleY );
1493         fhRefRecoAngleZ->Fill( angleZ );
1494         fhRefRecoNHits->Fill( mc.NHits() );
1495       } else {
1496         fhRefNotRecoX->Fill( mcX );
1497         fhRefNotRecoY->Fill( mcY );
1498         fhRefNotRecoZ->Fill( mcZ );
1499         fhRefNotRecoP->Fill( mc.P() );
1500         fhRefNotRecoPt->Fill( mc.Pt() );
1501         fhRefNotRecoAngleY->Fill( angleY );
1502         fhRefNotRecoAngleZ->Fill( angleZ );
1503         fhRefNotRecoNHits->Fill( mc.NHits() );
1504       }
1505     }
1506   }
1507
1508   fStatGBNRecTot += nRecTot;
1509   fStatGBNRecOut += nRecOut;
1510   fStatGBNGhost  += nGhost;
1511   fStatGBNMCAll  += nMCAll;
1512   fStatGBNRecAll  += nRecAll;
1513   fStatGBNClonesAll  += nClonesAll;
1514   fStatGBNMCRef  += nMCRef;
1515   fStatGBNRecRef  += nRecRef;
1516   fStatGBNClonesRef  += nClonesRef;
1517 }
1518
1519
1520
1521 void AliHLTTPCCAPerformance::ClusterPerformance()
1522 {
1523   // performance calculation for input clusters
1524
1525   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1526
1527   // distribution of cluster errors
1528
1529   for ( int iSlice = 0; iSlice < hlt.NSlices(); iSlice++ ) {
1530     const AliHLTTPCCAClusterData &data = hlt.ClusterData( iSlice );
1531     for ( int i = 0; i < data.NumberOfClusters(); i++ ) {
1532       AliHLTTPCCAHitLabel &l = fHitLabels[data.Id( i )];
1533       int nmc = 0;
1534       for ( int il = 0; il < 3; il++ ) if ( l.fLab[il] >= 0 ) nmc++;
1535       if ( nmc == 1 ) fhHitShared->Fill( data.RowNumber( i ), 0 );
1536       else if ( nmc > 1 ) fhHitShared->Fill( data.RowNumber( i ), 1 );
1537     }
1538   }
1539
1540   // cluster pulls
1541
1542   if ( !fDoClusterPulls || fNMCPoints <= 0 ) return;
1543
1544   // sort mc points
1545   if ( 1 ) {
1546     for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1547       AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1548       mc.SetNMCPoints( 0 );
1549     }
1550     sort( fMCPoints, fMCPoints + fNMCPoints, AliHLTTPCCAMCPoint::Compare );
1551
1552     for ( int ip = 0; ip < fNMCPoints; ip++ ) {
1553       AliHLTTPCCAMCPoint &p = fMCPoints[ip];
1554       AliHLTTPCCAMCTrack &t = fMCTracks[p.TrackID()];
1555       if ( t.NMCPoints() == 0 ) t.SetFirstMCPointID( ip );
1556       t.SetNMCPoints( t.NMCPoints() + 1 );
1557     }
1558   }
1559
1560   for ( int iSlice = 0; iSlice < hlt.NSlices(); iSlice++ ) {
1561
1562     const AliHLTTPCCAClusterData &data = hlt.ClusterData( iSlice );
1563
1564     for ( int ic = 0; ic < data.NumberOfClusters(); ic++ ) {
1565
1566       const AliHLTTPCCAHitLabel &l = fHitLabels[data.Id( ic )];
1567
1568       if ( l.fLab[0] < 0 || l.fLab[0] >= fNMCTracks
1569            || l.fLab[1] >= 0 || l.fLab[2] >= 0       ) continue;
1570
1571       int lab = l.fLab[0];
1572
1573       AliHLTTPCCAMCTrack &mc = fMCTracks[lab];
1574
1575       double x0 = data.X( ic );
1576       double y0 = data.Y( ic );
1577       double z0 = data.Z( ic );
1578
1579       if ( fabs( x0 ) < 1.e-4 ) continue;
1580       if ( mc.Pt() < .05 ) continue;
1581
1582       int ip1 = -1, ip2 = -1;
1583       double d1 = 1.e20, d2 = 1.e20;
1584
1585       AliHLTTPCCAMCPoint *pStart = lower_bound( fMCPoints + mc.FirstMCPointID(), fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), iSlice,  AliHLTTPCCAMCPoint::CompareSlice );
1586
1587       pStart = lower_bound( pStart, fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), x0 - 2.,  AliHLTTPCCAMCPoint::CompareX );
1588
1589       for ( int ip = ( pStart - fMCPoints ) - mc.FirstMCPointID(); ip < mc.NMCPoints(); ip++ ) {
1590         AliHLTTPCCAMCPoint &p = fMCPoints[mc.FirstMCPointID() + ip];
1591         if ( p.ISlice() != iSlice ) break;
1592         double dx = p.Sx() - x0;
1593         double dy = p.Sy() - y0;
1594         double dz = p.Sz() - z0;
1595         double d = dx * dx + dy * dy + dz * dz;
1596         if ( d > 9. ) continue;
1597         if ( dx <= 0  && dx > -2. ) {
1598           if ( fabs( dx ) < d1 ) {
1599             d1 = fabs( dx );
1600             ip1 = ip;
1601           }
1602         } else if ( dx > .2 ) {
1603           if ( dx >= 2. ) break;
1604           if ( fabs( dx ) < d2 ) {
1605             d2 = fabs( dx );
1606             ip2 = ip;
1607           }
1608         }
1609       }
1610
1611       if ( ip1 < 0 || ip2 < 0 ) continue;
1612
1613       AliHLTTPCCAMCPoint &p1 = fMCPoints[mc.FirstMCPointID() + ip1];
1614       AliHLTTPCCAMCPoint &p2 = fMCPoints[mc.FirstMCPointID() + ip2];
1615       double dx = p2.Sx() - p1.Sx();
1616       double dy = p2.Sy() - p1.Sy();
1617       double dz = p2.Sz() - p1.Sz();
1618       double sx = x0;
1619       double sy = p1.Sy() + dy / dx * ( sx - p1.Sx() );
1620       double sz = p1.Sz() + dz / dx * ( sx - p1.Sx() );
1621
1622       float errY, errZ;
1623       {
1624         AliHLTTPCCATrackParam t;
1625         double s = 1. / TMath::Sqrt( dx * dx + dy * dy );
1626         t.SetZ( sz );
1627         t.SetSinPhi( dy * s );
1628         t.SetSignCosPhi( dx );
1629         t.SetDzDs( dz * s );
1630         //hlt.SliceTracker( 0 ).GetErrors2( data.RowNumber( ic ), t, errY, errZ );
1631                 hlt.Param(0).GetClusterErrors2( data.RowNumber( ic ), t.GetZ(), t.SinPhi(), t.GetCosPhi(), t.DzDs(), errY, errZ );
1632         errY = TMath::Sqrt( errY );
1633         errZ = TMath::Sqrt( errZ );
1634       }
1635       fhHitErrY->Fill( errY );
1636       fhHitErrZ->Fill( errZ );
1637       fhHitResY->Fill( y0 - sy );
1638       fhHitResZ->Fill( z0 - sz );
1639       fhHitPullY->Fill( ( y0 - sy ) / errY );
1640       fhHitPullZ->Fill( ( z0 - sz ) / errZ );
1641       if ( mc.Pt() >= 1. ) {
1642         fhHitResY1->Fill( y0 - sy );
1643         fhHitResZ1->Fill( z0 - sz );
1644         fhHitPullY1->Fill( ( y0 - sy ) / errY );
1645         fhHitPullZ1->Fill( ( z0 - sz ) / errZ );
1646       }
1647     }
1648   }
1649 }
1650
1651
1652 void AliHLTTPCCAPerformance::SmearClustersMC()
1653 {
1654   // smear clusters with gaussian using MC info
1655
1656   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1657
1658   // cluster pulls
1659
1660   if ( fNMCPoints <= 0 ) return;
1661
1662   // sort mc points
1663   {
1664     for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1665       AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1666       mc.SetNMCPoints( 0 );
1667     }
1668     sort( fMCPoints, fMCPoints + fNMCPoints, AliHLTTPCCAMCPoint::Compare );
1669
1670     for ( int ip = 0; ip < fNMCPoints; ip++ ) {
1671       AliHLTTPCCAMCPoint &p = fMCPoints[ip];
1672       AliHLTTPCCAMCTrack &t = fMCTracks[p.TrackID()];
1673       if ( t.NMCPoints() == 0 ) t.SetFirstMCPointID( ip );
1674       t.SetNMCPoints( t.NMCPoints() + 1 );
1675     }
1676   }
1677
1678   for ( int iSlice = 0; iSlice < hlt.NSlices(); iSlice++ ) {
1679
1680     AliHLTTPCCAClusterData &data = hlt.ClusterData( iSlice );
1681
1682     for ( int ic = 0; ic < data.NumberOfClusters(); ic++ ) {
1683
1684       double x0 = data.X( ic );
1685       double y0 = data.Y( ic );
1686       double z0 = data.Z( ic );
1687       int row0 = data.RowNumber( ic );
1688
1689       AliHLTTPCCAClusterData::Data *cdata = data.GetClusterData( ic );
1690       cdata->fX = 0;
1691       cdata->fY = 0;
1692       cdata->fZ = 0;
1693
1694       const AliHLTTPCCAHitLabel &l = fHitLabels[data.Id( ic )];
1695
1696       if ( l.fLab[0] < 0 || l.fLab[0] >= fNMCTracks ) continue;
1697
1698       int lab = l.fLab[0];
1699
1700       AliHLTTPCCAMCTrack &mc = fMCTracks[lab];
1701
1702       int ip1 = -1, ip2 = -1;
1703       double d1 = 1.e20, d2 = 1.e20;
1704
1705       AliHLTTPCCAMCPoint *pStart = lower_bound( fMCPoints + mc.FirstMCPointID(), fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), iSlice,  AliHLTTPCCAMCPoint::CompareSlice );
1706
1707       pStart = lower_bound( pStart, fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), x0 - 2.,  AliHLTTPCCAMCPoint::CompareX );
1708
1709       for ( int ip = ( pStart - fMCPoints ) - mc.FirstMCPointID(); ip < mc.NMCPoints(); ip++ ) {
1710         AliHLTTPCCAMCPoint &p = fMCPoints[mc.FirstMCPointID() + ip];
1711         if ( p.ISlice() != iSlice ) break;
1712         double dx = p.Sx() - x0;
1713         double dy = p.Sy() - y0;
1714         double dz = p.Sz() - z0;
1715         double d = dx * dx + dy * dy + dz * dz;
1716         if ( d > 9. ) continue;
1717         if ( dx <= 0  && dx > -2. ) {
1718           if ( fabs( dx ) < d1 ) {
1719             d1 = fabs( dx );
1720             ip1 = ip;
1721           }
1722         } else if ( dx > .2 ) {
1723           if ( dx >= 2. ) break;
1724           if ( fabs( dx ) < d2 ) {
1725             d2 = fabs( dx );
1726             ip2 = ip;
1727           }
1728         }
1729       }
1730
1731       if ( ip1 < 0 || ip2 < 0 ) continue;
1732
1733       AliHLTTPCCAMCPoint &p1 = fMCPoints[mc.FirstMCPointID() + ip1];
1734       AliHLTTPCCAMCPoint &p2 = fMCPoints[mc.FirstMCPointID() + ip2];
1735       double dx = p2.Sx() - p1.Sx();
1736       double dy = p2.Sy() - p1.Sy();
1737       double dz = p2.Sz() - p1.Sz();
1738       double sx = x0;
1739       double sy = p1.Sy() + dy / dx * ( sx - p1.Sx() );
1740       double sz = p1.Sz() + dz / dx * ( sx - p1.Sx() );
1741
1742       float errY, errZ;
1743       {
1744         AliHLTTPCCATrackParam t;
1745         double s = 1. / TMath::Sqrt( dx * dx + dy * dy );
1746         t.SetZ( sz );
1747         t.SetSinPhi( dy * s );
1748         t.SetSignCosPhi( dx );
1749         t.SetDzDs( dz * s );
1750         //hlt.SliceTracker( 0 ).GetErrors2( row0, t, errY, errZ );
1751                 hlt.Param(0).GetClusterErrors2( row0, t.GetZ(), t.SinPhi(), t.GetCosPhi(), t.DzDs(), errY, errZ );
1752         errY = TMath::Sqrt( errY );
1753         errZ = TMath::Sqrt( errZ );
1754       }
1755
1756       cdata->fX = x0;
1757       cdata->fY = gRandom->Gaus( sy, errY );
1758       cdata->fZ = gRandom->Gaus( sz, errZ );
1759     }
1760   }
1761 }
1762
1763
1764 void AliHLTTPCCAPerformance::Performance( fstream *StatFile )
1765 {
1766   // main routine for performance calculation
1767
1768   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1769
1770   //SG!!!
1771   /*
1772   fStatNEvents=0;
1773     fStatNRecTot=0;
1774     fStatNRecOut=0;
1775     fStatNGhost=0;
1776     fStatNMCAll=0;
1777     fStatNRecAll=0;
1778     fStatNClonesAll=0;
1779     fStatNMCRef=0;
1780     fStatNRecRef=0;
1781     fStatNClonesRef=0;
1782   */
1783   fStatNEvents++;
1784   for ( int islice = 0; islice < hlt.NSlices(); islice++ ) {
1785     SliceTrackletPerformance( islice, 0 );
1786     SliceTrackCandPerformance( islice, 0 );
1787     SlicePerformance( islice, 0 );
1788   }
1789
1790   MergerPerformance();
1791   //ClusterPerformance();
1792
1793   {
1794     cout << "\nSlice Track Seed performance: \n" << endl;
1795     cout << " N tracks : "
1796          << fStatNMCAll / fStatNEvents << " mc all, "
1797          << fStatSeedNMCRef / fStatNEvents << " mc ref, "
1798          << fStatSeedNRecTot / fStatNEvents << " rec total, "
1799          << fStatSeedNRecAll / fStatNEvents << " rec all, "
1800          << fStatSeedNClonesAll / fStatNEvents << " clones all, "
1801          << fStatSeedNRecRef / fStatNEvents << " rec ref, "
1802          << fStatSeedNClonesRef / fStatNEvents << " clones ref, "
1803          << fStatSeedNRecOut / fStatNEvents << " out, "
1804          << fStatSeedNGhost / fStatNEvents << " ghost" << endl;
1805
1806     int nRecExtr = fStatSeedNRecAll - fStatSeedNRecRef;
1807     int nMCExtr = fStatNMCAll - fStatNMCRef;
1808     int nClonesExtr = fStatSeedNClonesAll - fStatSeedNClonesRef;
1809
1810     double dRecTot = ( fStatSeedNRecTot > 0 ) ? fStatSeedNRecTot : 1;
1811     double dMCAll = ( fStatNMCAll > 0 ) ? fStatNMCAll : 1;
1812     double dMCRef = ( fStatNMCRef > 0 ) ? fStatNMCRef : 1;
1813     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1814     double dRecAll = ( fStatSeedNRecAll + fStatSeedNClonesAll > 0 ) ? fStatSeedNRecAll + fStatSeedNClonesAll : 1;
1815     double dRecRef = ( fStatSeedNRecRef + fStatSeedNClonesRef > 0 ) ? fStatSeedNRecRef + fStatSeedNClonesRef : 1;
1816     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1817
1818     cout << " EffRef = " << fStatSeedNRecRef / dMCRef
1819          << ", CloneRef = " << fStatSeedNClonesRef / dRecRef << endl;
1820     cout << " EffExtra = " << nRecExtr / dMCExtr
1821          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1822     cout << " EffAll = " << fStatSeedNRecAll / dMCAll
1823          << ", CloneAll = " << fStatSeedNClonesAll / dRecAll << endl;
1824     cout << " Out = " << fStatSeedNRecOut / dRecTot
1825          << ", Ghost = " << fStatSeedNGhost / dRecTot << endl;
1826   }
1827
1828   {
1829     cout << "\nSlice Track candidate performance: \n" << endl;
1830     cout << " N tracks : "
1831          << fStatNMCAll / fStatNEvents << " mc all, "
1832          << fStatCandNMCRef / fStatNEvents << " mc ref, "
1833          << fStatCandNRecTot / fStatNEvents << " rec total, "
1834          << fStatCandNRecAll / fStatNEvents << " rec all, "
1835          << fStatCandNClonesAll / fStatNEvents << " clones all, "
1836          << fStatCandNRecRef / fStatNEvents << " rec ref, "
1837          << fStatCandNClonesRef / fStatNEvents << " clones ref, "
1838          << fStatCandNRecOut / fStatNEvents << " out, "
1839          << fStatCandNGhost / fStatNEvents << " ghost" << endl;
1840
1841     int nRecExtr = fStatCandNRecAll - fStatCandNRecRef;
1842     int nMCExtr = fStatNMCAll - fStatNMCRef;
1843     int nClonesExtr = fStatCandNClonesAll - fStatCandNClonesRef;
1844
1845     double dRecTot = ( fStatCandNRecTot > 0 ) ? fStatCandNRecTot : 1;
1846     double dMCAll = ( fStatNMCAll > 0 ) ? fStatNMCAll : 1;
1847     double dMCRef = ( fStatNMCRef > 0 ) ? fStatNMCRef : 1;
1848     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1849     double dRecAll = ( fStatCandNRecAll + fStatCandNClonesAll > 0 ) ? fStatCandNRecAll + fStatCandNClonesAll : 1;
1850     double dRecRef = ( fStatCandNRecRef + fStatCandNClonesRef > 0 ) ? fStatCandNRecRef + fStatCandNClonesRef : 1;
1851     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1852
1853     cout << " EffRef = " << fStatCandNRecRef / dMCRef
1854          << ", CloneRef = " << fStatCandNClonesRef / dRecRef << endl;
1855     cout << " EffExtra = " << nRecExtr / dMCExtr
1856          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1857     cout << " EffAll = " << fStatCandNRecAll / dMCAll
1858          << ", CloneAll = " << fStatCandNClonesAll / dRecAll << endl;
1859     cout << " Out = " << fStatCandNRecOut / dRecTot
1860          << ", Ghost = " << fStatCandNGhost / dRecTot << endl;
1861   }
1862
1863   {
1864     cout << "\nSlice tracker performance: \n" << endl;
1865     cout << " N tracks : "
1866          << fStatNMCAll / fStatNEvents << " mc all, "
1867          << fStatNMCRef / fStatNEvents << " mc ref, "
1868          << fStatNRecTot / fStatNEvents << " rec total, "
1869          << fStatNRecAll / fStatNEvents << " rec all, "
1870          << fStatNClonesAll / fStatNEvents << " clones all, "
1871          << fStatNRecRef / fStatNEvents << " rec ref, "
1872          << fStatNClonesRef / fStatNEvents << " clones ref, "
1873          << fStatNRecOut / fStatNEvents << " out, "
1874          << fStatNGhost / fStatNEvents << " ghost" << endl;
1875
1876     int nRecExtr = fStatNRecAll - fStatNRecRef;
1877     int nMCExtr = fStatNMCAll - fStatNMCRef;
1878     int nClonesExtr = fStatNClonesAll - fStatNClonesRef;
1879
1880     double dRecTot = ( fStatNRecTot > 0 ) ? fStatNRecTot : 1;
1881     double dMCAll = ( fStatNMCAll > 0 ) ? fStatNMCAll : 1;
1882     double dMCRef = ( fStatNMCRef > 0 ) ? fStatNMCRef : 1;
1883     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1884     double dRecAll = ( fStatNRecAll + fStatNClonesAll > 0 ) ? fStatNRecAll + fStatNClonesAll : 1;
1885     double dRecRef = ( fStatNRecRef + fStatNClonesRef > 0 ) ? fStatNRecRef + fStatNClonesRef : 1;
1886     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1887
1888     cout << " EffRef = " << fStatNRecRef / dMCRef
1889          << ", CloneRef = " << fStatNClonesRef / dRecRef << endl;
1890     cout << " EffExtra = " << nRecExtr / dMCExtr
1891          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1892     cout << " EffAll = " << fStatNRecAll / dMCAll
1893          << ", CloneAll = " << fStatNClonesAll / dRecAll << endl;
1894     cout << " Out = " << fStatNRecOut / dRecTot
1895          << ", Ghost = " << fStatNGhost / dRecTot << endl;
1896     cout << " Time = " << hlt.StatTime( 0 ) / hlt.StatNEvents()*1.e3 << " msec/event " << endl;
1897     cout << " Local timers = "
1898          << hlt.StatTime( 1 ) / hlt.StatNEvents()*1.e3 << " "
1899          << hlt.StatTime( 2 ) / hlt.StatNEvents()*1.e3 << " "
1900          << hlt.StatTime( 3 ) / hlt.StatNEvents()*1.e3 << " "
1901          << hlt.StatTime( 4 ) / hlt.StatNEvents()*1.e3 << " "
1902          << hlt.StatTime( 5 ) / hlt.StatNEvents()*1.e3 << " "
1903          << hlt.StatTime( 6 ) / hlt.StatNEvents()*1.e3 << " "
1904          << hlt.StatTime( 7 ) / hlt.StatNEvents()*1.e3 << " "
1905          << hlt.StatTime( 8 ) / hlt.StatNEvents()*1.e3 << " "
1906          << " msec/event " << endl;
1907   }
1908
1909
1910   {
1911     cout << "\nGlobal tracker performance for " << fStatNEvents << " events: \n" << endl;
1912     cout << " N tracks : "
1913          << fStatGBNMCAll << " mc all, "
1914          << fStatGBNMCRef << " mc ref, "
1915          << fStatGBNRecTot << " rec total, "
1916          << fStatGBNRecAll << " rec all, "
1917          << fStatGBNClonesAll << " clones all, "
1918          << fStatGBNRecRef << " rec ref, "
1919          << fStatGBNClonesRef << " clones ref, "
1920          << fStatGBNRecOut << " out, "
1921          << fStatGBNGhost << " ghost" << endl;
1922     cout << " N tracks average : "
1923          << fStatGBNMCAll / fStatNEvents << " mc all, "
1924          << fStatGBNMCRef / fStatNEvents << " mc ref, "
1925          << fStatGBNRecTot / fStatNEvents << " rec total, "
1926          << fStatGBNRecAll / fStatNEvents << " rec all, "
1927          << fStatGBNClonesAll / fStatNEvents << " clones all, "
1928          << fStatGBNRecRef / fStatNEvents << " rec ref, "
1929          << fStatGBNClonesRef / fStatNEvents << " clones ref, "
1930          << fStatGBNRecOut / fStatNEvents << " out, "
1931          << fStatGBNGhost / fStatNEvents << " ghost" << endl;
1932
1933     int nRecExtr = fStatGBNRecAll - fStatGBNRecRef;
1934     int nMCExtr = fStatGBNMCAll - fStatGBNMCRef;
1935     int nClonesExtr = fStatGBNClonesAll - fStatGBNClonesRef;
1936
1937     double dRecTot = ( fStatGBNRecTot > 0 ) ? fStatGBNRecTot : 1;
1938     double dMCAll = ( fStatGBNMCAll > 0 ) ? fStatGBNMCAll : 1;
1939     double dMCRef = ( fStatGBNMCRef > 0 ) ? fStatGBNMCRef : 1;
1940     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1941     double dRecAll = ( fStatGBNRecAll + fStatGBNClonesAll > 0 ) ? fStatGBNRecAll + fStatGBNClonesAll : 1;
1942     double dRecRef = ( fStatGBNRecRef + fStatGBNClonesRef > 0 ) ? fStatGBNRecRef + fStatGBNClonesRef : 1;
1943     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1944
1945     cout << " EffRef = " << fStatGBNRecRef / dMCRef
1946          << ", CloneRef = " << fStatGBNClonesRef / dRecRef << endl;
1947     cout << " EffExtra = " << nRecExtr / dMCExtr
1948          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1949     cout << " EffAll = " << fStatGBNRecAll / dMCAll
1950          << ", CloneAll = " << fStatGBNClonesAll / dRecAll << endl;
1951     cout << " Out = " << fStatGBNRecOut / dRecTot
1952          << ", Ghost = " << fStatGBNGhost / dRecTot << endl;
1953     cout << " Time = " << ( hlt.StatTime( 0 ) + hlt.StatTime( 9 ) ) / hlt.StatNEvents()*1.e3 << " msec/event " << endl;
1954     cout << " Local timers: " << endl;
1955     cout << " slice tracker " << hlt.StatTime( 0 ) / hlt.StatNEvents()*1.e3 << ": "
1956          << hlt.StatTime( 1 ) / hlt.StatNEvents()*1.e3 << " "
1957          << hlt.StatTime( 2 ) / hlt.StatNEvents()*1.e3 << " "
1958          << hlt.StatTime( 3 ) / hlt.StatNEvents()*1.e3 << " "
1959          << hlt.StatTime( 4 ) / hlt.StatNEvents()*1.e3 << " "
1960          << hlt.StatTime( 5 ) / hlt.StatNEvents()*1.e3 << "["
1961          << hlt.StatTime( 6 ) / hlt.StatNEvents()*1.e3 << "/"
1962          << hlt.StatTime( 7 ) / hlt.StatNEvents()*1.e3 << "] "
1963          << hlt.StatTime( 8 ) / hlt.StatNEvents()*1.e3
1964          << " msec/event " << endl;
1965     cout << " GB merger " << hlt.StatTime( 9 ) / hlt.StatNEvents()*1.e3 << ": "
1966          << hlt.StatTime( 10 ) / hlt.StatNEvents()*1.e3 << ", "
1967          << hlt.StatTime( 11 ) / hlt.StatNEvents()*1.e3 << ", "
1968          << hlt.StatTime( 12 ) / hlt.StatNEvents()*1.e3 << " "
1969          << " msec/event " << endl;
1970
1971     if ( StatFile && StatFile->is_open() ) {
1972       fstream &out = *StatFile;
1973
1974       //out<<"\nGlobal tracker performance for "<<fStatNEvents<<" events: \n"<<endl;
1975       //out<<" N tracks : "
1976       //<<fStatGBNMCAll/fStatNEvents<<" mc all, "
1977       //<<fStatGBNMCRef/fStatNEvents<<" mc ref, "
1978       // <<fStatGBNRecTot/fStatNEvents<<" rec total, "
1979       // <<fStatGBNRecAll/fStatNEvents<<" rec all, "
1980       // <<fStatGBNClonesAll/fStatNEvents<<" clones all, "
1981       // <<fStatGBNRecRef/fStatNEvents<<" rec ref, "
1982       // <<fStatGBNClonesRef/fStatNEvents<<" clones ref, "
1983       // <<fStatGBNRecOut/fStatNEvents<<" out, "
1984       // <<fStatGBNGhost/fStatNEvents<<" ghost"<<endl;
1985       fStatTime += hlt.StatTime( 0 );
1986       double timeHz = 0;
1987       if ( fStatTime > 1.e-4 ) timeHz = 1. / fStatTime * fStatNEvents;
1988
1989       out << "<table border>" << endl;
1990       out << "<tr>" << endl;
1991       out << "<td>      </td> <td align=center> RefSet </td> <td align=center> AllSet </td> <td align=center> ExtraSet </td>" << endl;
1992       out << "</tr>" << endl;
1993       out << "<tr>" << endl;
1994       out << "<td>Efficiency</td> <td align=center>" << fStatGBNRecRef / dMCRef
1995       << "</td> <td align=center>" << fStatGBNRecAll / dMCAll
1996       << "</td> <td align=center>" << nRecExtr / dMCExtr
1997       << "</td>" << endl;
1998       out << "</tr>" << endl;
1999       out << "<tr> " << endl;
2000       out << "<td>Clone</td>      <td align=center>" << fStatGBNClonesRef / dRecRef
2001       << "</td> <td align=center>" << fStatGBNClonesAll / dRecAll
2002       << "</td> <td align=center>" << nClonesExtr / dRecExtr
2003       << "</td>" << endl;
2004       out << "</tr>" << endl;
2005       out << "<tr> " << endl;
2006       out << "<td>Ghost</td>      <td colspan=3 align=center>" << fStatGBNGhost / dRecTot
2007       << "</td>" << endl;
2008       out << "</tr>" << endl;
2009       out << "<tr> " << endl;
2010       out << "<td>Time</td>      <td colspan=3 align=center>" << timeHz
2011       << " ev/s</td>" << endl;
2012       out << "</tr>" << endl;
2013       out << "<tr> " << endl;
2014       out << "<td>N Events</td>      <td colspan=3 align=center>" << fStatNEvents
2015       << "</td>" << endl;
2016       out << "</tr>" << endl;
2017       out << "</table>" << endl;
2018     }
2019
2020   }
2021
2022   WriteHistos();
2023 }
2024
2025
2026 void AliHLTTPCCAPerformance::WriteMCEvent( ostream &out ) const
2027 {
2028   // write MC information to the file
2029   out << fNMCTracks << endl;
2030   for ( int it = 0; it < fNMCTracks; it++ ) {
2031     AliHLTTPCCAMCTrack &t = fMCTracks[it];
2032     out << it << " ";
2033     out << t.PDG() << endl;
2034     for ( int i = 0; i < 7; i++ ) out << t.Par()[i] << " ";
2035     out << endl << "    ";
2036     for ( int i = 0; i < 7; i++ ) out << t.TPCPar()[i] << " ";
2037     out << endl << "    ";
2038     out << t.P() << " ";
2039     out << t.Pt() << " ";
2040     out << t.NMCPoints() << " ";
2041     out << t.FirstMCPointID() << " ";
2042     out << t.NHits() << " ";
2043     out << t.NReconstructed() << " ";
2044     out << t.Set() << " ";
2045     out << t.NTurns() << endl;
2046   }
2047
2048   out << fNHits << endl;
2049   for ( int ih = 0; ih < fNHits; ih++ ) {
2050     AliHLTTPCCAHitLabel &l = fHitLabels[ih];
2051     out << l.fLab[0] << " " << l.fLab[1] << " " << l.fLab[2] << endl;
2052   }
2053 }
2054
2055 void AliHLTTPCCAPerformance::WriteMCPoints( ostream &out ) const
2056 {
2057   // write Mc points to the file
2058   out << fNMCPoints << endl;
2059   for ( int ip = 0; ip < fNMCPoints; ip++ ) {
2060     AliHLTTPCCAMCPoint &p = fMCPoints[ip];
2061     out << p.X() << " ";
2062     out << p.Y() << " ";
2063     out << p.Z() << " ";
2064     out << p.Sx() << " ";
2065     out << p.Sy() << " ";
2066     out << p.Sz() << " ";
2067     out << p.Time() << " ";
2068     out << p.ISlice() << " ";
2069     out << p.TrackID() << endl;
2070   }
2071 }
2072
2073 void AliHLTTPCCAPerformance::ReadMCEvent( istream &in )
2074 {
2075   // read mc info from the file
2076   StartEvent();
2077   if ( fMCTracks ) delete[] fMCTracks;
2078   fMCTracks = 0;
2079   fNMCTracks = 0;
2080   if ( fHitLabels ) delete[] fHitLabels;
2081   fHitLabels = 0;
2082   fNHits = 0;
2083   if ( fMCPoints ) delete[] fMCPoints;
2084   fMCPoints = 0;
2085   fNMCPoints = 0;
2086
2087   in >> fNMCTracks;
2088   fMCTracks = new AliHLTTPCCAMCTrack[fNMCTracks];
2089   for ( int it = 0; it < fNMCTracks; it++ ) {
2090     AliHLTTPCCAMCTrack &t = fMCTracks[it];
2091     int j;
2092     float f;
2093     in >> j;
2094     in >> j; t.SetPDG( j );
2095     for ( int i = 0; i < 7; i++ ) { in >> f; t.SetPar( i, f );}
2096     for ( int i = 0; i < 7; i++ ) { in >> f; t.SetTPCPar( i, f );}
2097     in >> f; t.SetP( f );
2098     in >> f; t.SetPt( f );
2099     in >> j; t.SetNHits( j );
2100     in >> j; t.SetNMCPoints( j );
2101     in >> j; t.SetFirstMCPointID( j );
2102     in >> j; t.SetNReconstructed( j );
2103     in >> j; t.SetSet( j );
2104     in >> j; t.SetNTurns( j );
2105   }
2106
2107   in >> fNHits;
2108   fHitLabels = new AliHLTTPCCAHitLabel[fNHits];
2109   for ( int ih = 0; ih < fNHits; ih++ ) {
2110     AliHLTTPCCAHitLabel &l = fHitLabels[ih];
2111     in >> l.fLab[0] >> l.fLab[1] >> l.fLab[2];
2112   }
2113 }
2114
2115 void AliHLTTPCCAPerformance::ReadMCPoints( istream &in )
2116 {
2117   // read mc points from the file
2118   if ( fMCPoints ) delete[] fMCPoints;
2119   fMCPoints = 0;
2120   fNMCPoints = 0;
2121
2122   in >> fNMCPoints;
2123   fMCPoints = new AliHLTTPCCAMCPoint[fNMCPoints];
2124   for ( int ip = 0; ip < fNMCPoints; ip++ ) {
2125     AliHLTTPCCAMCPoint &p = fMCPoints[ip];
2126     float f;
2127     int i;
2128     in >> f;
2129     p.SetX( f );
2130     in >> f;
2131     p.SetY( f );
2132     in >> f;
2133     p.SetZ( f );
2134     in >> f;
2135     p.SetSx( f );
2136     in >> f;
2137     p.SetSy( f );
2138     in >> f;
2139     p.SetSz( f );
2140     in >> f;
2141     p.SetTime( f );
2142     in >> i;
2143     p.SetISlice( i );
2144     in >> i;
2145     p.SetTrackID( i );
2146   }
2147 }