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