]> git.uio.no Git - u/mrichter/AliRoot.git/blob - ANALYSIS/AliESDtrackCuts.cxx
Adding possibility to cut on length over which the track is sampled in the active...
[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 <AliMultiplicity.h>
24 #include <AliLog.h>
25
26 #include <TTree.h>
27 #include <TCanvas.h>
28 #include <TDirectory.h>
29 #include <TH2F.h>
30 #include <TF1.h>
31 #include <TBits.h>
32
33 //____________________________________________________________________
34 ClassImp(AliESDtrackCuts)
35
36 // Cut names
37 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
38  "require TPC refit",
39  "require TPC standalone",
40  "require ITS refit",
41  "n clusters TPC",
42  "n clusters ITS",
43  "#Chi^{2}/cluster TPC",
44  "#Chi^{2}/cluster ITS",
45  "cov 11",
46  "cov 22",
47  "cov 33",
48  "cov 44",
49  "cov 55",
50  "trk-to-vtx",
51  "trk-to-vtx failed",
52  "kink daughters",
53  "p",
54  "p_{T}",
55  "p_{x}",
56  "p_{y}",
57  "p_{z}",
58  "eta",
59  "y",
60  "trk-to-vtx max dca 2D absolute",
61  "trk-to-vtx max dca xy absolute",
62  "trk-to-vtx max dca z absolute",
63  "trk-to-vtx min dca 2D absolute",
64  "trk-to-vtx min dca xy absolute",
65  "trk-to-vtx min dca z absolute",
66  "SPD cluster requirement",
67  "SDD cluster requirement",
68  "SSD cluster requirement",
69  "require ITS stand-alone",
70  "rel 1/pt uncertainty",
71  "TPC n shared clusters",
72  "TPC rel shared clusters",
73  "require ITS Pid",
74  "n crossed rows TPC",
75  "n crossed rows / n findable clusters",
76  "missing ITS points",
77  "#Chi^{2} TPC constrained vs. global",
78  "require TOF out",
79  "TOF Distance cut",
80  "min length in active volume TPC"
81 };
82
83 AliESDtrackCuts* AliESDtrackCuts::fgMultEstTrackCuts[AliESDtrackCuts::kNMultEstTrackCuts] = { 0, 0, 0, 0 };
84 Char_t AliESDtrackCuts::fgBeamTypeFlag = -1;
85
86 //____________________________________________________________________
87 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
88   fCutMinNClusterTPC(0),
89   fCutMinNClusterITS(0),
90   fCutMinNCrossedRowsTPC(0),
91   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
92   f1CutMinNClustersTPCPtDep(0x0),
93   fCutMaxPtDepNClustersTPC(0),
94   fCutMinLengthActiveVolumeTPC(0),
95   fCutMaxChi2PerClusterTPC(0),
96   fCutMaxChi2PerClusterITS(0),
97   fCutMaxChi2TPCConstrainedVsGlobal(0),
98   fCutMaxChi2TPCConstrainedVsGlobalVertexType(kVertexTracks | kVertexSPD),
99   fCutMaxMissingITSPoints(0),
100   fCutMaxC11(0),
101   fCutMaxC22(0),
102   fCutMaxC33(0),
103   fCutMaxC44(0),
104   fCutMaxC55(0),
105   fCutMaxRel1PtUncertainty(0),
106   fCutAcceptKinkDaughters(0),
107   fCutAcceptSharedTPCClusters(0),
108   fCutMaxFractionSharedTPCClusters(0),
109   fCutRequireTPCRefit(0),
110   fCutRequireTPCStandAlone(0),
111   fCutRequireITSRefit(0), 
112   fCutRequireITSPid(0),
113   fCutRequireITSStandAlone(0),
114   fCutRequireITSpureSA(0),
115   fCutNsigmaToVertex(0),
116   fCutSigmaToVertexRequired(0),
117   fCutMaxDCAToVertexXY(0),
118   fCutMaxDCAToVertexZ(0),
119   fCutMinDCAToVertexXY(0),
120   fCutMinDCAToVertexZ(0),
121   fCutMaxDCAToVertexXYPtDep(""),
122   fCutMaxDCAToVertexZPtDep(""),
123   fCutMinDCAToVertexXYPtDep(""),
124   fCutMinDCAToVertexZPtDep(""),
125   f1CutMaxDCAToVertexXYPtDep(0x0),
126   f1CutMaxDCAToVertexZPtDep(0x0),
127   f1CutMinDCAToVertexXYPtDep(0x0),
128   f1CutMinDCAToVertexZPtDep(0x0),
129   fCutDCAToVertex2D(0),
130   fPMin(0),
131   fPMax(0),
132   fPtMin(0),
133   fPtMax(0),
134   fPxMin(0),
135   fPxMax(0),
136   fPyMin(0),
137   fPyMax(0),
138   fPzMin(0),
139   fPzMax(0),
140   fEtaMin(0),
141   fEtaMax(0),
142   fRapMin(0),
143   fRapMax(0),
144   fCutRequireTOFout(kFALSE),
145   fFlagCutTOFdistance(kFALSE),
146   fCutTOFdistance(3.),
147   fHistogramsOn(0),
148   ffDTheoretical(0),
149   fhCutStatistics(0),         
150   fhCutCorrelation(0)
151 {
152   //
153   // constructor
154   //
155
156   Init();
157
158   //##############################################################################
159   // setting default cuts
160   SetMinNClustersTPC();
161   SetMinNClustersITS();
162   SetMinNCrossedRowsTPC();
163   SetMinRatioCrossedRowsOverFindableClustersTPC();
164   SetMaxChi2PerClusterTPC();
165   SetMaxChi2PerClusterITS();
166   SetMaxChi2TPCConstrainedGlobal();
167   SetMaxChi2TPCConstrainedGlobalVertexType();
168   SetMaxNOfMissingITSPoints();
169   SetMaxCovDiagonalElements();
170   SetMaxRel1PtUncertainty();
171   SetRequireTPCRefit();
172   SetRequireTPCStandAlone();
173   SetRequireITSRefit();
174   SetRequireITSPid(kFALSE);
175   SetRequireITSStandAlone(kFALSE);
176   SetRequireITSPureStandAlone(kFALSE);
177   SetAcceptKinkDaughters();
178   SetAcceptSharedTPCClusters();
179   SetMaxFractionSharedTPCClusters();
180   SetMaxNsigmaToVertex();
181   SetMaxDCAToVertexXY();
182   SetMaxDCAToVertexZ();
183   SetDCAToVertex2D();
184   SetMinDCAToVertexXY();
185   SetMinDCAToVertexZ();
186   SetPRange();
187   SetPtRange();
188   SetPxRange();
189   SetPyRange();
190   SetPzRange();
191   SetEtaRange();
192   SetRapRange();
193   SetClusterRequirementITS(kSPD);
194   SetClusterRequirementITS(kSDD);
195   SetClusterRequirementITS(kSSD);
196
197   SetHistogramsOn();
198 }
199
200 //_____________________________________________________________________________
201 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
202   fCutMinNClusterTPC(0),
203   fCutMinNClusterITS(0),
204   fCutMinNCrossedRowsTPC(0),
205   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
206   f1CutMinNClustersTPCPtDep(0x0),
207   fCutMaxPtDepNClustersTPC(0),
208   fCutMinLengthActiveVolumeTPC(0),
209   fCutMaxChi2PerClusterTPC(0),
210   fCutMaxChi2PerClusterITS(0),
211   fCutMaxChi2TPCConstrainedVsGlobal(0),
212   fCutMaxChi2TPCConstrainedVsGlobalVertexType(kVertexTracks | kVertexSPD),
213   fCutMaxMissingITSPoints(0),
214   fCutMaxC11(0),
215   fCutMaxC22(0),
216   fCutMaxC33(0),
217   fCutMaxC44(0),
218   fCutMaxC55(0),
219   fCutMaxRel1PtUncertainty(0),
220   fCutAcceptKinkDaughters(0),
221   fCutAcceptSharedTPCClusters(0),
222   fCutMaxFractionSharedTPCClusters(0),
223   fCutRequireTPCRefit(0),
224   fCutRequireTPCStandAlone(0),
225   fCutRequireITSRefit(0),
226   fCutRequireITSPid(0),
227   fCutRequireITSStandAlone(0),
228   fCutRequireITSpureSA(0),
229   fCutNsigmaToVertex(0),
230   fCutSigmaToVertexRequired(0),
231   fCutMaxDCAToVertexXY(0),
232   fCutMaxDCAToVertexZ(0),
233   fCutMinDCAToVertexXY(0),
234   fCutMinDCAToVertexZ(0),
235   fCutMaxDCAToVertexXYPtDep(""),
236   fCutMaxDCAToVertexZPtDep(""),
237   fCutMinDCAToVertexXYPtDep(""),
238   fCutMinDCAToVertexZPtDep(""),
239   f1CutMaxDCAToVertexXYPtDep(0x0),
240   f1CutMaxDCAToVertexZPtDep(0x0),
241   f1CutMinDCAToVertexXYPtDep(0x0),
242   f1CutMinDCAToVertexZPtDep(0x0),
243   fCutDCAToVertex2D(0),
244   fPMin(0),
245   fPMax(0),
246   fPtMin(0),
247   fPtMax(0),
248   fPxMin(0),
249   fPxMax(0),
250   fPyMin(0),
251   fPyMax(0),
252   fPzMin(0),
253   fPzMax(0),
254   fEtaMin(0),
255   fEtaMax(0),
256   fRapMin(0),
257   fRapMax(0),
258   fCutRequireTOFout(kFALSE),
259   fFlagCutTOFdistance(kFALSE),
260   fCutTOFdistance(3.),
261   fHistogramsOn(0),
262   ffDTheoretical(0),
263   fhCutStatistics(0),
264   fhCutCorrelation(0)
265 {
266   //
267   // copy constructor
268   //
269
270   ((AliESDtrackCuts &) c).Copy(*this);
271 }
272
273 AliESDtrackCuts::~AliESDtrackCuts()
274 {
275   //
276   // destructor
277   //
278
279   for (Int_t i=0; i<2; i++) {
280     
281     if (fhNClustersITS[i])
282       delete fhNClustersITS[i];            
283     if (fhNClustersTPC[i])
284       delete fhNClustersTPC[i];                
285     if (fhNSharedClustersTPC[i])
286       delete fhNSharedClustersTPC[i];                
287     if (fhNCrossedRowsTPC[i])
288       delete fhNCrossedRowsTPC[i];                
289     if (fhRatioCrossedRowsOverFindableClustersTPC[i])
290       delete fhRatioCrossedRowsOverFindableClustersTPC[i];                
291     if (fhChi2PerClusterITS[i])
292       delete fhChi2PerClusterITS[i];       
293     if (fhChi2PerClusterTPC[i])
294       delete fhChi2PerClusterTPC[i];    
295     if (fhChi2TPCConstrainedVsGlobal[i])
296       delete fhChi2TPCConstrainedVsGlobal[i];
297     if(fhNClustersForITSPID[i])
298       delete fhNClustersForITSPID[i];
299     if(fhNMissingITSPoints[i])
300       delete fhNMissingITSPoints[i];
301     if (fhC11[i])
302       delete fhC11[i];                     
303     if (fhC22[i])
304       delete fhC22[i];                     
305     if (fhC33[i])
306       delete fhC33[i];                     
307     if (fhC44[i])
308       delete fhC44[i];                     
309     if (fhC55[i])
310       delete fhC55[i];
311
312     if (fhRel1PtUncertainty[i])
313       delete fhRel1PtUncertainty[i];
314     
315     if (fhDXY[i])
316       delete fhDXY[i];                     
317     if (fhDZ[i])
318       delete fhDZ[i];
319     if (fhDXYDZ[i])
320       delete fhDXYDZ[i];
321     if (fhDXYvsDZ[i])
322       delete fhDXYvsDZ[i];
323
324     if (fhDXYNormalized[i])
325       delete fhDXYNormalized[i];           
326     if (fhDZNormalized[i])
327       delete fhDZNormalized[i];
328     if (fhDXYvsDZNormalized[i])
329       delete fhDXYvsDZNormalized[i];
330     if (fhNSigmaToVertex[i])
331       delete fhNSigmaToVertex[i];
332     if (fhPt[i])
333       delete fhPt[i];
334     if (fhEta[i])
335       delete fhEta[i];
336     if (fhTOFdistance[i])
337       delete fhTOFdistance[i];
338   }
339
340   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
341   f1CutMaxDCAToVertexXYPtDep = 0;
342   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
343   f1CutMaxDCAToVertexZPtDep = 0;
344   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
345   f1CutMinDCAToVertexXYPtDep = 0;
346   if(f1CutMinDCAToVertexZPtDep)delete  f1CutMinDCAToVertexZPtDep; 
347   f1CutMinDCAToVertexZPtDep = 0;
348
349
350   if (ffDTheoretical)
351     delete ffDTheoretical;
352
353   if (fhCutStatistics)
354     delete fhCutStatistics;             
355   if (fhCutCorrelation)
356     delete fhCutCorrelation;    
357
358   if(f1CutMinNClustersTPCPtDep)
359     delete f1CutMinNClustersTPCPtDep;
360         
361 }
362
363 void AliESDtrackCuts::Init()
364 {
365   //
366   // sets everything to zero
367   //
368
369   fCutMinNClusterTPC = 0;
370   fCutMinNClusterITS = 0;
371
372   fCutMaxChi2PerClusterTPC = 0;
373   fCutMaxChi2PerClusterITS = 0;
374   fCutMaxChi2TPCConstrainedVsGlobal = 0;
375   fCutMaxChi2TPCConstrainedVsGlobalVertexType = kVertexTracks | kVertexSPD;
376   fCutMaxMissingITSPoints  = 0;
377   
378   for (Int_t i = 0; i < 3; i++)
379         fCutClusterRequirementITS[i] = kOff;
380
381   fCutMaxC11 = 0;
382   fCutMaxC22 = 0;
383   fCutMaxC33 = 0;
384   fCutMaxC44 = 0;
385   fCutMaxC55 = 0;
386   
387   fCutMaxRel1PtUncertainty = 0;
388
389   fCutAcceptKinkDaughters = 0;
390   fCutAcceptSharedTPCClusters = 0;
391   fCutMaxFractionSharedTPCClusters = 0;
392   fCutRequireTPCRefit = 0;
393   fCutRequireTPCStandAlone = 0;
394   fCutRequireITSRefit = 0;
395   fCutRequireITSPid = 0;
396   fCutRequireITSStandAlone = 0;
397   fCutRequireITSpureSA = 0;
398
399   fCutNsigmaToVertex = 0;
400   fCutSigmaToVertexRequired = 0;
401   fCutMaxDCAToVertexXY = 0;
402   fCutMaxDCAToVertexZ = 0;
403   fCutDCAToVertex2D = 0;
404   fCutMinDCAToVertexXY = 0;
405   fCutMinDCAToVertexZ = 0;
406   fCutMaxDCAToVertexXYPtDep = "";
407   fCutMaxDCAToVertexZPtDep = "";
408   fCutMinDCAToVertexXYPtDep = "";
409   fCutMinDCAToVertexZPtDep = "";
410
411   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
412   f1CutMaxDCAToVertexXYPtDep = 0;
413   if( f1CutMaxDCAToVertexXYPtDep) delete  f1CutMaxDCAToVertexXYPtDep;
414   f1CutMaxDCAToVertexXYPtDep = 0;
415   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
416   f1CutMaxDCAToVertexZPtDep = 0;
417   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
418   f1CutMinDCAToVertexXYPtDep = 0;
419   if(f1CutMinDCAToVertexZPtDep)delete f1CutMinDCAToVertexZPtDep;
420   f1CutMinDCAToVertexZPtDep = 0;
421
422   
423   fPMin = 0;
424   fPMax = 0;
425   fPtMin = 0;
426   fPtMax = 0;
427   fPxMin = 0;
428   fPxMax = 0;
429   fPyMin = 0;
430   fPyMax = 0;
431   fPzMin = 0;
432   fPzMax = 0;
433   fEtaMin = 0;
434   fEtaMax = 0;
435   fRapMin = 0;
436   fRapMax = 0;
437
438   fHistogramsOn = kFALSE;
439
440   for (Int_t i=0; i<2; ++i)
441   {
442     fhNClustersITS[i] = 0;
443     fhNClustersTPC[i] = 0;
444     fhNSharedClustersTPC[i] = 0;
445     fhNCrossedRowsTPC[i] = 0;
446     fhRatioCrossedRowsOverFindableClustersTPC[i] = 0;
447
448     fhChi2PerClusterITS[i] = 0;
449     fhChi2PerClusterTPC[i] = 0;
450     fhChi2TPCConstrainedVsGlobal[i] = 0;
451     fhNClustersForITSPID[i] = 0;
452     fhNMissingITSPoints[i] = 0;
453
454     fhC11[i] = 0;
455     fhC22[i] = 0;
456     fhC33[i] = 0;
457     fhC44[i] = 0;
458     fhC55[i] = 0;
459
460     fhRel1PtUncertainty[i] = 0;
461
462     fhDXY[i] = 0;
463     fhDZ[i] = 0;
464     fhDXYDZ[i] = 0;
465     fhDXYvsDZ[i] = 0;
466
467     fhDXYNormalized[i] = 0;
468     fhDZNormalized[i] = 0;
469     fhDXYvsDZNormalized[i] = 0;
470     fhNSigmaToVertex[i] = 0;
471     
472     fhPt[i] = 0;
473     fhEta[i] = 0;
474     fhTOFdistance[i] = 0;
475   }
476   ffDTheoretical = 0;
477
478   fhCutStatistics = 0;
479   fhCutCorrelation = 0;
480 }
481
482 //_____________________________________________________________________________
483 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
484 {
485   //
486   // Assignment operator
487   //
488
489   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
490   return *this;
491 }
492
493 //_____________________________________________________________________________
494 void AliESDtrackCuts::Copy(TObject &c) const
495 {
496   //
497   // Copy function
498   //
499
500   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
501
502   target.Init();
503
504   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
505   target.fCutMinNClusterITS = fCutMinNClusterITS;
506   target.fCutMinNCrossedRowsTPC = fCutMinNCrossedRowsTPC;
507   target.fCutMinRatioCrossedRowsOverFindableClustersTPC = fCutMinRatioCrossedRowsOverFindableClustersTPC;
508   if(f1CutMinNClustersTPCPtDep){
509     target.f1CutMinNClustersTPCPtDep = (TFormula*) f1CutMinNClustersTPCPtDep->Clone("f1CutMinNClustersTPCPtDep");
510   }
511   target.fCutMaxPtDepNClustersTPC =   fCutMaxPtDepNClustersTPC;
512   target.fCutMinLengthActiveVolumeTPC = fCutMinLengthActiveVolumeTPC;
513   
514   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
515   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
516   target.fCutMaxChi2TPCConstrainedVsGlobal = fCutMaxChi2TPCConstrainedVsGlobal;
517   target.fCutMaxChi2TPCConstrainedVsGlobalVertexType = fCutMaxChi2TPCConstrainedVsGlobalVertexType;
518   target.fCutMaxMissingITSPoints = fCutMaxMissingITSPoints;
519
520   for (Int_t i = 0; i < 3; i++)
521     target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
522
523   target.fCutMaxC11 = fCutMaxC11;
524   target.fCutMaxC22 = fCutMaxC22;
525   target.fCutMaxC33 = fCutMaxC33;
526   target.fCutMaxC44 = fCutMaxC44;
527   target.fCutMaxC55 = fCutMaxC55;
528
529   target.fCutMaxRel1PtUncertainty = fCutMaxRel1PtUncertainty;
530
531   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
532   target.fCutAcceptSharedTPCClusters = fCutAcceptSharedTPCClusters;
533   target.fCutMaxFractionSharedTPCClusters = fCutMaxFractionSharedTPCClusters;
534   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
535   target.fCutRequireTPCStandAlone = fCutRequireTPCStandAlone;
536   target.fCutRequireITSRefit = fCutRequireITSRefit;
537   target.fCutRequireITSPid = fCutRequireITSPid;
538   target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
539   target.fCutRequireITSpureSA = fCutRequireITSpureSA;
540
541   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
542   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
543   target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
544   target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
545   target.fCutDCAToVertex2D = fCutDCAToVertex2D;
546   target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
547   target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
548
549   target.fCutMaxDCAToVertexXYPtDep = fCutMaxDCAToVertexXYPtDep;
550   if(fCutMaxDCAToVertexXYPtDep.Length()>0)target.SetMaxDCAToVertexXYPtDep(fCutMaxDCAToVertexXYPtDep.Data());
551
552   target.fCutMaxDCAToVertexZPtDep = fCutMaxDCAToVertexZPtDep;
553   if(fCutMaxDCAToVertexZPtDep.Length()>0)target.SetMaxDCAToVertexZPtDep(fCutMaxDCAToVertexZPtDep.Data());
554
555   target.fCutMinDCAToVertexXYPtDep = fCutMinDCAToVertexXYPtDep;
556   if(fCutMinDCAToVertexXYPtDep.Length()>0)target.SetMinDCAToVertexXYPtDep(fCutMinDCAToVertexXYPtDep.Data());
557
558   target.fCutMinDCAToVertexZPtDep = fCutMinDCAToVertexZPtDep;
559   if(fCutMinDCAToVertexZPtDep.Length()>0)target.SetMinDCAToVertexZPtDep(fCutMinDCAToVertexZPtDep.Data());
560
561   target.fPMin = fPMin;
562   target.fPMax = fPMax;
563   target.fPtMin = fPtMin;
564   target.fPtMax = fPtMax;
565   target.fPxMin = fPxMin;
566   target.fPxMax = fPxMax;
567   target.fPyMin = fPyMin;
568   target.fPyMax = fPyMax;
569   target.fPzMin = fPzMin;
570   target.fPzMax = fPzMax;
571   target.fEtaMin = fEtaMin;
572   target.fEtaMax = fEtaMax;
573   target.fRapMin = fRapMin;
574   target.fRapMax = fRapMax;
575
576   target.fFlagCutTOFdistance = fFlagCutTOFdistance;
577   target.fCutTOFdistance = fCutTOFdistance;
578   target.fCutRequireTOFout = fCutRequireTOFout;
579
580   target.fHistogramsOn = fHistogramsOn;
581
582   for (Int_t i=0; i<2; ++i)
583   {
584     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
585     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
586     if (fhNSharedClustersTPC[i]) target.fhNSharedClustersTPC[i] = (TH1F*) fhNSharedClustersTPC[i]->Clone();
587     if (fhNCrossedRowsTPC[i]) target.fhNCrossedRowsTPC[i] = (TH1F*) fhNCrossedRowsTPC[i]->Clone();
588     if (fhRatioCrossedRowsOverFindableClustersTPC[i]) target.fhRatioCrossedRowsOverFindableClustersTPC[i] = (TH1F*) fhRatioCrossedRowsOverFindableClustersTPC[i]->Clone();
589
590     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
591     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
592     if (fhChi2TPCConstrainedVsGlobal[i]) target.fhChi2TPCConstrainedVsGlobal[i] = (TH1F*) fhChi2TPCConstrainedVsGlobal[i]->Clone();
593     if (fhNClustersForITSPID[i]) target.fhNClustersForITSPID[i] = (TH1F*) fhNClustersForITSPID[i]->Clone();
594     if (fhNMissingITSPoints[i]) target.fhNMissingITSPoints[i] = (TH1F*) fhNMissingITSPoints[i]->Clone();
595
596     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
597     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
598     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
599     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
600     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
601
602     if (fhRel1PtUncertainty[i]) target.fhRel1PtUncertainty[i] = (TH1F*) fhRel1PtUncertainty[i]->Clone();
603
604     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
605     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
606     if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
607     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
608
609     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
610     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
611     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
612     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
613     
614     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
615     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
616     if (fhTOFdistance[i]) target.fhTOFdistance[i] = (TH2F*) fhTOFdistance[i]->Clone();
617   }
618   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
619
620   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
621   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
622
623   TNamed::Copy(c);
624 }
625
626 //_____________________________________________________________________________
627 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
628   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
629   // Returns the number of merged objects (including this)
630   if (!list)
631     return 0;
632   if (list->IsEmpty())
633     return 1;
634   if (!fHistogramsOn)
635     return 0;
636   TIterator* iter = list->MakeIterator();
637   TObject* obj;
638
639   // collection of measured and generated histograms
640   Int_t count = 0;
641   while ((obj = iter->Next())) {
642
643     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
644     if (entry == 0)
645       continue;
646
647     if (!entry->fHistogramsOn)
648       continue;
649
650     for (Int_t i=0; i<2; i++) {
651       
652       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );      
653       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     ); 
654       if (fhNSharedClustersTPC[i])
655         fhNSharedClustersTPC[i]      ->Add(entry->fhNSharedClustersTPC[i]     ); 
656       if (fhNCrossedRowsTPC[i])
657         fhNCrossedRowsTPC[i]   ->Add(entry->fhNCrossedRowsTPC[i]     ); 
658       if (fhRatioCrossedRowsOverFindableClustersTPC[i])
659         fhRatioCrossedRowsOverFindableClustersTPC[i]      ->Add(entry->fhRatioCrossedRowsOverFindableClustersTPC[i]     ); 
660                                                                     
661       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]); 
662       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]); 
663       if (fhChi2TPCConstrainedVsGlobal[i])
664         fhChi2TPCConstrainedVsGlobal[i]->Add(entry->fhChi2TPCConstrainedVsGlobal[i]);
665       if (fhNClustersForITSPID[i])
666         fhNClustersForITSPID[i]->Add(entry->fhNClustersForITSPID[i]);
667       if (fhNMissingITSPoints[i])
668         fhNMissingITSPoints[i] ->Add(entry->fhNMissingITSPoints[i]);
669
670       fhC11[i]               ->Add(entry->fhC11[i]              ); 
671       fhC22[i]               ->Add(entry->fhC22[i]              ); 
672       fhC33[i]               ->Add(entry->fhC33[i]              ); 
673       fhC44[i]               ->Add(entry->fhC44[i]              ); 
674       fhC55[i]               ->Add(entry->fhC55[i]              );
675
676       fhRel1PtUncertainty[i] ->Add(entry->fhRel1PtUncertainty[i]);
677                                                                     
678       fhDXY[i]               ->Add(entry->fhDXY[i]              ); 
679       fhDZ[i]                ->Add(entry->fhDZ[i]               ); 
680       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
681       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
682
683       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
684       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
685       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
686       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]); 
687
688       fhPt[i]                ->Add(entry->fhPt[i]); 
689       fhEta[i]               ->Add(entry->fhEta[i]); 
690       fhTOFdistance[i]       ->Add(entry->fhTOFdistance[i]); 
691     }      
692
693     fhCutStatistics  ->Add(entry->fhCutStatistics);        
694     fhCutCorrelation ->Add(entry->fhCutCorrelation);      
695
696     count++;
697   }
698   return count+1;
699 }
700
701 void AliESDtrackCuts::SetMinNClustersTPCPtDep(TFormula *f1, Float_t ptmax) 
702 {
703   //
704   // Sets the pT dependent NClustersTPC cut
705   //
706
707   if(f1){ 
708     delete f1CutMinNClustersTPCPtDep;
709     f1CutMinNClustersTPCPtDep = (TFormula*)f1->Clone("f1CutMinNClustersTPCPtDep"); 
710   }
711   fCutMaxPtDepNClustersTPC=ptmax; 
712 }
713
714 //____________________________________________________________________
715 AliESDtrackCuts* AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
716 {
717   // creates an AliESDtrackCuts object and fills it with standard (pre data-taking) values for TPC-only cuts
718   
719   AliInfoClass("Creating track cuts for TPC-only.");
720   
721   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
722   
723   esdTrackCuts->SetMinNClustersTPC(50);
724   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
725   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
726   
727   esdTrackCuts->SetMaxDCAToVertexZ(3.2);
728   esdTrackCuts->SetMaxDCAToVertexXY(2.4);
729   esdTrackCuts->SetDCAToVertex2D(kTRUE);
730   
731   return esdTrackCuts;
732 }
733
734 //____________________________________________________________________
735 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
736 {
737   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2009 data
738   
739   AliInfoClass("Creating track cuts for ITS+TPC (2009 definition).");
740   
741   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
742
743   // TPC  
744   esdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin
745   esdTrackCuts->SetMinNClustersTPC(70);
746   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
747   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
748   esdTrackCuts->SetRequireTPCRefit(kTRUE);
749   // ITS
750   esdTrackCuts->SetRequireITSRefit(kTRUE);
751   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
752                                          AliESDtrackCuts::kAny);
753   if(selPrimaries) {
754     // 7*(0.0050+0.0060/pt^0.9)
755     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");
756     esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
757   }
758   esdTrackCuts->SetMaxDCAToVertexZ(1.e6);
759   esdTrackCuts->SetDCAToVertex2D(kFALSE);
760   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
761   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
762   
763   esdTrackCuts->SetMaxChi2PerClusterITS(36);
764   
765   return esdTrackCuts;
766 }
767
768 //____________________________________________________________________
769 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(Bool_t selPrimaries, Int_t clusterCut)
770 {
771   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2011 data  
772   // if clusterCut = 1, the cut on the number of clusters is replaced by
773   // a cut on the number of crossed rows and on the ration crossed
774   // rows/findable clusters
775
776   AliInfoClass("Creating track cuts for ITS+TPC (2011 definition).");
777   
778   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
779
780   // TPC  
781   if(clusterCut == 0)  esdTrackCuts->SetMinNClustersTPC(50);
782   else if (clusterCut == 1) {
783     esdTrackCuts->SetMinNCrossedRowsTPC(70);
784     esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
785   }
786   else {
787     AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut));
788     esdTrackCuts->SetMinNClustersTPC(50);
789   }
790   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
791   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
792   esdTrackCuts->SetRequireTPCRefit(kTRUE);
793   // ITS
794   esdTrackCuts->SetRequireITSRefit(kTRUE);
795   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
796                                          AliESDtrackCuts::kAny);
797   if(selPrimaries) {
798     // 7*(0.0015+0.0050/pt^1.1)
799     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
800     esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
801   }
802   esdTrackCuts->SetMaxDCAToVertexZ(2);
803   esdTrackCuts->SetDCAToVertex2D(kFALSE);
804   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
805   
806   esdTrackCuts->SetMaxChi2PerClusterITS(36);
807
808   return esdTrackCuts;
809 }
810
811 //____________________________________________________________________
812 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(Bool_t selPrimaries,Int_t clusterCut)
813 {
814   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2010 data  
815   // if clusterCut = 1, the cut on the number of clusters is replaced by
816   // a cut on the number of crossed rows and on the ration crossed
817   // rows/findable clusters
818
819   AliInfoClass("Creating track cuts for ITS+TPC (2010 definition).");
820   
821   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
822
823   // TPC  
824   if(clusterCut == 0)  esdTrackCuts->SetMinNClustersTPC(70);
825   else if (clusterCut == 1) {
826     esdTrackCuts->SetMinNCrossedRowsTPC(70);
827     esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
828   }
829   else {
830     AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut));
831     esdTrackCuts->SetMinNClustersTPC(70);
832   }
833   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
834   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
835   esdTrackCuts->SetRequireTPCRefit(kTRUE);
836   // ITS
837   esdTrackCuts->SetRequireITSRefit(kTRUE);
838   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
839                                          AliESDtrackCuts::kAny);
840   if(selPrimaries) {
841     // 7*(0.0026+0.0050/pt^1.01)
842     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
843     esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
844   }
845   esdTrackCuts->SetMaxDCAToVertexZ(2);
846   esdTrackCuts->SetDCAToVertex2D(kFALSE);
847   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
848   
849   esdTrackCuts->SetMaxChi2PerClusterITS(36);
850
851   return esdTrackCuts;
852 }
853
854 //____________________________________________________________________
855 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
856 {
857   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
858   
859   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
860   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
861   esdTrackCuts->SetRequireITSRefit(kTRUE); 
862   esdTrackCuts->SetMinNClustersITS(4);
863   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
864                                          AliESDtrackCuts::kAny);
865   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
866
867   if(selPrimaries) {
868     // 7*(0.0085+0.0026/pt^1.55)
869     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
870   }
871   if(useForPid){
872     esdTrackCuts->SetRequireITSPid(kTRUE);
873   }
874   return esdTrackCuts;
875 }
876
877 //____________________________________________________________________
878 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
879 {
880   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks - pp 2010
881   
882   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
883   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
884   esdTrackCuts->SetRequireITSRefit(kTRUE); 
885   esdTrackCuts->SetMinNClustersITS(4);
886   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
887                                          AliESDtrackCuts::kAny);
888   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
889
890   if(selPrimaries) {
891     // 7*(0.0033+0.0045/pt^1.3)
892     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
893   }
894   if(useForPid){
895     esdTrackCuts->SetRequireITSPid(kTRUE);
896   }
897   return esdTrackCuts;
898 }
899
900 //____________________________________________________________________
901 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
902 {
903   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
904   
905   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
906   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
907   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
908   esdTrackCuts->SetRequireITSRefit(kTRUE); 
909   esdTrackCuts->SetMinNClustersITS(4);
910   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
911                                          AliESDtrackCuts::kAny);
912   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
913
914   if(selPrimaries) {
915     // 7*(0.0085+0.0026/pt^1.55)
916     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
917   }
918   if(useForPid){
919     esdTrackCuts->SetRequireITSPid(kTRUE);
920   }
921   return esdTrackCuts;
922 }
923
924 //____________________________________________________________________
925 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
926 {
927   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks --pp 2010
928   
929   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
930   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
931   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
932   esdTrackCuts->SetRequireITSRefit(kTRUE); 
933   esdTrackCuts->SetMinNClustersITS(4);
934   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
935                                          AliESDtrackCuts::kAny);
936   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
937
938   if(selPrimaries) {
939     // 7*(0.0033+0.0045/pt^1.3)
940     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
941   }
942   if(useForPid){
943     esdTrackCuts->SetRequireITSPid(kTRUE);
944   }
945   return esdTrackCuts;
946 }
947
948 //____________________________________________________________________
949 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(Bool_t selPrimaries, Bool_t useForPid)
950 {
951   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks -- PbPb 2010
952   
953   AliESDtrackCuts* esdTrackCuts = GetStandardITSSATrackCuts2010(selPrimaries, useForPid);
954   esdTrackCuts->SetMaxNOfMissingITSPoints(1);
955
956   return esdTrackCuts;
957 }
958 //____________________________________________________________________
959
960 AliESDtrackCuts* AliESDtrackCuts::GetStandardV0DaughterCuts()
961 {
962   // creates a AliESDtrackCuts object and fills it with standard cuts for V0 daughters
963   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
964   esdTrackCuts->SetRequireTPCRefit(kTRUE);
965   esdTrackCuts->SetMinNClustersTPC(70);
966   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
967   return esdTrackCuts;
968 }
969
970 //____________________________________________________________________
971 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, Bool_t tpcOnly)
972 {
973   // Gets reference multiplicity following the standard cuts and a defined fiducial volume
974   // tpcOnly = kTRUE -> consider TPC-only tracks
975   //         = kFALSE -> consider global tracks
976   //
977   // DEPRECATED Use GetReferenceMultiplicity with the enum as second argument instead
978   
979   if (!tpcOnly)
980   {
981     AliErrorClass("Not implemented for global tracks!");
982     return -1;
983   }
984   
985   static AliESDtrackCuts* esdTrackCuts = 0;
986   if (!esdTrackCuts)
987   {
988     esdTrackCuts = GetStandardTPCOnlyTrackCuts();
989     esdTrackCuts->SetEtaRange(-0.8, 0.8);
990     esdTrackCuts->SetPtRange(0.15);
991   }
992   
993   Int_t nTracks = esdTrackCuts->CountAcceptedTracks(esd);
994   
995   return nTracks;
996 }
997
998 //____________________________________________________________________
999 Float_t AliESDtrackCuts::GetSigmaToVertex(const AliESDtrack* const esdTrack)
1000 {
1001   // Calculates the number of sigma to the vertex.
1002
1003   Float_t b[2];
1004   Float_t bRes[2];
1005   Float_t bCov[3];
1006   esdTrack->GetImpactParameters(b,bCov);
1007   
1008   if (bCov[0]<=0 || bCov[2]<=0) {
1009     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
1010     bCov[0]=0; bCov[2]=0;
1011   }
1012   bRes[0] = TMath::Sqrt(bCov[0]);
1013   bRes[1] = TMath::Sqrt(bCov[2]);
1014
1015   // -----------------------------------
1016   // How to get to a n-sigma cut?
1017   //
1018   // The accumulated statistics from 0 to d is
1019   //
1020   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1021   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1022   //
1023   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1024   // Can this be expressed in a different way?
1025
1026   if (bRes[0] == 0 || bRes[1] ==0)
1027     return -1;
1028
1029   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1030
1031   // work around precision problem
1032   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1033   // 1e-15 corresponds to nsigma ~ 7.7
1034   if (TMath::Exp(-d * d / 2) < 1e-15)
1035     return 1000;
1036
1037   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1038   return nSigma;
1039 }
1040
1041 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
1042 {
1043   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
1044
1045   tree->SetBranchStatus("fTracks.fFlags", 1);
1046   tree->SetBranchStatus("fTracks.fITSncls", 1);
1047   tree->SetBranchStatus("fTracks.fTPCncls", 1);
1048   tree->SetBranchStatus("fTracks.fITSchi2", 1);
1049   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
1050   tree->SetBranchStatus("fTracks.fC*", 1);
1051   tree->SetBranchStatus("fTracks.fD", 1);
1052   tree->SetBranchStatus("fTracks.fZ", 1);
1053   tree->SetBranchStatus("fTracks.fCdd", 1);
1054   tree->SetBranchStatus("fTracks.fCdz", 1);
1055   tree->SetBranchStatus("fTracks.fCzz", 1);
1056   tree->SetBranchStatus("fTracks.fP*", 1);
1057   tree->SetBranchStatus("fTracks.fR*", 1);
1058   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
1059 }
1060
1061 //____________________________________________________________________
1062 Bool_t AliESDtrackCuts::AcceptTrack(const AliESDtrack* esdTrack) 
1063 {
1064   // 
1065   // figure out if the tracks survives all the track cuts defined
1066   //
1067   // the different quality parameter and kinematic values are first
1068   // retrieved from the track. then it is found out what cuts the
1069   // track did not survive and finally the cuts are imposed.
1070
1071   // this function needs the following branches:
1072   // fTracks.fFlags
1073   // fTracks.fITSncls
1074   // fTracks.fTPCncls
1075   // fTracks.fITSchi2
1076   // fTracks.fTPCchi2
1077   // fTracks.fC   //GetExternalCovariance
1078   // fTracks.fD   //GetImpactParameters
1079   // fTracks.fZ   //GetImpactParameters
1080   // fTracks.fCdd //GetImpactParameters
1081   // fTracks.fCdz //GetImpactParameters
1082   // fTracks.fCzz //GetImpactParameters
1083   // fTracks.fP   //GetPxPyPz
1084   // fTracks.fR   //GetMass
1085   // fTracks.fP   //GetMass
1086   // fTracks.fKinkIndexes
1087   //
1088   // esdEvent is only required for the MaxChi2TPCConstrainedVsGlobal
1089
1090   UInt_t status = esdTrack->GetStatus();
1091
1092   // getting quality parameters from the ESD track
1093   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1094   Int_t nClustersTPC = -1;
1095   if(fCutRequireTPCStandAlone) {
1096     nClustersTPC = esdTrack->GetTPCNclsIter1();
1097   }
1098   else {
1099     nClustersTPC = esdTrack->GetTPCclusters(0);
1100   }
1101   
1102   //Pt dependent NClusters Cut
1103   if(f1CutMinNClustersTPCPtDep) {
1104     if(esdTrack->Pt()<fCutMaxPtDepNClustersTPC)
1105       fCutMinNClusterTPC = f1CutMinNClustersTPCPtDep->Eval(esdTrack->Pt());
1106     else
1107       fCutMinNClusterTPC = f1CutMinNClustersTPCPtDep->Eval(fCutMaxPtDepNClustersTPC);
1108   }
1109
1110   Float_t nCrossedRowsTPC = esdTrack->GetTPCCrossedRows();
1111   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1112   if (esdTrack->GetTPCNclsF()>0) {
1113     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / esdTrack->GetTPCNclsF();
1114   }
1115   
1116   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
1117   Float_t fracClustersTPCShared = -1.;
1118
1119   Float_t chi2PerClusterITS = -1;
1120   Float_t chi2PerClusterTPC = -1;
1121   if (nClustersITS!=0)
1122     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1123   if (nClustersTPC!=0) {
1124     if(fCutRequireTPCStandAlone) {
1125       chi2PerClusterTPC = esdTrack->GetTPCchi2Iter1()/Float_t(nClustersTPC);
1126     } else {
1127       chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1128     }
1129     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
1130   }
1131
1132   Double_t extCov[15];
1133   esdTrack->GetExternalCovariance(extCov);
1134
1135   Float_t b[2];
1136   Float_t bCov[3];
1137   esdTrack->GetImpactParameters(b,bCov);
1138   if (bCov[0]<=0 || bCov[2]<=0) {
1139     AliDebug(1, "Estimated b resolution lower or equal zero!");
1140     bCov[0]=0; bCov[2]=0;
1141   }
1142
1143
1144   // set pt-dependent DCA cuts, if requested
1145   SetPtDepDCACuts(esdTrack->Pt());
1146
1147
1148   Float_t dcaToVertexXY = b[0];
1149   Float_t dcaToVertexZ = b[1];
1150
1151   Float_t dcaToVertex = -1;
1152   
1153   if (fCutDCAToVertex2D)
1154   {
1155     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1156   }
1157   else
1158     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1159     
1160   // getting the kinematic variables of the track
1161   // (assuming the mass is known)
1162   Double_t p[3];
1163   esdTrack->GetPxPyPz(p);
1164
1165   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
1166   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
1167   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
1168
1169   //y-eta related calculations
1170   Float_t eta = -100.;
1171   Float_t y   = -100.;
1172   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
1173     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
1174   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
1175     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
1176     
1177   if (extCov[14] < 0) 
1178   {
1179     AliWarning(Form("GetSigma1Pt2() returns negative value for external covariance matrix element fC[14]: %f. Corrupted track information, track will not be accepted!", extCov[14]));
1180     return kFALSE;
1181   }
1182   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
1183   
1184   //########################################################################
1185   // cut the track?
1186   
1187   Bool_t cuts[kNCuts];
1188   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1189   
1190   // track quality cuts
1191   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1192     cuts[0]=kTRUE;
1193   if (fCutRequireTPCStandAlone && (status&AliESDtrack::kTPCin)==0)
1194     cuts[1]=kTRUE;
1195   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1196     cuts[2]=kTRUE;
1197   if (nClustersTPC<fCutMinNClusterTPC)
1198     cuts[3]=kTRUE;
1199   if (nClustersITS<fCutMinNClusterITS) 
1200     cuts[4]=kTRUE;
1201   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
1202     cuts[5]=kTRUE; 
1203   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
1204     cuts[6]=kTRUE;
1205   if (extCov[0]  > fCutMaxC11) 
1206     cuts[7]=kTRUE;  
1207   if (extCov[2]  > fCutMaxC22) 
1208     cuts[8]=kTRUE;  
1209   if (extCov[5]  > fCutMaxC33) 
1210     cuts[9]=kTRUE;  
1211   if (extCov[9]  > fCutMaxC44) 
1212     cuts[10]=kTRUE;  
1213   if (extCov[14]  > fCutMaxC55) 
1214     cuts[11]=kTRUE;  
1215   
1216   // cut 12 and 13 see below
1217   
1218   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1219     cuts[14]=kTRUE;
1220   // track kinematics cut
1221   if((momentum < fPMin) || (momentum > fPMax)) 
1222     cuts[15]=kTRUE;
1223   if((pt < fPtMin) || (pt > fPtMax)) 
1224     cuts[16] = kTRUE;
1225   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
1226     cuts[17] = kTRUE;
1227   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
1228     cuts[18] = kTRUE;
1229   if((p[2] < fPzMin) || (p[2] > fPzMax))
1230     cuts[19] = kTRUE;
1231   if((eta < fEtaMin) || (eta > fEtaMax))
1232     cuts[20] = kTRUE;
1233   if((y < fRapMin) || (y > fRapMax)) 
1234     cuts[21] = kTRUE;
1235   if (fCutDCAToVertex2D && dcaToVertex > 1)
1236     cuts[22] = kTRUE;
1237   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1238     cuts[23] = kTRUE;
1239   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1240     cuts[24] = kTRUE;
1241   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
1242     cuts[25] = kTRUE;
1243   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
1244     cuts[26] = kTRUE;
1245   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
1246     cuts[27] = kTRUE;
1247   
1248   for (Int_t i = 0; i < 3; i++) {
1249     if(!(esdTrack->GetStatus()&AliESDtrack::kITSupg)) { // current ITS
1250       cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
1251     } else { // upgraded ITS (7 layers)
1252       // at the moment, for L012 the layers 12 are considered together
1253       if(i==0) { // L012
1254         cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(0), (esdTrack->HasPointOnITSLayer(1))&(esdTrack->HasPointOnITSLayer(2)));
1255       } else { // L34 or L56
1256         cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2+1), esdTrack->HasPointOnITSLayer(i*2+2));
1257       }
1258     }
1259   }
1260   
1261   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
1262     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
1263       // TPC tracks
1264       cuts[31] = kTRUE; 
1265     }else{
1266       // ITS standalone tracks
1267       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
1268         if(status & AliESDtrack::kITSpureSA) cuts[31] = kTRUE;
1269       }else if(fCutRequireITSpureSA){
1270         if(!(status & AliESDtrack::kITSpureSA)) cuts[31] = kTRUE;
1271       }
1272     }
1273   }
1274
1275   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
1276      cuts[32] = kTRUE;
1277
1278   if (!fCutAcceptSharedTPCClusters && nClustersTPCShared!=0)
1279     cuts[33] = kTRUE;
1280
1281   if (fracClustersTPCShared > fCutMaxFractionSharedTPCClusters)
1282     cuts[34] = kTRUE;  
1283
1284   Int_t nITSPointsForPid=0;
1285   UChar_t clumap=esdTrack->GetITSClusterMap();
1286   for(Int_t i=2; i<6; i++){
1287     if(clumap&(1<<i)) ++nITSPointsForPid;
1288   }
1289   if(fCutRequireITSPid && nITSPointsForPid<3) cuts[35] = kTRUE;
1290   
1291
1292   if (nCrossedRowsTPC<fCutMinNCrossedRowsTPC)
1293     cuts[36]=kTRUE;
1294   if (ratioCrossedRowsOverFindableClustersTPC<fCutMinRatioCrossedRowsOverFindableClustersTPC) 
1295     cuts[37]=kTRUE;
1296
1297   Int_t nMissITSpts=0;
1298   Int_t idet,statusLay;
1299   Float_t xloc,zloc;
1300   for(Int_t iLay=0; iLay<6; iLay++){
1301     Bool_t retc=esdTrack->GetITSModuleIndexInfo(iLay,idet,statusLay,xloc,zloc);
1302     if(retc && statusLay==5) ++nMissITSpts;
1303   }
1304   if(nMissITSpts>fCutMaxMissingITSPoints) cuts[38] = kTRUE;
1305   
1306   //kTOFout
1307   if (fCutRequireTOFout && (status&AliESDtrack::kTOFout)==0)
1308     cuts[40]=kTRUE;
1309
1310   // TOF signal Dz cut
1311   Float_t dxTOF = esdTrack->GetTOFsignalDx();
1312   Float_t dzTOF = esdTrack->GetTOFsignalDz();
1313   if (fFlagCutTOFdistance && (esdTrack->GetStatus() & AliESDtrack::kTOFout) == AliESDtrack::kTOFout){ // applying the TOF distance cut only if requested, and only on tracks that reached the TOF and where associated with a TOF hit 
1314           if (fgBeamTypeFlag < 0) {  // the check on the beam type was not done yet
1315                   const AliESDEvent* event = esdTrack->GetESDEvent();
1316                   if (event){
1317                           TString beamTypeESD = event->GetBeamType();
1318                           AliDebug(2,Form("Beam type from ESD event = %s",beamTypeESD.Data()));
1319                           if (beamTypeESD.CompareTo("A-A",TString::kIgnoreCase) == 0){ // we are in PbPb collisions --> fgBeamTypeFlag will be set to 1, to apply the cut on TOF signal Dz
1320                                   fgBeamTypeFlag = 1;
1321                           }
1322                           else { // we are NOT in PbPb collisions --> fgBeamTypeFlag will be set to 0, to NOT apply the cu6 on TOF signal Dz
1323                                   fgBeamTypeFlag = 0;
1324                           }                               
1325                   }
1326                   else{
1327                           AliFatal("Beam type not available, but it is needed to apply the TOF cut!");
1328                   }
1329           }
1330
1331           if (fgBeamTypeFlag == 1){ // we are in PbPb collisions --> apply the cut on TOF signal Dz
1332                   Float_t radiusTOF = TMath::Sqrt(dxTOF*dxTOF + dzTOF*dzTOF);
1333                   AliDebug(3,Form("TOF check (with fCutTOFdistance = %f) --> dx = %f, dz = %f, radius = %f", fCutTOFdistance, dxTOF, dzTOF, radiusTOF));
1334                   if (radiusTOF > fCutTOFdistance){
1335                           AliDebug(2, Form("************* the radius is outside the range! %f > %f, the track will be skipped", radiusTOF, fCutTOFdistance));
1336                           cuts[41] = kTRUE;
1337                   }
1338           }
1339   }
1340   
1341   Bool_t cut=kFALSE;
1342   for (Int_t i=0; i<kNCuts; i++) 
1343     if (cuts[i]) {cut = kTRUE;}
1344
1345   // for performance evaluate the CPU intensive cuts only when the others have passed, and when they are requested
1346   Double_t chi2TPCConstrainedVsGlobal = -2;
1347   Float_t nSigmaToVertex = -2;
1348   if (!cut)
1349   {
1350     // getting the track to vertex parameters
1351     if (fCutSigmaToVertexRequired)
1352     {
1353       nSigmaToVertex = GetSigmaToVertex(esdTrack);
1354       if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
1355       {
1356         cuts[12] = kTRUE;
1357         cut = kTRUE;
1358       }
1359       // if n sigma could not be calculated
1360       if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
1361       {
1362         cuts[13] = kTRUE;
1363         cut = kTRUE;
1364       }
1365     }
1366       
1367     // max chi2 TPC constrained vs global track only if track passed the other cut
1368     if (fCutMaxChi2TPCConstrainedVsGlobal < 1e9)
1369     {
1370       const AliESDEvent* esdEvent = esdTrack->GetESDEvent();
1371       
1372       if (!esdEvent)
1373         AliFatal("fCutMaxChi2TPCConstrainedVsGlobal set but ESD event not set in AliESDTrack. Use AliESDTrack::SetESDEvent before calling AliESDtrackCuts.");
1374       
1375       // get vertex
1376       const AliESDVertex* vertex = 0;
1377       if (fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTracks)
1378         vertex = esdEvent->GetPrimaryVertexTracks();
1379       
1380       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexSPD)
1381         vertex = esdEvent->GetPrimaryVertexSPD();
1382         
1383       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTPC)
1384         vertex = esdEvent->GetPrimaryVertexTPC();
1385
1386       if (vertex->GetStatus())
1387         chi2TPCConstrainedVsGlobal = esdTrack->GetChi2TPCConstrainedVsGlobal(vertex);
1388       
1389       if (chi2TPCConstrainedVsGlobal < 0 || chi2TPCConstrainedVsGlobal > fCutMaxChi2TPCConstrainedVsGlobal)
1390       {
1391         cuts[39] = kTRUE;
1392         cut = kTRUE;
1393       }
1394     }
1395
1396     // max length in active volume
1397     Float_t lengthInActiveZoneTPC = -1;
1398     if (fCutMinLengthActiveVolumeTPC > 1.) { // do the calculation only if needed to save cpu-time
1399       if (esdTrack->GetESDEvent()) {
1400         lengthInActiveZoneTPC = esdTrack->GetLengthInActiveZone(1, 1.8, 220, esdTrack->GetESDEvent()->GetMagneticField()); 
1401         //
1402         if (lengthInActiveZoneTPC < fCutMinLengthActiveVolumeTPC ) {
1403           cuts[42] = kTRUE;
1404           cut = kTRUE;
1405         }
1406       }
1407     }
1408
1409     
1410   }
1411
1412   //########################################################################
1413   // filling histograms
1414   if (fHistogramsOn) {
1415     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1416     if (cut)
1417       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1418     
1419     for (Int_t i=0; i<kNCuts; i++) {
1420       if (fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i]) < 1)
1421         AliFatal(Form("Inconsistency! Cut %d with name %s not found", i, fgkCutNames[i]));
1422     
1423       if (cuts[i])
1424         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1425
1426       for (Int_t j=i; j<kNCuts; j++) {
1427         if (cuts[i] && cuts[j]) {
1428           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1429           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1430           fhCutCorrelation->Fill(xC, yC);
1431         }
1432       }
1433     }
1434   }
1435   
1436   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1437   // the code is not in a function due to too many local variables that would need to be passed
1438
1439   for (Int_t id = 0; id < 2; id++)
1440   {
1441     // id = 0 --> before cut
1442     // id = 1 --> after cut
1443
1444     if (fHistogramsOn)
1445     {
1446       fhNClustersITS[id]->Fill(nClustersITS);
1447       fhNClustersTPC[id]->Fill(nClustersTPC);
1448       fhNSharedClustersTPC[id]->Fill(nClustersTPCShared);
1449       fhNCrossedRowsTPC[id]->Fill(nCrossedRowsTPC);
1450       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
1451       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1452       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1453       fhChi2TPCConstrainedVsGlobal[id]->Fill(chi2TPCConstrainedVsGlobal);
1454       fhNClustersForITSPID[id]->Fill(nITSPointsForPid);
1455       fhNMissingITSPoints[id]->Fill(nMissITSpts);
1456
1457       fhC11[id]->Fill(extCov[0]);
1458       fhC22[id]->Fill(extCov[2]);
1459       fhC33[id]->Fill(extCov[5]);
1460       fhC44[id]->Fill(extCov[9]);
1461       fhC55[id]->Fill(extCov[14]);
1462
1463       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1464
1465       fhPt[id]->Fill(pt);
1466       fhEta[id]->Fill(eta);
1467       fhTOFdistance[id]->Fill(dxTOF, dzTOF);
1468
1469       Float_t bRes[2];
1470       bRes[0] = TMath::Sqrt(bCov[0]);
1471       bRes[1] = TMath::Sqrt(bCov[2]);
1472
1473       fhDZ[id]->Fill(b[1]);
1474       fhDXY[id]->Fill(b[0]);
1475       fhDXYDZ[id]->Fill(dcaToVertex);
1476       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1477
1478       if (bRes[0]!=0 && bRes[1]!=0) {
1479         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1480         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1481         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1482         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1483       }
1484     }
1485
1486     // cut the track
1487     if (cut)
1488       return kFALSE;
1489   }
1490
1491   return kTRUE;
1492 }
1493
1494 //____________________________________________________________________
1495 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1496 {
1497   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1498   
1499   switch (req)
1500   {
1501         case kOff:        return kTRUE;
1502         case kNone:       return !clusterL1 && !clusterL2;
1503         case kAny:        return clusterL1 || clusterL2;
1504         case kFirst:      return clusterL1;
1505         case kOnlyFirst:  return clusterL1 && !clusterL2;
1506         case kSecond:     return clusterL2;
1507         case kOnlySecond: return clusterL2 && !clusterL1;
1508         case kBoth:       return clusterL1 && clusterL2;
1509   }
1510   
1511   return kFALSE;
1512 }
1513
1514 //____________________________________________________________________
1515 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(const AliESDEvent* esd, Int_t iTrack)
1516 {
1517   // Utility function to create a TPC only track from the given esd track
1518   // 
1519   // IMPORTANT: The track has to be deleted by the user
1520   //
1521   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1522   // there are only missing propagations here that are needed for old data
1523   // this function will therefore become obsolete
1524   //
1525   // adapted from code provided by CKB
1526
1527   if (!esd->GetPrimaryVertexTPC())
1528     return 0; // No TPC vertex no TPC tracks
1529
1530   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1531     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1532
1533   AliESDtrack* track = esd->GetTrack(iTrack);
1534   if (!track)
1535     return 0;
1536
1537   AliESDtrack *tpcTrack = new AliESDtrack();
1538
1539   // only true if we have a tpc track
1540   if (!track->FillTPCOnlyTrack(*tpcTrack))
1541   {
1542     delete tpcTrack;
1543     return 0;
1544   }
1545
1546   return tpcTrack;
1547 }
1548
1549 //____________________________________________________________________
1550 TObjArray* AliESDtrackCuts::GetAcceptedTracks(const AliESDEvent* esd, Bool_t bTPC)
1551 {
1552   //
1553   // returns an array of all tracks that pass the cuts
1554   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1555   // tracks that pass the cut
1556   //
1557   // NOTE: List has to be deleted by the user
1558
1559   TObjArray* acceptedTracks = new TObjArray();
1560
1561   // loop over esd tracks
1562   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1563     if(bTPC){
1564       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1565       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1566
1567       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1568       if (!tpcTrack)
1569         continue;
1570
1571       if (AcceptTrack(tpcTrack)) {
1572         acceptedTracks->Add(tpcTrack);
1573       }
1574       else
1575         delete tpcTrack;
1576     }
1577     else
1578     {
1579       AliESDtrack* track = esd->GetTrack(iTrack);
1580       if(AcceptTrack(track))
1581         acceptedTracks->Add(track);
1582     }
1583   } 
1584   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1585   return acceptedTracks;
1586 }
1587
1588 //____________________________________________________________________
1589 Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
1590 {
1591   //
1592   // returns an the number of tracks that pass the cuts
1593   //
1594
1595   Int_t count = 0;
1596
1597   // loop over esd tracks
1598   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1599     AliESDtrack* track = esd->GetTrack(iTrack);
1600     if (AcceptTrack(track))
1601       count++;
1602   }
1603
1604   return count;
1605 }
1606
1607 //____________________________________________________________________
1608  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1609    // 
1610    // diagnostics histograms are defined
1611    // 
1612
1613    fHistogramsOn=kTRUE;
1614    
1615    Bool_t oldStatus = TH1::AddDirectoryStatus();
1616    TH1::AddDirectory(kFALSE);
1617    
1618    //###################################################################################
1619    // defining histograms
1620
1621    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1622
1623    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1624    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1625
1626    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1627   
1628    for (Int_t i=0; i<kNCuts; i++) {
1629      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1630      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1631      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1632    }
1633
1634   fhCutStatistics  ->SetLineColor(color);
1635   fhCutCorrelation ->SetLineColor(color);
1636   fhCutStatistics  ->SetLineWidth(2);
1637   fhCutCorrelation ->SetLineWidth(2);
1638
1639   for (Int_t i=0; i<2; i++) {
1640     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1641     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1642     fhNSharedClustersTPC[i]  = new TH1F("nSharedClustersTPC"     ,"",165,-0.5,164.5);
1643     fhNCrossedRowsTPC[i]     = new TH1F("nCrossedRowsTPC"  ,"",165,-0.5,164.5);
1644     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
1645     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1646     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1647     fhChi2TPCConstrainedVsGlobal[i]   = new TH1F("chi2TPCConstrainedVsGlobal","",600,-2,50);
1648     fhNClustersForITSPID[i]  = new TH1F("nPointsForITSpid","",5,-0.5,4.5);
1649     fhNMissingITSPoints[i]   = new TH1F("nMissingITSClusters","",7,-0.5,6.5);
1650
1651     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1652     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1653     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1654     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1655     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1656
1657     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1658
1659     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1660     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1661     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1662     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1663
1664     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1665     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1666     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1667
1668     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1669
1670     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1671     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1672     fhTOFdistance[i]         = new TH2F("TOFdistance"     ,"TOF distance;dx (cm};dz (cm)", 150, -15, 15, 150, -15, 15);
1673     
1674     fhNClustersITS[i]->SetTitle("n ITS clusters");
1675     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1676     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
1677     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1678     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1679     fhChi2TPCConstrainedVsGlobal[i]->SetTitle("#Chi^{2} TPC constrained track vs global track");
1680     fhNClustersForITSPID[i]->SetTitle("n ITS points for PID");
1681     fhNMissingITSPoints[i]->SetTitle("n ITS layers with missing cluster");
1682
1683     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1684     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1685     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1686     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1687     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1688
1689     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1690
1691     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1692     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1693     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1694     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1695     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1696
1697     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1698     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1699     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1700     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1701     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1702
1703     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1704     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1705     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
1706     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1707     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1708     fhChi2TPCConstrainedVsGlobal[i]->SetLineColor(color);   fhChi2TPCConstrainedVsGlobal[i]->SetLineWidth(2);
1709     fhNClustersForITSPID[i]->SetLineColor(color);  fhNClustersForITSPID[i]->SetLineWidth(2);
1710     fhNMissingITSPoints[i]->SetLineColor(color);   fhNMissingITSPoints[i]->SetLineWidth(2);
1711
1712     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1713     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1714     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1715     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1716     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1717
1718     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1719
1720     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1721     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1722     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1723
1724     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1725     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1726     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2); 
1727   }
1728
1729   // The number of sigmas to the vertex is per definition gaussian
1730   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1731   ffDTheoretical->SetParameter(0,1);
1732
1733   TH1::AddDirectory(oldStatus);
1734 }
1735
1736 //____________________________________________________________________
1737 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1738 {
1739   //
1740   // loads the histograms from a file
1741   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1742   //
1743
1744   if (!dir)
1745     dir = GetName();
1746
1747   if (!gDirectory->cd(dir))
1748     return kFALSE;
1749
1750   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1751
1752   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1753   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1754
1755   for (Int_t i=0; i<2; i++) {
1756     if (i==0)
1757     {
1758       gDirectory->cd("before_cuts");
1759     }
1760     else
1761       gDirectory->cd("after_cuts");
1762
1763     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1764     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1765     fhNSharedClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSharedClustersTPC"     ));
1766     fhNCrossedRowsTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("nCrossedRowsTPC"  ));
1767     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
1768     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1769     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1770     fhChi2TPCConstrainedVsGlobal[i] = dynamic_cast<TH1F*> (gDirectory->Get("fhChi2TPCConstrainedVsGlobal"));
1771     fhNClustersForITSPID[i] = dynamic_cast<TH1F*> (gDirectory->Get("nPointsForITSpid"));
1772     fhNMissingITSPoints[i] = dynamic_cast<TH1F*> (gDirectory->Get("nMissingITSClusters"));
1773
1774     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1775     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1776     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1777     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1778     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1779
1780     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1781
1782     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1783     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1784     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1785     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1786
1787     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1788     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1789     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1790     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1791
1792     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1793     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1794     fhTOFdistance[i] = dynamic_cast<TH2F*> (gDirectory->Get("TOFdistance"));
1795
1796     gDirectory->cd("../");
1797   }
1798
1799   gDirectory->cd("..");
1800
1801   return kTRUE;
1802 }
1803
1804 //____________________________________________________________________
1805 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1806   //
1807   // saves the histograms in a directory (dir)
1808   //
1809
1810   if (!fHistogramsOn) {
1811     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1812     return;
1813   }
1814
1815   if (!dir)
1816     dir = GetName();
1817
1818   gDirectory->mkdir(dir);
1819   gDirectory->cd(dir);
1820
1821   gDirectory->mkdir("before_cuts");
1822   gDirectory->mkdir("after_cuts");
1823  
1824   // a factor of 2 is needed since n sigma is positive
1825   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1826   ffDTheoretical->Write("nSigmaToVertexTheory");
1827
1828   fhCutStatistics->Write();
1829   fhCutCorrelation->Write();
1830
1831   for (Int_t i=0; i<2; i++) {
1832     if (i==0)
1833       gDirectory->cd("before_cuts");
1834     else
1835       gDirectory->cd("after_cuts");
1836
1837     fhNClustersITS[i]        ->Write();
1838     fhNClustersTPC[i]        ->Write();
1839     fhNSharedClustersTPC[i]  ->Write();
1840     fhNCrossedRowsTPC[i]     ->Write();
1841     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
1842     fhChi2PerClusterITS[i]   ->Write();
1843     fhChi2PerClusterTPC[i]   ->Write();
1844     fhChi2TPCConstrainedVsGlobal[i]   ->Write();
1845     fhNClustersForITSPID[i]  ->Write();
1846     fhNMissingITSPoints[i]   ->Write();
1847
1848     fhC11[i]                 ->Write();
1849     fhC22[i]                 ->Write();
1850     fhC33[i]                 ->Write();
1851     fhC44[i]                 ->Write();
1852     fhC55[i]                 ->Write();
1853
1854     fhRel1PtUncertainty[i]   ->Write();
1855
1856     fhDXY[i]                 ->Write();
1857     fhDZ[i]                  ->Write();
1858     fhDXYDZ[i]               ->Write();
1859     fhDXYvsDZ[i]             ->Write();
1860
1861     fhDXYNormalized[i]       ->Write();
1862     fhDZNormalized[i]        ->Write();
1863     fhDXYvsDZNormalized[i]   ->Write();
1864     fhNSigmaToVertex[i]      ->Write();
1865
1866     fhPt[i]                  ->Write();
1867     fhEta[i]                 ->Write();
1868     fhTOFdistance[i]         ->Write();
1869     
1870     gDirectory->cd("../");
1871   }
1872
1873   gDirectory->cd("../");
1874 }
1875
1876 //____________________________________________________________________
1877 void AliESDtrackCuts::DrawHistograms()
1878 {
1879   // draws some histograms
1880
1881   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1882   canvas1->Divide(2, 2);
1883
1884   canvas1->cd(1);
1885   fhNClustersTPC[0]->SetStats(kFALSE);
1886   fhNClustersTPC[0]->Draw();
1887
1888   canvas1->cd(2);
1889   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1890   fhChi2PerClusterTPC[0]->Draw();
1891
1892   canvas1->cd(3);
1893   fhNSigmaToVertex[0]->SetStats(kFALSE);
1894   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1895   fhNSigmaToVertex[0]->Draw();
1896
1897   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1898
1899   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1900   canvas2->Divide(3, 2);
1901
1902   canvas2->cd(1);
1903   fhC11[0]->SetStats(kFALSE);
1904   gPad->SetLogy();
1905   fhC11[0]->Draw();
1906
1907   canvas2->cd(2);
1908   fhC22[0]->SetStats(kFALSE);
1909   gPad->SetLogy();
1910   fhC22[0]->Draw();
1911
1912   canvas2->cd(3);
1913   fhC33[0]->SetStats(kFALSE);
1914   gPad->SetLogy();
1915   fhC33[0]->Draw();
1916
1917   canvas2->cd(4);
1918   fhC44[0]->SetStats(kFALSE);
1919   gPad->SetLogy();
1920   fhC44[0]->Draw();
1921
1922   canvas2->cd(5);
1923   fhC55[0]->SetStats(kFALSE);
1924   gPad->SetLogy();
1925   fhC55[0]->Draw();
1926
1927   canvas2->cd(6);
1928   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1929   gPad->SetLogy();
1930   fhRel1PtUncertainty[0]->Draw();
1931
1932   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1933
1934   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1935   canvas3->Divide(3, 2);
1936
1937   canvas3->cd(1);
1938   fhDXY[0]->SetStats(kFALSE);
1939   gPad->SetLogy();
1940   fhDXY[0]->Draw();
1941
1942   canvas3->cd(2);
1943   fhDZ[0]->SetStats(kFALSE);
1944   gPad->SetLogy();
1945   fhDZ[0]->Draw();
1946
1947   canvas3->cd(3);
1948   fhDXYvsDZ[0]->SetStats(kFALSE);
1949   gPad->SetLogz();
1950   gPad->SetRightMargin(0.15);
1951   fhDXYvsDZ[0]->Draw("COLZ");
1952
1953   canvas3->cd(4);
1954   fhDXYNormalized[0]->SetStats(kFALSE);
1955   gPad->SetLogy();
1956   fhDXYNormalized[0]->Draw();
1957
1958   canvas3->cd(5);
1959   fhDZNormalized[0]->SetStats(kFALSE);
1960   gPad->SetLogy();
1961   fhDZNormalized[0]->Draw();
1962
1963   canvas3->cd(6);
1964   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1965   gPad->SetLogz();
1966   gPad->SetRightMargin(0.15);
1967   fhDXYvsDZNormalized[0]->Draw("COLZ");
1968
1969   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1970
1971   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1972   canvas4->Divide(2, 1);
1973
1974   canvas4->cd(1);
1975   fhCutStatistics->SetStats(kFALSE);
1976   fhCutStatistics->LabelsOption("v");
1977   gPad->SetBottomMargin(0.3);
1978   fhCutStatistics->Draw();
1979
1980   canvas4->cd(2);
1981   fhCutCorrelation->SetStats(kFALSE);
1982   fhCutCorrelation->LabelsOption("v");
1983   gPad->SetBottomMargin(0.3);
1984   gPad->SetLeftMargin(0.3);
1985   fhCutCorrelation->Draw("COLZ");
1986
1987   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1988
1989   /*canvas->cd(1);
1990   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1991   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1992
1993   canvas->cd(2);
1994   fhNClustersTPC[0]->SetStats(kFALSE);
1995   fhNClustersTPC[0]->DrawCopy();
1996
1997   canvas->cd(3);
1998   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1999   fhChi2PerClusterITS[0]->DrawCopy();
2000   fhChi2PerClusterITS[1]->SetLineColor(2);
2001   fhChi2PerClusterITS[1]->DrawCopy("SAME");
2002
2003   canvas->cd(4);
2004   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
2005   fhChi2PerClusterTPC[0]->DrawCopy();
2006   fhChi2PerClusterTPC[1]->SetLineColor(2);
2007   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
2008 }
2009 //--------------------------------------------------------------------------
2010 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
2011   //
2012   // set the pt-dependent DCA cuts
2013   //
2014
2015   if(f1CutMaxDCAToVertexXYPtDep) {
2016      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
2017   }
2018
2019   if(f1CutMaxDCAToVertexZPtDep) {
2020     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
2021   }
2022
2023   if(f1CutMinDCAToVertexXYPtDep) {
2024     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
2025   }
2026
2027   if(f1CutMinDCAToVertexZPtDep) {
2028     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
2029   }
2030
2031
2032   return;
2033 }
2034
2035
2036
2037 //--------------------------------------------------------------------------
2038 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
2039   //
2040   // Check the correctness of the string syntax
2041   //
2042   Bool_t retval=kTRUE;
2043
2044   if(!dist.Contains("pt")) {
2045     if(print) AliError("string must contain \"pt\"");
2046     retval= kFALSE;
2047   } 
2048   return retval;
2049 }
2050
2051  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
2052
2053    if(f1CutMaxDCAToVertexXYPtDep){
2054      delete f1CutMaxDCAToVertexXYPtDep;
2055      // resetiing both
2056      f1CutMaxDCAToVertexXYPtDep = 0;
2057      fCutMaxDCAToVertexXYPtDep = "";
2058    }
2059    if(!CheckPtDepDCA(dist,kTRUE)){
2060      return;
2061    }  
2062    fCutMaxDCAToVertexXYPtDep = dist;
2063    TString tmp(dist);
2064    tmp.ReplaceAll("pt","x");
2065    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
2066  
2067 }
2068
2069  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
2070
2071
2072    if(f1CutMaxDCAToVertexZPtDep){
2073      delete f1CutMaxDCAToVertexZPtDep;
2074      // resetiing both
2075      f1CutMaxDCAToVertexZPtDep = 0;
2076      fCutMaxDCAToVertexZPtDep = "";
2077    }
2078    if(!CheckPtDepDCA(dist,kTRUE))return;
2079      
2080    fCutMaxDCAToVertexZPtDep = dist;
2081    TString tmp(dist);
2082    tmp.ReplaceAll("pt","x");
2083    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
2084
2085    
2086 }
2087
2088
2089  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
2090
2091
2092    if(f1CutMinDCAToVertexXYPtDep){
2093      delete f1CutMinDCAToVertexXYPtDep;
2094      // resetiing both
2095      f1CutMinDCAToVertexXYPtDep = 0;
2096      fCutMinDCAToVertexXYPtDep = "";
2097    }
2098    if(!CheckPtDepDCA(dist,kTRUE))return;
2099
2100    fCutMinDCAToVertexXYPtDep = dist;
2101    TString tmp(dist);
2102    tmp.ReplaceAll("pt","x");
2103    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
2104
2105 }
2106
2107
2108  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
2109
2110    
2111
2112    if(f1CutMinDCAToVertexZPtDep){
2113      delete f1CutMinDCAToVertexZPtDep;
2114      // resetiing both
2115      f1CutMinDCAToVertexZPtDep = 0;
2116      fCutMinDCAToVertexZPtDep = "";
2117    }
2118    if(!CheckPtDepDCA(dist,kTRUE))return;
2119    fCutMinDCAToVertexZPtDep = dist;
2120    TString tmp(dist);
2121    tmp.ReplaceAll("pt","x");
2122    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
2123 }
2124
2125 AliESDtrackCuts* AliESDtrackCuts::GetMultEstTrackCuts(MultEstTrackCuts cut)
2126 {
2127   // returns the multiplicity estimator track cuts objects to allow for user configuration
2128   // upon first call the objects are created
2129   //
2130   // the cut defined here correspond to GetStandardITSTPCTrackCuts2010 (apart from the one for without SPD)
2131   
2132   if (!fgMultEstTrackCuts[kMultEstTrackCutGlobal])
2133   {
2134     // quality cut on ITS+TPC tracks
2135     fgMultEstTrackCuts[kMultEstTrackCutGlobal] = new AliESDtrackCuts();
2136     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMinNClustersTPC(70);
2137     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMaxChi2PerClusterTPC(4);
2138     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetAcceptKinkDaughters(kFALSE);
2139     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireTPCRefit(kTRUE);
2140     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireITSRefit(kTRUE);
2141     //multiplicity underestimate if we use global tracks with |eta| > 0.9
2142     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetEtaRange(-0.9, 0.9);
2143     
2144     // quality cut on ITS_SA tracks (complementary to ITS+TPC)
2145     fgMultEstTrackCuts[kMultEstTrackCutITSSA] = new AliESDtrackCuts();
2146     fgMultEstTrackCuts[kMultEstTrackCutITSSA]->SetRequireITSRefit(kTRUE);
2147     
2148     // primary selection for tracks with SPD hits
2149     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD] = new AliESDtrackCuts();
2150     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
2151     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
2152     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexZ(2);
2153     
2154     // primary selection for tracks w/o SPD hits
2155     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD] = new AliESDtrackCuts();
2156     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
2157     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
2158     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexZ(2);
2159   }
2160   
2161   return fgMultEstTrackCuts[cut];
2162 }
2163
2164 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, MultEstTrackType trackType, Float_t etaRange) 
2165 {
2166   // Get multiplicity estimate based on TPC/ITS tracks and tracklets
2167   // Adapted for AliESDtrackCuts from a version developed by: Ruben Shahoyan, Anton Alkin, Arvinder Palaha
2168   //
2169   // Returns a negative value if no reliable estimate can be provided:
2170   //   -1 SPD vertex missing
2171   //   -2 SPD VertexerZ dispersion too large
2172   //   -3 Track vertex missing (not checked for kTracklets)
2173   //   -4 Distance between SPD and track vertices too large (not checked for kTracklets)
2174   //
2175   // WARNING This functions does not cut on the z vtx. Depending on the eta range requested, you need to restrict your z vertex range!
2176   //
2177   // Strategy for combined estimators
2178   //   1. Count global tracks and flag them
2179   //   2. Count ITSSA as complementaries for ITSTPC+ or as main tracks
2180   //   3. Count the complementary SPD tracklets
2181   
2182   const AliESDVertex* vertices[2];
2183   vertices[0] = esd->GetPrimaryVertexSPD();
2184   vertices[1] = esd->GetPrimaryVertexTracks();
2185   
2186   if (!vertices[0]->GetStatus())
2187   {
2188     AliDebugClass(AliLog::kDebug, "No SPD vertex. Not able to make a reliable multiplicity estimate.");
2189     return -1;
2190   }
2191   
2192   if (vertices[0]->IsFromVertexerZ() && vertices[0]->GetDispersion() > 0.02)
2193   {
2194     AliDebugClass(AliLog::kDebug, "Vertexer z dispersion > 0.02. Not able to make a reliable multiplicity estimate.");
2195     return -2;
2196   }
2197   
2198   Int_t multiplicityEstimate = 0;
2199   
2200   // SPD tracklet-only estimate
2201   if (trackType == kTracklets)
2202   {
2203     const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2204     for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i) 
2205     {
2206       if (TMath::Abs(spdmult->GetEta(i)) > etaRange) 
2207         continue; // eta selection for tracklets
2208       multiplicityEstimate++;
2209     }
2210     return multiplicityEstimate;
2211   }
2212
2213   if (!vertices[1]->GetStatus())
2214   {
2215     AliDebugClass(AliLog::kDebug, "No track vertex. Not able to make a reliable multiplicity estimate.");
2216     return -3;
2217   }
2218
2219   // TODO value of displacement to be studied
2220   const Float_t maxDisplacement = 0.5;
2221   //check for displaced vertices
2222   Double_t displacement = TMath::Abs(vertices[0]->GetZ() - vertices[1]->GetZ());
2223   if (displacement > maxDisplacement) 
2224   {
2225     AliDebugClass(AliLog::kDebug, Form("Displaced vertices %f > %f",displacement,maxDisplacement));
2226     return -4;
2227   }
2228   
2229   // update eta range in track cuts
2230   GetMultEstTrackCuts(kMultEstTrackCutITSSA)->SetEtaRange(-etaRange, etaRange);
2231   GetMultEstTrackCuts(kMultEstTrackCutDCAwSPD)->SetEtaRange(-etaRange, etaRange);
2232   GetMultEstTrackCuts(kMultEstTrackCutDCAwoSPD)->SetEtaRange(-etaRange, etaRange);
2233   
2234   //*******************************************************************************************************
2235   //set counters to initial zeros
2236   Int_t tracksITSTPC = 0;     //number of global tracks for a given event
2237   Int_t tracksITSSA = 0;      //number of ITS standalone tracks for a given event
2238   Int_t tracksITSTPCSA_complementary = 0; //number of ITS standalone tracks complementary to TPC for a given event
2239   Int_t trackletsITSTPC_complementary = 0;//number of SPD tracklets complementary to global/ITSSA tracks for a given events
2240   Int_t trackletsITSSA_complementary = 0; //number of SPD tracklets complementary to ITSSA tracks for a given event
2241
2242   const Int_t nESDTracks = esd->GetNumberOfTracks();
2243
2244   // flags for secondary and rejected tracks
2245   const Int_t kRejBit = BIT(15); // set this bit in global tracks if it is rejected by a cut
2246   const Int_t kSecBit = BIT(16); // set this bit in global tracks if it is secondary according to a cut
2247   
2248   for(Int_t itracks=0; itracks < nESDTracks; itracks++) {
2249     esd->GetTrack(itracks)->ResetBit(kSecBit|kRejBit); //reset bits used for flagging secondaries and rejected tracks in case they were changed before this analysis
2250   }
2251   const Int_t maxid = nESDTracks+1; // used to define bool array for check multiple associations of tracklets to one track. array starts at 0.
2252   
2253   // bit mask for esd tracks, to check if multiple tracklets are associated to it
2254   TBits globalBits(maxid), pureITSBits(maxid);
2255   // why labels are used with the data? RS
2256   //  Bool_t globalBits[maxid], pureITSBits[maxid];
2257   //  for(Int_t i=0; i<maxid; i++){ // set all bools to false
2258   //      globalBits[i]=kFALSE;
2259   //      pureITSBits[i]=kFALSE;
2260   //  }
2261
2262   //*******************************************************************************************************
2263   // get multiplicity from global tracks
2264   for(Int_t itracks = 0; itracks < nESDTracks; itracks++) { // flag the tracks
2265     AliESDtrack* track = esd->GetTrack(itracks);
2266
2267     // if track is a secondary from a V0, flag as a secondary
2268     if (track->IsOn(AliESDtrack::kMultInV0)) {
2269       track->SetBit(kSecBit);
2270       continue;
2271     }
2272     /* done via proper DCA cut
2273     //secondary?
2274     if (track->IsOn(AliESDtrack::kMultSec)) {
2275       track->SetBit(kSecBit);
2276       continue;
2277     }
2278     */
2279     // check tracks with ITS part
2280     //*******************************************************************************************************
2281     if (track->IsOn(AliESDtrack::kITSin) && !track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSTPC) { // track has ITS part but is not an ITS_SA
2282       //*******************************************************************************************************
2283       // TPC+ITS
2284       if (track->IsOn(AliESDtrack::kTPCin)) {  // Global track, has ITS and TPC contributions
2285           if (fgMultEstTrackCuts[kMultEstTrackCutGlobal]->AcceptTrack(track)) { // good ITSTPC track
2286             if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2287               tracksITSTPC++; //global track counted
2288               globalBits.SetBitNumber(itracks);
2289             }
2290             else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2291           }
2292           else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2293       }
2294       //*******************************************************************************************************
2295       // ITS complementary
2296       else if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITS complementary track
2297         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2298           tracksITSTPCSA_complementary++;
2299           globalBits.SetBitNumber(itracks);
2300         }
2301         else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2302       }
2303       else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2304     }
2305     //*******************************************************************************************************
2306     // check tracks from ITS_SA_PURE
2307     if (track->IsOn(AliESDtrack::kITSin) && track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSSA){
2308       if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITSSA track
2309         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2310           tracksITSSA++;
2311           pureITSBits.SetBitNumber(itracks);
2312         }
2313         else track->SetBit(kRejBit);
2314       }
2315       else track->SetBit(kRejBit);
2316     }
2317   }//ESD tracks counted
2318
2319   //*******************************************************************************************************
2320   // get multiplicity from ITS tracklets to complement TPC+ITS, and ITSpureSA
2321   const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2322   for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i) {
2323     if (TMath::Abs(spdmult->GetEta(i)) > etaRange) continue; // eta selection for tracklets
2324     
2325     // if counting tracks+tracklets, check if clusters were already used in tracks
2326     Int_t id1,id2,id3,id4;
2327     spdmult->GetTrackletTrackIDs(i,0,id1,id2); // references for eventual Global/ITS_SA tracks
2328     AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
2329     spdmult->GetTrackletTrackIDs(i,1,id3,id4); // references for eventual ITS_SA_pure tracks
2330     AliESDtrack* tr3 = id3>=0 ? esd->GetTrack(id3) : 0;
2331
2332     // are both clusters from the same tracks? If not, skip the tracklet (shouldn't change things much)
2333     if ((id1!=id2 && id1>=0 && id2>=0) || (id3!=id4 && id3>=0 && id4>=0)) continue;
2334
2335     Bool_t bUsedInGlobal = (id1 != -1) ? globalBits.TestBitNumber(id1) : 0;// has associated global track been associated to a previous tracklet?
2336     Bool_t bUsedInPureITS = (id3 != -1) ? pureITSBits.TestBitNumber(id3) : 0;// has associated pure ITS track been associated to a previous tracklet?
2337     //*******************************************************************************************************
2338     if (trackType == kTrackletsITSTPC) {
2339       // count tracklets towards global+complementary tracks
2340       if (  (tr1 && !tr1->TestBit(kSecBit)) &&   // reject as secondary
2341         (tr1 &&  tr1->TestBit(kRejBit)) ) {  // count tracklet as bad quality track
2342           if(!bUsedInGlobal){
2343             ++trackletsITSTPC_complementary;
2344             if(id1>0) globalBits.SetBitNumber(id1); // mark global track linked to this tracklet as "associated"
2345           }
2346       }
2347       else if(id1<0) {
2348         ++trackletsITSTPC_complementary; // if no associated track, count the tracklet
2349       }
2350     } else {
2351       // count tracklets towards ITS_SA_pure tracks
2352       if (  (tr3 && !tr3->TestBit(kSecBit)) &&   // reject as secondary
2353         (tr3 &&  tr3->TestBit(kRejBit)) ) { // count tracklet as bad quality track
2354         if(!bUsedInPureITS) {
2355           ++trackletsITSSA_complementary;
2356           if(id3>0) pureITSBits.SetBitNumber(id3); // mark global track linked to this tracklet as "associated"
2357         }
2358       }
2359       else if(id3<0) {
2360         ++trackletsITSSA_complementary; // if no associated track, count the tracklet
2361       }
2362     }
2363   }
2364   
2365   //*******************************************************************************************************
2366   if (trackType == kTrackletsITSTPC)
2367     multiplicityEstimate = tracksITSTPC + tracksITSTPCSA_complementary + trackletsITSTPC_complementary;
2368   else
2369     multiplicityEstimate = tracksITSSA + trackletsITSSA_complementary;
2370
2371   return multiplicityEstimate;
2372 }
2373
2374 //____________________________________________________________________
2375 void AliESDtrackCuts::SetRequireStandardTOFmatchCuts(){
2376
2377         // setting the TOF cuts flags (kTOFout = TOF matching distance) to true, to include the selection on the standard TOF matching 
2378
2379         SetRequireTOFout(kTRUE);
2380         SetFlagCutTOFdistance(kTRUE);
2381         SetCutTOFdistance(3.);
2382
2383 }
2384