Added AddTrackParams() method for convenience + some comments
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDtrackCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliESDtrackCuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
17
18 #include "AliESDtrackCuts.h"
19
20 #include <AliESDtrack.h>
21 #include <AliESDVertex.h>
22 #include <AliESDEvent.h>
23 #include <AliLog.h>
24
25 #include <TTree.h>
26 #include <TCanvas.h>
27 #include <TDirectory.h>
28 #include <TH2F.h>
29 #include <TF1.h>
30
31 //____________________________________________________________________
32 ClassImp(AliESDtrackCuts)
33
34 // Cut names
35 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
36  "require TPC refit",
37  "require ITS refit",
38  "n clusters TPC",
39  "n clusters ITS",
40  "#Chi^{2}/cluster TPC",
41  "#Chi^{2}/cluster ITS",
42  "cov 11",
43  "cov 22",
44  "cov 33",
45  "cov 44",
46  "cov 55",
47  "trk-to-vtx",
48  "trk-to-vtx failed",
49  "kink daughters",
50  "p",
51  "p_{T}",
52  "p_{x}",
53  "p_{y}",
54  "p_{z}",
55  "eta",
56  "y",
57  "trk-to-vtx max dca 2D absolute",
58  "trk-to-vtx max dca xy absolute",
59  "trk-to-vtx max dca z absolute",
60  "trk-to-vtx min dca 2D absolute",
61  "trk-to-vtx min dca xy absolute",
62  "trk-to-vtx min dca z absolute",
63  "SPD cluster requirement",
64  "SDD cluster requirement",
65  "SSD cluster requirement",
66  "require ITS stand-alone"
67 };
68
69 //____________________________________________________________________
70 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
71   fCutMinNClusterTPC(0),
72   fCutMinNClusterITS(0),
73   fCutMaxChi2PerClusterTPC(0),
74   fCutMaxChi2PerClusterITS(0),
75   fCutMaxC11(0),
76   fCutMaxC22(0),
77   fCutMaxC33(0),
78   fCutMaxC44(0),
79   fCutMaxC55(0),
80   fCutAcceptKinkDaughters(0),
81   fCutRequireTPCRefit(0),
82   fCutRequireITSRefit(0),
83   fCutRequireITSStandAlone(0),
84   fCutNsigmaToVertex(0),
85   fCutSigmaToVertexRequired(0),
86   fCutMaxDCAToVertexXY(0),
87   fCutMaxDCAToVertexZ(0),
88   fCutMinDCAToVertexXY(0),
89   fCutMinDCAToVertexZ(0),
90   fCutDCAToVertex2D(0),
91   fPMin(0),
92   fPMax(0),
93   fPtMin(0),
94   fPtMax(0),
95   fPxMin(0),
96   fPxMax(0),
97   fPyMin(0),
98   fPyMax(0),
99   fPzMin(0),
100   fPzMax(0),
101   fEtaMin(0),
102   fEtaMax(0),
103   fRapMin(0),
104   fRapMax(0),
105   fHistogramsOn(0),
106   ffDTheoretical(0),
107   fhCutStatistics(0),         
108   fhCutCorrelation(0)
109 {
110   //
111   // constructor
112   //
113
114   Init();
115
116   //##############################################################################
117   // setting default cuts
118   SetMinNClustersTPC();
119   SetMinNClustersITS();
120   SetMaxChi2PerClusterTPC();
121   SetMaxChi2PerClusterITS();                                
122   SetMaxCovDiagonalElements();                                      
123   SetRequireTPCRefit();
124   SetRequireITSRefit();
125   SetRequireITSStandAlone(kFALSE);
126   SetAcceptKinkDaughters();
127   SetMaxNsigmaToVertex();
128   SetMaxDCAToVertexXY();
129   SetMaxDCAToVertexZ();
130   SetDCAToVertex2D();
131   SetMinDCAToVertexXY();
132   SetMinDCAToVertexZ();
133   SetPRange();
134   SetPtRange();
135   SetPxRange();
136   SetPyRange();
137   SetPzRange();
138   SetEtaRange();
139   SetRapRange();
140   SetClusterRequirementITS(kSPD);
141   SetClusterRequirementITS(kSDD);
142   SetClusterRequirementITS(kSSD);
143
144   SetHistogramsOn();
145 }
146
147 //_____________________________________________________________________________
148 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
149   fCutMinNClusterTPC(0),
150   fCutMinNClusterITS(0),
151   fCutMaxChi2PerClusterTPC(0),
152   fCutMaxChi2PerClusterITS(0),
153   fCutMaxC11(0),
154   fCutMaxC22(0),
155   fCutMaxC33(0),
156   fCutMaxC44(0),
157   fCutMaxC55(0),
158   fCutAcceptKinkDaughters(0),
159   fCutRequireTPCRefit(0),
160   fCutRequireITSRefit(0),
161   fCutRequireITSStandAlone(0),
162   fCutNsigmaToVertex(0),
163   fCutSigmaToVertexRequired(0),
164   fCutMaxDCAToVertexXY(0),
165   fCutMaxDCAToVertexZ(0),
166   fCutMinDCAToVertexXY(0),
167   fCutMinDCAToVertexZ(0),
168   fCutDCAToVertex2D(0),
169   fPMin(0),
170   fPMax(0),
171   fPtMin(0),
172   fPtMax(0),
173   fPxMin(0),
174   fPxMax(0),
175   fPyMin(0),
176   fPyMax(0),
177   fPzMin(0),
178   fPzMax(0),
179   fEtaMin(0),
180   fEtaMax(0),
181   fRapMin(0),
182   fRapMax(0),
183   fHistogramsOn(0),
184   ffDTheoretical(0),                                 
185   fhCutStatistics(0),         
186   fhCutCorrelation(0)
187 {
188   //
189   // copy constructor
190   //
191
192   ((AliESDtrackCuts &) c).Copy(*this);
193 }
194
195 AliESDtrackCuts::~AliESDtrackCuts()
196 {
197   //
198   // destructor
199   //
200
201   for (Int_t i=0; i<2; i++) {
202     
203     if (fhNClustersITS[i])
204       delete fhNClustersITS[i];            
205     if (fhNClustersTPC[i])
206       delete fhNClustersTPC[i];                
207     if (fhChi2PerClusterITS[i])
208       delete fhChi2PerClusterITS[i];       
209     if (fhChi2PerClusterTPC[i])
210       delete fhChi2PerClusterTPC[i];       
211     if (fhC11[i])
212       delete fhC11[i];                     
213     if (fhC22[i])
214       delete fhC22[i];                     
215     if (fhC33[i])
216       delete fhC33[i];                     
217     if (fhC44[i])
218       delete fhC44[i];                     
219     if (fhC55[i])
220     delete fhC55[i];                     
221     
222     if (fhDXY[i])
223       delete fhDXY[i];                     
224     if (fhDZ[i])
225       delete fhDZ[i];
226     if (fhDXYDZ[i])
227       delete fhDXYDZ[i];
228     if (fhDXYvsDZ[i])
229       delete fhDXYvsDZ[i];
230
231     if (fhDXYNormalized[i])
232       delete fhDXYNormalized[i];           
233     if (fhDZNormalized[i])
234       delete fhDZNormalized[i];
235     if (fhDXYvsDZNormalized[i])
236       delete fhDXYvsDZNormalized[i];
237     if (fhNSigmaToVertex[i])
238       delete fhNSigmaToVertex[i];
239     if (fhPt[i])
240       delete fhPt[i];
241     if (fhEta[i])
242       delete fhEta[i];
243   }
244
245   if (ffDTheoretical)
246     delete ffDTheoretical;
247
248   if (fhCutStatistics)
249     delete fhCutStatistics;             
250   if (fhCutCorrelation)
251     delete fhCutCorrelation;            
252 }
253
254 void AliESDtrackCuts::Init()
255 {
256   //
257   // sets everything to zero
258   //
259
260   fCutMinNClusterTPC = 0;
261   fCutMinNClusterITS = 0;
262
263   fCutMaxChi2PerClusterTPC = 0;
264   fCutMaxChi2PerClusterITS = 0;
265   
266   for (Int_t i = 0; i < 3; i++)
267         fCutClusterRequirementITS[i] = kOff;
268
269   fCutMaxC11 = 0;
270   fCutMaxC22 = 0;
271   fCutMaxC33 = 0;
272   fCutMaxC44 = 0;
273   fCutMaxC55 = 0;
274
275   fCutAcceptKinkDaughters = 0;
276   fCutRequireTPCRefit = 0;
277   fCutRequireITSRefit = 0;
278   fCutRequireITSStandAlone = 0;
279
280   fCutNsigmaToVertex = 0;
281   fCutSigmaToVertexRequired = 0;
282   fCutMaxDCAToVertexXY = 0;
283   fCutMaxDCAToVertexZ = 0;
284   fCutDCAToVertex2D = 0;
285   fCutMinDCAToVertexXY = 0;
286   fCutMinDCAToVertexZ = 0;
287
288   
289   fPMin = 0;
290   fPMax = 0;
291   fPtMin = 0;
292   fPtMax = 0;
293   fPxMin = 0;
294   fPxMax = 0;
295   fPyMin = 0;
296   fPyMax = 0;
297   fPzMin = 0;
298   fPzMax = 0;
299   fEtaMin = 0;
300   fEtaMax = 0;
301   fRapMin = 0;
302   fRapMax = 0;
303
304   fHistogramsOn = kFALSE;
305
306   for (Int_t i=0; i<2; ++i)
307   {
308     fhNClustersITS[i] = 0;
309     fhNClustersTPC[i] = 0;
310
311     fhChi2PerClusterITS[i] = 0;
312     fhChi2PerClusterTPC[i] = 0;
313
314     fhC11[i] = 0;
315     fhC22[i] = 0;
316     fhC33[i] = 0;
317     fhC44[i] = 0;
318     fhC55[i] = 0;
319
320     fhDXY[i] = 0;
321     fhDZ[i] = 0;
322     fhDXYDZ[i] = 0;
323     fhDXYvsDZ[i] = 0;
324
325     fhDXYNormalized[i] = 0;
326     fhDZNormalized[i] = 0;
327     fhDXYvsDZNormalized[i] = 0;
328     fhNSigmaToVertex[i] = 0;
329     
330     fhPt[i] = 0;
331     fhEta[i] = 0;
332   }
333   ffDTheoretical = 0;
334
335   fhCutStatistics = 0;
336   fhCutCorrelation = 0;
337 }
338
339 //_____________________________________________________________________________
340 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
341 {
342   //
343   // Assignment operator
344   //
345
346   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
347   return *this;
348 }
349
350 //_____________________________________________________________________________
351 void AliESDtrackCuts::Copy(TObject &c) const
352 {
353   //
354   // Copy function
355   //
356
357   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
358
359   target.Init();
360
361   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
362   target.fCutMinNClusterITS = fCutMinNClusterITS;
363
364   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
365   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
366
367   for (Int_t i = 0; i < 3; i++)
368     target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
369
370   target.fCutMaxC11 = fCutMaxC11;
371   target.fCutMaxC22 = fCutMaxC22;
372   target.fCutMaxC33 = fCutMaxC33;
373   target.fCutMaxC44 = fCutMaxC44;
374   target.fCutMaxC55 = fCutMaxC55;
375
376   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
377   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
378   target.fCutRequireITSRefit = fCutRequireITSRefit;
379   target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
380
381   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
382   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
383   target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
384   target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
385   target.fCutDCAToVertex2D = fCutDCAToVertex2D;
386   target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
387   target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
388
389   target.fPMin = fPMin;
390   target.fPMax = fPMax;
391   target.fPtMin = fPtMin;
392   target.fPtMax = fPtMax;
393   target.fPxMin = fPxMin;
394   target.fPxMax = fPxMax;
395   target.fPyMin = fPyMin;
396   target.fPyMax = fPyMax;
397   target.fPzMin = fPzMin;
398   target.fPzMax = fPzMax;
399   target.fEtaMin = fEtaMin;
400   target.fEtaMax = fEtaMax;
401   target.fRapMin = fRapMin;
402   target.fRapMax = fRapMax;
403
404   target.fHistogramsOn = fHistogramsOn;
405
406   for (Int_t i=0; i<2; ++i)
407   {
408     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
409     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
410
411     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
412     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
413
414     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
415     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
416     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
417     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
418     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
419
420     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
421     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
422     if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
423     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
424
425     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
426     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
427     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
428     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
429     
430     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
431     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
432   }
433   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
434
435   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
436   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
437
438   TNamed::Copy(c);
439 }
440
441 //_____________________________________________________________________________
442 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
443   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
444   // Returns the number of merged objects (including this)
445   if (!list)
446     return 0;
447   if (list->IsEmpty())
448     return 1;
449   if (!fHistogramsOn)
450     return 0;
451   TIterator* iter = list->MakeIterator();
452   TObject* obj;
453
454   // collection of measured and generated histograms
455   Int_t count = 0;
456   while ((obj = iter->Next())) {
457
458     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
459     if (entry == 0)
460       continue;
461
462     if (!entry->fHistogramsOn)
463       continue;
464
465     for (Int_t i=0; i<2; i++) {
466       
467       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );      
468       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     ); 
469                                                                     
470       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]); 
471       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]); 
472                                                                     
473       fhC11[i]               ->Add(entry->fhC11[i]              ); 
474       fhC22[i]               ->Add(entry->fhC22[i]              ); 
475       fhC33[i]               ->Add(entry->fhC33[i]              ); 
476       fhC44[i]               ->Add(entry->fhC44[i]              ); 
477       fhC55[i]               ->Add(entry->fhC55[i]              ); 
478                                                                     
479       fhDXY[i]               ->Add(entry->fhDXY[i]              ); 
480       fhDZ[i]                ->Add(entry->fhDZ[i]               ); 
481       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
482       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
483
484       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
485       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
486       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
487       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]); 
488
489       fhPt[i]                ->Add(entry->fhPt[i]); 
490       fhEta[i]               ->Add(entry->fhEta[i]); 
491     }      
492
493     fhCutStatistics  ->Add(entry->fhCutStatistics);        
494     fhCutCorrelation ->Add(entry->fhCutCorrelation);      
495
496     count++;
497   }
498   return count+1;
499 }
500
501
502 //____________________________________________________________________
503 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
504 {
505   // Calculates the number of sigma to the vertex.
506
507   Float_t b[2];
508   Float_t bRes[2];
509   Float_t bCov[3];
510   esdTrack->GetImpactParameters(b,bCov);
511   
512   if (bCov[0]<=0 || bCov[2]<=0) {
513     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
514     bCov[0]=0; bCov[2]=0;
515   }
516   bRes[0] = TMath::Sqrt(bCov[0]);
517   bRes[1] = TMath::Sqrt(bCov[2]);
518
519   // -----------------------------------
520   // How to get to a n-sigma cut?
521   //
522   // The accumulated statistics from 0 to d is
523   //
524   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
525   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
526   //
527   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
528   // Can this be expressed in a different way?
529
530   if (bRes[0] == 0 || bRes[1] ==0)
531     return -1;
532
533   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
534
535   // work around precision problem
536   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
537   // 1e-15 corresponds to nsigma ~ 7.7
538   if (TMath::Exp(-d * d / 2) < 1e-15)
539     return 1000;
540
541   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
542   return nSigma;
543 }
544
545 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
546 {
547   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
548
549   tree->SetBranchStatus("fTracks.fFlags", 1);
550   tree->SetBranchStatus("fTracks.fITSncls", 1);
551   tree->SetBranchStatus("fTracks.fTPCncls", 1);
552   tree->SetBranchStatus("fTracks.fITSchi2", 1);
553   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
554   tree->SetBranchStatus("fTracks.fC*", 1);
555   tree->SetBranchStatus("fTracks.fD", 1);
556   tree->SetBranchStatus("fTracks.fZ", 1);
557   tree->SetBranchStatus("fTracks.fCdd", 1);
558   tree->SetBranchStatus("fTracks.fCdz", 1);
559   tree->SetBranchStatus("fTracks.fCzz", 1);
560   tree->SetBranchStatus("fTracks.fP*", 1);
561   tree->SetBranchStatus("fTracks.fR*", 1);
562   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
563 }
564
565 //____________________________________________________________________
566 Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) 
567 {
568   // 
569   // figure out if the tracks survives all the track cuts defined
570   //
571   // the different quality parameter and kinematic values are first
572   // retrieved from the track. then it is found out what cuts the
573   // track did not survive and finally the cuts are imposed.
574
575   // this function needs the following branches:
576   // fTracks.fFlags
577   // fTracks.fITSncls
578   // fTracks.fTPCncls
579   // fTracks.fITSchi2
580   // fTracks.fTPCchi2
581   // fTracks.fC   //GetExternalCovariance
582   // fTracks.fD   //GetImpactParameters
583   // fTracks.fZ   //GetImpactParameters
584   // fTracks.fCdd //GetImpactParameters
585   // fTracks.fCdz //GetImpactParameters
586   // fTracks.fCzz //GetImpactParameters
587   // fTracks.fP   //GetPxPyPz
588   // fTracks.fR   //GetMass
589   // fTracks.fP   //GetMass
590   // fTracks.fKinkIndexes
591
592   UInt_t status = esdTrack->GetStatus();
593
594   // getting quality parameters from the ESD track
595   Int_t nClustersITS = esdTrack->GetITSclusters(0);
596   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
597
598   Float_t chi2PerClusterITS = -1;
599   Float_t chi2PerClusterTPC = -1;
600   if (nClustersITS!=0)
601     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
602   if (nClustersTPC!=0)
603     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
604   Double_t extCov[15];
605   esdTrack->GetExternalCovariance(extCov);
606
607   // getting the track to vertex parameters
608   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
609       
610   Float_t b[2];
611   Float_t bCov[3];
612   esdTrack->GetImpactParameters(b,bCov);
613   if (bCov[0]<=0 || bCov[2]<=0) {
614     AliDebug(1, "Estimated b resolution lower or equal zero!");
615     bCov[0]=0; bCov[2]=0;
616   }
617
618   Float_t dcaToVertexXY = b[0];
619   Float_t dcaToVertexZ = b[1];
620
621   Float_t dcaToVertex = -1;
622   
623   if (fCutDCAToVertex2D)
624   {
625     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
626   }
627   else
628     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
629     
630   // getting the kinematic variables of the track
631   // (assuming the mass is known)
632   Double_t p[3];
633   esdTrack->GetPxPyPz(p);
634
635   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
636   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
637   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
638
639
640   //y-eta related calculations
641   Float_t eta = -100.;
642   Float_t y   = -100.;
643   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
644     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
645   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
646     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
647
648   
649   //########################################################################
650   // cut the track?
651   
652   Bool_t cuts[kNCuts];
653   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
654   
655   // track quality cuts
656   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
657     cuts[0]=kTRUE;
658   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
659     cuts[1]=kTRUE;
660   if (nClustersTPC<fCutMinNClusterTPC)
661     cuts[2]=kTRUE;
662   if (nClustersITS<fCutMinNClusterITS) 
663     cuts[3]=kTRUE;
664   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
665     cuts[4]=kTRUE; 
666   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
667     cuts[5]=kTRUE;
668   if (extCov[0]  > fCutMaxC11) 
669     cuts[6]=kTRUE;  
670   if (extCov[2]  > fCutMaxC22) 
671     cuts[7]=kTRUE;  
672   if (extCov[5]  > fCutMaxC33) 
673     cuts[8]=kTRUE;  
674   if (extCov[9]  > fCutMaxC44) 
675     cuts[9]=kTRUE;  
676   if (extCov[14]  > fCutMaxC55) 
677     cuts[10]=kTRUE;  
678   if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
679     cuts[11] = kTRUE;
680   // if n sigma could not be calculated
681   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
682     cuts[12]=kTRUE;
683   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
684     cuts[13]=kTRUE;
685   // track kinematics cut
686   if((momentum < fPMin) || (momentum > fPMax)) 
687     cuts[14]=kTRUE;
688   if((pt < fPtMin) || (pt > fPtMax)) 
689     cuts[15] = kTRUE;
690   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
691     cuts[16] = kTRUE;
692   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
693     cuts[17] = kTRUE;
694   if((p[2] < fPzMin) || (p[2] > fPzMax))
695     cuts[18] = kTRUE;
696   if((eta < fEtaMin) || (eta > fEtaMax))
697     cuts[19] = kTRUE;
698   if((y < fRapMin) || (y > fRapMax)) 
699     cuts[20] = kTRUE;
700   if (fCutDCAToVertex2D && dcaToVertex > 1)
701     cuts[21] = kTRUE;
702   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
703     cuts[22] = kTRUE;
704   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
705     cuts[23] = kTRUE;
706   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
707     cuts[24] = kTRUE;
708   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
709     cuts[25] = kTRUE;
710   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
711     cuts[26] = kTRUE;
712   
713   for (Int_t i = 0; i < 3; i++)
714     cuts[27+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
715   
716   if (fCutRequireITSStandAlone && ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)))
717     cuts[30]=kTRUE;
718   
719   Bool_t cut=kFALSE;
720   for (Int_t i=0; i<kNCuts; i++) 
721     if (cuts[i]) {cut = kTRUE;}
722
723
724
725   //########################################################################
726   // filling histograms
727   if (fHistogramsOn) {
728     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
729     if (cut)
730       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
731     
732     for (Int_t i=0; i<kNCuts; i++) {
733       if (cuts[i])
734         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
735
736       for (Int_t j=i; j<kNCuts; j++) {
737         if (cuts[i] && cuts[j]) {
738           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
739           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
740           fhCutCorrelation->Fill(xC, yC);
741         }
742       }
743     }
744   }
745
746   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
747   // the code is not in a function due to too many local variables that would need to be passed
748
749   for (Int_t id = 0; id < 2; id++)
750   {
751     // id = 0 --> before cut
752     // id = 1 --> after cut
753
754     if (fHistogramsOn)
755     {
756       fhNClustersITS[id]->Fill(nClustersITS);
757       fhNClustersTPC[id]->Fill(nClustersTPC);
758       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
759       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
760
761       fhC11[id]->Fill(extCov[0]);
762       fhC22[id]->Fill(extCov[2]);
763       fhC33[id]->Fill(extCov[5]);
764       fhC44[id]->Fill(extCov[9]);
765       fhC55[id]->Fill(extCov[14]);
766
767       fhPt[id]->Fill(pt);
768       fhEta[id]->Fill(eta);
769
770       Float_t bRes[2];
771       bRes[0] = TMath::Sqrt(bCov[0]);
772       bRes[1] = TMath::Sqrt(bCov[2]);
773
774       fhDZ[id]->Fill(b[1]);
775       fhDXY[id]->Fill(b[0]);
776       fhDXYDZ[id]->Fill(dcaToVertex);
777       fhDXYvsDZ[id]->Fill(b[1],b[0]);
778
779       if (bRes[0]!=0 && bRes[1]!=0) {
780         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
781         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
782         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
783         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
784       }
785     }
786
787     // cut the track
788     if (cut)
789       return kFALSE;
790   }
791
792   return kTRUE;
793 }
794
795 //____________________________________________________________________
796 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
797 {
798   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
799   
800   switch (req)
801   {
802         case kOff:        return kTRUE;
803         case kNone:       return !clusterL1 && !clusterL2;
804         case kAny:        return clusterL1 || clusterL2;
805         case kFirst:      return clusterL1;
806         case kOnlyFirst:  return clusterL1 && !clusterL2;
807         case kSecond:     return clusterL2;
808         case kOnlySecond: return clusterL2 && !clusterL1;
809         case kBoth:       return clusterL1 && clusterL2;
810   }
811   
812   return kFALSE;
813 }
814
815 //____________________________________________________________________
816 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
817 {
818   // creates a TPC only track from the given esd track
819   // the track has to be deleted by the user
820   //
821   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
822   // there are only missing propagations here that are needed for old data
823   // this function will therefore become obsolete
824   //
825   // adapted from code provided by CKB
826
827   if (!esd->GetPrimaryVertexTPC())
828     return 0; // No TPC vertex no TPC tracks
829
830   if(!esd->GetPrimaryVertexTPC()->GetStatus())
831     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
832
833   AliESDtrack* track = esd->GetTrack(iTrack);
834   if (!track)
835     return 0;
836
837   AliESDtrack *tpcTrack = new AliESDtrack();
838
839   // This should have been done during the reconstruction
840   // fixed by Juri in r26675
841   // but recalculate for older data CKB
842   Float_t p[2],cov[3];
843   track->GetImpactParametersTPC(p,cov);
844   if(p[0]==0&&p[1]==0)
845     track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
846   // BKC
847
848   // only true if we have a tpc track
849   if (!track->FillTPCOnlyTrack(*tpcTrack))
850   {
851     delete tpcTrack;
852     return 0;
853   }
854
855   // propagate to Vertex
856   // not needed for normal reconstructed ESDs...
857   // Double_t pTPC[2],covTPC[3];
858   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
859
860   return tpcTrack;
861 }
862
863 //____________________________________________________________________
864 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
865 {
866   //
867   // returns an array of all tracks that pass the cuts
868   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
869   // tracks that pass the cut
870
871   TObjArray* acceptedTracks = new TObjArray();
872
873   // loop over esd tracks
874   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
875     if(bTPC){
876       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
877       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
878
879       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
880       if (!tpcTrack)
881         continue;
882
883       if (AcceptTrack(tpcTrack)) {
884         acceptedTracks->Add(tpcTrack);
885       }
886       else
887         delete tpcTrack;
888     }
889     else
890     {
891       AliESDtrack* track = esd->GetTrack(iTrack);
892       if(AcceptTrack(track))
893         acceptedTracks->Add(track);
894     }
895   } 
896   if(bTPC)acceptedTracks->SetOwner(kTRUE);
897   return acceptedTracks;
898 }
899
900 //____________________________________________________________________
901 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
902 {
903   //
904   // returns an the number of tracks that pass the cuts
905   //
906
907   Int_t count = 0;
908
909   // loop over esd tracks
910   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
911     AliESDtrack* track = esd->GetTrack(iTrack);
912     if (AcceptTrack(track))
913       count++;
914   }
915
916   return count;
917 }
918
919 //____________________________________________________________________
920  void AliESDtrackCuts::DefineHistograms(Int_t color) {
921    // 
922    // diagnostics histograms are defined
923    // 
924
925    fHistogramsOn=kTRUE;
926    
927    Bool_t oldStatus = TH1::AddDirectoryStatus();
928    TH1::AddDirectory(kFALSE);
929    
930    //###################################################################################
931    // defining histograms
932
933    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
934
935    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
936    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
937
938    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
939   
940    for (Int_t i=0; i<kNCuts; i++) {
941      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
942      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
943      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
944    }
945
946   fhCutStatistics  ->SetLineColor(color);
947   fhCutCorrelation ->SetLineColor(color);
948   fhCutStatistics  ->SetLineWidth(2);
949   fhCutCorrelation ->SetLineWidth(2);
950
951   for (Int_t i=0; i<2; i++) {
952     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
953     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
954     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
955     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
956
957     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
958     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
959     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
960     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
961     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
962
963     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
964     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
965     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
966     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
967
968     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
969     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
970     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
971
972     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
973
974     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
975     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
976     
977     fhNClustersITS[i]->SetTitle("n ITS clusters");
978     fhNClustersTPC[i]->SetTitle("n TPC clusters");
979     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
980     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
981
982     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
983     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
984     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
985     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
986     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
987
988     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
989     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
990     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
991     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
992     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
993
994     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
995     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
996     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
997     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
998     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
999
1000     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1001     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1002     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1003     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1004
1005     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1006     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1007     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1008     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1009     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1010
1011     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1012     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1013     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1014
1015     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1016     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1017     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2); 
1018   }
1019
1020   // The number of sigmas to the vertex is per definition gaussian
1021   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1022   ffDTheoretical->SetParameter(0,1);
1023
1024   TH1::AddDirectory(oldStatus);
1025 }
1026
1027 //____________________________________________________________________
1028 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1029 {
1030   //
1031   // loads the histograms from a file
1032   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1033   //
1034
1035   if (!dir)
1036     dir = GetName();
1037
1038   if (!gDirectory->cd(dir))
1039     return kFALSE;
1040
1041   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1042
1043   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1044   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1045
1046   for (Int_t i=0; i<2; i++) {
1047     if (i==0)
1048     {
1049       gDirectory->cd("before_cuts");
1050     }
1051     else
1052       gDirectory->cd("after_cuts");
1053
1054     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1055     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1056     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1057     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1058
1059     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1060     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1061     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1062     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1063     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1064
1065     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1066     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1067     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1068     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1069
1070     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1071     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1072     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1073     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1074
1075     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1076     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1077
1078     gDirectory->cd("../");
1079   }
1080
1081   gDirectory->cd("..");
1082
1083   return kTRUE;
1084 }
1085
1086 //____________________________________________________________________
1087 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1088   //
1089   // saves the histograms in a directory (dir)
1090   //
1091
1092   if (!fHistogramsOn) {
1093     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1094     return;
1095   }
1096
1097   if (!dir)
1098     dir = GetName();
1099
1100   gDirectory->mkdir(dir);
1101   gDirectory->cd(dir);
1102
1103   gDirectory->mkdir("before_cuts");
1104   gDirectory->mkdir("after_cuts");
1105  
1106   // a factor of 2 is needed since n sigma is positive
1107   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1108   ffDTheoretical->Write("nSigmaToVertexTheory");
1109
1110   fhCutStatistics->Write();
1111   fhCutCorrelation->Write();
1112
1113   for (Int_t i=0; i<2; i++) {
1114     if (i==0)
1115       gDirectory->cd("before_cuts");
1116     else
1117       gDirectory->cd("after_cuts");
1118
1119     fhNClustersITS[i]        ->Write();
1120     fhNClustersTPC[i]        ->Write();
1121     fhChi2PerClusterITS[i]   ->Write();
1122     fhChi2PerClusterTPC[i]   ->Write();
1123
1124     fhC11[i]                 ->Write();
1125     fhC22[i]                 ->Write();
1126     fhC33[i]                 ->Write();
1127     fhC44[i]                 ->Write();
1128     fhC55[i]                 ->Write();
1129
1130     fhDXY[i]                 ->Write();
1131     fhDZ[i]                  ->Write();
1132     fhDXYDZ[i]               ->Write();
1133     fhDXYvsDZ[i]             ->Write();
1134
1135     fhDXYNormalized[i]       ->Write();
1136     fhDZNormalized[i]        ->Write();
1137     fhDXYvsDZNormalized[i]   ->Write();
1138     fhNSigmaToVertex[i]      ->Write();
1139
1140     fhPt[i]                  ->Write();
1141     fhEta[i]                 ->Write();
1142     
1143     gDirectory->cd("../");
1144   }
1145
1146   gDirectory->cd("../");
1147 }
1148
1149 //____________________________________________________________________
1150 void AliESDtrackCuts::DrawHistograms()
1151 {
1152   // draws some histograms
1153
1154   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1155   canvas1->Divide(2, 2);
1156
1157   canvas1->cd(1);
1158   fhNClustersTPC[0]->SetStats(kFALSE);
1159   fhNClustersTPC[0]->Draw();
1160
1161   canvas1->cd(2);
1162   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1163   fhChi2PerClusterTPC[0]->Draw();
1164
1165   canvas1->cd(3);
1166   fhNSigmaToVertex[0]->SetStats(kFALSE);
1167   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1168   fhNSigmaToVertex[0]->Draw();
1169
1170   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1171
1172   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1173   canvas2->Divide(3, 2);
1174
1175   canvas2->cd(1);
1176   fhC11[0]->SetStats(kFALSE);
1177   gPad->SetLogy();
1178   fhC11[0]->Draw();
1179
1180   canvas2->cd(2);
1181   fhC22[0]->SetStats(kFALSE);
1182   gPad->SetLogy();
1183   fhC22[0]->Draw();
1184
1185   canvas2->cd(3);
1186   fhC33[0]->SetStats(kFALSE);
1187   gPad->SetLogy();
1188   fhC33[0]->Draw();
1189
1190   canvas2->cd(4);
1191   fhC44[0]->SetStats(kFALSE);
1192   gPad->SetLogy();
1193   fhC44[0]->Draw();
1194
1195   canvas2->cd(5);
1196   fhC55[0]->SetStats(kFALSE);
1197   gPad->SetLogy();
1198   fhC55[0]->Draw();
1199
1200   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1201
1202   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1203   canvas3->Divide(3, 2);
1204
1205   canvas3->cd(1);
1206   fhDXY[0]->SetStats(kFALSE);
1207   gPad->SetLogy();
1208   fhDXY[0]->Draw();
1209
1210   canvas3->cd(2);
1211   fhDZ[0]->SetStats(kFALSE);
1212   gPad->SetLogy();
1213   fhDZ[0]->Draw();
1214
1215   canvas3->cd(3);
1216   fhDXYvsDZ[0]->SetStats(kFALSE);
1217   gPad->SetLogz();
1218   gPad->SetRightMargin(0.15);
1219   fhDXYvsDZ[0]->Draw("COLZ");
1220
1221   canvas3->cd(4);
1222   fhDXYNormalized[0]->SetStats(kFALSE);
1223   gPad->SetLogy();
1224   fhDXYNormalized[0]->Draw();
1225
1226   canvas3->cd(5);
1227   fhDZNormalized[0]->SetStats(kFALSE);
1228   gPad->SetLogy();
1229   fhDZNormalized[0]->Draw();
1230
1231   canvas3->cd(6);
1232   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1233   gPad->SetLogz();
1234   gPad->SetRightMargin(0.15);
1235   fhDXYvsDZNormalized[0]->Draw("COLZ");
1236
1237   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1238
1239   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1240   canvas4->Divide(2, 1);
1241
1242   canvas4->cd(1);
1243   fhCutStatistics->SetStats(kFALSE);
1244   fhCutStatistics->LabelsOption("v");
1245   gPad->SetBottomMargin(0.3);
1246   fhCutStatistics->Draw();
1247
1248   canvas4->cd(2);
1249   fhCutCorrelation->SetStats(kFALSE);
1250   fhCutCorrelation->LabelsOption("v");
1251   gPad->SetBottomMargin(0.3);
1252   gPad->SetLeftMargin(0.3);
1253   fhCutCorrelation->Draw("COLZ");
1254
1255   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1256
1257   /*canvas->cd(1);
1258   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1259   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1260
1261   canvas->cd(2);
1262   fhNClustersTPC[0]->SetStats(kFALSE);
1263   fhNClustersTPC[0]->DrawCopy();
1264
1265   canvas->cd(3);
1266   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1267   fhChi2PerClusterITS[0]->DrawCopy();
1268   fhChi2PerClusterITS[1]->SetLineColor(2);
1269   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1270
1271   canvas->cd(4);
1272   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1273   fhChi2PerClusterTPC[0]->DrawCopy();
1274   fhChi2PerClusterTPC[1]->SetLineColor(2);
1275   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1276 }
1277
1278 Float_t AliESDtrackCuts::GetMinNsigmaToVertex() const
1279 {
1280   // deprecated, please use GetMaxNsigmaToVertex
1281
1282   Printf("WARNING: AliESDtrackCuts::GetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use GetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
1283
1284   return GetMaxNsigmaToVertex();
1285 }
1286
1287 void AliESDtrackCuts::SetMinNsigmaToVertex(Float_t sigma)
1288 {
1289   // deprecated, will be removed in next release
1290
1291   Printf("WARNING: AliESDtrackCuts::SetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use SetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
1292   
1293   SetMaxNsigmaToVertex(sigma);
1294 }
1295
1296 void AliESDtrackCuts::SetAcceptKingDaughters(Bool_t b)
1297 {
1298   // deprecated, will be removed in next release
1299
1300   Printf("WARNING: AliESDtrackCuts::SetAcceptKingDaughters is DEPRECATED and will be removed in the next release. Please use SetAcceptKinkDaughters instead. Renaming was done to improve code readability.");
1301   
1302   SetAcceptKinkDaughters(b);
1303 }
1304
1305 Bool_t AliESDtrackCuts::GetAcceptKingDaughters() const
1306 {
1307   // deprecated, will be removed in next release
1308
1309   Printf("WARNING: AliESDtrackCuts::GetAcceptKingDaughters is DEPRECATED and will be removed in the next release. Please use GetAcceptKinkDaughters instead. Renaming was done to improve code readability.");
1310   
1311   return GetAcceptKinkDaughters();
1312 }