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