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