adding AliESDtrackCuts.
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx
1 #include "AliESDtrackCuts.h"
2
3 #include <Riostream.h>
4
5 //____________________________________________________________________
6 ClassImp(AliESDtrackCuts);
7
8 //____________________________________________________________________
9 AliESDtrackCuts::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 //____________________________________________________________________
61 Bool_t 
62 AliESDtrackCuts::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 //____________________________________________________________________
71 Bool_t 
72 AliESDtrackCuts::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 //____________________________________________________________________
80 Bool_t 
81 AliESDtrackCuts::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 //____________________________________________________________________
270 TObjArray*
271 AliESDtrackCuts::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 //____________________________________________________________________
287 void 
288 AliESDtrackCuts::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 //____________________________________________________________________
401 void 
402 AliESDtrackCuts::Print() {
403
404   AliInfo("AliESDtrackCuts...");
405 }
406
407
408 //____________________________________________________________________
409 void 
410 AliESDtrackCuts::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