Several fixes in the histogram definition needed for compilation and coding conv.
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / ESDtrackQualityCuts.cxx
CommitLineData
75ec0f41 1#include "ESDtrackQualityCuts.h"
2
3#include <Riostream.h>
4
5//____________________________________________________________________
6ClassImp(ESDtrackQualityCuts);
7
8//____________________________________________________________________
9ESDtrackQualityCuts::ESDtrackQualityCuts() {
10
11 //##############################################################################
12 // setting default cuts
13 SetMinNClustersTPC();
14 SetMinNClustersITS();
15 SetMaxChi2PerClusterTPC();
16 SetMaxChi2PerClusterITS();
17 SetMaxCovDiagonalElements();
18 SetRequireTPCRefit();
19 SetRequireITSRefit();
20 SetAcceptKingDaughters();
21 SetMinNsigmaToVertex();
22 SetRequireSigmaToVertex();
23
24 SetHistogramsOn();
25
26 // set the cut names
27 fCutNames[0] = "require TPC refit";
28 fCutNames[1] = "require ITS refit";
29 fCutNames[2] = "n clusters TPC";
30 fCutNames[3] = "n clusters ITS";
31 fCutNames[4] = "#Chi^{2}/clusters TPC";
32 fCutNames[5] = "#Chi^{2}/clusters ITS";
33 fCutNames[6] = "cov 11";
34 fCutNames[7] = "cov 22";
35 fCutNames[8] = "cov 33";
36 fCutNames[9] = "cov 44";
37 fCutNames[10] = "cov 55";
38 fCutNames[11] = "trk-to-vtx";
39 fCutNames[12] = "trk-to-vtx failed";
40 fCutNames[13] = "kink daughters";
41}
42
43//____________________________________________________________________
44Bool_t
45ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) {
46
47 // relate the track to the new vertex
48 esdTrack->RelateToVertex(esdVtx, field, 999999);
49
50 return AcceptTrack(esdTrack);
51}
52
53//____________________________________________________________________
54Bool_t
55ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field) {
56
57 AliESDVertex* esdVtx = new AliESDVertex(vtx, vtx_res,"new vertex");
58 esdTrack->RelateToVertex(esdVtx, field, 1e99);
59 return AcceptTrack(esdTrack);
60}
61
62//____________________________________________________________________
63Bool_t
64ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack) {
65
66 UInt_t status = esdTrack->GetStatus();
67
68 // getting variables from ESD
69 Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
70 Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
71
72 Float_t chi2PerClusterITS = -1;
73 Float_t chi2PerClusterTPC = -1;
74 if (nClustersITS!=0)
75 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
76 if (nClustersTPC!=0)
77 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
78
79 Double_t extCov[15];
80 esdTrack->GetExternalCovariance(extCov);
81
82 Float_t b[2];
83 Float_t bRes[2];
84 Float_t bCov[3];
85 esdTrack->GetImpactParameters(b,bCov);
86 if (bCov[0]<=0 || bCov[2]<=0) {
87 AliDebug(1, "Estimated b resolution zero!");
88 bCov[0]=0; bCov[1]=0;
89 }
90 bRes[0] = TMath::Sqrt(bCov[0]);
91 bRes[1] = TMath::Sqrt(bCov[2]);
92
93 //
94 // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
95 //
96 // this is not correct - it will not give n sigma!!!
97 //
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
104 //########################################################################
105 // cut the track?
106
107 Bool_t cuts[fNCuts];
108 for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
109
110 if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
111 cuts[0]=kTRUE;
112 if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
113 cuts[1]=kTRUE;
114 if (nClustersTPC<fCut_MinNClusterTPC)
115 cuts[2]=kTRUE;
116 if (nClustersITS<fCut_MinNClusterITS)
117 cuts[3]=kTRUE;
118 if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC)
119 cuts[4]=kTRUE;
120 if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS)
121 cuts[5]=kTRUE;
122 if (extCov[0] > fCut_MaxC11)
123 cuts[6]=kTRUE;
124 if (extCov[2] > fCut_MaxC22)
125 cuts[7]=kTRUE;
126 if (extCov[5] > fCut_MaxC33)
127 cuts[8]=kTRUE;
128 if (extCov[9] > fCut_MaxC44)
129 cuts[9]=kTRUE;
130 if (extCov[14] > fCut_MaxC55)
131 cuts[10]=kTRUE;
132 if (nSigmaToVertex > fCut_NsigmaToVertex)
133 cuts[11] = kTRUE;
134 // if n sigma could not be calculated
135 if (nSigmaToVertex==-1 && fCut_SigmaToVertexRequired)
136 cuts[12]=kTRUE;
137 if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
138 cuts[13]=kTRUE;
139
140 Bool_t cut=kFALSE;
141 for (Int_t i=0; i<fNCuts; i++)
142 if (cuts[i]) cut = kTRUE;
143
144 //########################################################################
145 // filling histograms
146 if (fHistogramsOn) {
147 hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
148
149 if (cut)
150 hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
151
152 for (Int_t i=0; i<fNCuts; i++) {
153 if (cuts[i])
154 hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
155
156 for (Int_t j=i; j<fNCuts; j++) {
157 if (cuts[i] && cuts[j]) {
158 Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
159 Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
160 hCutCorrelation->Fill(x,y);
161 }
162 }
163 }
164
165
166 hNClustersITS[0]->Fill(nClustersITS);
167 hNClustersTPC[0]->Fill(nClustersTPC);
168 hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
169 hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
170
171 hC11[0]->Fill(extCov[0]);
172 hC22[0]->Fill(extCov[2]);
173 hC33[0]->Fill(extCov[5]);
174 hC44[0]->Fill(extCov[9]);
175 hC55[0]->Fill(extCov[14]);
176
177 hDZ[0]->Fill(b[1]);
178 hDXY[0]->Fill(b[0]);
179 hDXYvsDZ[0]->Fill(b[1],b[0]);
180
181 if (bRes[0]!=0 && bRes[1]!=0) {
182 hDZNormalized[0]->Fill(b[1]/bRes[1]);
183 hDXYNormalized[0]->Fill(b[0]/bRes[0]);
184 hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
185 }
186 }
187
188 //########################################################################
189 // cut the track!
190 if (cut) return kFALSE;
191
192 //########################################################################
193 // filling histograms after cut
194 if (fHistogramsOn) {
195 hNClustersITS[1]->Fill(nClustersITS);
196 hNClustersTPC[1]->Fill(nClustersTPC);
197 hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
198 hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
199
200 hC11[1]->Fill(extCov[0]);
201 hC22[1]->Fill(extCov[2]);
202 hC33[1]->Fill(extCov[5]);
203 hC44[1]->Fill(extCov[9]);
204 hC55[1]->Fill(extCov[14]);
205
206 hDZ[1]->Fill(b[1]);
207 hDXY[1]->Fill(b[0]);
208 hDXYvsDZ[1]->Fill(b[1],b[0]);
209
210 hDZNormalized[1]->Fill(b[1]/bRes[1]);
211 hDXYNormalized[1]->Fill(b[0]/bRes[0]);
212 hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
213 }
214
215 return kTRUE;
216}
217
218//____________________________________________________________________
219void
220ESDtrackQualityCuts::DefineHistograms(Int_t color) {
221
222 fHistogramsOn=kTRUE;
223
224 //###################################################################################
225 // defining histograms
226
227 hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
228
229 hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
230 hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
231
232 hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
233
234 for (Int_t i=0; i<fNCuts; i++) {
235 hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
236 hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
237 hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
238 }
239
240 hCutStatistics ->SetLineColor(color);
241 hCutCorrelation ->SetLineColor(color);
242 hCutStatistics ->SetLineWidth(2);
243 hCutCorrelation ->SetLineWidth(2);
244
245
246 hNClustersITS = new TH1F*[2];
247 hNClustersTPC = new TH1F*[2];
248 hChi2PerClusterITS = new TH1F*[2];
249 hChi2PerClusterTPC = new TH1F*[2];
250
251 hC11 = new TH1F*[2];
252 hC22 = new TH1F*[2];
253 hC33 = new TH1F*[2];
254 hC44 = new TH1F*[2];
255 hC55 = new TH1F*[2];
256
257 hDXY = new TH1F*[2];
258 hDZ = new TH1F*[2];
259 hDXYvsDZ = new TH2F*[2];
260
261 hDXYNormalized = new TH1F*[2];
262 hDZNormalized = new TH1F*[2];
263 hDXYvsDZNormalized = new TH2F*[2];
264
265
266
267 Char_t str[256];
268 for (Int_t i=0; i<2; i++) {
269 if (i==0) sprintf(str,"");
270 else sprintf(str,"_cut");
271
272 hNClustersITS[i] = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
273 hNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
274 hChi2PerClusterITS[i] = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
275 hChi2PerClusterTPC[i] = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
276
277 hC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
278 hC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
279 hC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
280 hC44[i] = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
281 hC55[i] = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
282
283 hDXY[i] = new TH1F(Form("dXY%s",str),"",500,-10,10);
284 hDZ[i] = new TH1F(Form("dZ%s",str),"",500,-10,10);
285 hDXYvsDZ[i] = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
286
287 hDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
288 hDZNormalized[i] = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
289 hDXYvsDZNormalized[i] = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
290
291
292 hNClustersITS[i] ->SetXTitle("n ITS clusters");
293 hNClustersTPC[i] ->SetXTitle("n TPC clusters");
294 hChi2PerClusterITS[i] ->SetXTitle("#Chi^{2} per ITS cluster");
295 hChi2PerClusterTPC[i] ->SetXTitle("#Chi^{2} per TPC cluster");
296
297 hC11[i] ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
298 hC22[i] ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
299 hC33[i] ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
300 hC44[i] ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
301 hC55[i] ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
302
303 hDXY[i] ->SetXTitle("transverse impact parameter");
304 hDZ[i] ->SetXTitle("longitudinal impact parameter");
305 hDXYvsDZ[i] ->SetXTitle("longitudinal impact parameter");
306 hDXYvsDZ[i] ->SetYTitle("transverse impact parameter");
307
308 hDXYNormalized[i] ->SetXTitle("normalized trans impact par");
309 hDZNormalized[i] ->SetXTitle("normalized long impact par");
310 hDXYvsDZNormalized[i] ->SetXTitle("normalized long impact par");
311 hDXYvsDZNormalized[i] ->SetYTitle("normalized trans impact par");
312
313 hNClustersITS[i] ->SetLineColor(color); hNClustersITS[i] ->SetLineWidth(2);
314 hNClustersTPC[i] ->SetLineColor(color); hNClustersTPC[i] ->SetLineWidth(2);
315 hChi2PerClusterITS[i] ->SetLineColor(color); hChi2PerClusterITS[i] ->SetLineWidth(2);
316 hChi2PerClusterTPC[i] ->SetLineColor(color); hChi2PerClusterTPC[i] ->SetLineWidth(2);
317
318 hC11[i] ->SetLineColor(color); hC11[i] ->SetLineWidth(2);
319 hC22[i] ->SetLineColor(color); hC22[i] ->SetLineWidth(2);
320 hC33[i] ->SetLineColor(color); hC33[i] ->SetLineWidth(2);
321 hC44[i] ->SetLineColor(color); hC44[i] ->SetLineWidth(2);
322 hC55[i] ->SetLineColor(color); hC55[i] ->SetLineWidth(2);
323
324 hDXY[i] ->SetLineColor(color); hDXY[i] ->SetLineWidth(2);
325 hDZ[i] ->SetLineColor(color); hDZ[i] ->SetLineWidth(2);
326
327 hDXYNormalized[i] ->SetLineColor(color); hDXYNormalized[i] ->SetLineWidth(2);
328 hDZNormalized[i] ->SetLineColor(color); hDZNormalized[i] ->SetLineWidth(2);
329
330 }
331}
332
333//____________________________________________________________________
334void
335ESDtrackQualityCuts::SaveHistograms(Char_t* dir) {
336
337 if (!fHistogramsOn) {
338 AliDebug(0, "Histograms not on - cannot save histograms!!!");
339 return;
340 }
341
342 gDirectory->mkdir(dir);
343 gDirectory->cd(dir);
344
345 gDirectory->mkdir("before_cuts");
346 gDirectory->mkdir("after_cuts");
347
348 hCutStatistics->Write();
349 hCutCorrelation->Write();
350
351 for (Int_t i=0; i<2; i++) {
352 if (i==0)
353 gDirectory->cd("before_cuts");
354 else
355 gDirectory->cd("after_cuts");
356
357 hNClustersITS[i] ->Write();
358 hNClustersTPC[i] ->Write();
359 hChi2PerClusterITS[i] ->Write();
360 hChi2PerClusterTPC[i] ->Write();
361
362 hC11[i] ->Write();
363 hC22[i] ->Write();
364 hC33[i] ->Write();
365 hC44[i] ->Write();
366 hC55[i] ->Write();
367
368 hDXY[i] ->Write();
369 hDZ[i] ->Write();
370 hDXYvsDZ[i] ->Write();
371
372 hDXYNormalized[i] ->Write();
373 hDZNormalized[i] ->Write();
374 hDXYvsDZNormalized[i] ->Write();
375
376 gDirectory->cd("../");
377 }
378
379 gDirectory->cd("../");
380}
381
382
383