1c8896633a5a8f0a0c376c86e354480a079b4d0f
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx
1 #include "AliESDtrackCuts.h"
2
3
4 #include <AliESDtrack.h>
5 #include <AliESD.h>
6 #include <AliLog.h>
7
8 //____________________________________________________________________
9 ClassImp(AliESDtrackCuts)
10
11 // Cut names
12 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
13  "require TPC refit",
14  "require ITS refit",
15  "n clusters TPC",
16  "n clusters ITS",
17  "#Chi^{2}/clusters TPC",
18  "#Chi^{2}/clusters ITS",
19  "cov 11",
20  "cov 22",
21  "cov 33",
22  "cov 44",
23  "cov 55",
24  "trk-to-vtx",
25  "trk-to-vtx failed",
26  "kink daughters",
27  "p",
28  "p_{T}",
29  "p_{x}",
30  "p_{y}",
31  "p_{z}",
32  "y",
33  "eta"
34 };
35
36 //____________________________________________________________________
37 AliESDtrackCuts::AliESDtrackCuts() : TObject(),
38   fCutMinNClusterTPC(0),
39   fCutMinNClusterITS(0),
40   fCutMaxChi2PerClusterTPC(0),
41   fCutMaxChi2PerClusterITS(0),
42   fCutMaxC11(0),
43   fCutMaxC22(0),
44   fCutMaxC33(0),
45   fCutMaxC44(0),
46   fCutMaxC55(0),
47   fCutAcceptKinkDaughters(0),
48   fCutRequireTPCRefit(0),
49   fCutRequireITSRefit(0),
50   fCutNsigmaToVertex(0),
51   fCutSigmaToVertexRequired(0),
52   fPMin(0),
53   fPMax(0),
54   fPtMin(0),
55   fPtMax(0),
56   fPxMin(0),
57   fPxMax(0),
58   fPyMin(0),
59   fPyMax(0),
60   fPzMin(0),
61   fPzMax(0),
62   fEtaMin(0),
63   fEtaMax(0),
64   fRapMin(0),
65   fRapMax(0),
66   fHistogramsOn(0),
67   fhCutStatistics(0),         
68   fhCutCorrelation(0)
69 {
70   //
71   // constructor
72   //
73
74   Init();
75
76   //##############################################################################
77   // setting default cuts
78   SetMinNClustersTPC();
79   SetMinNClustersITS();
80   SetMaxChi2PerClusterTPC();
81   SetMaxChi2PerClusterITS();                                
82   SetMaxCovDiagonalElements();                                      
83   SetRequireTPCRefit();
84   SetRequireITSRefit();
85   SetAcceptKingDaughters();
86   SetMinNsigmaToVertex();
87   SetRequireSigmaToVertex();
88   SetPRange();
89   SetPtRange();
90   SetPxRange();
91   SetPyRange();
92   SetPzRange();
93   SetEtaRange();
94   SetRapRange();
95
96   SetHistogramsOn();
97 }
98
99 //_____________________________________________________________________________
100 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : TObject(c),
101   fCutMinNClusterTPC(0),
102   fCutMinNClusterITS(0),
103   fCutMaxChi2PerClusterTPC(0),
104   fCutMaxChi2PerClusterITS(0),
105   fCutMaxC11(0),
106   fCutMaxC22(0),
107   fCutMaxC33(0),
108   fCutMaxC44(0),
109   fCutMaxC55(0),
110   fCutAcceptKinkDaughters(0),
111   fCutRequireTPCRefit(0),
112   fCutRequireITSRefit(0),
113   fCutNsigmaToVertex(0),
114   fCutSigmaToVertexRequired(0),
115   fPMin(0),
116   fPMax(0),
117   fPtMin(0),
118   fPtMax(0),
119   fPxMin(0),
120   fPxMax(0),
121   fPyMin(0),
122   fPyMax(0),
123   fPzMin(0),
124   fPzMax(0),
125   fEtaMin(0),
126   fEtaMax(0),
127   fRapMin(0),
128   fRapMax(0),
129   fHistogramsOn(0),
130   fhCutStatistics(0),         
131   fhCutCorrelation(0)
132 {
133   //
134   // copy constructor
135   //
136
137   ((AliESDtrackCuts &) c).Copy(*this);
138 }
139
140 AliESDtrackCuts::~AliESDtrackCuts()
141 {
142   //
143   // destructor
144   //
145
146   // ## TODO to be implemented
147 }
148
149 void AliESDtrackCuts::Init()
150 {
151   //
152   // sets everything to zero
153   //
154
155   fCutMinNClusterTPC = 0;
156   fCutMinNClusterITS = 0;
157
158   fCutMaxChi2PerClusterTPC = 0;
159   fCutMaxChi2PerClusterITS = 0;
160
161   fCutMaxC11 = 0;
162   fCutMaxC22 = 0;
163   fCutMaxC33 = 0;
164   fCutMaxC44 = 0;
165   fCutMaxC55 = 0;
166
167   fCutAcceptKinkDaughters = 0;
168   fCutRequireTPCRefit = 0;
169   fCutRequireITSRefit = 0;
170
171   fCutNsigmaToVertex = 0;
172   fCutSigmaToVertexRequired = 0;
173
174   fPMin = 0;
175   fPMax = 0;
176   fPtMin = 0;
177   fPtMax = 0;
178   fPxMin = 0;
179   fPxMax = 0;
180   fPyMin = 0;
181   fPyMax = 0;
182   fPzMin = 0;
183   fPzMax = 0;
184   fEtaMin = 0;
185   fEtaMax = 0;
186   fRapMin = 0;
187   fRapMax = 0;
188
189   fHistogramsOn = kFALSE;
190
191   for (Int_t i=0; i<2; ++i)
192   {
193     fhNClustersITS[i] = 0;
194     fhNClustersTPC[i] = 0;
195
196     fhChi2PerClusterITS[i] = 0;
197     fhChi2PerClusterTPC[i] = 0;
198
199     fhC11[i] = 0;
200     fhC22[i] = 0;
201     fhC33[i] = 0;
202     fhC44[i] = 0;
203     fhC55[i] = 0;
204
205     fhDXY[i] = 0;
206     fhDZ[i] = 0;
207     fhDXYvsDZ[i] = 0;
208
209     fhDXYNormalized[i] = 0;
210     fhDZNormalized[i] = 0;
211     fhDXYvsDZNormalized[i] = 0;
212   }
213
214   fhCutStatistics = 0;
215   fhCutCorrelation = 0;
216 }
217
218 //_____________________________________________________________________________
219 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
220 {
221   //
222   // Assignment operator
223   //
224
225   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
226   return *this;
227 }
228
229 //_____________________________________________________________________________
230 void AliESDtrackCuts::Copy(TObject &c) const
231 {
232   //
233   // Copy function
234   //
235
236   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
237
238   target.Init();
239
240   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
241   target.fCutMinNClusterITS = fCutMinNClusterITS;
242
243   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
244   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
245
246   target.fCutMaxC11 = fCutMaxC11;
247   target.fCutMaxC22 = fCutMaxC22;
248   target.fCutMaxC33 = fCutMaxC33;
249   target.fCutMaxC44 = fCutMaxC44;
250   target.fCutMaxC55 = fCutMaxC55;
251
252   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
253   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
254   target.fCutRequireITSRefit = fCutRequireITSRefit;
255
256   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
257   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
258
259   target.fPMin = fPMin;
260   target.fPMax = fPMax;
261   target.fPtMin = fPtMin;
262   target.fPtMax = fPtMax;
263   target.fPxMin = fPxMin;
264   target.fPxMax = fPxMax;
265   target.fPyMin = fPyMin;
266   target.fPyMax = fPyMax;
267   target.fPzMin = fPzMin;
268   target.fPzMax = fPzMax;
269   target.fEtaMin = fEtaMin;
270   target.fEtaMax = fEtaMax;
271   target.fRapMin = fRapMin;
272   target.fRapMax = fRapMax;
273
274   target.fHistogramsOn = fHistogramsOn;
275
276   for (Int_t i=0; i<2; ++i)
277   {
278     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
279     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
280
281     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
282     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
283
284     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
285     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
286     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
287     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
288     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
289
290     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
291     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
292     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
293
294     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
295     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
296     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
297   }
298
299   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
300   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
301
302   TObject::Copy(c);
303 }
304
305 //____________________________________________________________________
306 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
307 {
308   //
309
310   Float_t b[2];
311   Float_t bRes[2];
312   Float_t bCov[3];
313   esdTrack->GetImpactParameters(b,bCov);
314   if (bCov[0]<=0 || bCov[2]<=0) {
315     AliDebug(1, "Estimated b resolution lower or equal zero!");
316     bCov[0]=0; bCov[2]=0;
317   }
318   bRes[0] = TMath::Sqrt(bCov[0]);
319   bRes[1] = TMath::Sqrt(bCov[2]);
320
321   // -----------------------------------
322   // How to get to a n-sigma cut?
323   //
324   // The accumulated statistics from 0 to d is
325   //
326   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
327   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
328   //
329   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
330   // Can this be expressed in a different way?
331   //
332   //
333   // FIX: I don't think this is correct!!! Keeping d as n_sigma for now...
334
335   if (bRes[0] == 0 || bRes[1] ==0)
336     return -1;
337
338   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
339
340   // stupid rounding problem screws up everything:
341   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
342   if (TMath::Exp(-d * d / 2) < 1e-10)
343     return 1000;
344
345   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
346   return d;
347 }
348
349 //____________________________________________________________________
350 Bool_t
351 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
352   // 
353   // figure out if the tracks survives all the track cuts defined
354   //
355   // the different quality parameter and kinematic values are first
356   // retrieved from the track. then it is found out what cuts the
357   // track did not survive and finally the cuts are imposed.
358
359
360
361   UInt_t status = esdTrack->GetStatus();
362
363   // dummy array
364   Int_t  fIdxInt[200];
365
366   // getting quality parameters from the ESD track
367   Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
368   Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
369   
370
371
372   Float_t chi2PerClusterITS = -1;
373   Float_t chi2PerClusterTPC = -1;
374   if (nClustersITS!=0)
375     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
376   if (nClustersTPC!=0)
377     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);  
378
379   Double_t extCov[15];
380   esdTrack->GetExternalCovariance(extCov);
381
382   // getting the track to vertex parameters
383   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
384
385   // getting the kinematic variables of the track
386   // (assuming the mass is known)
387   Double_t p[3];
388   esdTrack->GetPxPyPz(p);
389   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
390   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
391   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
392
393
394   //y-eta related calculations
395   Float_t eta = -100.;
396   Float_t y   = -100.;
397   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
398     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
399   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
400     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
401
402   
403   //########################################################################
404   // cut the track?
405   
406   Bool_t cuts[kNCuts];
407   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
408   
409   // track quality cuts
410   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
411     cuts[0]=kTRUE;
412   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
413     cuts[1]=kTRUE;
414   if (nClustersTPC<fCutMinNClusterTPC)
415     cuts[2]=kTRUE;
416   if (nClustersITS<fCutMinNClusterITS) 
417     cuts[3]=kTRUE;
418   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
419     cuts[4]=kTRUE; 
420   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
421     cuts[5]=kTRUE;
422   if (extCov[0]  > fCutMaxC11) 
423     cuts[6]=kTRUE;  
424   if (extCov[2]  > fCutMaxC22) 
425     cuts[7]=kTRUE;  
426   if (extCov[5]  > fCutMaxC33) 
427     cuts[8]=kTRUE;  
428   if (extCov[9]  > fCutMaxC44) 
429     cuts[9]=kTRUE;  
430   if (extCov[14]  > fCutMaxC55) 
431     cuts[10]=kTRUE;  
432   if (nSigmaToVertex > fCutNsigmaToVertex)
433     cuts[11] = kTRUE;
434   // if n sigma could not be calculated
435   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
436     cuts[12]=kTRUE;
437   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) 
438     cuts[13]=kTRUE;
439   // track kinematics cut
440   if((momentum < fPMin) || (momentum > fPMax)) 
441     cuts[14]=kTRUE;
442   if((pt < fPtMin) || (pt > fPtMax)) 
443     cuts[15] = kTRUE;
444   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
445     cuts[16] = kTRUE;
446   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
447     cuts[17] = kTRUE;
448   if((p[2] < fPzMin) || (p[2] > fPzMax))
449     cuts[18] = kTRUE;
450   if((eta < fEtaMin) || (eta > fEtaMax)) 
451     cuts[19] = kTRUE;
452   if((y < fRapMin) || (y > fRapMax)) 
453     cuts[20] = kTRUE;
454
455   Bool_t cut=kFALSE;
456   for (Int_t i=0; i<kNCuts; i++) 
457     if (cuts[i]) cut = kTRUE;
458   
459   //########################################################################
460   // filling histograms
461   if (fHistogramsOn) {
462     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
463     
464     if (cut)
465       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
466     
467     for (Int_t i=0; i<kNCuts; i++) {
468       if (cuts[i])
469         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
470       
471       for (Int_t j=i; j<kNCuts; j++) {
472         if (cuts[i] && cuts[j]) {
473           Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
474           Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
475           fhCutCorrelation->Fill(x,y);
476         }
477       }
478     }
479     
480
481     fhNClustersITS[0]->Fill(nClustersITS);
482     fhNClustersTPC[0]->Fill(nClustersTPC);
483     fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
484     fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
485
486     fhC11[0]->Fill(extCov[0]);
487     fhC22[0]->Fill(extCov[2]);
488     fhC33[0]->Fill(extCov[5]);
489     fhC44[0]->Fill(extCov[9]);
490     fhC55[0]->Fill(extCov[14]);
491
492     Float_t b[2];
493     Float_t bRes[2];
494     Float_t bCov[3];
495     esdTrack->GetImpactParameters(b,bCov);
496     if (bCov[0]<=0 || bCov[2]<=0) {
497       AliDebug(1, "Estimated b resolution lower or equal zero!");
498       bCov[0]=0; bCov[2]=0;
499     }
500     bRes[0] = TMath::Sqrt(bCov[0]);
501     bRes[1] = TMath::Sqrt(bCov[2]);
502
503     fhDZ[0]->Fill(b[1]);
504     fhDXY[0]->Fill(b[0]);
505     fhDXYvsDZ[0]->Fill(b[1],b[0]);
506
507     if (bRes[0]!=0 && bRes[1]!=0) {
508       fhDZNormalized[0]->Fill(b[1]/bRes[1]);
509       fhDXYNormalized[0]->Fill(b[0]/bRes[0]);
510       fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
511     }
512   }
513
514   //########################################################################
515   // cut the track!
516   if (cut) return kFALSE;
517
518   //########################################################################
519   // filling histograms after cut
520   if (fHistogramsOn) {
521     fhNClustersITS[1]->Fill(nClustersITS);
522     fhNClustersTPC[1]->Fill(nClustersTPC);
523     fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
524     fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
525
526     fhC11[1]->Fill(extCov[0]);
527     fhC22[1]->Fill(extCov[2]);
528     fhC33[1]->Fill(extCov[5]);
529     fhC44[1]->Fill(extCov[9]);
530     fhC55[1]->Fill(extCov[14]);
531
532     Float_t b[2];
533     Float_t bRes[2];
534     Float_t bCov[3];
535     esdTrack->GetImpactParameters(b,bCov);
536     if (bCov[0]<=0 || bCov[2]<=0) {
537       AliDebug(1, "Estimated b resolution lower or equal zero!");
538       bCov[0]=0; bCov[2]=0;
539     }
540     bRes[0] = TMath::Sqrt(bCov[0]);
541     bRes[1] = TMath::Sqrt(bCov[2]);
542
543     fhDZ[1]->Fill(b[1]);
544     fhDXY[1]->Fill(b[0]);
545     fhDXYvsDZ[1]->Fill(b[1],b[0]);
546
547     if (bRes[0]!=0 && bRes[1]!=0)
548     {
549       fhDZNormalized[1]->Fill(b[1]/bRes[1]);
550       fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
551       fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
552     }
553   }
554
555   return kTRUE;
556 }
557
558 //____________________________________________________________________
559 TObjArray*
560 AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
561 {
562   //
563   // returns an array of all tracks that pass the cuts
564   //
565
566   TObjArray* acceptedTracks = new TObjArray();
567
568   // loop over esd tracks
569   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
570     AliESDtrack* track = esd->GetTrack(iTrack);
571
572     if (AcceptTrack(track))
573       acceptedTracks->Add(track);
574   }
575
576   return acceptedTracks;
577 }
578
579 //____________________________________________________________________
580 Int_t
581 AliESDtrackCuts::CountAcceptedTracks(AliESD* esd)
582 {
583   //
584   // returns an the number of tracks that pass the cuts
585   //
586
587   Int_t count = 0;
588
589   // loop over esd tracks
590   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
591     AliESDtrack* track = esd->GetTrack(iTrack);
592
593     if (AcceptTrack(track))
594       count++;
595   }
596
597   return count;
598 }
599
600 //____________________________________________________________________
601  void AliESDtrackCuts::DefineHistograms(Int_t color) {
602    // 
603    // diagnostics histograms are defined
604    // 
605
606    fHistogramsOn=kTRUE;
607
608    //###################################################################################
609    // defining histograms
610
611    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
612
613    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
614    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
615
616    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
617   
618    for (Int_t i=0; i<kNCuts; i++) {
619      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
620      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
621      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
622    } 
623
624   fhCutStatistics  ->SetLineColor(color);
625   fhCutCorrelation ->SetLineColor(color);
626   fhCutStatistics  ->SetLineWidth(2);
627   fhCutCorrelation ->SetLineWidth(2);
628
629   Char_t str[256];
630   for (Int_t i=0; i<2; i++) {
631     if (i==0) sprintf(str," ");
632     else sprintf(str,"_cut");
633
634     fhNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
635     fhNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
636     fhChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
637     fhChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
638
639     fhC11[i]                 = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
640     fhC22[i]                 = new  TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
641     fhC33[i]                 = new  TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
642     fhC44[i]                 = new  TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
643     fhC55[i]                 = new  TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
644
645     fhDXY[i]                 = new  TH1F(Form("dXY%s",str),"",500,-10,10);
646     fhDZ[i]                  = new  TH1F(Form("dZ%s",str),"",500,-10,10);
647     fhDXYvsDZ[i]             = new  TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
648
649     fhDXYNormalized[i]       = new  TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
650     fhDZNormalized[i]        = new  TH1F(Form("dZNormalized%s",str),"",500,-10,10);
651     fhDXYvsDZNormalized[i]   = new  TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
652
653
654     fhNClustersITS[i]->SetXTitle("n ITS clusters");
655     fhNClustersTPC[i]->SetXTitle("n TPC clusters");
656     fhChi2PerClusterITS[i]->SetXTitle("#Chi^{2} per ITS cluster");
657     fhChi2PerClusterTPC[i]->SetXTitle("#Chi^{2} per TPC cluster");
658
659     fhC11[i]->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
660     fhC22[i]->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
661     fhC33[i]->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
662     fhC44[i]->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
663     fhC55[i]->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
664
665     fhDXY[i]->SetXTitle("transverse impact parameter");
666     fhDZ[i]->SetXTitle("longitudinal impact parameter");
667     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter");
668     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter");
669
670     fhDXYNormalized[i]->SetXTitle("normalized trans impact par");
671     fhDZNormalized[i]->SetXTitle("normalized long impact par");
672     fhDXYvsDZNormalized[i]->SetXTitle("normalized long impact par"); 
673     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par");
674
675     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
676     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
677     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
678     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
679
680     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
681     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
682     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
683     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
684     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
685
686     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
687     fhDZ[i]->SetLineColor(color);   fhDZ[i]->SetLineWidth(2);
688
689     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
690     fhDZNormalized[i]->SetLineColor(color);   fhDZNormalized[i]->SetLineWidth(2);
691   }
692 }
693
694 //____________________________________________________________________
695 void 
696 AliESDtrackCuts::Print(const Option_t*) const {
697   //
698   // print method - still to be implemented
699   //
700
701   AliInfo("AliESDtrackCuts...");
702 }
703
704
705 //____________________________________________________________________
706 void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
707   // 
708   // saves the histograms in a directory (dir)
709   // 
710
711   
712   if (!fHistogramsOn) {
713     AliDebug(0, "Histograms not on - cannot save histograms!!!");
714     return;
715   }
716
717   gDirectory->mkdir(dir);
718   gDirectory->cd(dir);
719
720   gDirectory->mkdir("before_cuts");
721   gDirectory->mkdir("after_cuts");
722  
723   fhCutStatistics->Write();
724   fhCutCorrelation->Write();
725
726   for (Int_t i=0; i<2; i++) {
727     if (i==0)
728       gDirectory->cd("before_cuts");
729     else
730       gDirectory->cd("after_cuts");
731     
732     fhNClustersITS[i]        ->Write();
733     fhNClustersTPC[i]        ->Write();
734     fhChi2PerClusterITS[i]   ->Write();
735     fhChi2PerClusterTPC[i]   ->Write();
736     
737     fhC11[i]                 ->Write();
738     fhC22[i]                 ->Write();
739     fhC33[i]                 ->Write();
740     fhC44[i]                 ->Write();
741     fhC55[i]                 ->Write();
742
743     fhDXY[i]                 ->Write();
744     fhDZ[i]                  ->Write();
745     fhDXYvsDZ[i]             ->Write();
746     
747     fhDXYNormalized[i]       ->Write();
748     fhDZNormalized[i]        ->Write();
749     fhDXYvsDZNormalized[i]   ->Write();
750
751     gDirectory->cd("../");
752   }
753
754   gDirectory->cd("../");
755 }
756
757
758