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().SliceTracker( iSlice ).Param().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         int ih = t.RowHit( irow );
1039         if ( ih < 0 ) continue;
1040         int index = firstSliceHit + slice.HitInputID( slice.Row( irow ), ih );
1041         AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
1042         if ( l.fLab[0] >= 0 ) lb[nla++] = l.fLab[0];
1043         if ( l.fLab[1] >= 0 ) lb[nla++] = l.fLab[1];
1044         if ( l.fLab[2] >= 0 ) lb[nla++] = l.fLab[2];
1045         nHits++;
1046       }
1047       if ( nHits < 10 ) continue;
1048
1049       sort( lb, lb + nla );
1050       int labmax = -1, labcur = -1, lmax = 0, lcurr = 0;
1051       for ( int i = 0; i < nla; i++ ) {
1052         if ( lb[i] != labcur ) {
1053           if ( labcur >= 0 && lmax < lcurr ) {
1054             lmax = lcurr;
1055             labmax = labcur;
1056           }
1057           labcur = lb[i];
1058           lcurr = 0;
1059         }
1060         lcurr++;
1061       }
1062       if ( labcur >= 0 && lmax < lcurr ) {
1063         lmax = lcurr;
1064         labmax = labcur;
1065       }
1066       lmax = 0;
1067       for ( int irow = firstRow; irow <= lastRow; irow++ ) {
1068         int ih = t.RowHit( irow );
1069         if ( ih < 0 ) continue;
1070         int index = firstSliceHit + slice.HitInputID( slice.Row( irow ), ih );
1071         AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
1072         if ( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax
1073            ) lmax++;
1074       }
1075       traLabels[itr] = labmax;
1076       traPurity[itr] = ( ( nHits > 0 ) ? double( lmax ) / double( nHits ) : 0 );
1077     }
1078   }
1079
1080   nRecTot += traN;
1081
1082   for ( int itr = 0; itr < traN; itr++ ) {
1083     if ( traPurity[itr] < .9 || traLabels[itr] < 0 || traLabels[itr] >= fNMCTracks ) {
1084       nGhost++;
1085       continue;
1086     }
1087
1088     AliHLTTPCCAMCTrack &mc = fMCTracks[traLabels[itr]];
1089     mc.SetNReconstructed( mc.NReconstructed() + 1 );
1090     if ( mc.Set() == 0 ) nRecOut++;
1091     else {
1092       if ( mc.NReconstructed() == 1 ) nRecAll++;
1093       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
1094       if ( mc.Set() == 2 ) {
1095         if ( mc.NReconstructed() == 1 ) nRecRef++;
1096         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
1097       }
1098     }
1099   }
1100
1101   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1102     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1103     if ( mc.Set() > 0 ) fhCandEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1104   }
1105
1106   if ( traLabels ) delete[] traLabels;
1107   if ( traPurity ) delete[] traPurity;
1108
1109   fStatCandNRecTot += nRecTot;
1110   fStatCandNRecOut += nRecOut;
1111   fStatCandNGhost  += nGhost;
1112   fStatCandNMCAll  += nMCAll;
1113   fStatCandNRecAll  += nRecAll;
1114   fStatCandNClonesAll  += nClonesAll;
1115   fStatCandNMCRef  += nMCRef;
1116   fStatCandNRecRef  += nRecRef;
1117   fStatCandNClonesRef  += nClonesRef;
1118
1119   if ( nMCAll == 0 ) return;
1120
1121   if ( PrintFlag ) {
1122     cout << "Track candidate performance for slice " << iSlice << " : " << endl;
1123     cout << " N tracks : "
1124          << nMCAll << " mc all, "
1125          << nMCRef << " mc ref, "
1126          << nRecTot << " rec total, "
1127          << nRecAll << " rec all, "
1128          << nClonesAll << " clones all, "
1129          << nRecRef << " rec ref, "
1130          << nClonesRef << " clones ref, "
1131          << nRecOut << " out, "
1132          << nGhost << " ghost" << endl;
1133
1134     int nRecExtr = nRecAll - nRecRef;
1135     int nMCExtr = nMCAll - nMCRef;
1136     int nClonesExtr = nClonesAll - nClonesRef;
1137
1138     double dRecTot = ( nRecTot > 0 ) ? nRecTot : 1;
1139     double dMCAll = ( nMCAll > 0 ) ? nMCAll : 1;
1140     double dMCRef = ( nMCRef > 0 ) ? nMCRef : 1;
1141     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1142     double dRecAll = ( nRecAll + nClonesAll > 0 ) ? nRecAll + nClonesAll : 1;
1143     double dRecRef = ( nRecRef + nClonesRef > 0 ) ? nRecRef + nClonesRef : 1;
1144     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1145
1146     cout << " EffRef = ";
1147     if ( nMCRef > 0 ) cout << nRecRef / dMCRef; else cout << "_";
1148     cout << ", CloneRef = ";
1149     if ( nRecRef > 0 ) cout << nClonesRef / dRecRef; else cout << "_";
1150     cout << endl;
1151     cout << " EffExtra = ";
1152     if ( nMCExtr > 0 ) cout << nRecExtr / dMCExtr; else cout << "_";
1153     cout << ", CloneExtra = ";
1154     if ( nRecExtr > 0 ) cout << nClonesExtr / dRecExtr; else cout << "_";
1155     cout << endl;
1156     cout << " EffAll = ";
1157     if ( nMCAll > 0 ) cout << nRecAll / dMCAll; else cout << "_";
1158     cout << ", CloneAll = ";
1159     if ( nRecAll > 0 ) cout << nClonesAll / dRecAll; else cout << "_";
1160     cout << endl;
1161     cout << " Out = ";
1162     if ( nRecTot > 0 ) cout << nRecOut / dRecTot; else cout << "_";
1163     cout << ", Ghost = ";
1164     if ( nRecTot > 0 ) cout << nGhost / dRecTot; else cout << "_";
1165     cout << endl;
1166   }
1167 #endif
1168 }
1169
1170
1171
1172 void AliHLTTPCCAPerformance::SlicePerformance( int iSlice, bool PrintFlag )
1173 {
1174   //* calculate slice tracker performance
1175
1176   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1177
1178   int nRecTot = 0, nGhost = 0, nRecOut = 0;
1179   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
1180   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
1181   const AliHLTTPCCATracker &tracker = hlt.SliceTracker( iSlice );
1182
1183   // Select reconstructable MC tracks
1184
1185   {
1186     for ( int imc = 0; imc < fNMCTracks; imc++ ) fMCTracks[imc].SetNHits( 0 );
1187
1188     for ( int ih = 0; ih < tracker.ClusterData()->NumberOfClusters(); ih++ ) {
1189       int id = tracker.ClusterData()->Id( ih );
1190       if ( id < 0 || id > fNHits ) break;
1191       AliHLTTPCCAHitLabel &l = fHitLabels[id];
1192       if ( l.fLab[0] >= 0 ) fMCTracks[l.fLab[0]].SetNHits( fMCTracks[l.fLab[0]].NHits() + 1 );
1193       if ( l.fLab[1] >= 0 ) fMCTracks[l.fLab[1]].SetNHits( fMCTracks[l.fLab[1]].NHits() + 1 );
1194       if ( l.fLab[2] >= 0 ) fMCTracks[l.fLab[2]].SetNHits( fMCTracks[l.fLab[2]].NHits() + 1 );
1195     }
1196
1197     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
1198       AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
1199       mc.SetSet( 0 );
1200       mc.SetNReconstructed( 0 );
1201       mc.SetNTurns( 1 );
1202       if ( mc.NHits() >=  30 && mc.P() >= .05 ) {
1203         mc.SetSet( 1 );
1204         nMCAll++;
1205         if ( mc.NHits() >=  30 && mc.P() >= 1. ) {
1206           mc.SetSet( 2 );
1207           nMCRef++;
1208         }
1209       }
1210     }
1211   }
1212
1213   if ( !tracker.Output() ) return;
1214
1215   const AliHLTTPCCASliceOutput &output = *tracker.Output();
1216
1217   int traN = output.NTracks();
1218
1219   nRecTot += traN;
1220
1221   for ( int itr = 0; itr < traN; itr++ ) {
1222
1223     const AliHLTTPCCASliceTrack &tCA = output.Track( itr );
1224     std::vector<int> clusterIDs;
1225     for ( int i = 0; i < tCA.NClusters(); i++ ) {
1226       clusterIDs.push_back( output.ClusterId( tCA.FirstClusterRef() + i ) );
1227     }
1228     int label;
1229     float purity;
1230     GetMCLabel( clusterIDs, label, purity );
1231
1232     if ( purity < .9 || label < 0 || label >= fNMCTracks ) {
1233       nGhost++;
1234       continue;
1235     }
1236
1237     AliHLTTPCCAMCTrack &mc = fMCTracks[label];
1238     mc.SetNReconstructed( mc.NReconstructed() + 1 );
1239     if ( mc.Set() == 0 ) nRecOut++;
1240     else {
1241       if ( mc.NReconstructed() == 1 ) nRecAll++;
1242       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
1243       if ( mc.Set() == 2 ) {
1244         if ( mc.NReconstructed() == 1 ) nRecRef++;
1245         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
1246       }
1247     }
1248
1249   }
1250
1251
1252   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1253     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1254     if ( mc.Set() > 0 ) fhEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1255   }
1256
1257
1258   fStatNRecTot += nRecTot;
1259   fStatNRecOut += nRecOut;
1260   fStatNGhost  += nGhost;
1261   fStatNMCAll  += nMCAll;
1262   fStatNRecAll  += nRecAll;
1263   fStatNClonesAll  += nClonesAll;
1264   fStatNMCRef  += nMCRef;
1265   fStatNRecRef  += nRecRef;
1266   fStatNClonesRef  += nClonesRef;
1267
1268   if ( nMCAll == 0 ) return;
1269
1270   if ( PrintFlag ) {
1271     cout << "Performance for slice " << iSlice << " : " << endl;
1272     cout << " N tracks : "
1273          << nMCAll << " mc all, "
1274          << nMCRef << " mc ref, "
1275          << nRecTot << " rec total, "
1276          << nRecAll << " rec all, "
1277          << nClonesAll << " clones all, "
1278          << nRecRef << " rec ref, "
1279          << nClonesRef << " clones ref, "
1280          << nRecOut << " out, "
1281          << nGhost << " ghost" << endl;
1282
1283     int nRecExtr = nRecAll - nRecRef;
1284     int nMCExtr = nMCAll - nMCRef;
1285     int nClonesExtr = nClonesAll - nClonesRef;
1286
1287     double dRecTot = ( nRecTot > 0 ) ? nRecTot : 1;
1288     double dMCAll = ( nMCAll > 0 ) ? nMCAll : 1;
1289     double dMCRef = ( nMCRef > 0 ) ? nMCRef : 1;
1290     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1291     double dRecAll = ( nRecAll + nClonesAll > 0 ) ? nRecAll + nClonesAll : 1;
1292     double dRecRef = ( nRecRef + nClonesRef > 0 ) ? nRecRef + nClonesRef : 1;
1293     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1294
1295     cout << " EffRef = ";
1296     if ( nMCRef > 0 ) cout << nRecRef / dMCRef; else cout << "_";
1297     cout << ", CloneRef = ";
1298     if ( nRecRef > 0 ) cout << nClonesRef / dRecRef; else cout << "_";
1299     cout << endl;
1300     cout << " EffExtra = ";
1301     if ( nMCExtr > 0 ) cout << nRecExtr / dMCExtr; else cout << "_";
1302     cout << ", CloneExtra = ";
1303     if ( nRecExtr > 0 ) cout << nClonesExtr / dRecExtr; else cout << "_";
1304     cout << endl;
1305     cout << " EffAll = ";
1306     if ( nMCAll > 0 ) cout << nRecAll / dMCAll; else cout << "_";
1307     cout << ", CloneAll = ";
1308     if ( nRecAll > 0 ) cout << nClonesAll / dRecAll; else cout << "_";
1309     cout << endl;
1310     cout << " Out = ";
1311     if ( nRecTot > 0 ) cout << nRecOut / dRecTot; else cout << "_";
1312     cout << ", Ghost = ";
1313     if ( nRecTot > 0 ) cout << nGhost / dRecTot; else cout << "_";
1314     cout << endl;
1315   }
1316 }
1317
1318
1319
1320 void AliHLTTPCCAPerformance::MergerPerformance()
1321 {
1322   // performance calculation for merged tracks
1323
1324   int nRecTot = 0, nGhost = 0, nRecOut = 0;
1325   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
1326   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
1327
1328   // Select reconstructable MC tracks
1329
1330   {
1331     for ( int imc = 0; imc < fNMCTracks; imc++ ) fMCTracks[imc].SetNHits( 0 );
1332
1333     for ( int ih = 0; ih < fNHits; ih++ ) {
1334       AliHLTTPCCAHitLabel &l = fHitLabels[ih];
1335       if ( l.fLab[0] >= 0 ) fMCTracks[l.fLab[0]].SetNHits( fMCTracks[l.fLab[0]].NHits() + 1 );
1336       if ( l.fLab[1] >= 0 ) fMCTracks[l.fLab[1]].SetNHits( fMCTracks[l.fLab[1]].NHits() + 1 );
1337       if ( l.fLab[2] >= 0 ) fMCTracks[l.fLab[2]].SetNHits( fMCTracks[l.fLab[2]].NHits() + 1 );
1338     }
1339
1340     for ( int imc = 0; imc < fNMCTracks; imc++ ) {
1341       AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
1342       mc.SetSet( 0 );
1343       mc.SetNReconstructed( 0 );
1344       mc.SetNTurns( 1 );
1345       if ( mc.NHits() >=  50 && mc.P() >= .05 ) {
1346         mc.SetSet( 1 );
1347         nMCAll++;
1348         if ( mc.P() >= 1. ) {
1349           mc.SetSet( 2 );
1350           nMCRef++;
1351         }
1352       }
1353     }
1354   }
1355
1356   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1357
1358   if ( !hlt.Merger().Output() ) return;
1359
1360   const AliHLTTPCCAMergerOutput &output = *( hlt.Merger().Output() );
1361
1362   int traN = output.NTracks();
1363
1364   nRecTot += traN;
1365
1366   for ( int itr = 0; itr < traN; itr++ ) {
1367
1368     const AliHLTTPCCAMergedTrack &tCA = output.Track( itr );
1369     std::vector<int> clusterIDs;
1370     for ( int i = 0; i < tCA.NClusters(); i++ ) {
1371       clusterIDs.push_back( output.ClusterId( tCA.FirstClusterRef() + i ) );
1372     }
1373     int label;
1374     float purity;
1375     GetMCLabel( clusterIDs, label, purity );
1376
1377     if ( purity < .9 || label < 0 || label >= fNMCTracks ) {
1378       nGhost++;
1379       continue;
1380     }
1381
1382     AliHLTTPCCAMCTrack &mc = fMCTracks[label];
1383     mc.SetNReconstructed( mc.NReconstructed() + 1 );
1384     if ( mc.Set() == 0 ) nRecOut++;
1385     else {
1386       if ( mc.NReconstructed() == 1 ) nRecAll++;
1387       else if ( mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
1388       if ( mc.Set() == 2 ) {
1389         if ( mc.NReconstructed() == 1 ) nRecRef++;
1390         else if ( mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
1391         fhTrackLengthRef->Fill( tCA.NClusters() / ( ( double ) mc.NHits() ) );
1392       }
1393     }
1394
1395     // track resolutions
1396     while ( mc.Set() == 2 && TMath::Abs( mc.TPCPar()[0] ) + TMath::Abs( mc.TPCPar()[1] ) > 1 ) {
1397
1398       if ( purity < .90 ) break;
1399       AliHLTTPCCATrackParam p = tCA.InnerParam();
1400       double cosA = TMath::Cos( tCA.InnerAlpha() );
1401       double sinA = TMath::Sin( tCA.InnerAlpha() );
1402       double mcX =  mc.TPCPar()[0] * cosA + mc.TPCPar()[1] * sinA;
1403       double mcY = -mc.TPCPar()[0] * sinA + mc.TPCPar()[1] * cosA;
1404       double mcZ =  mc.TPCPar()[2];
1405       double mcEx =  mc.TPCPar()[3] * cosA + mc.TPCPar()[4] * sinA;
1406       double mcEy = -mc.TPCPar()[3] * sinA + mc.TPCPar()[4] * cosA;
1407       double mcEz =  mc.TPCPar()[5];
1408       double mcEt = TMath::Sqrt( mcEx * mcEx + mcEy * mcEy );
1409       if ( TMath::Abs( mcEt ) < 1.e-4 ) break;
1410       double mcSinPhi = mcEy / mcEt;
1411       double mcDzDs   = mcEz / mcEt;
1412       double mcQPt = mc.TPCPar()[6] / mcEt;
1413       if ( TMath::Abs( mcQPt ) < 1.e-4 ) break;
1414       double mcPt = 1. / TMath::Abs( mcQPt );
1415
1416       if ( mcPt < 1. ) break;
1417
1418       if ( tCA.NClusters() <  50 ) break;
1419       if ( !p.TransportToXWithMaterial( mcX, hlt.Merger().SliceParam().GetBz( p ) ) ) break;
1420       if ( p.GetCosPhi()*mcEx < 0 ) { // change direction
1421         mcSinPhi = -mcSinPhi;
1422         mcDzDs = -mcDzDs;
1423         mcQPt = -mcQPt;
1424       }
1425
1426       double qPt = p.GetQPt();
1427       double pt = 100;
1428       if ( TMath::Abs( qPt ) > 1.e-4 ) pt = 1. / TMath::Abs( qPt );
1429
1430       fhResY->Fill( p.GetY() - mcY );
1431       fhResZ->Fill( p.GetZ() - mcZ );
1432       fhResSinPhi->Fill( p.GetSinPhi() - mcSinPhi );
1433       fhResDzDs->Fill( p.GetDzDs() - mcDzDs );
1434       fhResPt->Fill( ( pt - mcPt ) / mcPt );
1435
1436       if ( p.GetErr2Y() > 0 ) fhPullY->Fill( ( p.GetY() - mcY ) / TMath::Sqrt( p.GetErr2Y() ) );
1437       if ( p.GetErr2Z() > 0 ) fhPullZ->Fill( ( p.GetZ() - mcZ ) / TMath::Sqrt( p.GetErr2Z() ) );
1438
1439       if ( p.GetErr2SinPhi() > 0 ) fhPullSinPhi->Fill( ( p.GetSinPhi() - mcSinPhi ) / TMath::Sqrt( p.GetErr2SinPhi() ) );
1440       if ( p.GetErr2DzDs() > 0 ) fhPullDzDs->Fill( ( p.DzDs() - mcDzDs ) / TMath::Sqrt( p.GetErr2DzDs() ) );
1441       if ( p.GetErr2QPt() > 0 ) fhPullQPt->Fill( ( qPt - mcQPt ) / TMath::Sqrt( p.GetErr2QPt() ) );
1442       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 ) ) );
1443       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 ) ) );
1444
1445       break;
1446     } // end resolutions
1447
1448   }// end reco tracks
1449
1450
1451   for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1452     AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1453     if ( mc.Set() > 0 ) fhGBEffVsP->Fill( mc.P(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1454     if ( mc.Set() > 0 ) fhGBEffVsPt->Fill( mc.Pt(), ( mc.NReconstructed() > 0 ? 1 : 0 ) );
1455     if ( mc.Set() == 2 ) {
1456       const double *p = mc.TPCPar();
1457       double r = TMath::Sqrt( p[0] * p[0] + p[1] * p[1] );
1458       double cosA = p[0] / r;
1459       double sinA = p[1] / r;
1460
1461
1462       double phipos = TMath::Pi() + TMath::ATan2( -p[1], -p[0] );
1463       double alpha =  TMath::Pi() * ( 20 * ( ( ( ( int )( phipos * 180 / TMath::Pi() ) ) / 20 ) ) + 10 ) / 180.;
1464       cosA = TMath::Cos( alpha );
1465       sinA = TMath::Sin( alpha );
1466
1467       double mcX =  p[0] * cosA + p[1] * sinA;
1468       double mcY = -p[0] * sinA + p[1] * cosA;
1469       double mcZ =  p[2];
1470       double mcEx =  p[3] * cosA + p[4] * sinA;
1471       double mcEy = -p[3] * sinA + p[4] * cosA;
1472       double mcEz =  p[5];
1473       //double mcEt = TMath::Sqrt(mcEx*mcEx + mcEy*mcEy);
1474       double angleY = TMath::ATan2( mcEy, mcEx ) * 180. / TMath::Pi();
1475       double angleZ = TMath::ATan2( mcEz, mcEx ) * 180. / TMath::Pi();
1476
1477       if ( mc.NReconstructed() > 0 ) {
1478         fhRefRecoX->Fill( mcX );
1479         fhRefRecoY->Fill( mcY );
1480         fhRefRecoZ->Fill( mcZ );
1481         fhRefRecoP->Fill( mc.P() );
1482         fhRefRecoPt->Fill( mc.Pt() );
1483         fhRefRecoAngleY->Fill( angleY );
1484         fhRefRecoAngleZ->Fill( angleZ );
1485         fhRefRecoNHits->Fill( mc.NHits() );
1486       } else {
1487         fhRefNotRecoX->Fill( mcX );
1488         fhRefNotRecoY->Fill( mcY );
1489         fhRefNotRecoZ->Fill( mcZ );
1490         fhRefNotRecoP->Fill( mc.P() );
1491         fhRefNotRecoPt->Fill( mc.Pt() );
1492         fhRefNotRecoAngleY->Fill( angleY );
1493         fhRefNotRecoAngleZ->Fill( angleZ );
1494         fhRefNotRecoNHits->Fill( mc.NHits() );
1495       }
1496     }
1497   }
1498
1499   fStatGBNRecTot += nRecTot;
1500   fStatGBNRecOut += nRecOut;
1501   fStatGBNGhost  += nGhost;
1502   fStatGBNMCAll  += nMCAll;
1503   fStatGBNRecAll  += nRecAll;
1504   fStatGBNClonesAll  += nClonesAll;
1505   fStatGBNMCRef  += nMCRef;
1506   fStatGBNRecRef  += nRecRef;
1507   fStatGBNClonesRef  += nClonesRef;
1508 }
1509
1510
1511
1512 void AliHLTTPCCAPerformance::ClusterPerformance()
1513 {
1514   // performance calculation for input clusters
1515
1516   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1517
1518   // distribution of cluster errors
1519
1520   for ( int iSlice = 0; iSlice < hlt.NSlices(); iSlice++ ) {
1521     const AliHLTTPCCAClusterData &data = hlt.ClusterData( iSlice );
1522     for ( int i = 0; i < data.NumberOfClusters(); i++ ) {
1523       AliHLTTPCCAHitLabel &l = fHitLabels[data.Id( i )];
1524       int nmc = 0;
1525       for ( int il = 0; il < 3; il++ ) if ( l.fLab[il] >= 0 ) nmc++;
1526       if ( nmc == 1 ) fhHitShared->Fill( data.RowNumber( i ), 0 );
1527       else if ( nmc > 1 ) fhHitShared->Fill( data.RowNumber( i ), 1 );
1528     }
1529   }
1530
1531   // cluster pulls
1532
1533   if ( !fDoClusterPulls || fNMCPoints <= 0 ) return;
1534
1535   // sort mc points
1536   if ( 1 ) {
1537     for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1538       AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1539       mc.SetNMCPoints( 0 );
1540     }
1541     sort( fMCPoints, fMCPoints + fNMCPoints, AliHLTTPCCAMCPoint::Compare );
1542
1543     for ( int ip = 0; ip < fNMCPoints; ip++ ) {
1544       AliHLTTPCCAMCPoint &p = fMCPoints[ip];
1545       AliHLTTPCCAMCTrack &t = fMCTracks[p.TrackID()];
1546       if ( t.NMCPoints() == 0 ) t.SetFirstMCPointID( ip );
1547       t.SetNMCPoints( t.NMCPoints() + 1 );
1548     }
1549   }
1550
1551   for ( int iSlice = 0; iSlice < hlt.NSlices(); iSlice++ ) {
1552
1553     const AliHLTTPCCAClusterData &data = hlt.ClusterData( iSlice );
1554
1555     for ( int ic = 0; ic < data.NumberOfClusters(); ic++ ) {
1556
1557       const AliHLTTPCCAHitLabel &l = fHitLabels[data.Id( ic )];
1558
1559       if ( l.fLab[0] < 0 || l.fLab[0] >= fNMCTracks
1560            || l.fLab[1] >= 0 || l.fLab[2] >= 0       ) continue;
1561
1562       int lab = l.fLab[0];
1563
1564       AliHLTTPCCAMCTrack &mc = fMCTracks[lab];
1565
1566       double x0 = data.X( ic );
1567       double y0 = data.Y( ic );
1568       double z0 = data.Z( ic );
1569
1570       if ( fabs( x0 ) < 1.e-4 ) continue;
1571       if ( mc.Pt() < .05 ) continue;
1572
1573       int ip1 = -1, ip2 = -1;
1574       double d1 = 1.e20, d2 = 1.e20;
1575
1576       AliHLTTPCCAMCPoint *pStart = lower_bound( fMCPoints + mc.FirstMCPointID(), fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), iSlice,  AliHLTTPCCAMCPoint::CompareSlice );
1577
1578       pStart = lower_bound( pStart, fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), x0 - 2.,  AliHLTTPCCAMCPoint::CompareX );
1579
1580       for ( int ip = ( pStart - fMCPoints ) - mc.FirstMCPointID(); ip < mc.NMCPoints(); ip++ ) {
1581         AliHLTTPCCAMCPoint &p = fMCPoints[mc.FirstMCPointID() + ip];
1582         if ( p.ISlice() != iSlice ) break;
1583         double dx = p.Sx() - x0;
1584         double dy = p.Sy() - y0;
1585         double dz = p.Sz() - z0;
1586         double d = dx * dx + dy * dy + dz * dz;
1587         if ( d > 9. ) continue;
1588         if ( dx <= 0  && dx > -2. ) {
1589           if ( fabs( dx ) < d1 ) {
1590             d1 = fabs( dx );
1591             ip1 = ip;
1592           }
1593         } else if ( dx > .2 ) {
1594           if ( dx >= 2. ) break;
1595           if ( fabs( dx ) < d2 ) {
1596             d2 = fabs( dx );
1597             ip2 = ip;
1598           }
1599         }
1600       }
1601
1602       if ( ip1 < 0 || ip2 < 0 ) continue;
1603
1604       AliHLTTPCCAMCPoint &p1 = fMCPoints[mc.FirstMCPointID() + ip1];
1605       AliHLTTPCCAMCPoint &p2 = fMCPoints[mc.FirstMCPointID() + ip2];
1606       double dx = p2.Sx() - p1.Sx();
1607       double dy = p2.Sy() - p1.Sy();
1608       double dz = p2.Sz() - p1.Sz();
1609       double sx = x0;
1610       double sy = p1.Sy() + dy / dx * ( sx - p1.Sx() );
1611       double sz = p1.Sz() + dz / dx * ( sx - p1.Sx() );
1612
1613       float errY, errZ;
1614       {
1615         AliHLTTPCCATrackParam t;
1616         double s = 1. / TMath::Sqrt( dx * dx + dy * dy );
1617         t.SetZ( sz );
1618         t.SetSinPhi( dy * s );
1619         t.SetSignCosPhi( dx );
1620         t.SetDzDs( dz * s );
1621         hlt.SliceTracker( 0 ).GetErrors2( data.RowNumber( ic ), t, errY, errZ );
1622         errY = TMath::Sqrt( errY );
1623         errZ = TMath::Sqrt( errZ );
1624       }
1625       fhHitErrY->Fill( errY );
1626       fhHitErrZ->Fill( errZ );
1627       fhHitResY->Fill( y0 - sy );
1628       fhHitResZ->Fill( z0 - sz );
1629       fhHitPullY->Fill( ( y0 - sy ) / errY );
1630       fhHitPullZ->Fill( ( z0 - sz ) / errZ );
1631       if ( mc.Pt() >= 1. ) {
1632         fhHitResY1->Fill( y0 - sy );
1633         fhHitResZ1->Fill( z0 - sz );
1634         fhHitPullY1->Fill( ( y0 - sy ) / errY );
1635         fhHitPullZ1->Fill( ( z0 - sz ) / errZ );
1636       }
1637     }
1638   }
1639 }
1640
1641
1642 void AliHLTTPCCAPerformance::SmearClustersMC()
1643 {
1644   // smear clusters with gaussian using MC info
1645
1646   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1647
1648   // cluster pulls
1649
1650   if ( fNMCPoints <= 0 ) return;
1651
1652   // sort mc points
1653   {
1654     for ( int ipart = 0; ipart < fNMCTracks; ipart++ ) {
1655       AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
1656       mc.SetNMCPoints( 0 );
1657     }
1658     sort( fMCPoints, fMCPoints + fNMCPoints, AliHLTTPCCAMCPoint::Compare );
1659
1660     for ( int ip = 0; ip < fNMCPoints; ip++ ) {
1661       AliHLTTPCCAMCPoint &p = fMCPoints[ip];
1662       AliHLTTPCCAMCTrack &t = fMCTracks[p.TrackID()];
1663       if ( t.NMCPoints() == 0 ) t.SetFirstMCPointID( ip );
1664       t.SetNMCPoints( t.NMCPoints() + 1 );
1665     }
1666   }
1667
1668   for ( int iSlice = 0; iSlice < hlt.NSlices(); iSlice++ ) {
1669
1670     AliHLTTPCCAClusterData &data = hlt.ClusterData( iSlice );
1671
1672     for ( int ic = 0; ic < data.NumberOfClusters(); ic++ ) {
1673
1674       double x0 = data.X( ic );
1675       double y0 = data.Y( ic );
1676       double z0 = data.Z( ic );
1677       int row0 = data.RowNumber( ic );
1678
1679       AliHLTTPCCAClusterData::Data *cdata = data.GetClusterData( ic );
1680       cdata->fX = 0;
1681       cdata->fY = 0;
1682       cdata->fZ = 0;
1683
1684       const AliHLTTPCCAHitLabel &l = fHitLabels[data.Id( ic )];
1685
1686       if ( l.fLab[0] < 0 || l.fLab[0] >= fNMCTracks ) continue;
1687
1688       int lab = l.fLab[0];
1689
1690       AliHLTTPCCAMCTrack &mc = fMCTracks[lab];
1691
1692       int ip1 = -1, ip2 = -1;
1693       double d1 = 1.e20, d2 = 1.e20;
1694
1695       AliHLTTPCCAMCPoint *pStart = lower_bound( fMCPoints + mc.FirstMCPointID(), fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), iSlice,  AliHLTTPCCAMCPoint::CompareSlice );
1696
1697       pStart = lower_bound( pStart, fMCPoints + mc.FirstMCPointID() + mc.NMCPoints(), x0 - 2.,  AliHLTTPCCAMCPoint::CompareX );
1698
1699       for ( int ip = ( pStart - fMCPoints ) - mc.FirstMCPointID(); ip < mc.NMCPoints(); ip++ ) {
1700         AliHLTTPCCAMCPoint &p = fMCPoints[mc.FirstMCPointID() + ip];
1701         if ( p.ISlice() != iSlice ) break;
1702         double dx = p.Sx() - x0;
1703         double dy = p.Sy() - y0;
1704         double dz = p.Sz() - z0;
1705         double d = dx * dx + dy * dy + dz * dz;
1706         if ( d > 9. ) continue;
1707         if ( dx <= 0  && dx > -2. ) {
1708           if ( fabs( dx ) < d1 ) {
1709             d1 = fabs( dx );
1710             ip1 = ip;
1711           }
1712         } else if ( dx > .2 ) {
1713           if ( dx >= 2. ) break;
1714           if ( fabs( dx ) < d2 ) {
1715             d2 = fabs( dx );
1716             ip2 = ip;
1717           }
1718         }
1719       }
1720
1721       if ( ip1 < 0 || ip2 < 0 ) continue;
1722
1723       AliHLTTPCCAMCPoint &p1 = fMCPoints[mc.FirstMCPointID() + ip1];
1724       AliHLTTPCCAMCPoint &p2 = fMCPoints[mc.FirstMCPointID() + ip2];
1725       double dx = p2.Sx() - p1.Sx();
1726       double dy = p2.Sy() - p1.Sy();
1727       double dz = p2.Sz() - p1.Sz();
1728       double sx = x0;
1729       double sy = p1.Sy() + dy / dx * ( sx - p1.Sx() );
1730       double sz = p1.Sz() + dz / dx * ( sx - p1.Sx() );
1731
1732       float errY, errZ;
1733       {
1734         AliHLTTPCCATrackParam t;
1735         double s = 1. / TMath::Sqrt( dx * dx + dy * dy );
1736         t.SetZ( sz );
1737         t.SetSinPhi( dy * s );
1738         t.SetSignCosPhi( dx );
1739         t.SetDzDs( dz * s );
1740         hlt.SliceTracker( 0 ).GetErrors2( row0, t, errY, errZ );
1741         errY = TMath::Sqrt( errY );
1742         errZ = TMath::Sqrt( errZ );
1743       }
1744
1745       cdata->fX = x0;
1746       cdata->fY = gRandom->Gaus( sy, errY );
1747       cdata->fZ = gRandom->Gaus( sz, errZ );
1748     }
1749   }
1750 }
1751
1752
1753 void AliHLTTPCCAPerformance::Performance( fstream *StatFile )
1754 {
1755   // main routine for performance calculation
1756
1757   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
1758
1759   //SG!!!
1760   /*
1761   fStatNEvents=0;
1762     fStatNRecTot=0;
1763     fStatNRecOut=0;
1764     fStatNGhost=0;
1765     fStatNMCAll=0;
1766     fStatNRecAll=0;
1767     fStatNClonesAll=0;
1768     fStatNMCRef=0;
1769     fStatNRecRef=0;
1770     fStatNClonesRef=0;
1771   */
1772   fStatNEvents++;
1773   for ( int islice = 0; islice < hlt.NSlices(); islice++ ) {
1774     SliceTrackletPerformance( islice, 0 );
1775     SliceTrackCandPerformance( islice, 0 );
1776     SlicePerformance( islice, 0 );
1777   }
1778
1779   MergerPerformance();
1780   //ClusterPerformance();
1781
1782   {
1783     cout << "\nSlice Track Seed performance: \n" << endl;
1784     cout << " N tracks : "
1785          << fStatNMCAll / fStatNEvents << " mc all, "
1786          << fStatSeedNMCRef / fStatNEvents << " mc ref, "
1787          << fStatSeedNRecTot / fStatNEvents << " rec total, "
1788          << fStatSeedNRecAll / fStatNEvents << " rec all, "
1789          << fStatSeedNClonesAll / fStatNEvents << " clones all, "
1790          << fStatSeedNRecRef / fStatNEvents << " rec ref, "
1791          << fStatSeedNClonesRef / fStatNEvents << " clones ref, "
1792          << fStatSeedNRecOut / fStatNEvents << " out, "
1793          << fStatSeedNGhost / fStatNEvents << " ghost" << endl;
1794
1795     int nRecExtr = fStatSeedNRecAll - fStatSeedNRecRef;
1796     int nMCExtr = fStatNMCAll - fStatNMCRef;
1797     int nClonesExtr = fStatSeedNClonesAll - fStatSeedNClonesRef;
1798
1799     double dRecTot = ( fStatSeedNRecTot > 0 ) ? fStatSeedNRecTot : 1;
1800     double dMCAll = ( fStatNMCAll > 0 ) ? fStatNMCAll : 1;
1801     double dMCRef = ( fStatNMCRef > 0 ) ? fStatNMCRef : 1;
1802     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1803     double dRecAll = ( fStatSeedNRecAll + fStatSeedNClonesAll > 0 ) ? fStatSeedNRecAll + fStatSeedNClonesAll : 1;
1804     double dRecRef = ( fStatSeedNRecRef + fStatSeedNClonesRef > 0 ) ? fStatSeedNRecRef + fStatSeedNClonesRef : 1;
1805     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1806
1807     cout << " EffRef = " << fStatSeedNRecRef / dMCRef
1808          << ", CloneRef = " << fStatSeedNClonesRef / dRecRef << endl;
1809     cout << " EffExtra = " << nRecExtr / dMCExtr
1810          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1811     cout << " EffAll = " << fStatSeedNRecAll / dMCAll
1812          << ", CloneAll = " << fStatSeedNClonesAll / dRecAll << endl;
1813     cout << " Out = " << fStatSeedNRecOut / dRecTot
1814          << ", Ghost = " << fStatSeedNGhost / dRecTot << endl;
1815   }
1816
1817   {
1818     cout << "\nSlice Track candidate performance: \n" << endl;
1819     cout << " N tracks : "
1820          << fStatNMCAll / fStatNEvents << " mc all, "
1821          << fStatCandNMCRef / fStatNEvents << " mc ref, "
1822          << fStatCandNRecTot / fStatNEvents << " rec total, "
1823          << fStatCandNRecAll / fStatNEvents << " rec all, "
1824          << fStatCandNClonesAll / fStatNEvents << " clones all, "
1825          << fStatCandNRecRef / fStatNEvents << " rec ref, "
1826          << fStatCandNClonesRef / fStatNEvents << " clones ref, "
1827          << fStatCandNRecOut / fStatNEvents << " out, "
1828          << fStatCandNGhost / fStatNEvents << " ghost" << endl;
1829
1830     int nRecExtr = fStatCandNRecAll - fStatCandNRecRef;
1831     int nMCExtr = fStatNMCAll - fStatNMCRef;
1832     int nClonesExtr = fStatCandNClonesAll - fStatCandNClonesRef;
1833
1834     double dRecTot = ( fStatCandNRecTot > 0 ) ? fStatCandNRecTot : 1;
1835     double dMCAll = ( fStatNMCAll > 0 ) ? fStatNMCAll : 1;
1836     double dMCRef = ( fStatNMCRef > 0 ) ? fStatNMCRef : 1;
1837     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1838     double dRecAll = ( fStatCandNRecAll + fStatCandNClonesAll > 0 ) ? fStatCandNRecAll + fStatCandNClonesAll : 1;
1839     double dRecRef = ( fStatCandNRecRef + fStatCandNClonesRef > 0 ) ? fStatCandNRecRef + fStatCandNClonesRef : 1;
1840     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1841
1842     cout << " EffRef = " << fStatCandNRecRef / dMCRef
1843          << ", CloneRef = " << fStatCandNClonesRef / dRecRef << endl;
1844     cout << " EffExtra = " << nRecExtr / dMCExtr
1845          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1846     cout << " EffAll = " << fStatCandNRecAll / dMCAll
1847          << ", CloneAll = " << fStatCandNClonesAll / dRecAll << endl;
1848     cout << " Out = " << fStatCandNRecOut / dRecTot
1849          << ", Ghost = " << fStatCandNGhost / dRecTot << endl;
1850   }
1851
1852   {
1853     cout << "\nSlice tracker performance: \n" << endl;
1854     cout << " N tracks : "
1855          << fStatNMCAll / fStatNEvents << " mc all, "
1856          << fStatNMCRef / fStatNEvents << " mc ref, "
1857          << fStatNRecTot / fStatNEvents << " rec total, "
1858          << fStatNRecAll / fStatNEvents << " rec all, "
1859          << fStatNClonesAll / fStatNEvents << " clones all, "
1860          << fStatNRecRef / fStatNEvents << " rec ref, "
1861          << fStatNClonesRef / fStatNEvents << " clones ref, "
1862          << fStatNRecOut / fStatNEvents << " out, "
1863          << fStatNGhost / fStatNEvents << " ghost" << endl;
1864
1865     int nRecExtr = fStatNRecAll - fStatNRecRef;
1866     int nMCExtr = fStatNMCAll - fStatNMCRef;
1867     int nClonesExtr = fStatNClonesAll - fStatNClonesRef;
1868
1869     double dRecTot = ( fStatNRecTot > 0 ) ? fStatNRecTot : 1;
1870     double dMCAll = ( fStatNMCAll > 0 ) ? fStatNMCAll : 1;
1871     double dMCRef = ( fStatNMCRef > 0 ) ? fStatNMCRef : 1;
1872     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1873     double dRecAll = ( fStatNRecAll + fStatNClonesAll > 0 ) ? fStatNRecAll + fStatNClonesAll : 1;
1874     double dRecRef = ( fStatNRecRef + fStatNClonesRef > 0 ) ? fStatNRecRef + fStatNClonesRef : 1;
1875     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1876
1877     cout << " EffRef = " << fStatNRecRef / dMCRef
1878          << ", CloneRef = " << fStatNClonesRef / dRecRef << endl;
1879     cout << " EffExtra = " << nRecExtr / dMCExtr
1880          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1881     cout << " EffAll = " << fStatNRecAll / dMCAll
1882          << ", CloneAll = " << fStatNClonesAll / dRecAll << endl;
1883     cout << " Out = " << fStatNRecOut / dRecTot
1884          << ", Ghost = " << fStatNGhost / dRecTot << endl;
1885     cout << " Time = " << hlt.StatTime( 0 ) / hlt.StatNEvents()*1.e3 << " msec/event " << endl;
1886     cout << " Local timers = "
1887          << hlt.StatTime( 1 ) / hlt.StatNEvents()*1.e3 << " "
1888          << hlt.StatTime( 2 ) / hlt.StatNEvents()*1.e3 << " "
1889          << hlt.StatTime( 3 ) / hlt.StatNEvents()*1.e3 << " "
1890          << hlt.StatTime( 4 ) / hlt.StatNEvents()*1.e3 << " "
1891          << hlt.StatTime( 5 ) / hlt.StatNEvents()*1.e3 << " "
1892          << hlt.StatTime( 6 ) / hlt.StatNEvents()*1.e3 << " "
1893          << hlt.StatTime( 7 ) / hlt.StatNEvents()*1.e3 << " "
1894          << hlt.StatTime( 8 ) / hlt.StatNEvents()*1.e3 << " "
1895          << " msec/event " << endl;
1896   }
1897
1898
1899   {
1900     cout << "\nGlobal tracker performance for " << fStatNEvents << " events: \n" << endl;
1901     cout << " N tracks : "
1902          << fStatGBNMCAll << " mc all, "
1903          << fStatGBNMCRef << " mc ref, "
1904          << fStatGBNRecTot << " rec total, "
1905          << fStatGBNRecAll << " rec all, "
1906          << fStatGBNClonesAll << " clones all, "
1907          << fStatGBNRecRef << " rec ref, "
1908          << fStatGBNClonesRef << " clones ref, "
1909          << fStatGBNRecOut << " out, "
1910          << fStatGBNGhost << " ghost" << endl;
1911     cout << " N tracks average : "
1912          << fStatGBNMCAll / fStatNEvents << " mc all, "
1913          << fStatGBNMCRef / fStatNEvents << " mc ref, "
1914          << fStatGBNRecTot / fStatNEvents << " rec total, "
1915          << fStatGBNRecAll / fStatNEvents << " rec all, "
1916          << fStatGBNClonesAll / fStatNEvents << " clones all, "
1917          << fStatGBNRecRef / fStatNEvents << " rec ref, "
1918          << fStatGBNClonesRef / fStatNEvents << " clones ref, "
1919          << fStatGBNRecOut / fStatNEvents << " out, "
1920          << fStatGBNGhost / fStatNEvents << " ghost" << endl;
1921
1922     int nRecExtr = fStatGBNRecAll - fStatGBNRecRef;
1923     int nMCExtr = fStatGBNMCAll - fStatGBNMCRef;
1924     int nClonesExtr = fStatGBNClonesAll - fStatGBNClonesRef;
1925
1926     double dRecTot = ( fStatGBNRecTot > 0 ) ? fStatGBNRecTot : 1;
1927     double dMCAll = ( fStatGBNMCAll > 0 ) ? fStatGBNMCAll : 1;
1928     double dMCRef = ( fStatGBNMCRef > 0 ) ? fStatGBNMCRef : 1;
1929     double dMCExtr = ( nMCExtr > 0 ) ? nMCExtr : 1;
1930     double dRecAll = ( fStatGBNRecAll + fStatGBNClonesAll > 0 ) ? fStatGBNRecAll + fStatGBNClonesAll : 1;
1931     double dRecRef = ( fStatGBNRecRef + fStatGBNClonesRef > 0 ) ? fStatGBNRecRef + fStatGBNClonesRef : 1;
1932     double dRecExtr = ( nRecExtr + nClonesExtr > 0 ) ? nRecExtr + nClonesExtr : 1;
1933
1934     cout << " EffRef = " << fStatGBNRecRef / dMCRef
1935          << ", CloneRef = " << fStatGBNClonesRef / dRecRef << endl;
1936     cout << " EffExtra = " << nRecExtr / dMCExtr
1937          << ", CloneExtra = " << nClonesExtr / dRecExtr << endl;
1938     cout << " EffAll = " << fStatGBNRecAll / dMCAll
1939          << ", CloneAll = " << fStatGBNClonesAll / dRecAll << endl;
1940     cout << " Out = " << fStatGBNRecOut / dRecTot
1941          << ", Ghost = " << fStatGBNGhost / dRecTot << endl;
1942     cout << " Time = " << ( hlt.StatTime( 0 ) + hlt.StatTime( 9 ) ) / hlt.StatNEvents()*1.e3 << " msec/event " << endl;
1943     cout << " Local timers: " << endl;
1944     cout << " slice tracker " << hlt.StatTime( 0 ) / hlt.StatNEvents()*1.e3 << ": "
1945          << hlt.StatTime( 1 ) / hlt.StatNEvents()*1.e3 << " "
1946          << hlt.StatTime( 2 ) / hlt.StatNEvents()*1.e3 << " "
1947          << hlt.StatTime( 3 ) / hlt.StatNEvents()*1.e3 << " "
1948          << hlt.StatTime( 4 ) / hlt.StatNEvents()*1.e3 << " "
1949          << hlt.StatTime( 5 ) / hlt.StatNEvents()*1.e3 << "["
1950          << hlt.StatTime( 6 ) / hlt.StatNEvents()*1.e3 << "/"
1951          << hlt.StatTime( 7 ) / hlt.StatNEvents()*1.e3 << "] "
1952          << hlt.StatTime( 8 ) / hlt.StatNEvents()*1.e3
1953          << " msec/event " << endl;
1954     cout << " GB merger " << hlt.StatTime( 9 ) / hlt.StatNEvents()*1.e3 << ": "
1955          << hlt.StatTime( 10 ) / hlt.StatNEvents()*1.e3 << ", "
1956          << hlt.StatTime( 11 ) / hlt.StatNEvents()*1.e3 << ", "
1957          << hlt.StatTime( 12 ) / hlt.StatNEvents()*1.e3 << " "
1958          << " msec/event " << endl;
1959
1960     if ( StatFile && StatFile->is_open() ) {
1961       fstream &out = *StatFile;
1962
1963       //out<<"\nGlobal tracker performance for "<<fStatNEvents<<" events: \n"<<endl;
1964       //out<<" N tracks : "
1965       //<<fStatGBNMCAll/fStatNEvents<<" mc all, "
1966       //<<fStatGBNMCRef/fStatNEvents<<" mc ref, "
1967       // <<fStatGBNRecTot/fStatNEvents<<" rec total, "
1968       // <<fStatGBNRecAll/fStatNEvents<<" rec all, "
1969       // <<fStatGBNClonesAll/fStatNEvents<<" clones all, "
1970       // <<fStatGBNRecRef/fStatNEvents<<" rec ref, "
1971       // <<fStatGBNClonesRef/fStatNEvents<<" clones ref, "
1972       // <<fStatGBNRecOut/fStatNEvents<<" out, "
1973       // <<fStatGBNGhost/fStatNEvents<<" ghost"<<endl;
1974       fStatTime += hlt.StatTime( 0 );
1975       double timeHz = 0;
1976       if ( fStatTime > 1.e-4 ) timeHz = 1. / fStatTime * fStatNEvents;
1977
1978       out << "<table border>" << endl;
1979       out << "<tr>" << endl;
1980       out << "<td>      </td> <td align=center> RefSet </td> <td align=center> AllSet </td> <td align=center> ExtraSet </td>" << endl;
1981       out << "</tr>" << endl;
1982       out << "<tr>" << endl;
1983       out << "<td>Efficiency</td> <td align=center>" << fStatGBNRecRef / dMCRef
1984       << "</td> <td align=center>" << fStatGBNRecAll / dMCAll
1985       << "</td> <td align=center>" << nRecExtr / dMCExtr
1986       << "</td>" << endl;
1987       out << "</tr>" << endl;
1988       out << "<tr> " << endl;
1989       out << "<td>Clone</td>      <td align=center>" << fStatGBNClonesRef / dRecRef
1990       << "</td> <td align=center>" << fStatGBNClonesAll / dRecAll
1991       << "</td> <td align=center>" << nClonesExtr / dRecExtr
1992       << "</td>" << endl;
1993       out << "</tr>" << endl;
1994       out << "<tr> " << endl;
1995       out << "<td>Ghost</td>      <td colspan=3 align=center>" << fStatGBNGhost / dRecTot
1996       << "</td>" << endl;
1997       out << "</tr>" << endl;
1998       out << "<tr> " << endl;
1999       out << "<td>Time</td>      <td colspan=3 align=center>" << timeHz
2000       << " ev/s</td>" << endl;
2001       out << "</tr>" << endl;
2002       out << "<tr> " << endl;
2003       out << "<td>N Events</td>      <td colspan=3 align=center>" << fStatNEvents
2004       << "</td>" << endl;
2005       out << "</tr>" << endl;
2006       out << "</table>" << endl;
2007     }
2008
2009   }
2010
2011   WriteHistos();
2012 }
2013
2014
2015 void AliHLTTPCCAPerformance::WriteMCEvent( ostream &out ) const
2016 {
2017   // write MC information to the file
2018   out << fNMCTracks << endl;
2019   for ( int it = 0; it < fNMCTracks; it++ ) {
2020     AliHLTTPCCAMCTrack &t = fMCTracks[it];
2021     out << it << " ";
2022     out << t.PDG() << endl;
2023     for ( int i = 0; i < 7; i++ ) out << t.Par()[i] << " ";
2024     out << endl << "    ";
2025     for ( int i = 0; i < 7; i++ ) out << t.TPCPar()[i] << " ";
2026     out << endl << "    ";
2027     out << t.P() << " ";
2028     out << t.Pt() << " ";
2029     out << t.NMCPoints() << " ";
2030     out << t.FirstMCPointID() << " ";
2031     out << t.NHits() << " ";
2032     out << t.NReconstructed() << " ";
2033     out << t.Set() << " ";
2034     out << t.NTurns() << endl;
2035   }
2036
2037   out << fNHits << endl;
2038   for ( int ih = 0; ih < fNHits; ih++ ) {
2039     AliHLTTPCCAHitLabel &l = fHitLabels[ih];
2040     out << l.fLab[0] << " " << l.fLab[1] << " " << l.fLab[2] << endl;
2041   }
2042 }
2043
2044 void AliHLTTPCCAPerformance::WriteMCPoints( ostream &out ) const
2045 {
2046   // write Mc points to the file
2047   out << fNMCPoints << endl;
2048   for ( int ip = 0; ip < fNMCPoints; ip++ ) {
2049     AliHLTTPCCAMCPoint &p = fMCPoints[ip];
2050     out << p.X() << " ";
2051     out << p.Y() << " ";
2052     out << p.Z() << " ";
2053     out << p.Sx() << " ";
2054     out << p.Sy() << " ";
2055     out << p.Sz() << " ";
2056     out << p.Time() << " ";
2057     out << p.ISlice() << " ";
2058     out << p.TrackID() << endl;
2059   }
2060 }
2061
2062 void AliHLTTPCCAPerformance::ReadMCEvent( istream &in )
2063 {
2064   // read mc info from the file
2065   StartEvent();
2066   if ( fMCTracks ) delete[] fMCTracks;
2067   fMCTracks = 0;
2068   fNMCTracks = 0;
2069   if ( fHitLabels ) delete[] fHitLabels;
2070   fHitLabels = 0;
2071   fNHits = 0;
2072   if ( fMCPoints ) delete[] fMCPoints;
2073   fMCPoints = 0;
2074   fNMCPoints = 0;
2075
2076   in >> fNMCTracks;
2077   fMCTracks = new AliHLTTPCCAMCTrack[fNMCTracks];
2078   for ( int it = 0; it < fNMCTracks; it++ ) {
2079     AliHLTTPCCAMCTrack &t = fMCTracks[it];
2080     int j;
2081     float f;
2082     in >> j;
2083     in >> j; t.SetPDG( j );
2084     for ( int i = 0; i < 7; i++ ) { in >> f; t.SetPar( i, f );}
2085     for ( int i = 0; i < 7; i++ ) { in >> f; t.SetTPCPar( i, f );}
2086     in >> f; t.SetP( f );
2087     in >> f; t.SetPt( f );
2088     in >> j; t.SetNHits( j );
2089     in >> j; t.SetNMCPoints( j );
2090     in >> j; t.SetFirstMCPointID( j );
2091     in >> j; t.SetNReconstructed( j );
2092     in >> j; t.SetSet( j );
2093     in >> j; t.SetNTurns( j );
2094   }
2095
2096   in >> fNHits;
2097   fHitLabels = new AliHLTTPCCAHitLabel[fNHits];
2098   for ( int ih = 0; ih < fNHits; ih++ ) {
2099     AliHLTTPCCAHitLabel &l = fHitLabels[ih];
2100     in >> l.fLab[0] >> l.fLab[1] >> l.fLab[2];
2101   }
2102 }
2103
2104 void AliHLTTPCCAPerformance::ReadMCPoints( istream &in )
2105 {
2106   // read mc points from the file
2107   if ( fMCPoints ) delete[] fMCPoints;
2108   fMCPoints = 0;
2109   fNMCPoints = 0;
2110
2111   in >> fNMCPoints;
2112   fMCPoints = new AliHLTTPCCAMCPoint[fNMCPoints];
2113   for ( int ip = 0; ip < fNMCPoints; ip++ ) {
2114     AliHLTTPCCAMCPoint &p = fMCPoints[ip];
2115     float f;
2116     int i;
2117     in >> f;
2118     p.SetX( f );
2119     in >> f;
2120     p.SetY( f );
2121     in >> f;
2122     p.SetZ( f );
2123     in >> f;
2124     p.SetSx( f );
2125     in >> f;
2126     p.SetSy( f );
2127     in >> f;
2128     p.SetSz( f );
2129     in >> f;
2130     p.SetTime( f );
2131     in >> i;
2132     p.SetISlice( i );
2133     in >> i;
2134     p.SetTrackID( i );
2135   }
2136 }