output of the slice tracker optimised
[u/mrichter/AliRoot.git] / 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       int id, row;
1236       float x,y,z;
1237       tCA->Cluster(i).Get(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 }