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