Fixed warning
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.cxx
CommitLineData
979d80f1 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
34ClassImp(AliESDpidCuts)\r
35\r
36const Int_t AliESDpidCuts::kNcuts = 3;\r
37\r
38//_____________________________________________________________________\r
39AliESDpidCuts::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
61cfa442 54 memset(fCutTPCnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);\r
55 memset(fCutTOFnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);\r
979d80f1 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
63AliESDpidCuts::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
92AliESDpidCuts &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
102AliESDpidCuts::~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
120Bool_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
138void 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
61cfa442 163 memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
979d80f1 164 \r
165 AliESDpidCuts::Copy(c);\r
166}\r
167\r
168//_____________________________________________________________________\r
169Long64_t AliESDpidCuts::Merge(TCollection *coll){\r
170 //\r
171 // Merge Cut objects\r
172 //\r
2c881cf4 173 if(!coll) return 0;\r
174 if(coll->IsEmpty()) return 1;\r
979d80f1 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
201void 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
accb0197 217 snprintf(hname, 256, "fHclusterRatio%s", imode ? "After" : "Before");\r
218 snprintf(htitle, 256, "TPC cluster Ratio %s cuts;Ratio;Entries", imode ? "after" : "before");\r
979d80f1 219 fHclusterRatio[imode] = new TH1F(hname, htitle, 20, 0., 1.);\r
220 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
accb0197 221 snprintf(hname, 256, "fHnSigma%sTPC%s", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
222 snprintf(htitle, 256, "TPC sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
979d80f1 223 fHnSigmaTPC[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);\r
accb0197 224 snprintf(hname, 256, "fHnSigma%sTOF%s", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
225 snprintf(htitle, 256, "TOF sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
979d80f1 226 fHnSigmaTOF[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);\r
227 }\r
228 }\r
229}\r
230\r
231//_____________________________________________________________________\r
232Bool_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
2c881cf4 242 if(fTOFsigmaCutRequired && event == 0) {\r
979d80f1 243 AliError("No event pointer. Need event pointer for T0 for TOF cut");\r
2c881cf4 244 return (0);\r
245 }\r
979d80f1 246 Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;\r
247 if(fCutTPCclusterRatio > 0.){\r
248 SETBIT(cutRequired, kCutClusterRatioTPC);\r
249 if(clusterRatio >= fCutTPCclusterRatio) \r
250 SETBIT(cutFullfiled, kCutClusterRatioTPC);\r
251 }\r
252 // check TPC nSigma cut\r
253 Float_t nsigmaTPC[AliPID::kSPECIES], nsigma; // need all sigmas for QA plotting\r
254 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
255 nsigmaTPC[ispec] = nsigma = fESDpid->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
256 if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;\r
257 SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required\r
258 if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC); // Fullfiled for at least one species\r
259 }\r
260 // check TOF nSigma cut\r
261 Float_t nsigmaTOF[AliPID::kSPECIES]; // see above\r
262 Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set\r
263 Double_t times[AliPID::kSPECIES];\r
264 track->GetIntegratedTimes(times);\r
265 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
266 \r
2c881cf4 267 if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fESDpid->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec, event->GetT0());\r
979d80f1 268 if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;\r
269 SETBIT(cutRequired, kCutNsigmaTOF);\r
270 if(track->GetOuterParam()->P() >= fMinMomentumTOF){\r
271 if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);\r
272 }\r
273 }\r
274\r
275 // Fill Histograms before cuts\r
276 if(HasHistograms()){\r
277 fHclusterRatio[0]->Fill(clusterRatio);\r
278 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
279 fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);\r
280 if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);\r
281 }\r
282 }\r
283 if(cutRequired != cutFullfiled){\r
284 // Fill cut statistics\r
285 if(HasHistograms()){\r
286 for(Int_t icut = 0; icut < kNcuts; icut++){\r
287 if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){\r
288 // cut not fullfiled\r
289 fHcutStatistics->Fill(icut);\r
290 for(Int_t jcut = 0; jcut <= icut; jcut++)\r
291 if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);\r
292 }\r
293 }\r
294 }\r
295 return kFALSE; // At least one cut is not fullfiled\r
296 }\r
297\r
298 // Fill Histograms after cuts\r
299 if(HasHistograms()){\r
300 fHclusterRatio[1]->Fill(clusterRatio);\r
301 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
302 fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);\r
303 if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);\r
304 }\r
305 }\r
306\r
307 return kTRUE;\r
308}\r
309\r
310//_____________________________________________________________________\r
311void AliESDpidCuts::SaveHistograms(const Char_t * location){\r
312 //\r
313 // Save the histograms to a file\r
314 //\r
315 if(!HasHistograms()){\r
316 AliError("Histograms not on - Exiting");\r
317 return;\r
318 }\r
319 if(!location) location = GetName();\r
320 gDirectory->mkdir(location);\r
321 gDirectory->cd(location);\r
322 fHcutStatistics->Write();\r
323 fHcutCorrelation->Write();\r
324\r
325 gDirectory->mkdir("before_cuts");\r
326 gDirectory->mkdir("after_cuts");\r
327\r
328 gDirectory->cd("before_cuts");\r
329 fHclusterRatio[0]->Write();\r
330 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
331 fHnSigmaTPC[ispec][0]->Write();\r
332 fHnSigmaTOF[ispec][0]->Write();\r
333 }\r
334\r
335 gDirectory->cd("../after_cuts");\r
336 fHclusterRatio[1]->Write();\r
337 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
338 fHnSigmaTPC[ispec][1]->Write();\r
339 fHnSigmaTOF[ispec][1]->Write();\r
340 }\r
341\r
342 gDirectory->cd("..");\r
343}\r
344\r
345//_____________________________________________________________________\r
346void AliESDpidCuts::DrawHistograms(){\r
347 //\r
348 // Draw the Histograms\r
349 //\r
350 TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);\r
351 stat->cd();\r
352 fHcutStatistics->SetStats(kFALSE);\r
353 fHcutStatistics->Draw();\r
354 stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));\r
355\r
356 TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);\r
357 correl->cd();\r
358 fHcutCorrelation->SetStats(kFALSE);\r
359 fHcutCorrelation->Draw("colz");\r
360 correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));\r
361\r
362 TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);\r
363 cRatio->cd();\r
364 fHclusterRatio[0]->SetLineColor(kRed);\r
365 fHclusterRatio[0]->SetStats(kFALSE);\r
366 fHclusterRatio[0]->Draw();\r
367 fHclusterRatio[1]->SetLineColor(kBlue);\r
368 fHclusterRatio[1]->SetStats(kFALSE);\r
369 fHclusterRatio[1]->Draw("same");\r
370 cRatio->SaveAs(Form("%s_%s.gif", GetName(), cRatio->GetName()));\r
371\r
372 TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);\r
373 cNsigTPC->Divide(3,2);\r
374 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
375 cNsigTPC->cd(ispec + 1);\r
376 fHnSigmaTPC[ispec][0]->SetLineColor(kRed);\r
377 fHnSigmaTPC[ispec][0]->SetStats(kFALSE);\r
378 fHnSigmaTPC[ispec][0]->Draw();\r
379 fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);\r
380 fHnSigmaTPC[ispec][1]->SetStats(kFALSE);\r
381 fHnSigmaTPC[ispec][1]->Draw("same");\r
382 }\r
383 cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));\r
384\r
385 TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);\r
386 cNsigTOF->Divide(3,2);\r
387 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
388 cNsigTOF->cd(ispec + 1);\r
389 fHnSigmaTOF[ispec][0]->SetLineColor(kRed);\r
390 fHnSigmaTOF[ispec][0]->SetStats(kFALSE);\r
391 fHnSigmaTOF[ispec][0]->Draw();\r
392 fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);\r
393 fHnSigmaTOF[ispec][1]->SetStats(kFALSE);\r
394 fHnSigmaTOF[ispec][1]->Draw("same");\r
395 }\r
396 cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));\r
397}\r
398\r