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