]> 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
Fix for checking dependencies of a LEGO module
[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
250 AliError("PID Response not available");\r
251 return 0;\r
252 }\r
979d80f1 253 enum{\r
254 kCutClusterRatioTPC,\r
255 kCutNsigmaTPC,\r
256 kCutNsigmaTOF\r
257 };\r
258 Long64_t cutRequired=0, cutFullfiled = 0;\r
2c881cf4 259 if(fTOFsigmaCutRequired && event == 0) {\r
979d80f1 260 AliError("No event pointer. Need event pointer for T0 for TOF cut");\r
2c881cf4 261 return (0);\r
262 }\r
979d80f1 263 Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;\r
264 if(fCutTPCclusterRatio > 0.){\r
265 SETBIT(cutRequired, kCutClusterRatioTPC);\r
266 if(clusterRatio >= fCutTPCclusterRatio) \r
267 SETBIT(cutFullfiled, kCutClusterRatioTPC);\r
268 }\r
269 // check TPC nSigma cut\r
270 Float_t nsigmaTPC[AliPID::kSPECIES], nsigma; // need all sigmas for QA plotting\r
271 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
320386a4 272 nsigmaTPC[ispec] = nsigma = fPIDresponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
979d80f1 273 if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;\r
274 SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required\r
275 if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC); // Fullfiled for at least one species\r
276 }\r
277 // check TOF nSigma cut\r
278 Float_t nsigmaTOF[AliPID::kSPECIES]; // see above\r
279 Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set\r
280 Double_t times[AliPID::kSPECIES];\r
281 track->GetIntegratedTimes(times);\r
282 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
283 \r
320386a4 284 if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fPIDresponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);\r
979d80f1 285 if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;\r
286 SETBIT(cutRequired, kCutNsigmaTOF);\r
287 if(track->GetOuterParam()->P() >= fMinMomentumTOF){\r
288 if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);\r
289 }\r
290 }\r
291\r
292 // Fill Histograms before cuts\r
293 if(HasHistograms()){\r
294 fHclusterRatio[0]->Fill(clusterRatio);\r
295 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
296 fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);\r
297 if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);\r
298 }\r
299 }\r
300 if(cutRequired != cutFullfiled){\r
301 // Fill cut statistics\r
302 if(HasHistograms()){\r
303 for(Int_t icut = 0; icut < kNcuts; icut++){\r
304 if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){\r
305 // cut not fullfiled\r
306 fHcutStatistics->Fill(icut);\r
307 for(Int_t jcut = 0; jcut <= icut; jcut++)\r
308 if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);\r
309 }\r
310 }\r
311 }\r
312 return kFALSE; // At least one cut is not fullfiled\r
313 }\r
314\r
315 // Fill Histograms after cuts\r
316 if(HasHistograms()){\r
317 fHclusterRatio[1]->Fill(clusterRatio);\r
318 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
319 fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);\r
320 if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);\r
321 }\r
322 }\r
323\r
324 return kTRUE;\r
325}\r
326\r
327//_____________________________________________________________________\r
328void AliESDpidCuts::SaveHistograms(const Char_t * location){\r
329 //\r
330 // Save the histograms to a file\r
331 //\r
332 if(!HasHistograms()){\r
333 AliError("Histograms not on - Exiting");\r
334 return;\r
335 }\r
336 if(!location) location = GetName();\r
337 gDirectory->mkdir(location);\r
338 gDirectory->cd(location);\r
339 fHcutStatistics->Write();\r
340 fHcutCorrelation->Write();\r
341\r
342 gDirectory->mkdir("before_cuts");\r
343 gDirectory->mkdir("after_cuts");\r
344\r
345 gDirectory->cd("before_cuts");\r
346 fHclusterRatio[0]->Write();\r
347 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
348 fHnSigmaTPC[ispec][0]->Write();\r
349 fHnSigmaTOF[ispec][0]->Write();\r
350 }\r
351\r
352 gDirectory->cd("../after_cuts");\r
353 fHclusterRatio[1]->Write();\r
354 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
355 fHnSigmaTPC[ispec][1]->Write();\r
356 fHnSigmaTOF[ispec][1]->Write();\r
357 }\r
358\r
359 gDirectory->cd("..");\r
360}\r
361\r
362//_____________________________________________________________________\r
363void AliESDpidCuts::DrawHistograms(){\r
364 //\r
365 // Draw the Histograms\r
366 //\r
367 TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);\r
368 stat->cd();\r
369 fHcutStatistics->SetStats(kFALSE);\r
370 fHcutStatistics->Draw();\r
371 stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));\r
372\r
373 TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);\r
374 correl->cd();\r
375 fHcutCorrelation->SetStats(kFALSE);\r
376 fHcutCorrelation->Draw("colz");\r
377 correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));\r
378\r
379 TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);\r
380 cRatio->cd();\r
381 fHclusterRatio[0]->SetLineColor(kRed);\r
382 fHclusterRatio[0]->SetStats(kFALSE);\r
383 fHclusterRatio[0]->Draw();\r
384 fHclusterRatio[1]->SetLineColor(kBlue);\r
385 fHclusterRatio[1]->SetStats(kFALSE);\r
386 fHclusterRatio[1]->Draw("same");\r
387 cRatio->SaveAs(Form("%s_%s.gif", GetName(), cRatio->GetName()));\r
388\r
389 TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);\r
390 cNsigTPC->Divide(3,2);\r
391 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
392 cNsigTPC->cd(ispec + 1);\r
393 fHnSigmaTPC[ispec][0]->SetLineColor(kRed);\r
394 fHnSigmaTPC[ispec][0]->SetStats(kFALSE);\r
395 fHnSigmaTPC[ispec][0]->Draw();\r
396 fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);\r
397 fHnSigmaTPC[ispec][1]->SetStats(kFALSE);\r
398 fHnSigmaTPC[ispec][1]->Draw("same");\r
399 }\r
400 cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));\r
401\r
402 TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);\r
403 cNsigTOF->Divide(3,2);\r
404 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
405 cNsigTOF->cd(ispec + 1);\r
406 fHnSigmaTOF[ispec][0]->SetLineColor(kRed);\r
407 fHnSigmaTOF[ispec][0]->SetStats(kFALSE);\r
408 fHnSigmaTOF[ispec][0]->Draw();\r
409 fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);\r
410 fHnSigmaTOF[ispec][1]->SetStats(kFALSE);\r
411 fHnSigmaTOF[ispec][1]->Draw("same");\r
412 }\r
413 cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));\r
414}\r
415\r