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