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