b76fcba7f4539ad15fb71d43e8d83869b4962936
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / 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$ */
17
18 #include "AliESDtrackCuts.h"
19
20
21 #include <AliESDtrack.h>
22 #include <AliESD.h>
23 #include <AliLog.h>
24 #include <TTree.h>
25 #include <TDirectory.h>
26
27 //____________________________________________________________________
28 ClassImp(AliESDtrackCuts)
29
30 // Cut names
31 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
32  "require TPC refit",
33  "require ITS refit",
34  "n clusters TPC",
35  "n clusters ITS",
36  "#Chi^{2}/clusters TPC",
37  "#Chi^{2}/clusters ITS",
38  "cov 11",
39  "cov 22",
40  "cov 33",
41  "cov 44",
42  "cov 55",
43  "trk-to-vtx",
44  "trk-to-vtx failed",
45  "kink daughters",
46  "p",
47  "p_{T}",
48  "p_{x}",
49  "p_{y}",
50  "p_{z}",
51  "y",
52  "eta"
53 };
54
55 //____________________________________________________________________
56 AliESDtrackCuts::AliESDtrackCuts() : TNamed(),
57   fCutMinNClusterTPC(0),
58   fCutMinNClusterITS(0),
59   fCutMaxChi2PerClusterTPC(0),
60   fCutMaxChi2PerClusterITS(0),
61   fCutMaxC11(0),
62   fCutMaxC22(0),
63   fCutMaxC33(0),
64   fCutMaxC44(0),
65   fCutMaxC55(0),
66   fCutAcceptKinkDaughters(0),
67   fCutRequireTPCRefit(0),
68   fCutRequireITSRefit(0),
69   fCutNsigmaToVertex(0),
70   fCutSigmaToVertexRequired(0),
71   fPMin(0),
72   fPMax(0),
73   fPtMin(0),
74   fPtMax(0),
75   fPxMin(0),
76   fPxMax(0),
77   fPyMin(0),
78   fPyMax(0),
79   fPzMin(0),
80   fPzMax(0),
81   fEtaMin(0),
82   fEtaMax(0),
83   fRapMin(0),
84   fRapMax(0),
85   fHistogramsOn(0),
86   ffDTheoretical(0),                                 
87   fhCutStatistics(0),         
88   fhCutCorrelation(0)
89 {
90   //
91   // default constructor
92   //
93
94   Init();
95 }
96
97 //____________________________________________________________________
98 AliESDtrackCuts::AliESDtrackCuts(Char_t* name, Char_t* title) : TNamed(name,title),
99   fCutMinNClusterTPC(0),
100   fCutMinNClusterITS(0),
101   fCutMaxChi2PerClusterTPC(0),
102   fCutMaxChi2PerClusterITS(0),
103   fCutMaxC11(0),
104   fCutMaxC22(0),
105   fCutMaxC33(0),
106   fCutMaxC44(0),
107   fCutMaxC55(0),
108   fCutAcceptKinkDaughters(0),
109   fCutRequireTPCRefit(0),
110   fCutRequireITSRefit(0),
111   fCutNsigmaToVertex(0),
112   fCutSigmaToVertexRequired(0),
113   fPMin(0),
114   fPMax(0),
115   fPtMin(0),
116   fPtMax(0),
117   fPxMin(0),
118   fPxMax(0),
119   fPyMin(0),
120   fPyMax(0),
121   fPzMin(0),
122   fPzMax(0),
123   fEtaMin(0),
124   fEtaMax(0),
125   fRapMin(0),
126   fRapMax(0),
127   fHistogramsOn(0),
128   fhCutStatistics(0),         
129   fhCutCorrelation(0)
130 {
131   //
132   // constructor
133   //
134
135   Init();
136
137   //##############################################################################
138   // setting default cuts
139   SetMinNClustersTPC();
140   SetMinNClustersITS();
141   SetMaxChi2PerClusterTPC();
142   SetMaxChi2PerClusterITS();                                
143   SetMaxCovDiagonalElements();                                      
144   SetRequireTPCRefit();
145   SetRequireITSRefit();
146   SetAcceptKingDaughters();
147   SetMinNsigmaToVertex();
148   SetRequireSigmaToVertex();
149   SetPRange();
150   SetPtRange();
151   SetPxRange();
152   SetPyRange();
153   SetPzRange();
154   SetEtaRange();
155   SetRapRange();
156
157   SetHistogramsOn();
158 }
159
160 //_____________________________________________________________________________
161 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : TNamed(c),
162   fCutMinNClusterTPC(0),
163   fCutMinNClusterITS(0),
164   fCutMaxChi2PerClusterTPC(0),
165   fCutMaxChi2PerClusterITS(0),
166   fCutMaxC11(0),
167   fCutMaxC22(0),
168   fCutMaxC33(0),
169   fCutMaxC44(0),
170   fCutMaxC55(0),
171   fCutAcceptKinkDaughters(0),
172   fCutRequireTPCRefit(0),
173   fCutRequireITSRefit(0),
174   fCutNsigmaToVertex(0),
175   fCutSigmaToVertexRequired(0),
176   fPMin(0),
177   fPMax(0),
178   fPtMin(0),
179   fPtMax(0),
180   fPxMin(0),
181   fPxMax(0),
182   fPyMin(0),
183   fPyMax(0),
184   fPzMin(0),
185   fPzMax(0),
186   fEtaMin(0),
187   fEtaMax(0),
188   fRapMin(0),
189   fRapMax(0),
190   fHistogramsOn(0),
191   ffDTheoretical(0),                                 
192   fhCutStatistics(0),         
193   fhCutCorrelation(0)
194 {
195   //
196   // copy constructor
197   //
198
199   ((AliESDtrackCuts &) c).Copy(*this);
200 }
201
202 AliESDtrackCuts::~AliESDtrackCuts()
203 {
204   //
205   // destructor
206   //
207
208   for (Int_t i=0; i<2; i++) {
209     
210     if (fhNClustersITS[i])
211       delete fhNClustersITS[i];            
212     if (fhNClustersTPC[i])
213       delete fhNClustersTPC[i];                
214     if (fhChi2PerClusterITS[i])
215       delete fhChi2PerClusterITS[i];       
216     if (fhChi2PerClusterTPC[i])
217       delete fhChi2PerClusterTPC[i];       
218     if (fhC11[i])
219       delete fhC11[i];                     
220     if (fhC22[i])
221       delete fhC22[i];                     
222     if (fhC33[i])
223       delete fhC33[i];                     
224     if (fhC44[i])
225       delete fhC44[i];                     
226     if (fhC55[i])
227     delete fhC55[i];                     
228     
229     if (fhDXY[i])
230       delete fhDXY[i];                     
231     if (fhDZ[i])
232       delete fhDZ[i];                      
233     if (fhDXYvsDZ[i])
234       delete fhDXYvsDZ[i];                 
235     
236     if (fhDXYNormalized[i])
237       delete fhDXYNormalized[i];           
238     if (fhDZNormalized[i])
239       delete fhDZNormalized[i];            
240     if (fhDXYvsDZNormalized[i])
241       delete fhDXYvsDZNormalized[i];       
242     if (fhNSigmaToVertex[i])
243       delete fhNSigmaToVertex[i];
244   }
245   
246   if (ffDTheoretical)
247     delete ffDTheoretical;
248
249   if (fhCutStatistics)
250     delete fhCutStatistics;             
251   if (fhCutCorrelation)
252     delete fhCutCorrelation;            
253 }
254
255 void AliESDtrackCuts::Init()
256 {
257   //
258   // sets everything to zero
259   //
260
261   fCutMinNClusterTPC = 0;
262   fCutMinNClusterITS = 0;
263
264   fCutMaxChi2PerClusterTPC = 0;
265   fCutMaxChi2PerClusterITS = 0;
266
267   fCutMaxC11 = 0;
268   fCutMaxC22 = 0;
269   fCutMaxC33 = 0;
270   fCutMaxC44 = 0;
271   fCutMaxC55 = 0;
272
273   fCutAcceptKinkDaughters = 0;
274   fCutRequireTPCRefit = 0;
275   fCutRequireITSRefit = 0;
276
277   fCutNsigmaToVertex = 0;
278   fCutSigmaToVertexRequired = 0;
279
280   fPMin = 0;
281   fPMax = 0;
282   fPtMin = 0;
283   fPtMax = 0;
284   fPxMin = 0;
285   fPxMax = 0;
286   fPyMin = 0;
287   fPyMax = 0;
288   fPzMin = 0;
289   fPzMax = 0;
290   fEtaMin = 0;
291   fEtaMax = 0;
292   fRapMin = 0;
293   fRapMax = 0;
294
295   fHistogramsOn = kFALSE;
296
297   for (Int_t i=0; i<2; ++i)
298   {
299     fhNClustersITS[i] = 0;
300     fhNClustersTPC[i] = 0;
301
302     fhChi2PerClusterITS[i] = 0;
303     fhChi2PerClusterTPC[i] = 0;
304
305     fhC11[i] = 0;
306     fhC22[i] = 0;
307     fhC33[i] = 0;
308     fhC44[i] = 0;
309     fhC55[i] = 0;
310
311     fhDXY[i] = 0;
312     fhDZ[i] = 0;
313     fhDXYvsDZ[i] = 0;
314
315     fhDXYNormalized[i] = 0;
316     fhDZNormalized[i] = 0;
317     fhDXYvsDZNormalized[i] = 0;
318     fhNSigmaToVertex[i] = 0;
319   }
320   ffDTheoretical = 0;
321
322   fhCutStatistics = 0;
323   fhCutCorrelation = 0;
324 }
325
326 //_____________________________________________________________________________
327 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
328 {
329   //
330   // Assignment operator
331   //
332
333   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
334   return *this;
335 }
336
337 //_____________________________________________________________________________
338 void AliESDtrackCuts::Copy(TObject &c) const
339 {
340   //
341   // Copy function
342   //
343
344   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
345
346   target.Init();
347
348   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
349   target.fCutMinNClusterITS = fCutMinNClusterITS;
350
351   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
352   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
353
354   target.fCutMaxC11 = fCutMaxC11;
355   target.fCutMaxC22 = fCutMaxC22;
356   target.fCutMaxC33 = fCutMaxC33;
357   target.fCutMaxC44 = fCutMaxC44;
358   target.fCutMaxC55 = fCutMaxC55;
359
360   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
361   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
362   target.fCutRequireITSRefit = fCutRequireITSRefit;
363
364   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
365   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
366
367   target.fPMin = fPMin;
368   target.fPMax = fPMax;
369   target.fPtMin = fPtMin;
370   target.fPtMax = fPtMax;
371   target.fPxMin = fPxMin;
372   target.fPxMax = fPxMax;
373   target.fPyMin = fPyMin;
374   target.fPyMax = fPyMax;
375   target.fPzMin = fPzMin;
376   target.fPzMax = fPzMax;
377   target.fEtaMin = fEtaMin;
378   target.fEtaMax = fEtaMax;
379   target.fRapMin = fRapMin;
380   target.fRapMax = fRapMax;
381
382   target.fHistogramsOn = fHistogramsOn;
383
384   for (Int_t i=0; i<2; ++i)
385   {
386     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
387     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
388
389     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
390     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
391
392     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
393     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
394     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
395     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
396     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
397
398     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
399     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
400     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
401
402     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
403     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
404     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
405     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
406   }
407   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
408
409   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
410   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
411
412   TNamed::Copy(c);
413 }
414
415 //_____________________________________________________________________________
416 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
417   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
418   // Returns the number of merged objects (including this)
419
420   if (!list)
421     return 0;
422   
423   if (list->IsEmpty())
424     return 1;
425
426   if (!fHistogramsOn)
427     return 0;
428
429   TIterator* iter = list->MakeIterator();
430   TObject* obj;
431
432
433   // collection of measured and generated histograms
434   Int_t count = 0;
435   while ((obj = iter->Next())) {
436
437     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
438     if (entry == 0)
439       continue;
440
441     if (!entry->fHistogramsOn)
442       continue;
443     
444     for (Int_t i=0; i<2; i++) {
445       
446       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );      
447       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     ); 
448                                                                     
449       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]); 
450       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]); 
451                                                                     
452       fhC11[i]               ->Add(entry->fhC11[i]              ); 
453       fhC22[i]               ->Add(entry->fhC22[i]              ); 
454       fhC33[i]               ->Add(entry->fhC33[i]              ); 
455       fhC44[i]               ->Add(entry->fhC44[i]              ); 
456       fhC55[i]               ->Add(entry->fhC55[i]              ); 
457                                                                     
458       fhDXY[i]               ->Add(entry->fhDXY[i]              ); 
459       fhDZ[i]                ->Add(entry->fhDZ[i]               ); 
460       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          ); 
461                                                                     
462       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    ); 
463       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     ); 
464       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]); 
465       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]); 
466
467     }      
468
469     fhCutStatistics  ->Add(entry->fhCutStatistics);        
470     fhCutCorrelation ->Add(entry->fhCutCorrelation);      
471
472     count++;
473   }
474
475   return count+1;
476 }
477
478
479 //____________________________________________________________________
480 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
481 {
482   // Calculates the number of sigma to the vertex.
483
484   Float_t b[2];
485   Float_t bRes[2];
486   Float_t bCov[3];
487   esdTrack->GetImpactParameters(b,bCov);
488   if (bCov[0]<=0 || bCov[2]<=0) {
489     AliDebug(1, "Estimated b resolution lower or equal zero!");
490     bCov[0]=0; bCov[2]=0;
491   }
492   bRes[0] = TMath::Sqrt(bCov[0]);
493   bRes[1] = TMath::Sqrt(bCov[2]);
494
495   // -----------------------------------
496   // How to get to a n-sigma cut?
497   //
498   // The accumulated statistics from 0 to d is
499   //
500   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
501   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
502   //
503   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
504   // Can this be expressed in a different way?
505
506   if (bRes[0] == 0 || bRes[1] ==0)
507     return -1;
508
509   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
510
511   // stupid rounding problem screws up everything:
512   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
513   if (TMath::Exp(-d * d / 2) < 1e-10)
514     return 1000;
515
516   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
517   return d;
518 }
519
520 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
521 {
522   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
523
524   tree->SetBranchStatus("fTracks.fFlags", 1);
525   tree->SetBranchStatus("fTracks.fITSncls", 1);
526   tree->SetBranchStatus("fTracks.fTPCncls", 1);
527   tree->SetBranchStatus("fTracks.fITSchi2", 1);
528   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
529   tree->SetBranchStatus("fTracks.fC*", 1);
530   tree->SetBranchStatus("fTracks.fD", 1);
531   tree->SetBranchStatus("fTracks.fZ", 1);
532   tree->SetBranchStatus("fTracks.fCdd", 1);
533   tree->SetBranchStatus("fTracks.fCdz", 1);
534   tree->SetBranchStatus("fTracks.fCzz", 1);
535   tree->SetBranchStatus("fTracks.fP*", 1);
536   tree->SetBranchStatus("fTracks.fR*", 1);
537   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
538 }
539
540 //____________________________________________________________________
541 Bool_t
542 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
543   // 
544   // figure out if the tracks survives all the track cuts defined
545   //
546   // the different quality parameter and kinematic values are first
547   // retrieved from the track. then it is found out what cuts the
548   // track did not survive and finally the cuts are imposed.
549
550   // this function needs the following branches:
551   // fTracks.fFlags
552   // fTracks.fITSncls
553   // fTracks.fTPCncls
554   // fTracks.fITSchi2
555   // fTracks.fTPCchi2
556   // fTracks.fC   //GetExternalCovariance
557   // fTracks.fD   //GetImpactParameters
558   // fTracks.fZ   //GetImpactParameters
559   // fTracks.fCdd //GetImpactParameters
560   // fTracks.fCdz //GetImpactParameters
561   // fTracks.fCzz //GetImpactParameters
562   // fTracks.fP   //GetPxPyPz
563   // fTracks.fR   //GetMass
564   // fTracks.fP   //GetMass
565   // fTracks.fKinkIndexes
566
567   UInt_t status = esdTrack->GetStatus();
568
569   // dummy array
570   Int_t  fIdxInt[200];
571
572   // getting quality parameters from the ESD track
573   Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
574   Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
575   
576
577
578   Float_t chi2PerClusterITS = -1;
579   Float_t chi2PerClusterTPC = -1;
580   if (nClustersITS!=0)
581     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
582   if (nClustersTPC!=0)
583     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
584
585   Double_t extCov[15];
586   esdTrack->GetExternalCovariance(extCov);
587
588   // getting the track to vertex parameters
589   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
590
591   // getting the kinematic variables of the track
592   // (assuming the mass is known)
593   Double_t p[3];
594   esdTrack->GetPxPyPz(p);
595   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
596   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
597   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
598
599
600   //y-eta related calculations
601   Float_t eta = -100.;
602   Float_t y   = -100.;
603   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
604     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
605   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
606     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
607
608   
609   //########################################################################
610   // cut the track?
611   
612   Bool_t cuts[kNCuts];
613   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
614   
615   // track quality cuts
616   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
617     cuts[0]=kTRUE;
618   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
619     cuts[1]=kTRUE;
620   if (nClustersTPC<fCutMinNClusterTPC)
621     cuts[2]=kTRUE;
622   if (nClustersITS<fCutMinNClusterITS) 
623     cuts[3]=kTRUE;
624   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
625     cuts[4]=kTRUE; 
626   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
627     cuts[5]=kTRUE;
628   if (extCov[0]  > fCutMaxC11) 
629     cuts[6]=kTRUE;  
630   if (extCov[2]  > fCutMaxC22) 
631     cuts[7]=kTRUE;  
632   if (extCov[5]  > fCutMaxC33) 
633     cuts[8]=kTRUE;  
634   if (extCov[9]  > fCutMaxC44) 
635     cuts[9]=kTRUE;  
636   if (extCov[14]  > fCutMaxC55) 
637     cuts[10]=kTRUE;  
638   if (nSigmaToVertex > fCutNsigmaToVertex)
639     cuts[11] = kTRUE;
640   // if n sigma could not be calculated
641   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
642     cuts[12]=kTRUE;
643   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
644     cuts[13]=kTRUE;
645   // track kinematics cut
646   if((momentum < fPMin) || (momentum > fPMax)) 
647     cuts[14]=kTRUE;
648   if((pt < fPtMin) || (pt > fPtMax)) 
649     cuts[15] = kTRUE;
650   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
651     cuts[16] = kTRUE;
652   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
653     cuts[17] = kTRUE;
654   if((p[2] < fPzMin) || (p[2] > fPzMax))
655     cuts[18] = kTRUE;
656   if((eta < fEtaMin) || (eta > fEtaMax)) 
657     cuts[19] = kTRUE;
658   if((y < fRapMin) || (y > fRapMax)) 
659     cuts[20] = kTRUE;
660
661   Bool_t cut=kFALSE;
662   for (Int_t i=0; i<kNCuts; i++) 
663     if (cuts[i]) cut = kTRUE;
664   
665   //########################################################################
666   // filling histograms
667   if (fHistogramsOn) {
668     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
669     
670     if (cut)
671       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
672     
673     for (Int_t i=0; i<kNCuts; i++) {
674       if (cuts[i])
675         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
676       
677       for (Int_t j=i; j<kNCuts; j++) {
678         if (cuts[i] && cuts[j]) {
679           Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
680           Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
681           fhCutCorrelation->Fill(x,y);
682         }
683       }
684     }
685     
686
687     fhNClustersITS[0]->Fill(nClustersITS);
688     fhNClustersTPC[0]->Fill(nClustersTPC);
689     fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
690     fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
691
692     fhC11[0]->Fill(extCov[0]);
693     fhC22[0]->Fill(extCov[2]);
694     fhC33[0]->Fill(extCov[5]);
695     fhC44[0]->Fill(extCov[9]);
696     fhC55[0]->Fill(extCov[14]);
697
698     Float_t b[2];
699     Float_t bRes[2];
700     Float_t bCov[3];
701     esdTrack->GetImpactParameters(b,bCov);
702     if (bCov[0]<=0 || bCov[2]<=0) {
703       AliDebug(1, "Estimated b resolution lower or equal zero!");
704       bCov[0]=0; bCov[2]=0;
705     }
706     bRes[0] = TMath::Sqrt(bCov[0]);
707     bRes[1] = TMath::Sqrt(bCov[2]);
708
709     fhDZ[0]->Fill(b[1]);
710     fhDXY[0]->Fill(b[0]);
711     fhDXYvsDZ[0]->Fill(b[1],b[0]);
712
713     if (bRes[0]!=0 && bRes[1]!=0) {
714       fhDZNormalized[0]->Fill(b[1]/bRes[1]);
715       fhDXYNormalized[0]->Fill(b[0]/bRes[0]);
716       fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
717       fhNSigmaToVertex[0]->Fill(nSigmaToVertex);
718     }
719   }
720
721   //########################################################################
722   // cut the track!
723   if (cut) return kFALSE;
724
725   //########################################################################
726   // filling histograms after cut
727   if (fHistogramsOn) {
728     fhNClustersITS[1]->Fill(nClustersITS);
729     fhNClustersTPC[1]->Fill(nClustersTPC);
730     fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
731     fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
732
733     fhC11[1]->Fill(extCov[0]);
734     fhC22[1]->Fill(extCov[2]);
735     fhC33[1]->Fill(extCov[5]);
736     fhC44[1]->Fill(extCov[9]);
737     fhC55[1]->Fill(extCov[14]);
738
739     Float_t b[2];
740     Float_t bRes[2];
741     Float_t bCov[3];
742     esdTrack->GetImpactParameters(b,bCov);
743     if (bCov[0]<=0 || bCov[2]<=0) {
744       AliDebug(1, "Estimated b resolution lower or equal zero!");
745       bCov[0]=0; bCov[2]=0;
746     }
747     bRes[0] = TMath::Sqrt(bCov[0]);
748     bRes[1] = TMath::Sqrt(bCov[2]);
749
750     fhDZ[1]->Fill(b[1]);
751     fhDXY[1]->Fill(b[0]);
752     fhDXYvsDZ[1]->Fill(b[1],b[0]);
753
754     if (bRes[0]!=0 && bRes[1]!=0)
755     {
756       fhDZNormalized[1]->Fill(b[1]/bRes[1]);
757       fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
758       fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
759       fhNSigmaToVertex[1]->Fill(nSigmaToVertex);
760     }
761   }
762
763   return kTRUE;
764 }
765
766 //____________________________________________________________________
767 TObjArray*
768 AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
769 {
770   //
771   // returns an array of all tracks that pass the cuts
772   //
773
774   TObjArray* acceptedTracks = new TObjArray();
775
776   // loop over esd tracks
777   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
778     AliESDtrack* track = esd->GetTrack(iTrack);
779
780     if (AcceptTrack(track))
781       acceptedTracks->Add(track);
782   }
783
784   return acceptedTracks;
785 }
786
787 //____________________________________________________________________
788 Int_t
789 AliESDtrackCuts::CountAcceptedTracks(AliESD* esd)
790 {
791   //
792   // returns an the number of tracks that pass the cuts
793   //
794
795   Int_t count = 0;
796
797   // loop over esd tracks
798   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
799     AliESDtrack* track = esd->GetTrack(iTrack);
800
801     if (AcceptTrack(track))
802       count++;
803   }
804
805   return count;
806 }
807
808 //____________________________________________________________________
809  void AliESDtrackCuts::DefineHistograms(Int_t color) {
810    // 
811    // diagnostics histograms are defined
812    // 
813
814    fHistogramsOn=kTRUE;
815
816    //###################################################################################
817    // defining histograms
818
819    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
820
821    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
822    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
823
824    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
825   
826    for (Int_t i=0; i<kNCuts; i++) {
827      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
828      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
829      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
830    } 
831
832   fhCutStatistics  ->SetLineColor(color);
833   fhCutCorrelation ->SetLineColor(color);
834   fhCutStatistics  ->SetLineWidth(2);
835   fhCutCorrelation ->SetLineWidth(2);
836
837   Char_t str[256];
838   for (Int_t i=0; i<2; i++) {
839     if (i==0) sprintf(str," ");
840     else sprintf(str,"_cut");
841
842     fhNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
843     fhNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
844     fhChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
845     fhChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
846
847     fhC11[i]                 = new TH1F(Form("covMatrixDiagonal11%s",str),"",2000,0,20);
848     fhC22[i]                 = new  TH1F(Form("covMatrixDiagonal22%s",str),"",2000,0,20);
849     fhC33[i]                 = new  TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,1);
850     fhC44[i]                 = new  TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
851     fhC55[i]                 = new  TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
852
853     fhDXY[i]                 = new  TH1F(Form("dXY%s",str),"",500,-10,10);
854     fhDZ[i]                  = new  TH1F(Form("dZ%s",str),"",500,-10,10);
855     fhDXYvsDZ[i]             = new  TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
856
857     fhDXYNormalized[i]       = new  TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
858     fhDZNormalized[i]        = new  TH1F(Form("dZNormalized%s",str),"",500,-10,10);
859     fhDXYvsDZNormalized[i]   = new  TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
860
861     fhNSigmaToVertex[i]         = new  TH1F(Form("nSigmaToVertex%s",str),"",500,0,50);
862
863     fhNClustersITS[i]->SetXTitle("n ITS clusters");
864     fhNClustersTPC[i]->SetXTitle("n TPC clusters");
865     fhChi2PerClusterITS[i]->SetXTitle("#Chi^{2} per ITS cluster");
866     fhChi2PerClusterTPC[i]->SetXTitle("#Chi^{2} per TPC cluster");
867
868     fhC11[i]->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
869     fhC22[i]->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
870     fhC33[i]->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
871     fhC44[i]->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
872     fhC55[i]->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
873
874     fhDXY[i]->SetXTitle("transverse impact parameter");
875     fhDZ[i]->SetXTitle("longitudinal impact parameter");
876     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter");
877     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter");
878
879     fhDXYNormalized[i]->SetXTitle("normalized trans impact par");
880     fhDZNormalized[i]->SetXTitle("normalized long impact par");
881     fhDXYvsDZNormalized[i]->SetXTitle("normalized long impact par"); 
882     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par");
883     fhNSigmaToVertex[i]->SetXTitle("n #sigma to vertex");
884
885     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
886     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
887     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
888     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
889
890     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
891     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
892     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
893     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
894     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
895
896     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
897     fhDZ[i]->SetLineColor(color);   fhDZ[i]->SetLineWidth(2);
898
899     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
900     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
901     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2); 
902   }
903
904   // The number of sigmas to the vertex is per definition gaussian
905   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
906   ffDTheoretical->SetParameter(0,1);
907 }
908
909
910
911 //____________________________________________________________________
912 void 
913 AliESDtrackCuts::Print(const Option_t*) const {
914   //
915   // print method - still to be implemented
916   //
917
918   AliInfo("AliESDtrackCuts...");
919 }
920
921
922 //____________________________________________________________________
923 void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
924   // 
925   // saves the histograms in a directory (dir)
926   // 
927
928   
929   if (!fHistogramsOn) {
930     AliDebug(0, "Histograms not on - cannot save histograms!!!");
931     return;
932   }
933
934   gDirectory->mkdir(dir);
935   gDirectory->cd(dir);
936
937   gDirectory->mkdir("before_cuts");
938   gDirectory->mkdir("after_cuts");
939  
940   // a factor of 2 is needed since n sigma is positive
941   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
942   ffDTheoretical->Write("nSigmaToVertexTheory");
943
944   fhCutStatistics->Write();
945   fhCutCorrelation->Write();
946
947   for (Int_t i=0; i<2; i++) {
948     if (i==0)
949       gDirectory->cd("before_cuts");
950     else
951       gDirectory->cd("after_cuts");
952     
953     fhNClustersITS[i]        ->Write();
954     fhNClustersTPC[i]        ->Write();
955     fhChi2PerClusterITS[i]   ->Write();
956     fhChi2PerClusterTPC[i]   ->Write();
957     
958     fhC11[i]                 ->Write();
959     fhC22[i]                 ->Write();
960     fhC33[i]                 ->Write();
961     fhC44[i]                 ->Write();
962     fhC55[i]                 ->Write();
963
964     fhDXY[i]                 ->Write();
965     fhDZ[i]                  ->Write();
966     fhDXYvsDZ[i]             ->Write();
967     
968     fhDXYNormalized[i]       ->Write();
969     fhDZNormalized[i]        ->Write();
970     fhDXYvsDZNormalized[i]   ->Write();
971     fhNSigmaToVertex[i]      ->Write();
972
973     gDirectory->cd("../");
974   }
975
976   gDirectory->cd("../");
977 }
978
979
980