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