Several fixes in the histogram definition needed for compilation and coding conv.
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx
CommitLineData
0b75bef2 1#include "AliESDtrackCuts.h"
2
3#include <Riostream.h>
4
5//____________________________________________________________________
f58f1a93 6ClassImp(AliESDtrackCuts)
0b75bef2 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
0b75bef2 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) {
f58f1a93 179 fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
0b75bef2 180
181 if (cut)
f58f1a93 182 fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
0b75bef2 183
184 for (Int_t i=0; i<fNCuts; i++) {
185 if (cuts[i])
f58f1a93 186 fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
0b75bef2 187
188 for (Int_t j=i; j<fNCuts; j++) {
189 if (cuts[i] && cuts[j]) {
f58f1a93 190 Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
191 Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
192 fhCutCorrelation->Fill(x,y);
0b75bef2 193 }
194 }
195 }
196
197
f58f1a93 198 fhNClustersITS[0]->Fill(nClustersITS);
199 fhNClustersTPC[0]->Fill(nClustersTPC);
200 fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
201 fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
0b75bef2 202
f58f1a93 203 fhC11[0]->Fill(extCov[0]);
204 fhC22[0]->Fill(extCov[2]);
205 fhC33[0]->Fill(extCov[5]);
206 fhC44[0]->Fill(extCov[9]);
207 fhC55[0]->Fill(extCov[14]);
0b75bef2 208
f58f1a93 209 fhDZ[0]->Fill(b[1]);
210 fhDXY[0]->Fill(b[0]);
211 fhDXYvsDZ[0]->Fill(b[1],b[0]);
0b75bef2 212
213 if (bRes[0]!=0 && bRes[1]!=0) {
f58f1a93 214 fhDZNormalized[0]->Fill(b[1]/bRes[1]);
215 fhDXYNormalized[0]->Fill(b[0]/bRes[0]);
216 fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
0b75bef2 217 }
218 }
219
220 //########################################################################
221 // cut the track!
222 if (cut) return kFALSE;
223
224 //########################################################################
225 // filling histograms after cut
226 if (fHistogramsOn) {
f58f1a93 227 fhNClustersITS[1]->Fill(nClustersITS);
228 fhNClustersTPC[1]->Fill(nClustersTPC);
229 fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
230 fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
0b75bef2 231
f58f1a93 232 fhC11[1]->Fill(extCov[0]);
233 fhC22[1]->Fill(extCov[2]);
234 fhC33[1]->Fill(extCov[5]);
235 fhC44[1]->Fill(extCov[9]);
236 fhC55[1]->Fill(extCov[14]);
0b75bef2 237
f58f1a93 238 fhDZ[1]->Fill(b[1]);
239 fhDXY[1]->Fill(b[0]);
240 fhDXYvsDZ[1]->Fill(b[1],b[0]);
0b75bef2 241
f58f1a93 242 fhDZNormalized[1]->Fill(b[1]/bRes[1]);
243 fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
244 fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
0b75bef2 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//____________________________________________________________________
f58f1a93 268 void AliESDtrackCuts::DefineHistograms(Int_t color) {
0b75bef2 269
f58f1a93 270 fHistogramsOn=kTRUE;
0b75bef2 271
f58f1a93 272// //###################################################################################
273 // defining histograms
0b75bef2 274
f58f1a93 275 fhCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
0b75bef2 276
f58f1a93 277 fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
278 fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
0b75bef2 279
f58f1a93 280 fhCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
0b75bef2 281
f58f1a93 282 for (Int_t i=0; i<fNCuts; i++) {
283 fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
284 fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
285 fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
286 }
287
288 fhCutStatistics ->SetLineColor(color);
289 fhCutCorrelation ->SetLineColor(color);
290 fhCutStatistics ->SetLineWidth(2);
291 fhCutCorrelation ->SetLineWidth(2);
292
293
294 TH1F *fhNClustersITS = new TH1F[2];
295 TH1F *fhNClustersTPC = new TH1F[2];
296 TH1F *fhChi2PerClusterITS = new TH1F[2];
297 TH1F *fhChi2PerClusterTPC = new TH1F[2];
0b75bef2 298
f58f1a93 299 TH1F *fhC11 = new TH1F[2];
300 TH1F *fhC22 = new TH1F[2];
301 TH1F *fhC33 = new TH1F[2];
302 TH1F *fhC44 = new TH1F[2];
303 TH1F *fhC55 = new TH1F[2];
0b75bef2 304
f58f1a93 305 TH1F *fhDXY = new TH1F[2];
306 TH1F *fhDZ = new TH1F[2];
307 TH2F *fhDXYvsDZ = new TH2F[2];
0b75bef2 308
f58f1a93 309 TH1F *fhDXYNormalized = new TH1F[2];
310 TH1F *fhDZNormalized = new TH1F[2];
311 TH2F *fhDXYvsDZNormalized = new TH2F[2];
0b75bef2 312
313
314 Char_t str[256];
315 for (Int_t i=0; i<2; i++) {
316 if (i==0) sprintf(str," ");
317 else sprintf(str,"_cut");
318
f58f1a93 319 fhNClustersITS[i] = TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
320 fhNClustersTPC[i] = TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
321 fhChi2PerClusterITS[i] = TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
322 fhChi2PerClusterTPC[i] = TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
0b75bef2 323
f58f1a93 324 fhC11[i] = TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
325 fhC22[i] = TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
326 fhC33[i] = TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
327 fhC44[i] = TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
328 fhC55[i] = TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
0b75bef2 329
f58f1a93 330 fhDXY[i] = TH1F(Form("dXY%s",str),"",500,-10,10);
331 fhDZ[i] = TH1F(Form("dZ%s",str),"",500,-10,10);
332 fhDXYvsDZ[i] = TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
0b75bef2 333
f58f1a93 334 fhDXYNormalized[i] = TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
335 fhDZNormalized[i] = TH1F(Form("dZNormalized%s",str),"",500,-10,10);
336 fhDXYvsDZNormalized[i] = TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
0b75bef2 337
338
f58f1a93 339 fhNClustersITS[i].SetXTitle("n ITS clusters");
340 fhNClustersTPC[i].SetXTitle("n TPC clusters");
341 fhChi2PerClusterITS[i].SetXTitle("#Chi^{2} per ITS cluster");
342 fhChi2PerClusterTPC[i].SetXTitle("#Chi^{2} per TPC cluster");
0b75bef2 343
f58f1a93 344 fhC11[i].SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
345 fhC22[i].SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
346 fhC33[i].SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
347 fhC44[i].SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
348 fhC55[i].SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
0b75bef2 349
f58f1a93 350 fhDXY[i].SetXTitle("transverse impact parameter");
351 fhDZ[i].SetXTitle("longitudinal impact parameter");
352 fhDXYvsDZ[i].SetXTitle("longitudinal impact parameter");
353 fhDXYvsDZ[i].SetYTitle("transverse impact parameter");
354
355 fhDXYNormalized[i].SetXTitle("normalized trans impact par");
356 fhDZNormalized[i].SetXTitle("normalized long impact par");
357 fhDXYvsDZNormalized[i].SetXTitle("normalized long impact par");
358 fhDXYvsDZNormalized[i].SetYTitle("normalized trans impact par");
359
360 fhNClustersITS[i].SetLineColor(color); fhNClustersITS[i].SetLineWidth(2);
361 fhNClustersTPC[i].SetLineColor(color); fhNClustersTPC[i].SetLineWidth(2);
362 fhChi2PerClusterITS[i].SetLineColor(color); fhChi2PerClusterITS[i].SetLineWidth(2);
363 fhChi2PerClusterTPC[i].SetLineColor(color); fhChi2PerClusterTPC[i].SetLineWidth(2);
0b75bef2 364
f58f1a93 365 fhC11[i].SetLineColor(color); fhC11[i].SetLineWidth(2);
366 fhC22[i].SetLineColor(color); fhC22[i].SetLineWidth(2);
367 fhC33[i].SetLineColor(color); fhC33[i].SetLineWidth(2);
368 fhC44[i].SetLineColor(color); fhC44[i].SetLineWidth(2);
369 fhC55[i].SetLineColor(color); fhC55[i].SetLineWidth(2);
0b75bef2 370
f58f1a93 371 fhDXY[i].SetLineColor(color); fhDXY[i].SetLineWidth(2);
372 fhDZ[i].SetLineColor(color); fhDZ[i].SetLineWidth(2);
0b75bef2 373
f58f1a93 374 fhDXYNormalized[i].SetLineColor(color); fhDXYNormalized[i].SetLineWidth(2);
375 fhDZNormalized[i].SetLineColor(color); fhDZNormalized[i].SetLineWidth(2);
0b75bef2 376
377 }
378}
379
380//____________________________________________________________________
381void
49dc84d9 382AliESDtrackCuts::Print(const Option_t*) const {
0b75bef2 383
384 AliInfo("AliESDtrackCuts...");
385}
386
387
388//____________________________________________________________________
f58f1a93 389void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
0b75bef2 390
391 if (!fHistogramsOn) {
392 AliDebug(0, "Histograms not on - cannot save histograms!!!");
393 return;
394 }
395
396 gDirectory->mkdir(dir);
397 gDirectory->cd(dir);
398
399 gDirectory->mkdir("before_cuts");
400 gDirectory->mkdir("after_cuts");
401
f58f1a93 402 fhCutStatistics->Write();
403 fhCutCorrelation->Write();
0b75bef2 404
405 for (Int_t i=0; i<2; i++) {
406 if (i==0)
407 gDirectory->cd("before_cuts");
408 else
409 gDirectory->cd("after_cuts");
410
f58f1a93 411 fhNClustersITS[i] ->Write();
412 fhNClustersTPC[i] ->Write();
413 fhChi2PerClusterITS[i] ->Write();
414 fhChi2PerClusterTPC[i] ->Write();
0b75bef2 415
f58f1a93 416 fhC11[i] ->Write();
417 fhC22[i] ->Write();
418 fhC33[i] ->Write();
419 fhC44[i] ->Write();
420 fhC55[i] ->Write();
421
422 fhDXY[i] ->Write();
423 fhDZ[i] ->Write();
424 fhDXYvsDZ[i] ->Write();
0b75bef2 425
f58f1a93 426 fhDXYNormalized[i] ->Write();
427 fhDZNormalized[i] ->Write();
428 fhDXYvsDZNormalized[i] ->Write();
0b75bef2 429
430 gDirectory->cd("../");
431 }
432
433 gDirectory->cd("../");
434}
435
436
437