Minor improvement
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx
CommitLineData
0b75bef2 1#include "AliESDtrackCuts.h"
2
3#include <Riostream.h>
4
5//____________________________________________________________________
6ClassImp(AliESDtrackCuts);
7
8//____________________________________________________________________
9AliESDtrackCuts::AliESDtrackCuts() {
10
11 //##############################################################################
12 // setting default cuts
13
14 SetMinNClustersTPC();
15 SetMinNClustersITS();
16 SetMaxChi2PerClusterTPC();
17 SetMaxChi2PerClusterITS();
18 SetMaxCovDiagonalElements();
19 SetRequireTPCRefit();
20 SetRequireITSRefit();
21 SetAcceptKingDaughters();
22 SetMinNsigmaToVertex();
23 SetRequireSigmaToVertex();
24 SetPRange();
25 SetPtRange();
26 SetPxRange();
27 SetPyRange();
28 SetPzRange();
29 SetEtaRange();
30 SetRapRange();
31
32 SetHistogramsOn();
33
34 // set the cut names
35 fCutNames[0] = "require TPC refit";
36 fCutNames[1] = "require ITS refit";
37 fCutNames[2] = "n clusters TPC";
38 fCutNames[3] = "n clusters ITS";
39 fCutNames[4] = "#Chi^{2}/clusters TPC";
40 fCutNames[5] = "#Chi^{2}/clusters ITS";
41 fCutNames[6] = "cov 11";
42 fCutNames[7] = "cov 22";
43 fCutNames[8] = "cov 33";
44 fCutNames[9] = "cov 44";
45 fCutNames[10] = "cov 55";
46 fCutNames[11] = "trk-to-vtx";
47 fCutNames[12] = "trk-to-vtx failed";
48 fCutNames[13] = "kink daughters";
49
50 fCutNames[14] = "p";
51 fCutNames[15] = "p_{T}";
52 fCutNames[16] = "p_{x}";
53 fCutNames[17] = "p_{y}";
54 fCutNames[18] = "p_{z}";
55 fCutNames[19] = "y";
56 fCutNames[20] = "eta";
57
58}
59
0b75bef2 60//____________________________________________________________________
61Bool_t
62AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
63 //
64 // figure out if the tracks survives all the track cuts defined
65 //
66
67 UInt_t status = esdTrack->GetStatus();
68
69 // getting quality parameters from the ESD track
70 Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
71 Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
72
73 Float_t chi2PerClusterITS = -1;
74 Float_t chi2PerClusterTPC = -1;
75 if (nClustersITS!=0)
76 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
77 if (nClustersTPC!=0)
78 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
79
80 Double_t extCov[15];
81 esdTrack->GetExternalCovariance(extCov);
82
83 // getting the track to vertex parameters
84 Float_t b[2];
85 Float_t bRes[2];
86 Float_t bCov[3];
87 esdTrack->GetImpactParameters(b,bCov);
88 if (bCov[0]<=0 || bCov[2]<=0) {
89 AliDebug(1, "Estimated b resolution zero!");
90 bCov[0]=0; bCov[1]=0;
91 }
92 bRes[0] = TMath::Sqrt(bCov[0]);
93 bRes[1] = TMath::Sqrt(bCov[2]);
94
95 // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
96 //
97 // this is not correct - it will not give n sigma!!!
98 //
99 Float_t nSigmaToVertex = -1;
100 if (bRes[0]!=0 && bRes[1]!=0)
101 nSigmaToVertex = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
102
103 // getting the kinematic variables of the track
104 // (assuming the mass is known)
105 Double_t p[3];
106 esdTrack->GetPxPyPz(p);
107 Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
108 Float_t pt = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
109 Float_t energy = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
110
111 //y-eta related calculations
112 Float_t eta = -100.;
113 Float_t y = -100.;
114 if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
115 eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
116 if((energy != TMath::Abs(p[2]))&&(momentum != 0))
117 y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
118
119
120 //########################################################################
121 // cut the track?
122
123 Bool_t cuts[fNCuts];
124 for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
125
126 // track quality cuts
127 if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
128 cuts[0]=kTRUE;
129 if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
130 cuts[1]=kTRUE;
131 if (nClustersTPC<fCut_MinNClusterTPC)
132 cuts[2]=kTRUE;
133 if (nClustersITS<fCut_MinNClusterITS)
134 cuts[3]=kTRUE;
135 if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC)
136 cuts[4]=kTRUE;
137 if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS)
138 cuts[5]=kTRUE;
139 if (extCov[0] > fCut_MaxC11)
140 cuts[6]=kTRUE;
141 if (extCov[2] > fCut_MaxC22)
142 cuts[7]=kTRUE;
143 if (extCov[5] > fCut_MaxC33)
144 cuts[8]=kTRUE;
145 if (extCov[9] > fCut_MaxC44)
146 cuts[9]=kTRUE;
147 if (extCov[14] > fCut_MaxC55)
148 cuts[10]=kTRUE;
149 if (nSigmaToVertex > fCut_NsigmaToVertex)
150 cuts[11] = kTRUE;
151 // if n sigma could not be calculated
152 if (nSigmaToVertex<0 && fCut_SigmaToVertexRequired)
153 cuts[12]=kTRUE;
154 if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
155 cuts[13]=kTRUE;
156 // track kinematics cut
157 if((momentum < fPMin) || (momentum > fPMax))
158 cuts[14]=kTRUE;
159 if((pt < fPtMin) || (pt > fPtMax))
160 cuts[15] = kTRUE;
161 if((p[0] < fPxMin) || (p[0] > fPxMax))
162 cuts[16] = kTRUE;
163 if((p[1] < fPyMin) || (p[1] > fPyMax))
164 cuts[17] = kTRUE;
165 if((p[2] < fPzMin) || (p[2] > fPzMax))
166 cuts[18] = kTRUE;
167 if((eta < fEtaMin) || (eta > fEtaMax))
168 cuts[19] = kTRUE;
169 if((y < fRapMin) || (y > fRapMax))
170 cuts[20] = kTRUE;
171
172 Bool_t cut=kFALSE;
173 for (Int_t i=0; i<fNCuts; i++)
174 if (cuts[i]) cut = kTRUE;
175
176 //########################################################################
177 // filling histograms
178 if (fHistogramsOn) {
179 hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
180
181 if (cut)
182 hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
183
184 for (Int_t i=0; i<fNCuts; i++) {
185 if (cuts[i])
186 hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
187
188 for (Int_t j=i; j<fNCuts; j++) {
189 if (cuts[i] && cuts[j]) {
190 Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
191 Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
192 hCutCorrelation->Fill(x,y);
193 }
194 }
195 }
196
197
198 hNClustersITS[0]->Fill(nClustersITS);
199 hNClustersTPC[0]->Fill(nClustersTPC);
200 hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
201 hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
202
203 hC11[0]->Fill(extCov[0]);
204 hC22[0]->Fill(extCov[2]);
205 hC33[0]->Fill(extCov[5]);
206 hC44[0]->Fill(extCov[9]);
207 hC55[0]->Fill(extCov[14]);
208
209 hDZ[0]->Fill(b[1]);
210 hDXY[0]->Fill(b[0]);
211 hDXYvsDZ[0]->Fill(b[1],b[0]);
212
213 if (bRes[0]!=0 && bRes[1]!=0) {
214 hDZNormalized[0]->Fill(b[1]/bRes[1]);
215 hDXYNormalized[0]->Fill(b[0]/bRes[0]);
216 hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
217 }
218 }
219
220 //########################################################################
221 // cut the track!
222 if (cut) return kFALSE;
223
224 //########################################################################
225 // filling histograms after cut
226 if (fHistogramsOn) {
227 hNClustersITS[1]->Fill(nClustersITS);
228 hNClustersTPC[1]->Fill(nClustersTPC);
229 hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
230 hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
231
232 hC11[1]->Fill(extCov[0]);
233 hC22[1]->Fill(extCov[2]);
234 hC33[1]->Fill(extCov[5]);
235 hC44[1]->Fill(extCov[9]);
236 hC55[1]->Fill(extCov[14]);
237
238 hDZ[1]->Fill(b[1]);
239 hDXY[1]->Fill(b[0]);
240 hDXYvsDZ[1]->Fill(b[1],b[0]);
241
242 hDZNormalized[1]->Fill(b[1]/bRes[1]);
243 hDXYNormalized[1]->Fill(b[0]/bRes[0]);
244 hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
245 }
246
247 return kTRUE;
248}
249
250//____________________________________________________________________
251TObjArray*
252AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) {
253
254 // returns an array of all tracks that pass the cuts
255 fAcceptedTracks->Clear();
256
257 // loop over esd tracks
258 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
259 AliESDtrack* track = esd->GetTrack(iTrack);
260
261 if(AcceptTrack(track)) fAcceptedTracks->Add(track);
262 }
263
264 return fAcceptedTracks;
265}
266
267//____________________________________________________________________
268void
269AliESDtrackCuts::DefineHistograms(Int_t color) {
270
271 fHistogramsOn=kTRUE;
272
273 //###################################################################################
274 // defining histograms
275
276 hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
277
278 hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
279 hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
280
281 hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
282
283 for (Int_t i=0; i<fNCuts; i++) {
284 hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
285 hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
286 hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
287 }
288
289 hCutStatistics ->SetLineColor(color);
290 hCutCorrelation ->SetLineColor(color);
291 hCutStatistics ->SetLineWidth(2);
292 hCutCorrelation ->SetLineWidth(2);
293
294
295 hNClustersITS = new TH1F*[2];
296 hNClustersTPC = new TH1F*[2];
297 hChi2PerClusterITS = new TH1F*[2];
298 hChi2PerClusterTPC = new TH1F*[2];
299
300 hC11 = new TH1F*[2];
301 hC22 = new TH1F*[2];
302 hC33 = new TH1F*[2];
303 hC44 = new TH1F*[2];
304 hC55 = new TH1F*[2];
305
306 hDXY = new TH1F*[2];
307 hDZ = new TH1F*[2];
308 hDXYvsDZ = new TH2F*[2];
309
310 hDXYNormalized = new TH1F*[2];
311 hDZNormalized = new TH1F*[2];
312 hDXYvsDZNormalized = new TH2F*[2];
313
314
315 Char_t str[256];
316 for (Int_t i=0; i<2; i++) {
317 if (i==0) sprintf(str," ");
318 else sprintf(str,"_cut");
319
320 hNClustersITS[i] = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
321 hNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
322 hChi2PerClusterITS[i] = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
323 hChi2PerClusterTPC[i] = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
324
325 hC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
326 hC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
327 hC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
328 hC44[i] = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
329 hC55[i] = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
330
331 hDXY[i] = new TH1F(Form("dXY%s",str),"",500,-10,10);
332 hDZ[i] = new TH1F(Form("dZ%s",str),"",500,-10,10);
333 hDXYvsDZ[i] = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
334
335 hDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
336 hDZNormalized[i] = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
337 hDXYvsDZNormalized[i] = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
338
339
340 hNClustersITS[i] ->SetXTitle("n ITS clusters");
341 hNClustersTPC[i] ->SetXTitle("n TPC clusters");
342 hChi2PerClusterITS[i] ->SetXTitle("#Chi^{2} per ITS cluster");
343 hChi2PerClusterTPC[i] ->SetXTitle("#Chi^{2} per TPC cluster");
344
345 hC11[i] ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
346 hC22[i] ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
347 hC33[i] ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
348 hC44[i] ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
349 hC55[i] ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
350
351 hDXY[i] ->SetXTitle("transverse impact parameter");
352 hDZ[i] ->SetXTitle("longitudinal impact parameter");
353 hDXYvsDZ[i] ->SetXTitle("longitudinal impact parameter");
354 hDXYvsDZ[i] ->SetYTitle("transverse impact parameter");
355
356 hDXYNormalized[i] ->SetXTitle("normalized trans impact par");
357 hDZNormalized[i] ->SetXTitle("normalized long impact par");
358 hDXYvsDZNormalized[i] ->SetXTitle("normalized long impact par");
359 hDXYvsDZNormalized[i] ->SetYTitle("normalized trans impact par");
360
361 hNClustersITS[i] ->SetLineColor(color); hNClustersITS[i] ->SetLineWidth(2);
362 hNClustersTPC[i] ->SetLineColor(color); hNClustersTPC[i] ->SetLineWidth(2);
363 hChi2PerClusterITS[i] ->SetLineColor(color); hChi2PerClusterITS[i] ->SetLineWidth(2);
364 hChi2PerClusterTPC[i] ->SetLineColor(color); hChi2PerClusterTPC[i] ->SetLineWidth(2);
365
366 hC11[i] ->SetLineColor(color); hC11[i] ->SetLineWidth(2);
367 hC22[i] ->SetLineColor(color); hC22[i] ->SetLineWidth(2);
368 hC33[i] ->SetLineColor(color); hC33[i] ->SetLineWidth(2);
369 hC44[i] ->SetLineColor(color); hC44[i] ->SetLineWidth(2);
370 hC55[i] ->SetLineColor(color); hC55[i] ->SetLineWidth(2);
371
372 hDXY[i] ->SetLineColor(color); hDXY[i] ->SetLineWidth(2);
373 hDZ[i] ->SetLineColor(color); hDZ[i] ->SetLineWidth(2);
374
375 hDXYNormalized[i] ->SetLineColor(color); hDXYNormalized[i] ->SetLineWidth(2);
376 hDZNormalized[i] ->SetLineColor(color); hDZNormalized[i] ->SetLineWidth(2);
377
378 }
379}
380
381//____________________________________________________________________
382void
49dc84d9 383AliESDtrackCuts::Print(const Option_t*) const {
0b75bef2 384
385 AliInfo("AliESDtrackCuts...");
386}
387
388
389//____________________________________________________________________
390void
391AliESDtrackCuts::SaveHistograms(Char_t* dir) {
392
393 if (!fHistogramsOn) {
394 AliDebug(0, "Histograms not on - cannot save histograms!!!");
395 return;
396 }
397
398 gDirectory->mkdir(dir);
399 gDirectory->cd(dir);
400
401 gDirectory->mkdir("before_cuts");
402 gDirectory->mkdir("after_cuts");
403
404 hCutStatistics->Write();
405 hCutCorrelation->Write();
406
407 for (Int_t i=0; i<2; i++) {
408 if (i==0)
409 gDirectory->cd("before_cuts");
410 else
411 gDirectory->cd("after_cuts");
412
413 hNClustersITS[i] ->Write();
414 hNClustersTPC[i] ->Write();
415 hChi2PerClusterITS[i] ->Write();
416 hChi2PerClusterTPC[i] ->Write();
417
418 hC11[i] ->Write();
419 hC22[i] ->Write();
420 hC33[i] ->Write();
421 hC44[i] ->Write();
422 hC55[i] ->Write();
423
424 hDXY[i] ->Write();
425 hDZ[i] ->Write();
426 hDXYvsDZ[i] ->Write();
427
428 hDXYNormalized[i] ->Write();
429 hDZNormalized[i] ->Write();
430 hDXYvsDZNormalized[i] ->Write();
431
432 gDirectory->cd("../");
433 }
434
435 gDirectory->cd("../");
436}
437
438
439