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