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