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