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