Getter added (M. Fasel)
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 #include <TCanvas.h>\r
16 #include <TClass.h>\r
17 #include <TCollection.h>\r
18 #include <TDirectory.h>\r
19 #include <TH1F.h>\r
20 #include <TH1I.h>\r
21 #include <TH2I.h>\r
22 #include <TMath.h>\r
23 #include <TIterator.h>\r
24 #include <TString.h>\r
25 #include <TList.h>\r
26 \r
27 #include "AliESDtrack.h"\r
28 #include "AliESDEvent.h"\r
29 #include "AliLog.h"\r
30 #include "AliESDpid.h"\r
31 \r
32 #include "AliESDpidCuts.h"\r
33 \r
34 ClassImp(AliESDpidCuts)\r
35 \r
36 const Int_t AliESDpidCuts::kNcuts = 3;\r
37 \r
38 //_____________________________________________________________________\r
39 AliESDpidCuts::AliESDpidCuts(const Char_t *name, const Char_t *title):\r
40     AliAnalysisCuts(name, title)\r
41   , fESDpid(NULL)\r
42   , fTPCsigmaCutRequired(0)\r
43   , fTOFsigmaCutRequired(0)\r
44   , fCutTPCclusterRatio(0.)\r
45   , fMinMomentumTOF(0.5)\r
46   , fHcutStatistics(NULL)\r
47   , fHcutCorrelation(NULL)\r
48 {\r
49   //\r
50   // Default constructor\r
51   //\r
52   \r
53   fESDpid = new AliESDpid;\r
54   memset(fCutTPCnSigma, 0, sizeof(Float_t) * 2);\r
55   memset(fCutTOFnSigma, 0, sizeof(Float_t) * 2);\r
56 \r
57   memset(fHclusterRatio, 0, sizeof(TH1F *) * 2);\r
58   memset(fHnSigmaTPC, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);\r
59   memset(fHnSigmaTOF, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);\r
60 }\r
61 \r
62 //_____________________________________________________________________\r
63 AliESDpidCuts::AliESDpidCuts(const AliESDpidCuts &ref):\r
64     AliAnalysisCuts(ref)\r
65   , fESDpid(NULL)\r
66   , fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)\r
67   , fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)\r
68   , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)\r
69   , fMinMomentumTOF(ref.fMinMomentumTOF)\r
70   , fHcutStatistics(NULL)\r
71   , fHcutCorrelation(NULL)\r
72 {\r
73   //\r
74   // Copy constructor\r
75   //\r
76   fESDpid = new AliESDpid(*ref.fESDpid);\r
77   memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
78   memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
79   \r
80   if(ref.fHcutStatistics) fHcutStatistics = dynamic_cast<TH1I *>(ref.fHcutStatistics->Clone());\r
81   if(ref.fHcutCorrelation) fHcutCorrelation = dynamic_cast<TH2I *>(ref.fHcutCorrelation->Clone());\r
82   for(Int_t imode = 0; imode < 2; imode++){\r
83     if(ref.fHclusterRatio[imode]) fHclusterRatio[imode] = dynamic_cast<TH1F *>(ref.fHclusterRatio[imode]->Clone());\r
84     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
85       if(fHnSigmaTPC[ispec][imode]) fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());\r
86       if(fHnSigmaTOF[ispec][imode]) fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());\r
87     }\r
88   }\r
89 }\r
90 \r
91 //_____________________________________________________________________\r
92 AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){\r
93   //\r
94   // Assignment operator\r
95   //\r
96   if(this != &ref)\r
97     ref.Copy(*this);\r
98   return *this;\r
99 }\r
100 \r
101 //_____________________________________________________________________\r
102 AliESDpidCuts::~AliESDpidCuts(){\r
103   //\r
104   // Destructor\r
105   //\r
106   delete fESDpid;\r
107 \r
108   delete fHcutStatistics;\r
109   delete fHcutCorrelation;\r
110   for(Int_t imode = 0; imode < 2; imode++){\r
111     delete fHclusterRatio[imode];\r
112     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
113       delete fHnSigmaTPC[ispec][imode];\r
114       delete fHnSigmaTOF[ispec][imode];\r
115     }\r
116   }\r
117 }\r
118 \r
119 //_____________________________________________________________________\r
120 Bool_t AliESDpidCuts::IsSelected(TObject *obj){\r
121   //\r
122   // Select Track\r
123   // \r
124   AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);\r
125   if(!trk){\r
126     AliError("Provided object is not AliESDtrack!");\r
127     return kFALSE;\r
128   }\r
129   AliESDEvent* evt = trk->GetESDEvent();\r
130   if(!evt){\r
131     AliError("No AliESDEvent!");\r
132     return kFALSE;\r
133   }\r
134   return AcceptTrack(trk, evt);\r
135 }\r
136 \r
137 //_____________________________________________________________________\r
138 void AliESDpidCuts::Copy(TObject &c) const {\r
139   //\r
140   // Copy function\r
141   //\r
142   AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);\r
143 \r
144   target.fESDpid = new AliESDpid(*fESDpid);\r
145 \r
146   target.fCutTPCclusterRatio = fCutTPCclusterRatio;\r
147   target.fMinMomentumTOF = fMinMomentumTOF;\r
148 \r
149   target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;\r
150   target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;\r
151   \r
152   if(fHcutStatistics) target.fHcutStatistics = dynamic_cast<TH1I *>(fHcutStatistics->Clone());\r
153   if(fHcutCorrelation) target.fHcutCorrelation = dynamic_cast<TH2I *>(fHcutCorrelation->Clone());\r
154   for(Int_t imode = 0; imode < 2; imode++){\r
155     if(fHclusterRatio[imode]) target.fHclusterRatio[imode] = dynamic_cast<TH1F *>(fHclusterRatio[imode]->Clone());\r
156     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
157       if(fHnSigmaTPC[ispec][imode]) target.fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());\r
158       if(fHnSigmaTOF[ispec][imode]) target.fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTOF[ispec][imode]->Clone());\r
159     }\r
160   }\r
161  \r
162   memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
163   memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES);\r
164  \r
165   AliESDpidCuts::Copy(c);\r
166 }\r
167 \r
168 //_____________________________________________________________________\r
169 Long64_t AliESDpidCuts::Merge(TCollection *coll){\r
170   //\r
171   // Merge Cut objects\r
172   //\r
173   if(coll) return 0;\r
174   if(coll->IsEmpty()) return 1;\r
175   if(!HasHistograms())  return 0;\r
176   \r
177   TIterator *iter = coll->MakeIterator();\r
178   TObject *o = NULL; \r
179   AliESDpidCuts *ref = NULL;\r
180   Int_t counter = 0;\r
181   while((o = iter->Next())){\r
182     ref = dynamic_cast<AliESDpidCuts *>(o);\r
183     if(!ref) continue;\r
184     if(!ref->HasHistograms()) continue;\r
185 \r
186     fHcutStatistics->Add(ref->fHcutStatistics);\r
187     fHcutCorrelation->Add(ref->fHcutCorrelation);\r
188     for(Int_t imode = 0; imode < 2; imode++){\r
189       fHclusterRatio[imode]->Add(ref->fHclusterRatio[imode]);\r
190       for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
191         fHnSigmaTPC[ispec][imode]->Add(ref->fHnSigmaTPC[ispec][imode]);\r
192         fHnSigmaTOF[ispec][imode]->Add(ref->fHnSigmaTOF[ispec][imode]);\r
193       }\r
194     }\r
195     ++counter;\r
196   }\r
197   return ++counter;\r
198 }\r
199 \r
200 //_____________________________________________________________________\r
201 void AliESDpidCuts::DefineHistograms(Color_t color){\r
202   //\r
203   // Swich on QA and create the histograms\r
204   //\r
205   SetBit(kHasHistograms, kTRUE);\r
206   fHcutStatistics = new TH1I("fHcutStatistics", "Cut Statistics", kNcuts, 0, kNcuts);\r
207   fHcutStatistics->SetLineColor(color);\r
208   fHcutCorrelation = new TH2I("fHcutCorrelation", "Cut Correlation", kNcuts, 0, kNcuts, kNcuts, 0, kNcuts);\r
209   TString cutname[kNcuts] = {"TPCclusterRatio", "TPC sigma", "TOF sigma"};\r
210   for(Int_t icut = 0; icut < kNcuts; icut++){\r
211     fHcutStatistics->GetXaxis()->SetBinLabel(fHcutStatistics->GetXaxis()->GetFirst() + icut, cutname[icut].Data());\r
212     fHcutCorrelation->GetXaxis()->SetBinLabel(fHcutCorrelation->GetXaxis()->GetFirst() + icut, cutname[icut].Data());\r
213     fHcutCorrelation->GetYaxis()->SetBinLabel(fHcutCorrelation->GetYaxis()->GetFirst() + icut, cutname[icut].Data());\r
214   }\r
215   Char_t hname[256], htitle[256];\r
216   for(Int_t imode = 0; imode < 2; imode++){\r
217     sprintf(hname, "fHclusterRatio%s", imode ? "After" : "Before");\r
218     sprintf(htitle, "TPC cluster Ratio %s cuts;Ratio;Entries", imode ? "after" : "before");\r
219     fHclusterRatio[imode] = new TH1F(hname, htitle, 20, 0., 1.);\r
220     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
221       sprintf(hname, "fHnSigma%sTPC%s", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
222       sprintf(htitle, "TPC sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
223       fHnSigmaTPC[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);\r
224       sprintf(hname, "fHnSigma%sTOF%s", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
225       sprintf(htitle, "TOF sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
226       fHnSigmaTOF[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);\r
227     }\r
228   }\r
229 }\r
230 \r
231 //_____________________________________________________________________\r
232 Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack *track, const AliESDEvent *event){\r
233   //\r
234   // Check whether the tracks survived the cuts\r
235   //\r
236   enum{\r
237     kCutClusterRatioTPC,\r
238     kCutNsigmaTPC,\r
239     kCutNsigmaTOF\r
240   };\r
241   Long64_t cutRequired=0, cutFullfiled = 0;\r
242   if(fTOFsigmaCutRequired && event == 0) \r
243     AliError("No event pointer. Need event pointer for T0 for TOF cut");\r
244   Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;\r
245   if(fCutTPCclusterRatio > 0.){\r
246     SETBIT(cutRequired, kCutClusterRatioTPC);\r
247     if(clusterRatio >= fCutTPCclusterRatio) \r
248       SETBIT(cutFullfiled, kCutClusterRatioTPC);\r
249   }\r
250   // check TPC nSigma cut\r
251   Float_t nsigmaTPC[AliPID::kSPECIES], nsigma;   // need all sigmas for QA plotting\r
252   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
253     nsigmaTPC[ispec] = nsigma = fESDpid->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
254     if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;\r
255     SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required\r
256     if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC);    // Fullfiled for at least one species\r
257   }\r
258   // check TOF nSigma cut\r
259   Float_t nsigmaTOF[AliPID::kSPECIES];    // see above\r
260   Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set\r
261   Double_t times[AliPID::kSPECIES];\r
262   track->GetIntegratedTimes(times);\r
263   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
264     \r
265     if(hasTOFpid) nsigmaTOF[ispec] = nsigma = fESDpid->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec, event->GetT0());\r
266     if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;\r
267     SETBIT(cutRequired, kCutNsigmaTOF);\r
268     if(track->GetOuterParam()->P() >= fMinMomentumTOF){\r
269       if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);\r
270     }\r
271   }\r
272 \r
273   // Fill Histograms before cuts\r
274   if(HasHistograms()){\r
275     fHclusterRatio[0]->Fill(clusterRatio);\r
276     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
277       fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);\r
278       if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);\r
279     }\r
280   }\r
281   if(cutRequired != cutFullfiled){\r
282     // Fill cut statistics\r
283     if(HasHistograms()){\r
284       for(Int_t icut = 0; icut < kNcuts; icut++){\r
285         if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){\r
286           // cut not fullfiled\r
287           fHcutStatistics->Fill(icut);\r
288           for(Int_t jcut = 0; jcut <= icut; jcut++)\r
289             if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);\r
290         }\r
291       }\r
292     }\r
293     return kFALSE;    // At least one cut is not fullfiled\r
294   }\r
295 \r
296   // Fill Histograms after cuts\r
297   if(HasHistograms()){\r
298     fHclusterRatio[1]->Fill(clusterRatio);\r
299     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
300       fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);\r
301       if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);\r
302     }\r
303   }\r
304 \r
305   return kTRUE;\r
306 }\r
307 \r
308 //_____________________________________________________________________\r
309 void AliESDpidCuts::SaveHistograms(const Char_t * location){\r
310   //\r
311   // Save the histograms to a file\r
312   //\r
313   if(!HasHistograms()){\r
314     AliError("Histograms not on - Exiting");\r
315     return;\r
316   }\r
317   if(!location) location = GetName();\r
318   gDirectory->mkdir(location);\r
319   gDirectory->cd(location);\r
320   fHcutStatistics->Write();\r
321   fHcutCorrelation->Write();\r
322 \r
323   gDirectory->mkdir("before_cuts");\r
324   gDirectory->mkdir("after_cuts");\r
325 \r
326   gDirectory->cd("before_cuts");\r
327   fHclusterRatio[0]->Write();\r
328   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
329     fHnSigmaTPC[ispec][0]->Write();\r
330     fHnSigmaTOF[ispec][0]->Write();\r
331   }\r
332 \r
333   gDirectory->cd("../after_cuts");\r
334   fHclusterRatio[1]->Write();\r
335   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
336     fHnSigmaTPC[ispec][1]->Write();\r
337     fHnSigmaTOF[ispec][1]->Write();\r
338   }\r
339 \r
340   gDirectory->cd("..");\r
341 }\r
342 \r
343 //_____________________________________________________________________\r
344 void AliESDpidCuts::DrawHistograms(){\r
345   //\r
346   // Draw the Histograms\r
347   //\r
348   TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);\r
349   stat->cd();\r
350   fHcutStatistics->SetStats(kFALSE);\r
351   fHcutStatistics->Draw();\r
352   stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));\r
353 \r
354   TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);\r
355   correl->cd();\r
356   fHcutCorrelation->SetStats(kFALSE);\r
357   fHcutCorrelation->Draw("colz");\r
358   correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));\r
359 \r
360   TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);\r
361   cRatio->cd();\r
362   fHclusterRatio[0]->SetLineColor(kRed);\r
363   fHclusterRatio[0]->SetStats(kFALSE);\r
364   fHclusterRatio[0]->Draw();\r
365   fHclusterRatio[1]->SetLineColor(kBlue);\r
366   fHclusterRatio[1]->SetStats(kFALSE);\r
367   fHclusterRatio[1]->Draw("same");\r
368   cRatio->SaveAs(Form("%s_%s.gif",  GetName(), cRatio->GetName()));\r
369 \r
370   TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);\r
371   cNsigTPC->Divide(3,2);\r
372   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
373     cNsigTPC->cd(ispec + 1);\r
374     fHnSigmaTPC[ispec][0]->SetLineColor(kRed);\r
375     fHnSigmaTPC[ispec][0]->SetStats(kFALSE);\r
376     fHnSigmaTPC[ispec][0]->Draw();\r
377     fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);\r
378     fHnSigmaTPC[ispec][1]->SetStats(kFALSE);\r
379     fHnSigmaTPC[ispec][1]->Draw("same");\r
380   }\r
381   cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));\r
382 \r
383   TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);\r
384   cNsigTOF->Divide(3,2);\r
385   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
386     cNsigTOF->cd(ispec + 1);\r
387     fHnSigmaTOF[ispec][0]->SetLineColor(kRed);\r
388     fHnSigmaTOF[ispec][0]->SetStats(kFALSE);\r
389     fHnSigmaTOF[ispec][0]->Draw();\r
390     fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);\r
391     fHnSigmaTOF[ispec][1]->SetStats(kFALSE);\r
392     fHnSigmaTOF[ispec][1]->Draw("same");\r
393   }\r
394   cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));\r
395 }\r
396 \r