Several fixes in the DefineHistograms method needed for compilation
[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) {
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) {
179     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
180     
181     if (cut)
182       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
183     
184     for (Int_t i=0; i<fNCuts; i++) {
185       if (cuts[i])
186         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
187       
188       for (Int_t j=i; j<fNCuts; j++) {
189         if (cuts[i] && cuts[j]) {
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);
193         }
194       }
195     }
196     
197
198     fhNClustersITS[0]->Fill(nClustersITS);        
199     fhNClustersTPC[0]->Fill(nClustersTPC);        
200     fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
201     fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);   
202     
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]);                                  
208     
209     fhDZ[0]->Fill(b[1]);     
210     fhDXY[0]->Fill(b[0]);    
211     fhDXYvsDZ[0]->Fill(b[1],b[0]);
212
213     if (bRes[0]!=0 && bRes[1]!=0) {
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]);
217     }
218   }
219
220   //########################################################################  
221   // cut the track!
222   if (cut) return kFALSE;
223
224   //########################################################################  
225   // filling histograms after cut
226   if (fHistogramsOn) {
227     fhNClustersITS[1]->Fill(nClustersITS);        
228     fhNClustersTPC[1]->Fill(nClustersTPC);        
229     fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
230     fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);   
231     
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]);                                  
237     
238     fhDZ[1]->Fill(b[1]);     
239     fhDXY[1]->Fill(b[0]);    
240     fhDXYvsDZ[1]->Fill(b[1],b[0]);
241
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]);
245   }
246   
247   return kTRUE;
248 }
249
250 //____________________________________________________________________
251 TObjArray*
252 AliESDtrackCuts::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 //____________________________________________________________________
268  void AliESDtrackCuts::DefineHistograms(Int_t color) {
269
270    fHistogramsOn=kTRUE;
271
272 //   //###################################################################################
273    // defining histograms
274
275    fhCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
276
277    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
278    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
279
280    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
281   
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];
298                        
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];
304   
305   TH1F *fhDXY                 = new TH1F[2];
306   TH1F *fhDZ                   = new TH1F[2];
307   TH2F *fhDXYvsDZ              = new TH2F[2];
308
309   TH1F *fhDXYNormalized       = new TH1F[2];
310   TH1F *fhDZNormalized        = new TH1F[2];
311   TH2F *fhDXYvsDZNormalized   = new TH2F[2];
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
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);
323     
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);
329     
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);
333
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);
337
338
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"); 
343     
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]"); 
349    
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);
364                                                                                               
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);
370                                                                                               
371     fhDXY[i].SetLineColor(color);   fhDXY[i].SetLineWidth(2);
372     fhDZ[i].SetLineColor(color);   fhDZ[i].SetLineWidth(2);
373                                                      
374     fhDXYNormalized[i].SetLineColor(color);   fhDXYNormalized[i].SetLineWidth(2);
375     fhDZNormalized[i].SetLineColor(color);   fhDZNormalized[i].SetLineWidth(2);
376
377   }
378 }
379
380 //____________________________________________________________________
381 void 
382 AliESDtrackCuts::Print(const Option_t*) const {
383
384   AliInfo("AliESDtrackCuts...");
385 }
386
387
388 //____________________________________________________________________
389 void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
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  
402   fhCutStatistics->Write();
403   fhCutCorrelation->Write();
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     
411     fhNClustersITS[i]        ->Write();
412     fhNClustersTPC[i]        ->Write();
413     fhChi2PerClusterITS[i]   ->Write();
414     fhChi2PerClusterTPC[i]   ->Write();
415     
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();
425     
426     fhDXYNormalized[i]       ->Write();
427     fhDZNormalized[i]        ->Write();
428     fhDXYvsDZNormalized[i]   ->Write();
429
430     gDirectory->cd("../");
431   }
432
433   gDirectory->cd("../");
434 }
435
436
437