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