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