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