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