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