Coverity 18632
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskPIDCombined.cxx
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
38 const 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
47 ClassImp(AliAnalysisTaskPIDCombined)\r
48 \r
49 //_________________________________________________________________________________\r
50 AliAnalysisTaskPIDCombined::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
75 AliAnalysisTaskPIDCombined::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
102 void 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\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<AliPID::kSPECIES; ++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
232 void 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
312 void 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
321 void 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
331 Int_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