]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/pid/AliAnalysisTaskLambdaBayes.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added macro for iterative computationn of priors
[u/mrichter/AliRoot.git] / PWGPP / pid / AliAnalysisTaskLambdaBayes.cxx
CommitLineData
a8ad4709 1#include "AliAnalysisTaskLambdaBayes.h"
2
3// ROOT includes
4#include <TMath.h>
5
6// AliRoot includes
7#include "AliInputEventHandler.h"
8#include "AliAODEvent.h"
9#include "AliAODVertex.h"
10#include "AliAODTrack.h"
11#include "AliESDtrack.h"
12#include "AliCentrality.h"
13#include "AliVHeader.h"
14#include "AliAODVZERO.h"
15#include "TFile.h"
16#include "AliOADBContainer.h"
17#include "TH2F.h"
18#include "TF1.h"
19#include "AliGenHijingEventHeader.h"
20#include "AliMCEvent.h"
21#include "AliAODMCHeader.h"
22#include "AliAODMCParticle.h"
23#include "TChain.h"
24#include "AliESDtrackCuts.h"
25#include "AliESDVertex.h"
26#include "AliEventplane.h"
27#include "AliAnalysisManager.h"
28#include "TRandom.h"
29
30
31Float_t AliAnalysisTaskLambdaBayes::fPtLambdaMin[AliAnalysisTaskLambdaBayes::nPtBin] = {0.,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.};// ptmin bin
32Float_t AliAnalysisTaskLambdaBayes::fPtLambdaMax[AliAnalysisTaskLambdaBayes::nPtBin] = {0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.,10.};// ptmax bin
33
34// STL includes
35//#include <iostream>
36//using namespace std;
37
38ClassImp(AliAnalysisTaskLambdaBayes)
39//_____________________________________________________________________________
40AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes():
41 AliAnalysisTaskSE(),
42 fVtxCut(10.0), // cut on |vertex| < fVtxCut
43 fEtaCut(0.8), // cut on |eta| < fEtaCut
44 fMinPt(0.15), // cut on pt > fMinPt
45 fIsMC(kFALSE),
46 fQAsw(kFALSE),
47 fNcluster(70),
48 fFilterBit(1),
49 fList(new TList()),
50 fList2(new TList()),
51 fList3(new TList()),
52 fCentrality(-1),
53 fPsi(0),
54 fPtLambdaC(0.),
55 fPhiLambdaC(0.),
56 fEtaLambda(0.),
57 fMassV0(0.),
58 fPtKp(0.),
59 fPhiKp(0.),
60 fEtaKp(0.),
61 fPtKn(0.),
62 fPhiKn(0.),
63 fEtaKn(0.),
64 fPidKp(0),
65 fPidKn(0),
66 fTOFTPCsignal(0),
67 fCombsignal(0),
68 fCutsDaughter(NULL),
69 fPIDCombined(NULL),
70 fContPid(NULL),
71 fContPid2(NULL),
520899fb 72 fContUser(NULL),
73 fContUser2(NULL),
a8ad4709 74 fNLambda(0),
75 fNpPos(0),
76 fNpNeg(0),
77 fHmismTOF(0),
9c668341 78 fHchannelTOFdistr(0),
520899fb 79 fTypeCol(2),
8c1aafae 80 fPIDuserCut(NULL),
bc629e38 81 fToEP(kFALSE),
82 fSpeciesRef(4)
a8ad4709 83{
84 // Default constructor (should not be used)
85 fList->SetName("contLambdaBayes1");
86 fList2->SetName("contLambdaBayes2");
87 fList3->SetName("contLambdaBayes3");
88
89 fList->SetOwner(kTRUE);
90 fList2->SetOwner(kTRUE);
91 fList3->SetOwner(kTRUE);
92
93 TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
94 fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
95
96 TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
97 fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
f4c4ac9e 98
99 for(Int_t i=0;i < nCentrBin;i++){
100 fElTOF[i] = NULL;
101 fElTPC[i] = NULL;
102 fPiTOF[i] = NULL;
103 fPiTPC[i] = NULL;
104 fKaTOF[i] = NULL;
105 fKaTPC[i] = NULL;
106 fPrTOF[i] = NULL;
107 fPrTPC[i] = NULL;
108 }
109 for(Int_t i=0;i < 4;i++){
110 hMatching[i] = NULL;
111 hTracking[i] = NULL;
112 }
113
114 for(Int_t i=0;i < 1000;i++){
115 fPhiLambda[i] = 0.0;
116 fPtLambda[i] = 0.0;
117 fIPPos[i] = 0;
118 fIPNeg[i] = 0;
119 fIpP[i] = 0;
120 fIpN[i] = 0;
121 fMassLambda[i] = 0.0;
122 }
a8ad4709 123}
124
125//______________________________________________________________________________
126AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes(const char *name):
127 AliAnalysisTaskSE(name),
128 fVtxCut(10.0), // cut on |vertex| < fVtxCut
129 fEtaCut(0.8), // cut on |eta| < fEtaCut
130 fMinPt(0.15), // cut on pt > fMinPt
131 fIsMC(kFALSE),
132 fQAsw(kFALSE),
133 fNcluster(70),
134 fFilterBit(1),
135 fList(new TList()),
136 fList2(new TList()),
137 fList3(new TList()),
138 fCentrality(-1),
139 fPsi(0),
140 fPtLambdaC(0.),
141 fPhiLambdaC(0.),
142 fEtaLambda(0.),
143 fMassV0(0.),
144 fPtKp(0.),
145 fPhiKp(0.),
146 fEtaKp(0.),
147 fPtKn(0.),
148 fPhiKn(0.),
149 fEtaKn(0.),
150 fPidKp(0),
151 fPidKn(0),
152 fTOFTPCsignal(0),
153 fCombsignal(0),
154 fCutsDaughter(NULL),
155 fPIDCombined(NULL),
156 fContPid(NULL),
157 fContPid2(NULL),
520899fb 158 fContUser(NULL),
159 fContUser2(NULL),
a8ad4709 160 fNLambda(0),
161 fNpPos(0),
162 fNpNeg(0),
163 fHmismTOF(0),
9c668341 164 fHchannelTOFdistr(0),
520899fb 165 fTypeCol(2),
8c1aafae 166 fPIDuserCut(NULL),
bc629e38 167 fToEP(kFALSE),
168 fSpeciesRef(4)
a8ad4709 169{
170
171 DefineOutput(1, TList::Class());
172 DefineOutput(2, TList::Class());
173 DefineOutput(3, TList::Class());
174 DefineOutput(4, TList::Class());
175
176 // Output slot #1 writes into a TTree
177 fList->SetName("contLambdaBayes1");
178 fList2->SetName("contLambdaBayes2");
179 fList3->SetName("contLambdaBayes3");
180
181 fList->SetOwner(kTRUE);
182 fList2->SetOwner(kTRUE);
183 fList3->SetOwner(kTRUE);
184
185 TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
186 fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
187
188 TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
189 fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
f4c4ac9e 190
191 for(Int_t i=0;i < nCentrBin;i++){
192 fElTOF[i] = NULL;
193 fElTPC[i] = NULL;
194 fPiTOF[i] = NULL;
195 fPiTPC[i] = NULL;
196 fKaTOF[i] = NULL;
197 fKaTPC[i] = NULL;
198 fPrTOF[i] = NULL;
199 fPrTPC[i] = NULL;
200 }
201 for(Int_t i=0;i < 4;i++){
202 hMatching[i] = NULL;
203 hTracking[i] = NULL;
204 }
205
206 for(Int_t i=0;i < 1000;i++){
207 fPhiLambda[i] = 0.0;
208 fPtLambda[i] = 0.0;
209 fIPPos[i] = 0;
210 fIPNeg[i] = 0;
211 fIpP[i] = 0;
212 fIpN[i] = 0;
213 fMassLambda[i] = 0.0;
214 }
a8ad4709 215}
216//_____________________________________________________________________________
217AliAnalysisTaskLambdaBayes::~AliAnalysisTaskLambdaBayes()
218{
219
220}
221
222//______________________________________________________________________________
223void AliAnalysisTaskLambdaBayes::UserCreateOutputObjects()
224{
225
226 fPIDCombined=new AliPIDCombined;
227 fPIDCombined->SetDefaultTPCPriors();
228 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
229
230 Float_t invMmin = 1.115-0.005*8;
231 Float_t invMmax = 1.115+0.005*8;
232
233 const Int_t nBinPid = 14;
234
8c1aafae 235 Int_t binPid[nBinPid] = {1/*ptKs*/,1+7*(!fToEP)/*EtaPiP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1+4*(fToEP)/*DeltaPhi+*/,1/*DeltaPhi-*/,1+4*(fToEP)/*Psi*/};
a8ad4709 236
8c1aafae 237 Int_t binPid2[nBinPid] = {1/*ptKs*/,1+7*(!fToEP)/*EtaPiN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1+4*(fToEP)/*DeltaPhi-*/,1+4*(fToEP)/*Psi*/};
a8ad4709 238
239 fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
240 fContPid->SetTitleX("M_{#Lambda}");
241 fContPid->SetTitleY("centrality (%)");
242 fContPid->SetVarName(0,"p_{T}^{#Lambda}");
243 fContPid->SetVarRange(0,0,10);
244 fContPid->SetVarName(1,"#eta^{#Lambda}");
245 fContPid->SetVarRange(1,-0.8,0.8);
246 fContPid->SetVarName(2,"p_{T}^{Kp}");
247 fContPid->SetVarRange(2,0.3,4.3);
248 fContPid->SetVarName(3,"p_{T}^{Kn}");
249 fContPid->SetVarRange(3,0.3,4.3);
250 fContPid->SetVarName(4,"BayesProb^{Kp}");
251 fContPid->SetVarRange(4,0,1.);
252 fContPid->SetVarName(5,"BayesProb^{Kn}");
253 fContPid->SetVarRange(5,0,1.);
254 fContPid->SetVarName(6,"isTOFmatch^{Kp}");
255 fContPid->SetVarRange(6,-0.5,1.5);
256 fContPid->SetVarName(7,"isTOFmatch^{Kn}");
257 fContPid->SetVarRange(7,-0.5,1.5);
258 fContPid->SetVarName(8,"isLambdaTrue^{Kn}");
259 fContPid->SetVarRange(8,-0.5,1.5);
260 fContPid->SetVarName(9,"N#sigma^{Kp}");
261 fContPid->SetVarRange(9,1.25,6.25);
262 fContPid->SetVarName(10,"N#sigma^{Kn}");
263 fContPid->SetVarRange(10,1.25,6.25);
264 fContPid->SetVarName(11,"#Delta#phi^{Kp}");
265 fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
266 fContPid->SetVarName(12,"#Delta#phi^{Kn}");
267 fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
268 fContPid->SetVarName(13,"#Psi");
269 fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
270
271 fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
272 fContPid2->SetTitleX("M_{#Lambda}");
273 fContPid2->SetTitleY("centrality (%)");
274 fContPid2->SetVarName(0,"p_{T}^{#Lambda}");
275 fContPid2->SetVarRange(0,0,10);
276 fContPid2->SetVarName(1,"#eta^{#Lambda}");
277 fContPid2->SetVarRange(1,-0.8,0.8);
278 fContPid2->SetVarName(2,"p_{T}^{Kp}");
279 fContPid2->SetVarRange(2,0.3,4.3);
280 fContPid2->SetVarName(3,"p_{T}^{Kn}");
281 fContPid2->SetVarRange(3,0.3,4.3);
282 fContPid2->SetVarName(4,"BayesProb^{Kp}");
283 fContPid2->SetVarRange(4,0,1.);
284 fContPid2->SetVarName(5,"BayesProb^{Kn}");
285 fContPid2->SetVarRange(5,0,1.);
286 fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
287 fContPid2->SetVarRange(6,-0.5,1.5);
288 fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
289 fContPid2->SetVarRange(7,-0.5,1.5);
290 fContPid2->SetVarName(8,"isLambdaTrue^{Kn}");
291 fContPid2->SetVarRange(8,-0.5,1.5);
292 fContPid2->SetVarName(9,"N#sigma^{Kp}");
293 fContPid2->SetVarRange(9,1.25,6.25);
294 fContPid2->SetVarName(10,"N#sigma^{Kn}");
295 fContPid2->SetVarRange(10,1.25,6.25);
296 fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
297 fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
298 fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
299 fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
300 fContPid2->SetVarName(13,"#Psi");
301 fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
302
303
304
305 const Int_t nDETsignal = 100; // mass
306 Double_t binDETsignal[nDETsignal+1];
307 for(Int_t i=0;i<nDETsignal+1;i++){
308 binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
309 }
310 const Int_t nDETsignal2 = 10; // centrality
311 Double_t binDETsignal2[nDETsignal2+1];
312 for(Int_t i=0;i<nDETsignal2+1;i++){
313 binDETsignal2[i] = i*100./ nDETsignal2;
314 }
315 fContPid->AddSpecies("Lambda",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
316 fContPid2->AddSpecies("Lambda2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
317
318 fList->Add(fContPid);
319 fList->Add(fContPid2);
320
520899fb 321 const Int_t nBinUser = 6;
322 Int_t binUser[nBinUser] = {8/*Eta*/,20/*pt*/,2/*istrue*/,4/*whatSelection*/,1/*DeltaPhi*/,1/*Psi*/};
323 fContUser = new AliPIDperfContainer("contUserPID",nBinUser,binUser);
324 fContUser->SetTitleX("M_{#Lambda}");
325 fContUser->SetTitleY("centrality (%)");
326 fContUser->SetVarName(0,"#eta");
327 fContUser->SetVarRange(0,-0.8,0.8);
328 fContUser->SetVarName(1,"p_{T}");
329 fContUser->SetVarRange(1,0.3,4.3);
330 fContUser->SetVarName(2,"isLambdaTrue");
331 fContUser->SetVarRange(2,-0.5,1.5);
332 fContUser->SetVarName(3,"whatSelected"); // 0=no, 1=pi, 2=K, 3=p
333 fContUser->SetVarRange(3,-0.5,3.5);
334 fContUser->SetVarName(4,"#Delta#phi");
335 fContUser->SetVarRange(4,-TMath::Pi(),TMath::Pi());
336 fContUser->SetVarName(5,"#Psi");
337 fContUser->SetVarRange(5,-TMath::Pi()/2,TMath::Pi()/2);
338
339 fContUser2 = new AliPIDperfContainer("contUserPID2",nBinUser,binUser);
340 fContUser2->SetTitleX("M_{#Lambda}");
341 fContUser2->SetTitleY("centrality (%)");
342 fContUser2->SetVarName(0,"#eta");
343 fContUser2->SetVarRange(0,-0.8,0.8);
344 fContUser2->SetVarName(1,"p_{T}");
345 fContUser2->SetVarRange(1,0.3,4.3);
346 fContUser2->SetVarName(2,"isLambdaTrue");
347 fContUser2->SetVarRange(2,-0.5,1.5);
348 fContUser2->SetVarName(3,"whatSelected");
349 fContUser2->SetVarRange(3,-0.5,3.5);
350 fContUser2->SetVarName(4,"#Delta#phi");
351 fContUser2->SetVarRange(4,-TMath::Pi(),TMath::Pi());
352 fContUser2->SetVarName(5,"#Psi");
353 fContUser2->SetVarRange(5,-TMath::Pi()/2,TMath::Pi()/2);
354
355 fContUser->AddSpecies("Lambda",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
356 fContUser2->AddSpecies("Lambda2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
357
358 fList->Add(fContUser);
359 fList->Add(fContUser2);
360
a8ad4709 361 hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
362 hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
363 hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
364 hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
365
366 hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
367 hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
368 hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
369 hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
370
371 fList2->Add(hMatching[0]);
372 fList2->Add(hMatching[1]);
373 fList2->Add(hMatching[2]);
374 fList2->Add(hMatching[3]);
375 fList2->Add(hTracking[0]);
376 fList2->Add(hTracking[1]);
377 fList2->Add(hTracking[2]);
378 fList2->Add(hTracking[3]);
379
380 fTOFTPCsignal = new TH2F("hTOFTPCsignal","TOF-TPC signal for protons (1.2 < p_{T} < 1.3 GeV/#it{c}, cent. = 0-20);N#sigma_{TOF};N#sigma_{TPC}",100,-5,5,100,-5,5);
381 fList2->Add(fTOFTPCsignal);
382 fCombsignal = new TH1F("hCombsignal","Combined signal for protons (1.2 < p_{T} < 1.3 GeV/#it{c}, cent. = 0-20);N#sigma",100,0,5);
383 fList2->Add(fCombsignal);
384
385 // QA plots
386 char name[100],title[400];
387 for(Int_t i=0;i < nCentrBin;i++){
388 snprintf(name,100,"hPiTPCQA_%i",i);
389 snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
390 fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
391 fList3->Add(fPiTPC[i]);
392
393 snprintf(name,100,"hKaTPCQA_%i",i);
394 snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
395 fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
396 fList3->Add(fKaTPC[i]);
397
398 snprintf(name,100,"hPrTPCQA_%i",i);
399 snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
400 fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
401 fList3->Add(fPrTPC[i]);
402
403 snprintf(name,100,"hElTPCQA_%i",i);
404 snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
405 fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
406 fList3->Add(fElTPC[i]);
407
408 snprintf(name,100,"hPiTOFQA_%i",i);
409 snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
410 fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
411 fList3->Add(fPiTOF[i]);
412
413 snprintf(name,100,"hKaTOFQA_%i",i);
414 snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
415 fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
416 fList3->Add(fKaTOF[i]);
417
418 snprintf(name,100,"hPrTOFQA_%i",i);
419 snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
420 fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
421 fList3->Add(fPrTOF[i]);
422
423 snprintf(name,100,"hElTOFQA_%i",i);
424 snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
425 fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
426 fList3->Add(fElTOF[i]);
427
428 }
429
430 // Post output data.
431 PostData(1, fList);
432 PostData(2, fList2);
433 PostData(3, fList3);
434}
435
436//______________________________________________________________________________
437void AliAnalysisTaskLambdaBayes::UserExec(Option_t *)
438{
439 // Main loop
440 // Called for each event
441
442 fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
443 if(!fOutputAOD){
444 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
445 this->Dump();
446 return;
447 }
448
449 Float_t zvtx = GetVertex(fOutputAOD);
450
451
452
453 //Get the MC object
454 if(fIsMC){
455 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
456 if (!mcHeader) {
457 AliError("Could not find MC Header in AOD");
458 return;
459 }
460 }
461
462 if (TMath::Abs(zvtx) < fVtxCut) {
463
464 SelectLambda();
465
466 //Centrality
467 Float_t v0Centr = -10.;
468 Float_t trkCentr = -10.;
469 AliCentrality *centrality = fOutputAOD->GetCentrality();
470 if (centrality){
471 trkCentr = centrality->GetCentralityPercentile("V0M");
472 v0Centr = centrality->GetCentralityPercentile("TRK");
473 }
474
cd71f9f4 475 if(!fTypeCol){
476 v0Centr=100./(fOutputAOD->GetNumberOfTracks()/12.+1);
477 trkCentr=v0Centr;
478 }
479
9c668341 480 if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
a8ad4709 481 fCentrality = v0Centr;
482 Analyze(fOutputAOD); // Do analysis!!!!
483
484 }
485 }
486
487}
488
489//________________________________________________________________________
490void AliAnalysisTaskLambdaBayes::Analyze(AliAODEvent* aodEvent)
491{
492
493 Int_t ntrack = aodEvent->GetNumberOfTracks();
494
495 fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
496 fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
497 fPidKp=0,fPidKn=0;
498 fMassV0=-1;
499
500 TLorentzVector ksV;
501
502 Double_t px,py,pz,E;
503
504 Float_t invMmin = 1.115-0.005*8;
505 Float_t invMmax = 1.115+0.005*8;
506
507 Int_t icentr = 8;
508 if(fCentrality < 0) icentr = 8;
509 else if(fCentrality < 10) icentr = 0;
510 else if(fCentrality < 20) icentr = 1;
511 else if(fCentrality < 30) icentr = 2;
512 else if(fCentrality < 40) icentr = 3;
513 else if(fCentrality < 50) icentr = 4;
514 else if(fCentrality < 60) icentr = 5;
515 else if(fCentrality < 70) icentr = 6;
516 else if(fCentrality < 80) icentr = 7;
517
518 Float_t addMismatchForMC = 0.005;
faa6d35f 519 if(fCentrality < 50) addMismatchForMC += 0.005;
a8ad4709 520 if(fCentrality < 20) addMismatchForMC += 0.02;
521
faa6d35f 522 if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
9c668341 523
a8ad4709 524 fPsi = 0;
525 /* Compute TPC EP */
526 Double_t Qx2 = 0, Qy2 = 0;
527 Double_t Qx3 = 0, Qy3 = 0;
528 for(Int_t iT = 0; iT < ntrack; iT++) {
529 AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
530
531 if (!aodTrack){
532 continue;
533 }
534
535 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
536
537 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
538 continue;
539
540 Double_t b[2] = {-99., -99.};
541 Double_t bCov[3] = {-99., -99., -99.};
542 if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
543 continue;
544
545 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
546 continue;
547
548 Qx2 += TMath::Cos(2*aodTrack->Phi());
549 Qy2 += TMath::Sin(2*aodTrack->Phi());
550 Qx3 += TMath::Cos(3*aodTrack->Phi());
551 Qy3 += TMath::Sin(3*aodTrack->Phi());
552
553 }
554
555 fPsi = TMath::ATan2(Qy2, Qx2)/2.;
556
557 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
558 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
559 AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
560
cd71f9f4 561// PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
562// PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
563// PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
564// PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
a8ad4709 565
566 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
567
f4c4ac9e 568
569 Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
570 Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
a8ad4709 571 Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
bc629e38 572 Double_t nSigmaTPCRef,nSigmaTOFRef=6,nSigmaCombRef;
a8ad4709 573
574 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
575 TClonesArray *mcArray = NULL;
576 if (mcHeader)
577 mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
578
579 Int_t nmc = 0;
580 if(mcArray)
581 nmc = mcArray->GetEntries();
582
583 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
584 AliAODTrack* track = aodEvent->GetTrack(i);
585
586 AliAODMCParticle *mcp = NULL;
587 Int_t pdg = 0;
588
589 if (!track){
590 continue;
591 }
592
593 Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
594
595 Int_t label = -1;
596 if(mcArray){
597 label = track->GetLabel();
598 if(label != -1 && label < nmc){
599 label = TMath::Abs(label);
600 mcp = (AliAODMCParticle*)mcArray->At(label);
601 pdg = TMath::Abs(mcp->GetPdgCode());
602 }
603 else
604 label = -1;
605 }
606 else{
607 /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
608 }
609
610 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
611 hTracking[0]->Fill(track->Pt(),fCentrality);
612 if(pdg == 211)
613 hTracking[1]->Fill(track->Pt(),fCentrality);
614 else if(pdg == 321)
615 hTracking[2]->Fill(track->Pt(),fCentrality);
616 else if(pdg == 2212)
617 hTracking[3]->Fill(track->Pt(),fCentrality);
618 else if(! mcp){ // Fill matching histo with the prob
619 hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
620 hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
621 hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
622 }
623 }
624
625 if(!tofMatch) continue;
626
627 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
628 hMatching[0]->Fill(track->Pt(),fCentrality);
629 if(pdg == 211)
630 hMatching[1]->Fill(track->Pt(),fCentrality);
631 else if(pdg == 321)
632 hMatching[2]->Fill(track->Pt(),fCentrality);
633 else if(pdg == 2212)
634 hMatching[3]->Fill(track->Pt(),fCentrality);
635 else if(! mcp){ // Fill matching histo with the prob
636 hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
637 hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
638 hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
639 }
640 }
641 }
642
007df888 643// Int_t pdg1 = -1;
644// Int_t pdg2 = -1;
a8ad4709 645
646
647 // start analysis Lambda
648 for(Int_t i=0;i < ntrack;i++){ // loop on proton candidate tracks
649 AliAODTrack* KpTrack = aodEvent->GetTrack(i);
650
651 if (!KpTrack){
652 continue;
653 }
654
655 if(!(KpTrack->Pt() > 0.3 && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
656
657 Bool_t isLambda = (KpTrack->Charge() > 0);
658
659 nSigmaComb=5;
bc629e38 660 nSigmaCombRef=5;
e34b28fe 661 nSigmaTOF=5;
bc629e38 662 nSigmaTOFRef=5;
a8ad4709 663 fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
664 fPidKp=0;
665
666 UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
667
668 Double_t oldpP[10];
669 fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
670
671
a8ad4709 672 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
bc629e38 673 nSigmaTPCRef = PIDResponse->NumberOfSigmasTPC(KpTrack,(AliPID::EParticleType) fSpeciesRef);
a8ad4709 674
675 if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
676
677 Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
9c668341 678 /*
a8ad4709 679 if(mcArray){
680 Int_t labelK = TMath::Abs(KpTrack->GetLabel());
681 AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 682 pdg1 = TMath::Abs(mcp1->GetPdgCode());
a8ad4709 683 }
9c668341 684 */
a8ad4709 685
686 fPidKp = Int_t(probP[4]*100);
687
688 if(tofMatch1){
689 if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
690 // remove this tof hit
691 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
692 detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
693 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
694 fPidKp = Int_t(probP[4]*100);
695 tofMatch1=0;
696 }
697 else{
698 if(probP[4] > probP[3] && probP[4] > probP[2] && probP[4] > probP[0]) fPidKp += 128; // max prob
699
cd71f9f4 700 if(isLambda){
701 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
702 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
703 if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
704 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
705 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
706 if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
707 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
708 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
709 if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
710 }
a8ad4709 711 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
bc629e38 712 nSigmaTOFRef = PIDResponse->NumberOfSigmasTOF(KpTrack,(AliPID::EParticleType) fSpeciesRef);
cd71f9f4 713 if(isLambda){
714 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
715 if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
716 }
a8ad4709 717
a8ad4709 718 if(fIsMC){
719 Float_t mismAdd = addMismatchForMC;
720 if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
721
722 if(gRandom->Rndm() < mismAdd){
723 Float_t etaAbs = TMath::Abs(KpTrack->Eta());
724 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
725 channel = channel % 8736;
726 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
727
728 // generate random time
a16898fb 729 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
a8ad4709 730 Double_t times[10];
731 KpTrack->GetIntegratedTimes(times);
bc629e38 732 nSigmaTOF = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[4], AliPID::kProton);
733 nSigmaTOFRef = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[fSpeciesRef], (AliPID::EParticleType) fSpeciesRef);
a8ad4709 734 }
735 }
736
737 if(fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
63d1d548 738 nSigmaTOF = TMath::Abs(nSigmaTOF);
a8ad4709 739
740 if(nSigmaTOF < 2) fPidKp += 256;
741 else if(nSigmaTOF < 3) fPidKp += 512;
742 }
743 }
744
745 if(tofMatch1){
746 nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
bc629e38 747 nSigmaCombRef = TMath::Sqrt(nSigmaTOFRef*nSigmaTOFRef + nSigmaTPCRef*nSigmaTPCRef);
a8ad4709 748 if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2){
749 fCombsignal->Fill(nSigmaComb);
750 }
751 } else {
752 nSigmaComb = TMath::Abs(nSigmaTPC);
bc629e38 753 nSigmaCombRef = TMath::Abs(nSigmaTPCRef);
a8ad4709 754 }
755
756 // use sigmaTOF instead of sigmaComb
bc629e38 757 if(tofMatch1){
758 nSigmaComb = nSigmaTOF;
759 nSigmaCombRef = nSigmaTOFRef;
760 }
a8ad4709 761
762 if(nSigmaComb < 2) nSigmaComb = 2;
763 else if(nSigmaComb < 3) nSigmaComb = 3;
764 else if(nSigmaComb < 5) nSigmaComb = 4.99;
765 else nSigmaComb = 6;
766
bc629e38 767 if(nSigmaCombRef < 2) nSigmaCombRef = 2;
768 else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
769 else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
770 else nSigmaCombRef = 6;
771
a8ad4709 772 Int_t iks=-1;
773 for(Int_t k=0;k < fNLambda;k++){ // find the Lambda which contains the positive track
774 if(i == fIpP[k]) iks = k;
775 }
776
777 if(fPtKp > 4.299) fPtKp = 4.299;
778
779 if(iks > -1 && fIpN[iks] > -1){
780 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
781 Int_t j = fIpN[iks];
782 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
783
784 if (!KnTrack){
785 continue;
786 }
787
e34b28fe 788 nSigmaComb2 = 5;
789 nSigmaTOF2 = 5;
790
a8ad4709 791 if(!(KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
792
793 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
794 fPidKn=0;
795
796 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
797 Double_t oldpN[10];
798 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
799
800 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
bc629e38 801
a8ad4709 802 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
803
a8ad4709 804 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 805 /*
a8ad4709 806 if(mcArray){
807 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
808 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 809 pdg2 = TMath::Abs(mcp2->GetPdgCode());
810 }
811 */
a8ad4709 812
813 fPidKn = Int_t(probN[2]*100);
814
815 if(tofMatch2){
816 if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
817
818 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
819
bc629e38 820 if(fIsMC){
821 Float_t mismAdd = addMismatchForMC;
822 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
823
824 if(gRandom->Rndm() < mismAdd){
825 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
826 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
827 channel = channel % 8736;
828 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
829
830 // generate random time
831 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
832 Double_t times[10];
833 KnTrack->GetIntegratedTimes(times);
834 nSigmaTOF2 = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[4], AliPID::kProton);
835 }
836 }
837
838
839
a8ad4709 840 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
841
842 if(nSigmaTOF2 < 2) fPidKn += 256;
843 else if(nSigmaTOF2 < 3) fPidKn += 512;
844 }
845
846 px = KpTrack->Px() + KnTrack->Px();
847 py = KpTrack->Py() + KnTrack->Py();
848 pz = KpTrack->Pz() + KnTrack->Pz();
849 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 0.938e-01*0.938e-01);
850 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
851
852 ksV.SetPxPyPzE(px,py,pz,E);
853 fMassV0 = fMassLambda[iks];
854
855 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
856
857
858 fPtLambdaC = ksV.Pt();
859 fEtaLambda = ksV.Eta();
860 fPhiLambdaC = ksV.Phi();
861
862 if(tofMatch2){
863 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
864 } else {
865 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
866 }
867
868 // use sigmaTOF instead of sigmaComb
e34b28fe 869 if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
a8ad4709 870
871 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
872 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
873 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
874 else nSigmaComb2 = 6;
875
876 Bool_t isTrue = kFALSE;
877
878 if(mcArray){
879 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
880 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
881
882 if(labelP > -1 && labelN > -1){
883 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
884 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
885
886 Int_t mP = partP->GetMother();
887 Int_t mN = partN->GetMother();
888 if(mP == mN && mP > -1){
889 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
890 Int_t pdgM = partM->GetPdgCode();
891 if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
892 }
893 }
894
895 }
896
897 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
898 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
899
f8f53ae0 900 if(gRandom->Rndm() < 0.5){
901 deltaphi1 += TMath::Pi();
902 deltaphi2 += TMath::Pi();
903 }
904
a8ad4709 905 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
906 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
907 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
908 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
909
910 if(fPtKn > 4.299) fPtKn = 4.299;
911
912 Float_t xTOfill[] = {fPtLambdaC,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
913 Float_t xTOfill2[] = {fPtLambdaC,KpTrack->Eta(),fPtKn,fPtKp,(fPidKn%128)*0.01,(fPidKp%128)*0.01,tofMatch2,tofMatch1,isTrue,nSigmaComb2,nSigmaComb,deltaphi2,deltaphi1,fPsi};
914
915
916 Int_t ipt = 0;
917 while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
918 ipt++;
919 }
920 ipt--;
921 if(ipt < 0) ipt = 0; // just to be sure
922
923 if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
bc629e38 924 if(fSpeciesRef != 4){
925 xTOfill[4] = probP[fSpeciesRef];
926 xTOfill2[5] = probP[fSpeciesRef];
927
928 xTOfill[9] = nSigmaCombRef;
929 xTOfill2[10] = nSigmaCombRef;
930
931 }
932
a8ad4709 933 if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
934 else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
520899fb 935
936
937
938 if(fPIDuserCut){
939 Float_t xUser[] = {KpTrack->Eta(),fPtKp,isTrue,0,deltaphi1,fPsi};
940
941 if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled
942 xUser[3] = 1;
943 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
944 xUser[3] = 2;
945 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
946 xUser[3] = 3;
947 }
948
949 if(isLambda) fContUser->Fill(0,fMassV0,fCentrality,xUser);
950 else fContUser2->Fill(0,fMassV0,fCentrality,xUser);
951 }
a8ad4709 952 }
953
954
955 }
956 } // end analysi Lambda
957}
958
959//_____________________________________________________________________________
960Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
961{
962
963 Float_t zvtx = -999;
964
965 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
966 if (!vtxAOD)
967 return zvtx;
968 if(vtxAOD->GetNContributors()>0)
969 zvtx = vtxAOD->GetZ();
970
971 return zvtx;
972}
973//_____________________________________________________________________________
974void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
975{
976 // Terminate loop
977 Printf("Terminate()");
978}
979//=======================================================================
980void AliAnalysisTaskLambdaBayes::SelectLambda(){
981 fNLambda=0;
982 fNpPos=0;
983 fNpNeg=0;
984
985 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
986 AliAODv0 *myV0;
987 Double_t dMASS=0.0;
988 for (Int_t i=0; i!=nV0s; ++i) {
989 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
990 if(!myV0) continue;
991 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
992 Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
993 if(pass) {
994 if(pass==1) dMASS = myV0->MassLambda();
995 else dMASS = myV0->MassAntiLambda();
996
f4c4ac9e 997 Float_t massKs = myV0->MassK0Short();
998
999 if(TMath::Abs(dMASS-1.115)/0.005 < 8 && TMath::Abs(massKs - 0.497)/0.005 > 8){
a8ad4709 1000 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1001 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1002 if(iT->Charge()<0){
1003 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1004 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1005 }
1006 fPhiLambda[fNLambda] = myV0->Phi();
1007 fPtLambda[fNLambda] = myV0->Pt();
1008 fIpP[fNLambda] = FindDaugheterIndex(iT);
1009 fIpN[fNLambda] = FindDaugheterIndex(jT);
1010
1011 if(pass==2){
1012 Int_t itemp = fIpP[fNLambda];
1013 fIpP[fNLambda] = fIpN[fNLambda];
1014 fIpN[fNLambda] = itemp;
1015 }
1016
1017 fMassLambda[fNLambda] = dMASS;
1018 if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
1019 fNLambda++;
1020 }
1021
1022 }
1023 }
1024 }
1025}
1026
1027//=======================================================================
1028Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1029{
1030 if (myV0->GetOnFlyStatus() ) return 0;
1031
1032 //the following is needed in order to evualuate track-quality
1033 AliAODTrack *iT, *jT;
1034 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1035 Double_t pos[3],cov[6];
1036 vV0s->GetXYZ(pos);
1037 vV0s->GetCovarianceMatrix(cov);
1038 const AliESDVertex vESD(pos,cov,100.,100);
1039
1040 // TESTING CHARGE
1041 int iPos, iNeg;
1042 iT=(AliAODTrack*) myV0->GetDaughter(0);
1043 if(iT->Charge()>0) {
1044 iPos = 0; iNeg = 1;
1045 } else {
1046 iPos = 1; iNeg = 0;
1047 }
1048 // END OF TEST
1049
1050 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1051
1052 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1053
1054 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1055 if(!trkFlag) return 0;
1056 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1057 if(!trkFlag2) return 0;
1058
1059 Double_t pvertex[3];
1060 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1061 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1062 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1063
1064 Double_t dDL=myV0->DecayLengthV0( pvertex );
1065 if(dDL < 0.5 || dDL > 25) return 0;
1066
1067 Double_t dDCA=myV0->DcaV0Daughters();
1068 if(dDCA >0.5) return 0;
1069
1070 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1071 if(dCTP < -1) return 0;
1072
1073// AliESDtrack ieT( iT );
1074// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1075// if(TMath::Abs(dD0P) < 0]) return 0;
1076
1077// AliESDtrack jeT( jT );
1078// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1079// if(TMath::Abs(dD0M) < 0) return 0;
1080
1081// Double_t dD0D0=dD0P*dD0M;
1082// if(dD0D0>0) return 0;
1083
1084// Double_t dETA=myV0->Eta();
1085// if(dETA <-0.8) return 0;
1086// if(dETA >0.8) return 0;
1087
1088// Double_t dQT=myV0->PtArmV0();
1089// if(specie==0) if(dQT<???) return 0;
1090
1091 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1092 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1093
1094 if(specie==1 && dALPHA<0) return 2; // antilambda
1095 return 1; // Ks or lambda
1096}
1097//-------------------------------------------------------------------------
1098Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
1099 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1100
1101 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1102 AliAODTrack* track = fOutputAOD->GetTrack(i);
1103 if(track == trk) return i;
1104 }
1105
1106 printf("Daughter for %p not found\n",trk);
1107 return -1;
1108}
1109//-------------------------------------------------------------------------
1110Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
1111 if(!fIsMC) return 1; // used only on MC
1112
9c668341 1113 if( fTypeCol==2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 1114 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1115
1116 if(!(channel%20)) return 0; // 5% additional loss in MC
1117 }
1118
1119 return 1;
1120}