]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ANALYSIS/AliAnalysisTaskPIDCombined.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
including new Salvatore's patch
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskPIDCombined.cxx
... / ...
CommitLineData
1/*************************************************************************\r
2* Copyright(c) 1998-2009, 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\r
16///////////////////////////////////////////////////////////////////////////\r
17// //\r
18// Basic Analysis Task //\r
19// for PID Analysis //\r
20// //\r
21///////////////////////////////////////////////////////////////////////////\r
22\r
23#include <TChain.h>\r
24#include <TH1D.h>\r
25#include <TH2D.h>\r
26\r
27#include <AliVEvent.h>\r
28#include <AliInputEventHandler.h>\r
29#include <AliAnalysisManager.h>\r
30\r
31#include <AliLog.h>\r
32#include <AliPID.h>\r
33#include <AliPIDCombined.h>\r
34#include <AliPIDResponse.h>\r
35\r
36#include "AliAnalysisTaskPIDCombined.h"\r
37\r
38const char *AliAnalysisTaskPIDCombined::fgkBinMomDesc[AliAnalysisTaskPIDCombined::kPtBins] = {\r
39 " 0 <= p < 0.5 GeV/c",\r
40 " 0.5 <= p < 0.7 GeV/c",\r
41 " 0.7 <= p < 1.0 GeV/c",\r
42 " 1.0 <= p < 1.5 GeV/c",\r
43 " 1.5 <= p < 2.0 GeV/c",\r
44 " p >= 2.0 GeV/c"\r
45};\r
46\r
47ClassImp(AliAnalysisTaskPIDCombined)\r
48\r
49//_________________________________________________________________________________\r
50AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined() :\r
51 AliAnalysisTaskSE(),\r
52 fHistList(),\r
53 fProbTPCnSigma(),\r
54 fProbTOFnSigma(),\r
55 fProbTPCTOFnSigmaTPC(),\r
56 fProbTPC(),\r
57 fProbTOF(),\r
58 fProbTPCTOF(),\r
59 fPriors(),\r
60 fProbTPCTOFnSigTPCMom(),\r
61 fProbTPCnSigTPCMom(),\r
62 fProbTOFnSigTOFMom(),\r
63 fPriorsUsed(),\r
64 fPIDResponse(0x0),\r
65 fPIDCombined(0x0),\r
66 fTrackCuts(0x0),\r
67 fTrackFilter(0x0)\r
68{\r
69 //\r
70 // Constructor\r
71 //\r
72}\r
73\r
74//_________________________________________________________________________________\r
75AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined(const char *name) :\r
76 AliAnalysisTaskSE(name),\r
77 fHistList(),\r
78 fProbTPCnSigma(),\r
79 fProbTOFnSigma(),\r
80 fProbTPCTOFnSigmaTPC(),\r
81 fProbTPC(),\r
82 fProbTOF(),\r
83 fProbTPCTOF(),\r
84 fPriors(),\r
85 fProbTPCTOFnSigTPCMom(),\r
86 fProbTPCnSigTPCMom(),\r
87 fProbTOFnSigTOFMom(),\r
88 fPriorsUsed(),\r
89 fPIDResponse(0x0),\r
90 fPIDCombined(0x0),\r
91 fTrackCuts(0x0),\r
92 fTrackFilter(0x0)\r
93{\r
94 //\r
95 // Constructor\r
96 //\r
97 DefineInput(0,TChain::Class());\r
98 DefineOutput(1, TList::Class());\r
99}\r
100\r
101//_________________________________________________________________________________\r
102void AliAnalysisTaskPIDCombined::UserCreateOutputObjects()\r
103{\r
104 //\r
105 // Initialise the framework objects\r
106 //\r
107\r
108\r
109 // ------- track cuts\r
110 fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");\r
111 fTrackCuts->SetAcceptKinkDaughters(kFALSE);\r
112 fTrackCuts->SetMinNClustersTPC(80);\r
113 fTrackCuts->SetMaxChi2PerClusterTPC(4);\r
114 fTrackCuts->SetMaxDCAToVertexXY(3);\r
115 fTrackCuts->SetMaxDCAToVertexZ(3);\r
116 fTrackCuts->SetRequireTPCRefit(kTRUE);\r
117 fTrackCuts->SetRequireITSRefit(kTRUE);\r
118 fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
119 fTrackFilter = new AliAnalysisFilter("trackFilter");\r
120 fTrackFilter->AddCuts(fTrackCuts);\r
121\r
122\r
123\r
124 // ------- setup PIDCombined\r
125 fPIDCombined=new AliPIDCombined;\r
126 fPIDCombined->SetDefaultTPCPriors();\r
127 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);\r
128\r
129 // no light nuclei - no need to call it, this is default\r
130 // fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);\r
131\r
132\r
133 fHistList.Add(new TH1D("nEvents","Number of Evnets;Selection",2,0,2));\r
134\r
135 for (Int_t ispec=0; ispec<5; ++ispec){\r
136\r
137 \r
138 fProbTPC[ispec]=new TH2D(Form("prob%s_mom_TPC",AliPID::ParticleName(ispec)),\r
139 Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
140 100,0.,20.,50,0.,1.);\r
141 fHistList.Add(fProbTPC[ispec]);\r
142 fProbTPCnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TPC",AliPID::ParticleName(ispec)),\r
143 Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
144 20,-5.,5.,50,0.,1.);\r
145 fHistList.Add(fProbTPCnSigma[ispec]);\r
146\r
147 for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
148 fProbTPCnSigTPCMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
149 Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
150 20,-5.,5.,50,0.,1.);\r
151 fHistList.Add(fProbTPCnSigTPCMom[ibin][ispec]);\r
152 }\r
153\r
154\r
155\r
156 fProbTOF[ispec]=new TH2D(Form("prob%s_mom_TOF",AliPID::ParticleName(ispec)),\r
157 Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
158 100,0.,20.,50,0.,1.);\r
159 fHistList.Add(fProbTOF[ispec]);\r
160 fProbTOFnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TOF",AliPID::ParticleName(ispec)),\r
161 Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
162 20,-5.,5.,50,0.,1.);\r
163 fHistList.Add(fProbTOFnSigma[ispec]);\r
164 for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
165 fProbTOFnSigTOFMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TOF (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
166 Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
167 20,-5.,5.,50,0.,1.);\r
168 fHistList.Add(fProbTOFnSigTOFMom[ibin][ispec]);\r
169 }\r
170\r
171\r
172\r
173 fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),\r
174 Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
175 100,0.,20.,50,0.,1.);\r
176 fHistList.Add(fProbTPCTOF[ispec]);\r
177 fProbTPCTOFnSigmaTPC[ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC",AliPID::ParticleName(ispec)),\r
178 Form("%s TPCTOF probability vs. n#sigmaTPC;n#sigma;probability",AliPID::ParticleName(ispec)),\r
179 20,-5.,5.,50,0.,1.);\r
180 fHistList.Add(fProbTPCTOFnSigmaTPC[ispec]);\r
181 for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
182 fProbTPCTOFnSigTPCMom[ibin][ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
183 Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
184 20,-5.,5.,50,0.,1.);\r
185 fHistList.Add(fProbTPCTOFnSigTPCMom[ibin][ispec]);\r
186 }\r
187\r
188\r
189\r
190 // basic priors\r
191 fPriors[ispec]=new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),\r
192 Form("%s priors vs momentum",AliPID::ParticleName(ispec)),\r
193 100,0.,20.);\r
194 fHistList.Add(fPriors[ispec]);\r
195 switch (ispec) {\r
196 case AliPID::kElectron:\r
197 for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);\r
198 break;\r
199 case AliPID::kMuon:\r
200 for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);\r
201 break;\r
202 case AliPID::kPion:\r
203 for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);\r
204 break;\r
205 case AliPID::kKaon:\r
206 for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);\r
207 break;\r
208 case AliPID::kProton:\r
209 for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);\r
210 break;\r
211 default:\r
212 break;\r
213 }\r
214 fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);\r
215\r
216 // priors used\r
217 fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),\r
218 Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),\r
219 100,0.,20.,101,0,1.01); \r
220 fHistList.Add(fPriorsUsed[ispec]);\r
221 }\r
222\r
223\r
224\r
225 fHistList.SetOwner();\r
226 PostData(1,&fHistList);\r
227\r
228\r
229}\r
230\r
231//_________________________________________________________________________________\r
232void AliAnalysisTaskPIDCombined::UserExec(Option_t *)\r
233{\r
234 //\r
235 // Main loop. Called for every event\r
236 //\r
237 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();\r
238 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());\r
239 fPIDResponse=inputHandler->GetPIDResponse();\r
240 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
241\r
242 // Printf(" ---------------------- UserExec PID task ---------------------");\r
243 \r
244 FillHistogram("nEvents",0.);\r
245\r
246 AliVEvent *event=InputEvent();\r
247 AliVTrack *track=0x0;\r
248 Int_t ntracks=event->GetNumberOfTracks();\r
249\r
250 Double_t probTPC[AliPID::kSPECIES]={0.};\r
251 Double_t probTOF[AliPID::kSPECIES]={0.};\r
252 Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
253 \r
254 //loop over all tracks\r
255 for (Int_t itrack=0; itrack<ntracks; ++itrack){\r
256\r
257 track=(AliVTrack*)event->GetTrack(itrack);\r
258\r
259 if ( fTrackFilter->IsSelected(track) ) {\r
260\r
261 Double_t mom=track->GetTPCmomentum();\r
262 Int_t ibin=GetMomBin(mom);\r
263\r
264 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
265 UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
266 \r
267 if (detUsed != 0) { // TPC is available\r
268 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
269 Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
270 fProbTPC[ispec]->Fill(mom,probTPC[ispec]);\r
271 fProbTPCnSigma[ispec]->Fill(nSigmaTPC,probTPC[ispec]);\r
272 fProbTPCnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPC[ispec]);\r
273 }\r
274\r
275 // compute priors for TPC+TOF, even if we ask just TOF for PID\r
276 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
277 detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
278 Double_t priors[5]; // check priors used for TOF\r
279 fPIDCombined->GetPriors(track,priors,fPIDResponse,detUsed);\r
280 for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(track->Pt()),priors[ispec]);\r
281\r
282 if (detUsed != 0) { // TOF is available\r
283 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
284 Double_t nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);\r
285 fProbTOF[ispec]->Fill(mom,probTOF[ispec]);\r
286 fProbTOFnSigma[ispec]->Fill(nSigmaTOF,probTOF[ispec]);\r
287 fProbTOFnSigTOFMom[ibin][ispec]->Fill(nSigmaTOF,probTOF[ispec]);\r
288 }\r
289 }\r
290\r
291 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
292 detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
293 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {\r
294 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
295 Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
296 fProbTPCTOF[ispec]->Fill(mom,probTPCTOF[ispec]);\r
297 fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);\r
298 fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);\r
299 }\r
300 }\r
301\r
302 }\r
303\r
304\r
305 }\r
306 }\r
307\r
308 PostData(1, &fHistList);\r
309}\r
310\r
311//_________________________________________________________________________________\r
312void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t weight)\r
313{\r
314 //\r
315 // Fill 1D histogram by name\r
316 //\r
317 ((TH1*)fHistList.FindObject(name))->Fill(x,weight);\r
318}\r
319\r
320//_________________________________________________________________________________\r
321void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t y, Double_t weight)\r
322{\r
323 //\r
324 // Fill 2D histogram by name\r
325 //\r
326 ((TH2*)fHistList.FindObject(name))->Fill(x,y,weight);\r
327}\r
328\r
329\r
330//_________________________________________________________________________________\r
331Int_t AliAnalysisTaskPIDCombined::GetMomBin(Float_t mom)\r
332{\r
333 //\r
334 // Given momentum return histogram to be filled\r
335 //\r
336 if (mom>0. && mom < 0.5) return 0;\r
337 if (mom>=0.5 && mom < 0.7) return 1;\r
338 if (mom>=0.7 && mom < 1.0) return 2;\r
339 if (mom>=1.0 && mom < 1.5) return 3;\r
340 if (mom>=1.5 && mom < 2.0) return 4;\r
341 return kPtBins-1;\r
342}\r
343\r