]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliESDpidCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
renamed Cosmics folder to Cosmic
[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
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
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
163 memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES);\r
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
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
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
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
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
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
309void 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
344void 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