AliCentrality for ESD and AOD analysis
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDtrackCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliESDtrackCuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
17
18 #include "AliESDtrackCuts.h"
19
20 #include <AliESDtrack.h>
21 #include <AliESDVertex.h>
22 #include <AliESDEvent.h>
23 #include <AliLog.h>
24
25 #include <TTree.h>
26 #include <TCanvas.h>
27 #include <TDirectory.h>
28 #include <TH2F.h>
29 #include <TF1.h>
30
31 //____________________________________________________________________
32 ClassImp(AliESDtrackCuts)
33
34 // Cut names
35 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
36  "require TPC refit",
37  "require TPC standalone",
38  "require ITS refit",
39  "n clusters TPC",
40  "n clusters ITS",
41  "#Chi^{2}/cluster TPC",
42  "#Chi^{2}/cluster ITS",
43  "cov 11",
44  "cov 22",
45  "cov 33",
46  "cov 44",
47  "cov 55",
48  "trk-to-vtx",
49  "trk-to-vtx failed",
50  "kink daughters",
51  "p",
52  "p_{T}",
53  "p_{x}",
54  "p_{y}",
55  "p_{z}",
56  "eta",
57  "y",
58  "trk-to-vtx max dca 2D absolute",
59  "trk-to-vtx max dca xy absolute",
60  "trk-to-vtx max dca z absolute",
61  "trk-to-vtx min dca 2D absolute",
62  "trk-to-vtx min dca xy absolute",
63  "trk-to-vtx min dca z absolute",
64  "SPD cluster requirement",
65  "SDD cluster requirement",
66  "SSD cluster requirement",
67  "require ITS stand-alone",
68  "rel 1/pt uncertainty",
69  "require ITS Pid"
70 };
71
72 //____________________________________________________________________
73 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
74   fCutMinNClusterTPC(0),
75   fCutMinNClusterITS(0),
76   fCutMaxChi2PerClusterTPC(0),
77   fCutMaxChi2PerClusterITS(0),
78   fCutMaxC11(0),
79   fCutMaxC22(0),
80   fCutMaxC33(0),
81   fCutMaxC44(0),
82   fCutMaxC55(0),
83   fCutMaxRel1PtUncertainty(0),
84   fCutAcceptKinkDaughters(0),
85   fCutAcceptSharedTPCClusters(0),
86   fCutMaxFractionSharedTPCClusters(0),
87   fCutRequireTPCRefit(0),
88   fCutRequireTPCStandAlone(0),
89   fCutRequireITSRefit(0), 
90   fCutRequireITSPid(0),
91   fCutRequireITSStandAlone(0),
92   fCutRequireITSpureSA(0),
93   fCutNsigmaToVertex(0),
94   fCutSigmaToVertexRequired(0),
95   fCutMaxDCAToVertexXY(0),
96   fCutMaxDCAToVertexZ(0),
97   fCutMinDCAToVertexXY(0),
98   fCutMinDCAToVertexZ(0),
99   fCutMaxDCAToVertexXYPtDep(""),
100   fCutMaxDCAToVertexZPtDep(""),
101   fCutMinDCAToVertexXYPtDep(""),
102   fCutMinDCAToVertexZPtDep(""),
103   f1CutMaxDCAToVertexXYPtDep(0x0),
104   f1CutMaxDCAToVertexZPtDep(0x0),
105   f1CutMinDCAToVertexXYPtDep(0x0),
106   f1CutMinDCAToVertexZPtDep(0x0),
107   fCutDCAToVertex2D(0),
108   fPMin(0),
109   fPMax(0),
110   fPtMin(0),
111   fPtMax(0),
112   fPxMin(0),
113   fPxMax(0),
114   fPyMin(0),
115   fPyMax(0),
116   fPzMin(0),
117   fPzMax(0),
118   fEtaMin(0),
119   fEtaMax(0),
120   fRapMin(0),
121   fRapMax(0),
122   fHistogramsOn(0),
123   ffDTheoretical(0),
124   fhCutStatistics(0),         
125   fhCutCorrelation(0)
126 {
127   //
128   // constructor
129   //
130
131   Init();
132
133   //##############################################################################
134   // setting default cuts
135   SetMinNClustersTPC();
136   SetMinNClustersITS();
137   SetMaxChi2PerClusterTPC();
138   SetMaxChi2PerClusterITS();                                
139   SetMaxCovDiagonalElements();
140   SetMaxRel1PtUncertainty();
141   SetRequireTPCRefit();
142   SetRequireTPCStandAlone();
143   SetRequireITSRefit();
144   SetRequireITSPid(kFALSE);
145   SetRequireITSStandAlone(kFALSE);
146   SetRequireITSPureStandAlone(kFALSE);
147   SetAcceptKinkDaughters();
148   SetMaxNsigmaToVertex();
149   SetMaxDCAToVertexXY();
150   SetMaxDCAToVertexZ();
151   SetDCAToVertex2D();
152   SetMinDCAToVertexXY();
153   SetMinDCAToVertexZ();
154   SetPRange();
155   SetPtRange();
156   SetPxRange();
157   SetPyRange();
158   SetPzRange();
159   SetEtaRange();
160   SetRapRange();
161   SetClusterRequirementITS(kSPD);
162   SetClusterRequirementITS(kSDD);
163   SetClusterRequirementITS(kSSD);
164
165   SetHistogramsOn();
166 }
167
168 //_____________________________________________________________________________
169 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
170   fCutMinNClusterTPC(0),
171   fCutMinNClusterITS(0),
172   fCutMaxChi2PerClusterTPC(0),
173   fCutMaxChi2PerClusterITS(0),
174   fCutMaxC11(0),
175   fCutMaxC22(0),
176   fCutMaxC33(0),
177   fCutMaxC44(0),
178   fCutMaxC55(0),
179   fCutMaxRel1PtUncertainty(0),
180   fCutAcceptKinkDaughters(0),
181   fCutAcceptSharedTPCClusters(0),
182   fCutMaxFractionSharedTPCClusters(0),
183   fCutRequireTPCRefit(0),
184   fCutRequireTPCStandAlone(0),
185   fCutRequireITSRefit(0),
186   fCutRequireITSPid(0),
187   fCutRequireITSStandAlone(0),
188   fCutRequireITSpureSA(0),
189   fCutNsigmaToVertex(0),
190   fCutSigmaToVertexRequired(0),
191   fCutMaxDCAToVertexXY(0),
192   fCutMaxDCAToVertexZ(0),
193   fCutMinDCAToVertexXY(0),
194   fCutMinDCAToVertexZ(0),
195   fCutMaxDCAToVertexXYPtDep(""),
196   fCutMaxDCAToVertexZPtDep(""),
197   fCutMinDCAToVertexXYPtDep(""),
198   fCutMinDCAToVertexZPtDep(""),
199   f1CutMaxDCAToVertexXYPtDep(0x0),
200   f1CutMaxDCAToVertexZPtDep(0x0),
201   f1CutMinDCAToVertexXYPtDep(0x0),
202   f1CutMinDCAToVertexZPtDep(0x0),
203   fCutDCAToVertex2D(0),
204   fPMin(0),
205   fPMax(0),
206   fPtMin(0),
207   fPtMax(0),
208   fPxMin(0),
209   fPxMax(0),
210   fPyMin(0),
211   fPyMax(0),
212   fPzMin(0),
213   fPzMax(0),
214   fEtaMin(0),
215   fEtaMax(0),
216   fRapMin(0),
217   fRapMax(0),
218   fHistogramsOn(0),
219   ffDTheoretical(0),                                 
220   fhCutStatistics(0),         
221   fhCutCorrelation(0)
222 {
223   //
224   // copy constructor
225   //
226
227   ((AliESDtrackCuts &) c).Copy(*this);
228 }
229
230 AliESDtrackCuts::~AliESDtrackCuts()
231 {
232   //
233   // destructor
234   //
235
236   for (Int_t i=0; i<2; i++) {
237     
238     if (fhNClustersITS[i])
239       delete fhNClustersITS[i];            
240     if (fhNClustersTPC[i])
241       delete fhNClustersTPC[i];                
242     if (fhChi2PerClusterITS[i])
243       delete fhChi2PerClusterITS[i];       
244     if (fhChi2PerClusterTPC[i])
245       delete fhChi2PerClusterTPC[i];       
246     if (fhC11[i])
247       delete fhC11[i];                     
248     if (fhC22[i])
249       delete fhC22[i];                     
250     if (fhC33[i])
251       delete fhC33[i];                     
252     if (fhC44[i])
253       delete fhC44[i];                     
254     if (fhC55[i])
255       delete fhC55[i];
256
257     if (fhRel1PtUncertainty[i])
258       delete fhRel1PtUncertainty[i];
259     
260     if (fhDXY[i])
261       delete fhDXY[i];                     
262     if (fhDZ[i])
263       delete fhDZ[i];
264     if (fhDXYDZ[i])
265       delete fhDXYDZ[i];
266     if (fhDXYvsDZ[i])
267       delete fhDXYvsDZ[i];
268
269     if (fhDXYNormalized[i])
270       delete fhDXYNormalized[i];           
271     if (fhDZNormalized[i])
272       delete fhDZNormalized[i];
273     if (fhDXYvsDZNormalized[i])
274       delete fhDXYvsDZNormalized[i];
275     if (fhNSigmaToVertex[i])
276       delete fhNSigmaToVertex[i];
277     if (fhPt[i])
278       delete fhPt[i];
279     if (fhEta[i])
280       delete fhEta[i];
281   }
282
283   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
284   f1CutMaxDCAToVertexXYPtDep = 0;
285   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
286   f1CutMaxDCAToVertexZPtDep = 0;
287   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
288   f1CutMinDCAToVertexXYPtDep = 0;
289   if(f1CutMinDCAToVertexZPtDep)delete  f1CutMinDCAToVertexZPtDep; 
290   f1CutMinDCAToVertexZPtDep = 0;
291
292
293   if (ffDTheoretical)
294     delete ffDTheoretical;
295
296   if (fhCutStatistics)
297     delete fhCutStatistics;             
298   if (fhCutCorrelation)
299     delete fhCutCorrelation;            
300 }
301
302 void AliESDtrackCuts::Init()
303 {
304   //
305   // sets everything to zero
306   //
307
308   fCutMinNClusterTPC = 0;
309   fCutMinNClusterITS = 0;
310
311   fCutMaxChi2PerClusterTPC = 0;
312   fCutMaxChi2PerClusterITS = 0;
313   
314   for (Int_t i = 0; i < 3; i++)
315         fCutClusterRequirementITS[i] = kOff;
316
317   fCutMaxC11 = 0;
318   fCutMaxC22 = 0;
319   fCutMaxC33 = 0;
320   fCutMaxC44 = 0;
321   fCutMaxC55 = 0;
322   
323   fCutMaxRel1PtUncertainty = 0;
324
325   fCutAcceptKinkDaughters = 0;
326   fCutAcceptSharedTPCClusters = 0;
327   fCutMaxFractionSharedTPCClusters = 0;
328   fCutRequireTPCRefit = 0;
329   fCutRequireTPCStandAlone = 0;
330   fCutRequireITSRefit = 0;
331   fCutRequireITSPid = 0;
332   fCutRequireITSStandAlone = 0;
333   fCutRequireITSpureSA = 0;
334
335   fCutNsigmaToVertex = 0;
336   fCutSigmaToVertexRequired = 0;
337   fCutMaxDCAToVertexXY = 0;
338   fCutMaxDCAToVertexZ = 0;
339   fCutDCAToVertex2D = 0;
340   fCutMinDCAToVertexXY = 0;
341   fCutMinDCAToVertexZ = 0;
342   fCutMaxDCAToVertexXYPtDep = "";
343   fCutMaxDCAToVertexZPtDep = "";
344   fCutMinDCAToVertexXYPtDep = "";
345   fCutMinDCAToVertexZPtDep = "";
346
347   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
348   f1CutMaxDCAToVertexXYPtDep = 0;
349   if( f1CutMaxDCAToVertexXYPtDep) delete  f1CutMaxDCAToVertexXYPtDep;
350   f1CutMaxDCAToVertexXYPtDep = 0;
351   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
352   f1CutMaxDCAToVertexZPtDep = 0;
353   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
354   f1CutMinDCAToVertexXYPtDep = 0;
355   if(f1CutMinDCAToVertexZPtDep)delete f1CutMinDCAToVertexZPtDep;
356   f1CutMinDCAToVertexZPtDep = 0;
357
358   
359   fPMin = 0;
360   fPMax = 0;
361   fPtMin = 0;
362   fPtMax = 0;
363   fPxMin = 0;
364   fPxMax = 0;
365   fPyMin = 0;
366   fPyMax = 0;
367   fPzMin = 0;
368   fPzMax = 0;
369   fEtaMin = 0;
370   fEtaMax = 0;
371   fRapMin = 0;
372   fRapMax = 0;
373
374   fHistogramsOn = kFALSE;
375
376   for (Int_t i=0; i<2; ++i)
377   {
378     fhNClustersITS[i] = 0;
379     fhNClustersTPC[i] = 0;
380
381     fhChi2PerClusterITS[i] = 0;
382     fhChi2PerClusterTPC[i] = 0;
383
384     fhC11[i] = 0;
385     fhC22[i] = 0;
386     fhC33[i] = 0;
387     fhC44[i] = 0;
388     fhC55[i] = 0;
389
390     fhRel1PtUncertainty[i] = 0;
391
392     fhDXY[i] = 0;
393     fhDZ[i] = 0;
394     fhDXYDZ[i] = 0;
395     fhDXYvsDZ[i] = 0;
396
397     fhDXYNormalized[i] = 0;
398     fhDZNormalized[i] = 0;
399     fhDXYvsDZNormalized[i] = 0;
400     fhNSigmaToVertex[i] = 0;
401     
402     fhPt[i] = 0;
403     fhEta[i] = 0;
404   }
405   ffDTheoretical = 0;
406
407   fhCutStatistics = 0;
408   fhCutCorrelation = 0;
409 }
410
411 //_____________________________________________________________________________
412 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
413 {
414   //
415   // Assignment operator
416   //
417
418   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
419   return *this;
420 }
421
422 //_____________________________________________________________________________
423 void AliESDtrackCuts::Copy(TObject &c) const
424 {
425   //
426   // Copy function
427   //
428
429   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
430
431   target.Init();
432
433   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
434   target.fCutMinNClusterITS = fCutMinNClusterITS;
435
436   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
437   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
438
439   for (Int_t i = 0; i < 3; i++)
440     target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
441
442   target.fCutMaxC11 = fCutMaxC11;
443   target.fCutMaxC22 = fCutMaxC22;
444   target.fCutMaxC33 = fCutMaxC33;
445   target.fCutMaxC44 = fCutMaxC44;
446   target.fCutMaxC55 = fCutMaxC55;
447
448   target.fCutMaxRel1PtUncertainty = fCutMaxRel1PtUncertainty;
449
450   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
451   target.fCutAcceptSharedTPCClusters = fCutAcceptSharedTPCClusters;
452   target.fCutMaxFractionSharedTPCClusters = fCutMaxFractionSharedTPCClusters;
453   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
454   target.fCutRequireTPCStandAlone = fCutRequireTPCStandAlone;
455   target.fCutRequireITSRefit = fCutRequireITSRefit;
456   target.fCutRequireITSPid = fCutRequireITSPid;
457   target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
458   target.fCutRequireITSpureSA = fCutRequireITSpureSA;
459
460   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
461   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
462   target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
463   target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
464   target.fCutDCAToVertex2D = fCutDCAToVertex2D;
465   target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
466   target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
467
468   target.fCutMaxDCAToVertexXYPtDep = fCutMaxDCAToVertexXYPtDep;
469   target.SetMaxDCAToVertexXYPtDep(fCutMaxDCAToVertexXYPtDep.Data());
470
471   target.fCutMaxDCAToVertexZPtDep = fCutMaxDCAToVertexZPtDep;
472   target.SetMaxDCAToVertexZPtDep(fCutMaxDCAToVertexZPtDep.Data());
473
474   target.fCutMinDCAToVertexXYPtDep = fCutMinDCAToVertexXYPtDep;
475   target.SetMinDCAToVertexXYPtDep(fCutMinDCAToVertexXYPtDep.Data());
476
477   target.fCutMinDCAToVertexZPtDep = fCutMinDCAToVertexZPtDep;
478   target.SetMinDCAToVertexZPtDep(fCutMinDCAToVertexZPtDep.Data());
479
480   target.fPMin = fPMin;
481   target.fPMax = fPMax;
482   target.fPtMin = fPtMin;
483   target.fPtMax = fPtMax;
484   target.fPxMin = fPxMin;
485   target.fPxMax = fPxMax;
486   target.fPyMin = fPyMin;
487   target.fPyMax = fPyMax;
488   target.fPzMin = fPzMin;
489   target.fPzMax = fPzMax;
490   target.fEtaMin = fEtaMin;
491   target.fEtaMax = fEtaMax;
492   target.fRapMin = fRapMin;
493   target.fRapMax = fRapMax;
494
495   target.fHistogramsOn = fHistogramsOn;
496
497   for (Int_t i=0; i<2; ++i)
498   {
499     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
500     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
501
502     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
503     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
504
505     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
506     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
507     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
508     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
509     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
510
511     if (fhRel1PtUncertainty[i]) target.fhRel1PtUncertainty[i] = (TH1F*) fhRel1PtUncertainty[i]->Clone();
512
513     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
514     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
515     if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
516     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
517
518     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
519     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
520     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
521     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
522     
523     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
524     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
525   }
526   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
527
528   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
529   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
530
531   TNamed::Copy(c);
532 }
533
534 //_____________________________________________________________________________
535 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
536   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
537   // Returns the number of merged objects (including this)
538   if (!list)
539     return 0;
540   if (list->IsEmpty())
541     return 1;
542   if (!fHistogramsOn)
543     return 0;
544   TIterator* iter = list->MakeIterator();
545   TObject* obj;
546
547   // collection of measured and generated histograms
548   Int_t count = 0;
549   while ((obj = iter->Next())) {
550
551     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
552     if (entry == 0)
553       continue;
554
555     if (!entry->fHistogramsOn)
556       continue;
557
558     for (Int_t i=0; i<2; i++) {
559       
560       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );      
561       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     ); 
562                                                                     
563       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]); 
564       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]); 
565                                                                     
566       fhC11[i]               ->Add(entry->fhC11[i]              ); 
567       fhC22[i]               ->Add(entry->fhC22[i]              ); 
568       fhC33[i]               ->Add(entry->fhC33[i]              ); 
569       fhC44[i]               ->Add(entry->fhC44[i]              ); 
570       fhC55[i]               ->Add(entry->fhC55[i]              );
571
572       fhRel1PtUncertainty[i] ->Add(entry->fhRel1PtUncertainty[i]);
573                                                                     
574       fhDXY[i]               ->Add(entry->fhDXY[i]              ); 
575       fhDZ[i]                ->Add(entry->fhDZ[i]               ); 
576       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
577       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
578
579       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
580       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
581       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
582       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]); 
583
584       fhPt[i]                ->Add(entry->fhPt[i]); 
585       fhEta[i]               ->Add(entry->fhEta[i]); 
586     }      
587
588     fhCutStatistics  ->Add(entry->fhCutStatistics);        
589     fhCutCorrelation ->Add(entry->fhCutCorrelation);      
590
591     count++;
592   }
593   return count+1;
594 }
595
596 //____________________________________________________________________
597 AliESDtrackCuts* AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
598 {
599   // creates an AliESDtrackCuts object and fills it with standard (pre data-taking) values for TPC-only cuts
600   
601   Printf("AliESDtrackCuts::GetStandardTPCOnlyTrackCuts: Creating track cuts for TPC-only.");
602   
603   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
604   
605   esdTrackCuts->SetMinNClustersTPC(50);
606   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
607   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
608   
609   esdTrackCuts->SetMaxDCAToVertexZ(3.2);
610   esdTrackCuts->SetMaxDCAToVertexXY(2.4);
611   esdTrackCuts->SetDCAToVertex2D(kTRUE);
612   
613   return esdTrackCuts;
614 }
615
616 //____________________________________________________________________
617 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
618 {
619   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2009 data
620   
621   Printf("AliESDtrackCuts::GetStandardITSTPCTrackCuts: Creating track cuts for ITS+TPC.");
622   
623   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
624
625   // TPC  
626   esdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin
627   esdTrackCuts->SetMinNClustersTPC(70);
628   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
629   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
630   esdTrackCuts->SetRequireTPCRefit(kTRUE);
631   // ITS
632   esdTrackCuts->SetRequireITSRefit(kTRUE);
633   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
634                                          AliESDtrackCuts::kAny);
635   if(selPrimaries) {
636     // 7*(0.0050+0.0060/pt^0.9)
637     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");
638   }
639   esdTrackCuts->SetMaxDCAToVertexZ(1.e6);
640   esdTrackCuts->SetDCAToVertex2D(kFALSE);
641   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
642   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
643   
644   return esdTrackCuts;
645 }
646
647 //____________________________________________________________________
648 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(Bool_t selPrimaries)
649 {
650   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2010 data
651   
652   Printf("AliESDtrackCuts::GetStandardITSTPCTrackCuts: Creating track cuts for ITS+TPC.");
653   
654   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
655
656   // TPC  
657   esdTrackCuts->SetMinNClustersTPC(70);
658   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
659   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
660   esdTrackCuts->SetRequireTPCRefit(kTRUE);
661   // ITS
662   esdTrackCuts->SetRequireITSRefit(kTRUE);
663   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
664                                          AliESDtrackCuts::kAny);
665   if(selPrimaries) {
666     // 7*(0.0026+0.0050/pt^1.01)
667     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
668   }
669   esdTrackCuts->SetMaxDCAToVertexZ(2);
670   esdTrackCuts->SetDCAToVertex2D(kFALSE);
671   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
672   
673   return esdTrackCuts;
674 }
675
676 //____________________________________________________________________
677 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
678 {
679   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
680   
681   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
682   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
683   esdTrackCuts->SetRequireITSRefit(kTRUE); 
684   esdTrackCuts->SetMinNClustersITS(4);
685   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
686                                          AliESDtrackCuts::kAny);
687   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
688
689   if(selPrimaries) {
690     // 7*(0.0085+0.0026/pt^1.55)
691     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
692   }
693   if(useForPid){
694     esdTrackCuts->SetRequireITSPid(kTRUE);
695   }
696   return esdTrackCuts;
697 }
698
699 //____________________________________________________________________
700 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
701 {
702   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
703   
704   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
705   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
706   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
707   esdTrackCuts->SetRequireITSRefit(kTRUE); 
708   esdTrackCuts->SetMinNClustersITS(4);
709   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
710                                          AliESDtrackCuts::kAny);
711   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
712
713   if(selPrimaries) {
714     // 7*(0.0085+0.0026/pt^1.55)
715     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
716   }
717   if(useForPid){
718     esdTrackCuts->SetRequireITSPid(kTRUE);
719   }
720   return esdTrackCuts;
721 }
722
723
724 //____________________________________________________________________
725 Int_t AliESDtrackCuts::GetReferenceMultiplicity(AliESDEvent* esd, Bool_t tpcOnly)
726 {
727   // Gets reference multiplicity following the standard cuts and a defined fiducial volume
728   // tpcOnly = kTRUE -> consider TPC-only tracks
729   //         = kFALSE -> consider global tracks
730   
731   if (!tpcOnly)
732   {
733     Printf("AliESDtrackCuts::GetReferenceMultiplicity: Not implemented for global tracks!");
734     return -1;
735   }
736   
737   AliESDtrackCuts* esdTrackCuts = GetStandardTPCOnlyTrackCuts();
738   esdTrackCuts->SetEtaRange(-0.8, 0.8);
739   esdTrackCuts->SetPtRange(0.15);
740   
741   Int_t nTracks = esdTrackCuts->CountAcceptedTracks(esd);
742   
743   delete esdTrackCuts;
744   esdTrackCuts = 0;
745   
746   return nTracks;
747 }
748
749 //____________________________________________________________________
750 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
751 {
752   // Calculates the number of sigma to the vertex.
753
754   Float_t b[2];
755   Float_t bRes[2];
756   Float_t bCov[3];
757   esdTrack->GetImpactParameters(b,bCov);
758   
759   if (bCov[0]<=0 || bCov[2]<=0) {
760     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
761     bCov[0]=0; bCov[2]=0;
762   }
763   bRes[0] = TMath::Sqrt(bCov[0]);
764   bRes[1] = TMath::Sqrt(bCov[2]);
765
766   // -----------------------------------
767   // How to get to a n-sigma cut?
768   //
769   // The accumulated statistics from 0 to d is
770   //
771   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
772   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
773   //
774   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
775   // Can this be expressed in a different way?
776
777   if (bRes[0] == 0 || bRes[1] ==0)
778     return -1;
779
780   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
781
782   // work around precision problem
783   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
784   // 1e-15 corresponds to nsigma ~ 7.7
785   if (TMath::Exp(-d * d / 2) < 1e-15)
786     return 1000;
787
788   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
789   return nSigma;
790 }
791
792 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
793 {
794   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
795
796   tree->SetBranchStatus("fTracks.fFlags", 1);
797   tree->SetBranchStatus("fTracks.fITSncls", 1);
798   tree->SetBranchStatus("fTracks.fTPCncls", 1);
799   tree->SetBranchStatus("fTracks.fITSchi2", 1);
800   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
801   tree->SetBranchStatus("fTracks.fC*", 1);
802   tree->SetBranchStatus("fTracks.fD", 1);
803   tree->SetBranchStatus("fTracks.fZ", 1);
804   tree->SetBranchStatus("fTracks.fCdd", 1);
805   tree->SetBranchStatus("fTracks.fCdz", 1);
806   tree->SetBranchStatus("fTracks.fCzz", 1);
807   tree->SetBranchStatus("fTracks.fP*", 1);
808   tree->SetBranchStatus("fTracks.fR*", 1);
809   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
810 }
811
812 //____________________________________________________________________
813 Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) 
814 {
815   // 
816   // figure out if the tracks survives all the track cuts defined
817   //
818   // the different quality parameter and kinematic values are first
819   // retrieved from the track. then it is found out what cuts the
820   // track did not survive and finally the cuts are imposed.
821
822   // this function needs the following branches:
823   // fTracks.fFlags
824   // fTracks.fITSncls
825   // fTracks.fTPCncls
826   // fTracks.fITSchi2
827   // fTracks.fTPCchi2
828   // fTracks.fC   //GetExternalCovariance
829   // fTracks.fD   //GetImpactParameters
830   // fTracks.fZ   //GetImpactParameters
831   // fTracks.fCdd //GetImpactParameters
832   // fTracks.fCdz //GetImpactParameters
833   // fTracks.fCzz //GetImpactParameters
834   // fTracks.fP   //GetPxPyPz
835   // fTracks.fR   //GetMass
836   // fTracks.fP   //GetMass
837   // fTracks.fKinkIndexes
838
839   UInt_t status = esdTrack->GetStatus();
840
841   // getting quality parameters from the ESD track
842   Int_t nClustersITS = esdTrack->GetITSclusters(0);
843   Int_t nClustersTPC = -1;
844   if(fCutRequireTPCStandAlone) {
845     nClustersTPC = esdTrack->GetTPCNclsIter1();
846   }
847   else {
848     nClustersTPC = esdTrack->GetTPCclusters(0);
849   }
850
851   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
852   Float_t fracClustersTPCShared = -1.;
853
854   Float_t chi2PerClusterITS = -1;
855   Float_t chi2PerClusterTPC = -1;
856   if (nClustersITS!=0)
857     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
858   if (nClustersTPC!=0) {
859     if(fCutRequireTPCStandAlone) {
860       chi2PerClusterTPC = esdTrack->GetTPCchi2Iter1()/Float_t(nClustersTPC);
861     } else {
862       chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
863     }
864     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
865   }
866
867   Double_t extCov[15];
868   esdTrack->GetExternalCovariance(extCov);
869
870   // getting the track to vertex parameters
871   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
872       
873   Float_t b[2];
874   Float_t bCov[3];
875   esdTrack->GetImpactParameters(b,bCov);
876   if (bCov[0]<=0 || bCov[2]<=0) {
877     AliDebug(1, "Estimated b resolution lower or equal zero!");
878     bCov[0]=0; bCov[2]=0;
879   }
880
881
882   // set pt-dependent DCA cuts, if requested
883   SetPtDepDCACuts(esdTrack->Pt());
884
885
886   Float_t dcaToVertexXY = b[0];
887   Float_t dcaToVertexZ = b[1];
888
889   Float_t dcaToVertex = -1;
890   
891   if (fCutDCAToVertex2D)
892   {
893     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
894   }
895   else
896     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
897     
898   // getting the kinematic variables of the track
899   // (assuming the mass is known)
900   Double_t p[3];
901   esdTrack->GetPxPyPz(p);
902
903   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
904   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
905   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
906
907   //y-eta related calculations
908   Float_t eta = -100.;
909   Float_t y   = -100.;
910   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
911     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
912   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
913     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
914     
915   if (extCov[14] < 0) 
916   {
917     Printf("AliESDtrackCuts::AcceptTrack: WARNING: GetSigma1Pt2() returns negative value for external covariance matrix element fC[14]: %f. Corrupted track information, track will not be accepted!", extCov[14]);
918     return kFALSE;
919   }
920   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
921   
922   //########################################################################
923   // cut the track?
924   
925   Bool_t cuts[kNCuts];
926   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
927   
928   // track quality cuts
929   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
930     cuts[0]=kTRUE;
931   if (fCutRequireTPCStandAlone && (status&AliESDtrack::kTPCin)==0)
932     cuts[1]=kTRUE;
933   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
934     cuts[2]=kTRUE;
935   if (nClustersTPC<fCutMinNClusterTPC)
936     cuts[3]=kTRUE;
937   if (nClustersITS<fCutMinNClusterITS) 
938     cuts[4]=kTRUE;
939   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
940     cuts[5]=kTRUE; 
941   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
942     cuts[6]=kTRUE;
943   if (extCov[0]  > fCutMaxC11) 
944     cuts[7]=kTRUE;  
945   if (extCov[2]  > fCutMaxC22) 
946     cuts[8]=kTRUE;  
947   if (extCov[5]  > fCutMaxC33) 
948     cuts[9]=kTRUE;  
949   if (extCov[9]  > fCutMaxC44) 
950     cuts[10]=kTRUE;  
951   if (extCov[14]  > fCutMaxC55) 
952     cuts[11]=kTRUE;  
953   if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
954     cuts[12] = kTRUE;
955   // if n sigma could not be calculated
956   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
957     cuts[13]=kTRUE;
958   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
959     cuts[14]=kTRUE;
960   // track kinematics cut
961   if((momentum < fPMin) || (momentum > fPMax)) 
962     cuts[15]=kTRUE;
963   if((pt < fPtMin) || (pt > fPtMax)) 
964     cuts[16] = kTRUE;
965   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
966     cuts[17] = kTRUE;
967   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
968     cuts[18] = kTRUE;
969   if((p[2] < fPzMin) || (p[2] > fPzMax))
970     cuts[19] = kTRUE;
971   if((eta < fEtaMin) || (eta > fEtaMax))
972     cuts[20] = kTRUE;
973   if((y < fRapMin) || (y > fRapMax)) 
974     cuts[21] = kTRUE;
975   if (fCutDCAToVertex2D && dcaToVertex > 1)
976     cuts[22] = kTRUE;
977   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
978     cuts[23] = kTRUE;
979   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
980     cuts[24] = kTRUE;
981   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
982     cuts[25] = kTRUE;
983   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
984     cuts[26] = kTRUE;
985   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
986     cuts[27] = kTRUE;
987   
988   for (Int_t i = 0; i < 3; i++)
989     cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
990   
991   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
992     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
993       // TPC tracks
994       cuts[31] = kTRUE; 
995     }else{
996       // ITS standalone tracks
997       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
998         if(status & AliESDtrack::kITSpureSA) cuts[31] = kTRUE;
999       }else if(fCutRequireITSpureSA){
1000         if(!(status & AliESDtrack::kITSpureSA)) cuts[31] = kTRUE;
1001       }
1002     }
1003   }
1004
1005   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
1006      cuts[32] = kTRUE;
1007
1008   if (!fCutAcceptSharedTPCClusters && nClustersTPCShared!=0)
1009     cuts[33] = kTRUE;
1010
1011   if (fracClustersTPCShared > fCutMaxFractionSharedTPCClusters)
1012     cuts[34] = kTRUE;  
1013
1014   if(fCutRequireITSPid){
1015     UChar_t clumap=esdTrack->GetITSClusterMap();
1016     Int_t nPointsForPid=0;
1017     for(Int_t i=2; i<6; i++){
1018       if(clumap&(1<<i)) ++nPointsForPid;
1019     }
1020     if(nPointsForPid<3) cuts[35] = kTRUE;
1021   }
1022
1023   Bool_t cut=kFALSE;
1024   for (Int_t i=0; i<kNCuts; i++) 
1025     if (cuts[i]) {cut = kTRUE;}
1026
1027
1028   //########################################################################
1029   // filling histograms
1030   if (fHistogramsOn) {
1031     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1032     if (cut)
1033       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1034     
1035     for (Int_t i=0; i<kNCuts; i++) {
1036       if (cuts[i])
1037         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1038
1039       for (Int_t j=i; j<kNCuts; j++) {
1040         if (cuts[i] && cuts[j]) {
1041           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1042           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1043           fhCutCorrelation->Fill(xC, yC);
1044         }
1045       }
1046     }
1047   }
1048
1049   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1050   // the code is not in a function due to too many local variables that would need to be passed
1051
1052   for (Int_t id = 0; id < 2; id++)
1053   {
1054     // id = 0 --> before cut
1055     // id = 1 --> after cut
1056
1057     if (fHistogramsOn)
1058     {
1059       fhNClustersITS[id]->Fill(nClustersITS);
1060       fhNClustersTPC[id]->Fill(nClustersTPC);
1061       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1062       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1063
1064       fhC11[id]->Fill(extCov[0]);
1065       fhC22[id]->Fill(extCov[2]);
1066       fhC33[id]->Fill(extCov[5]);
1067       fhC44[id]->Fill(extCov[9]);
1068       fhC55[id]->Fill(extCov[14]);
1069
1070       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1071
1072       fhPt[id]->Fill(pt);
1073       fhEta[id]->Fill(eta);
1074
1075       Float_t bRes[2];
1076       bRes[0] = TMath::Sqrt(bCov[0]);
1077       bRes[1] = TMath::Sqrt(bCov[2]);
1078
1079       fhDZ[id]->Fill(b[1]);
1080       fhDXY[id]->Fill(b[0]);
1081       fhDXYDZ[id]->Fill(dcaToVertex);
1082       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1083
1084       if (bRes[0]!=0 && bRes[1]!=0) {
1085         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1086         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1087         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1088         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1089       }
1090     }
1091
1092     // cut the track
1093     if (cut)
1094       return kFALSE;
1095   }
1096
1097   return kTRUE;
1098 }
1099
1100 //____________________________________________________________________
1101 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1102 {
1103   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1104   
1105   switch (req)
1106   {
1107         case kOff:        return kTRUE;
1108         case kNone:       return !clusterL1 && !clusterL2;
1109         case kAny:        return clusterL1 || clusterL2;
1110         case kFirst:      return clusterL1;
1111         case kOnlyFirst:  return clusterL1 && !clusterL2;
1112         case kSecond:     return clusterL2;
1113         case kOnlySecond: return clusterL2 && !clusterL1;
1114         case kBoth:       return clusterL1 && clusterL2;
1115   }
1116   
1117   return kFALSE;
1118 }
1119
1120 //____________________________________________________________________
1121 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
1122 {
1123   
1124   // Utility function to 
1125   // create a TPC only track from the given esd track
1126   // 
1127   // IMPORTANT: The track has to be deleted by the user
1128   //
1129   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1130   // there are only missing propagations here that are needed for old data
1131   // this function will therefore become obsolete
1132   //
1133   // adapted from code provided by CKB
1134
1135   if (!esd->GetPrimaryVertexTPC())
1136     return 0; // No TPC vertex no TPC tracks
1137
1138   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1139     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1140
1141   AliESDtrack* track = esd->GetTrack(iTrack);
1142   if (!track)
1143     return 0;
1144
1145   AliESDtrack *tpcTrack = new AliESDtrack();
1146
1147   // only true if we have a tpc track
1148   if (!track->FillTPCOnlyTrack(*tpcTrack))
1149   {
1150     delete tpcTrack;
1151     return 0;
1152   }
1153
1154   // propagate to Vertex
1155   // not needed for normal reconstructed ESDs...
1156   // Double_t pTPC[2],covTPC[3];
1157   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
1158
1159   return tpcTrack;
1160 }
1161
1162 //____________________________________________________________________
1163 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
1164 {
1165   //
1166   // returns an array of all tracks that pass the cuts
1167   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1168   // tracks that pass the cut
1169
1170   TObjArray* acceptedTracks = new TObjArray();
1171
1172   // loop over esd tracks
1173   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1174     if(bTPC){
1175       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1176       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1177
1178       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1179       if (!tpcTrack)
1180         continue;
1181
1182       if (AcceptTrack(tpcTrack)) {
1183         acceptedTracks->Add(tpcTrack);
1184       }
1185       else
1186         delete tpcTrack;
1187     }
1188     else
1189     {
1190       AliESDtrack* track = esd->GetTrack(iTrack);
1191       if(AcceptTrack(track))
1192         acceptedTracks->Add(track);
1193     }
1194   } 
1195   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1196   return acceptedTracks;
1197 }
1198
1199 //____________________________________________________________________
1200 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
1201 {
1202   //
1203   // returns an the number of tracks that pass the cuts
1204   //
1205
1206   Int_t count = 0;
1207
1208   // loop over esd tracks
1209   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1210     AliESDtrack* track = esd->GetTrack(iTrack);
1211     if (AcceptTrack(track))
1212       count++;
1213   }
1214
1215   return count;
1216 }
1217
1218 //____________________________________________________________________
1219  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1220    // 
1221    // diagnostics histograms are defined
1222    // 
1223
1224    fHistogramsOn=kTRUE;
1225    
1226    Bool_t oldStatus = TH1::AddDirectoryStatus();
1227    TH1::AddDirectory(kFALSE);
1228    
1229    //###################################################################################
1230    // defining histograms
1231
1232    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1233
1234    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1235    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1236
1237    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1238   
1239    for (Int_t i=0; i<kNCuts; i++) {
1240      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1241      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1242      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1243    }
1244
1245   fhCutStatistics  ->SetLineColor(color);
1246   fhCutCorrelation ->SetLineColor(color);
1247   fhCutStatistics  ->SetLineWidth(2);
1248   fhCutCorrelation ->SetLineWidth(2);
1249
1250   for (Int_t i=0; i<2; i++) {
1251     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1252     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1253     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1254     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1255
1256     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1257     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1258     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1259     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1260     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1261
1262     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1263
1264     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1265     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1266     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1267     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1268
1269     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1270     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1271     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1272
1273     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1274
1275     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1276     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1277     
1278     fhNClustersITS[i]->SetTitle("n ITS clusters");
1279     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1280     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1281     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1282
1283     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1284     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1285     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1286     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1287     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1288
1289     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1290
1291     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1292     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1293     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1294     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1295     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1296
1297     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1298     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1299     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1300     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1301     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1302
1303     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1304     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1305     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1306     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1307
1308     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1309     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1310     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1311     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1312     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1313
1314     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1315
1316     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1317     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1318     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1319
1320     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1321     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1322     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2); 
1323   }
1324
1325   // The number of sigmas to the vertex is per definition gaussian
1326   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1327   ffDTheoretical->SetParameter(0,1);
1328
1329   TH1::AddDirectory(oldStatus);
1330 }
1331
1332 //____________________________________________________________________
1333 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1334 {
1335   //
1336   // loads the histograms from a file
1337   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1338   //
1339
1340   if (!dir)
1341     dir = GetName();
1342
1343   if (!gDirectory->cd(dir))
1344     return kFALSE;
1345
1346   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1347
1348   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1349   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1350
1351   for (Int_t i=0; i<2; i++) {
1352     if (i==0)
1353     {
1354       gDirectory->cd("before_cuts");
1355     }
1356     else
1357       gDirectory->cd("after_cuts");
1358
1359     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1360     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1361     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1362     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1363
1364     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1365     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1366     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1367     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1368     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1369
1370     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1371
1372     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1373     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1374     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1375     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1376
1377     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1378     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1379     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1380     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1381
1382     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1383     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1384
1385     gDirectory->cd("../");
1386   }
1387
1388   gDirectory->cd("..");
1389
1390   return kTRUE;
1391 }
1392
1393 //____________________________________________________________________
1394 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1395   //
1396   // saves the histograms in a directory (dir)
1397   //
1398
1399   if (!fHistogramsOn) {
1400     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1401     return;
1402   }
1403
1404   if (!dir)
1405     dir = GetName();
1406
1407   gDirectory->mkdir(dir);
1408   gDirectory->cd(dir);
1409
1410   gDirectory->mkdir("before_cuts");
1411   gDirectory->mkdir("after_cuts");
1412  
1413   // a factor of 2 is needed since n sigma is positive
1414   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1415   ffDTheoretical->Write("nSigmaToVertexTheory");
1416
1417   fhCutStatistics->Write();
1418   fhCutCorrelation->Write();
1419
1420   for (Int_t i=0; i<2; i++) {
1421     if (i==0)
1422       gDirectory->cd("before_cuts");
1423     else
1424       gDirectory->cd("after_cuts");
1425
1426     fhNClustersITS[i]        ->Write();
1427     fhNClustersTPC[i]        ->Write();
1428     fhChi2PerClusterITS[i]   ->Write();
1429     fhChi2PerClusterTPC[i]   ->Write();
1430
1431     fhC11[i]                 ->Write();
1432     fhC22[i]                 ->Write();
1433     fhC33[i]                 ->Write();
1434     fhC44[i]                 ->Write();
1435     fhC55[i]                 ->Write();
1436
1437     fhRel1PtUncertainty[i]   ->Write();
1438
1439     fhDXY[i]                 ->Write();
1440     fhDZ[i]                  ->Write();
1441     fhDXYDZ[i]               ->Write();
1442     fhDXYvsDZ[i]             ->Write();
1443
1444     fhDXYNormalized[i]       ->Write();
1445     fhDZNormalized[i]        ->Write();
1446     fhDXYvsDZNormalized[i]   ->Write();
1447     fhNSigmaToVertex[i]      ->Write();
1448
1449     fhPt[i]                  ->Write();
1450     fhEta[i]                 ->Write();
1451     
1452     gDirectory->cd("../");
1453   }
1454
1455   gDirectory->cd("../");
1456 }
1457
1458 //____________________________________________________________________
1459 void AliESDtrackCuts::DrawHistograms()
1460 {
1461   // draws some histograms
1462
1463   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1464   canvas1->Divide(2, 2);
1465
1466   canvas1->cd(1);
1467   fhNClustersTPC[0]->SetStats(kFALSE);
1468   fhNClustersTPC[0]->Draw();
1469
1470   canvas1->cd(2);
1471   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1472   fhChi2PerClusterTPC[0]->Draw();
1473
1474   canvas1->cd(3);
1475   fhNSigmaToVertex[0]->SetStats(kFALSE);
1476   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1477   fhNSigmaToVertex[0]->Draw();
1478
1479   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1480
1481   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1482   canvas2->Divide(3, 2);
1483
1484   canvas2->cd(1);
1485   fhC11[0]->SetStats(kFALSE);
1486   gPad->SetLogy();
1487   fhC11[0]->Draw();
1488
1489   canvas2->cd(2);
1490   fhC22[0]->SetStats(kFALSE);
1491   gPad->SetLogy();
1492   fhC22[0]->Draw();
1493
1494   canvas2->cd(3);
1495   fhC33[0]->SetStats(kFALSE);
1496   gPad->SetLogy();
1497   fhC33[0]->Draw();
1498
1499   canvas2->cd(4);
1500   fhC44[0]->SetStats(kFALSE);
1501   gPad->SetLogy();
1502   fhC44[0]->Draw();
1503
1504   canvas2->cd(5);
1505   fhC55[0]->SetStats(kFALSE);
1506   gPad->SetLogy();
1507   fhC55[0]->Draw();
1508
1509   canvas2->cd(6);
1510   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1511   gPad->SetLogy();
1512   fhRel1PtUncertainty[0]->Draw();
1513
1514   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1515
1516   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1517   canvas3->Divide(3, 2);
1518
1519   canvas3->cd(1);
1520   fhDXY[0]->SetStats(kFALSE);
1521   gPad->SetLogy();
1522   fhDXY[0]->Draw();
1523
1524   canvas3->cd(2);
1525   fhDZ[0]->SetStats(kFALSE);
1526   gPad->SetLogy();
1527   fhDZ[0]->Draw();
1528
1529   canvas3->cd(3);
1530   fhDXYvsDZ[0]->SetStats(kFALSE);
1531   gPad->SetLogz();
1532   gPad->SetRightMargin(0.15);
1533   fhDXYvsDZ[0]->Draw("COLZ");
1534
1535   canvas3->cd(4);
1536   fhDXYNormalized[0]->SetStats(kFALSE);
1537   gPad->SetLogy();
1538   fhDXYNormalized[0]->Draw();
1539
1540   canvas3->cd(5);
1541   fhDZNormalized[0]->SetStats(kFALSE);
1542   gPad->SetLogy();
1543   fhDZNormalized[0]->Draw();
1544
1545   canvas3->cd(6);
1546   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1547   gPad->SetLogz();
1548   gPad->SetRightMargin(0.15);
1549   fhDXYvsDZNormalized[0]->Draw("COLZ");
1550
1551   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1552
1553   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1554   canvas4->Divide(2, 1);
1555
1556   canvas4->cd(1);
1557   fhCutStatistics->SetStats(kFALSE);
1558   fhCutStatistics->LabelsOption("v");
1559   gPad->SetBottomMargin(0.3);
1560   fhCutStatistics->Draw();
1561
1562   canvas4->cd(2);
1563   fhCutCorrelation->SetStats(kFALSE);
1564   fhCutCorrelation->LabelsOption("v");
1565   gPad->SetBottomMargin(0.3);
1566   gPad->SetLeftMargin(0.3);
1567   fhCutCorrelation->Draw("COLZ");
1568
1569   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1570
1571   /*canvas->cd(1);
1572   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1573   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1574
1575   canvas->cd(2);
1576   fhNClustersTPC[0]->SetStats(kFALSE);
1577   fhNClustersTPC[0]->DrawCopy();
1578
1579   canvas->cd(3);
1580   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1581   fhChi2PerClusterITS[0]->DrawCopy();
1582   fhChi2PerClusterITS[1]->SetLineColor(2);
1583   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1584
1585   canvas->cd(4);
1586   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1587   fhChi2PerClusterTPC[0]->DrawCopy();
1588   fhChi2PerClusterTPC[1]->SetLineColor(2);
1589   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1590 }
1591 //--------------------------------------------------------------------------
1592 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
1593   //
1594   // set the pt-dependent DCA cuts
1595   //
1596
1597   if(f1CutMaxDCAToVertexXYPtDep) {
1598      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
1599   }
1600
1601   if(f1CutMaxDCAToVertexZPtDep) {
1602     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
1603   }
1604
1605   if(f1CutMinDCAToVertexXYPtDep) {
1606     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
1607   }
1608
1609   if(f1CutMinDCAToVertexZPtDep) {
1610     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
1611   }
1612
1613
1614   return;
1615 }
1616
1617
1618
1619 //--------------------------------------------------------------------------
1620 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
1621   //
1622   // Check the correctness of the string syntax
1623   //
1624   Bool_t retval=kTRUE;
1625
1626   if(!dist.Contains("pt")) {
1627     if(print) printf("AliESDtrackCuts::CheckPtDepDCA(): string must contain \"pt\"\n");
1628     retval= kFALSE;
1629   } 
1630   return retval;
1631 }
1632
1633  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
1634
1635    if(f1CutMaxDCAToVertexXYPtDep){
1636      delete f1CutMaxDCAToVertexXYPtDep;
1637      // resetiing both
1638      f1CutMaxDCAToVertexXYPtDep = 0;
1639      fCutMaxDCAToVertexXYPtDep = "";
1640    }
1641    if(!CheckPtDepDCA(dist,kTRUE)){
1642      return;
1643    }  
1644    fCutMaxDCAToVertexXYPtDep = dist;
1645    TString tmp(dist);
1646    tmp.ReplaceAll("pt","x");
1647    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
1648  
1649 }
1650
1651  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
1652
1653
1654    if(f1CutMaxDCAToVertexZPtDep){
1655      delete f1CutMaxDCAToVertexZPtDep;
1656      // resetiing both
1657      f1CutMaxDCAToVertexZPtDep = 0;
1658      fCutMaxDCAToVertexZPtDep = "";
1659    }
1660    if(!CheckPtDepDCA(dist,kTRUE))return;
1661      
1662    fCutMaxDCAToVertexZPtDep = dist;
1663    TString tmp(dist);
1664    tmp.ReplaceAll("pt","x");
1665    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
1666
1667    
1668 }
1669
1670
1671  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
1672
1673
1674    if(f1CutMinDCAToVertexXYPtDep){
1675      delete f1CutMinDCAToVertexXYPtDep;
1676      // resetiing both
1677      f1CutMinDCAToVertexXYPtDep = 0;
1678      fCutMinDCAToVertexXYPtDep = "";
1679    }
1680    if(!CheckPtDepDCA(dist,kTRUE))return;
1681
1682    fCutMinDCAToVertexXYPtDep = dist;
1683    TString tmp(dist);
1684    tmp.ReplaceAll("pt","x");
1685    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
1686
1687 }
1688
1689
1690  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
1691
1692    
1693
1694    if(f1CutMinDCAToVertexZPtDep){
1695      delete f1CutMinDCAToVertexZPtDep;
1696      // resetiing both
1697      f1CutMinDCAToVertexZPtDep = 0;
1698      fCutMinDCAToVertexZPtDep = "";
1699    }
1700    if(!CheckPtDepDCA(dist,kTRUE))return;
1701    fCutMinDCAToVertexZPtDep = dist;
1702    TString tmp(dist);
1703    tmp.ReplaceAll("pt","x");
1704    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
1705 }