Several fixes in the histogram definition needed for compilation and coding conv.
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / ESDtrackQualityCuts.cxx
1 #include "ESDtrackQualityCuts.h"
2
3 #include <Riostream.h>
4
5 //____________________________________________________________________
6 ClassImp(ESDtrackQualityCuts);
7
8 //____________________________________________________________________
9 ESDtrackQualityCuts::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 //____________________________________________________________________
44 Bool_t 
45 ESDtrackQualityCuts::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 //____________________________________________________________________
54 Bool_t 
55 ESDtrackQualityCuts::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 //____________________________________________________________________
63 Bool_t 
64 ESDtrackQualityCuts::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 //____________________________________________________________________
219 void 
220 ESDtrackQualityCuts::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 //____________________________________________________________________
334 void 
335 ESDtrackQualityCuts::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