]> 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
PWGPP-3, ATO-19 - Reduce size of the space points for alignment (Laser calibration...
[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;
fc9b31a7 730 Double_t times[AliPID::kSPECIESC];
731 KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
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
d82096df 757 nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
758
bc629e38 759 if(tofMatch1){
760 nSigmaComb = nSigmaTOF;
761 nSigmaCombRef = nSigmaTOFRef;
762 }
a8ad4709 763
764 if(nSigmaComb < 2) nSigmaComb = 2;
765 else if(nSigmaComb < 3) nSigmaComb = 3;
766 else if(nSigmaComb < 5) nSigmaComb = 4.99;
767 else nSigmaComb = 6;
768
bc629e38 769 if(nSigmaCombRef < 2) nSigmaCombRef = 2;
770 else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
771 else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
772 else nSigmaCombRef = 6;
773
a8ad4709 774 Int_t iks=-1;
775 for(Int_t k=0;k < fNLambda;k++){ // find the Lambda which contains the positive track
776 if(i == fIpP[k]) iks = k;
777 }
778
779 if(fPtKp > 4.299) fPtKp = 4.299;
780
781 if(iks > -1 && fIpN[iks] > -1){
782 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
783 Int_t j = fIpN[iks];
784 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
785
786 if (!KnTrack){
787 continue;
788 }
789
e34b28fe 790 nSigmaComb2 = 5;
791 nSigmaTOF2 = 5;
792
a8ad4709 793 if(!(KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
794
795 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
796 fPidKn=0;
797
798 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
799 Double_t oldpN[10];
800 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
801
802 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
bc629e38 803
a8ad4709 804 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
805
a8ad4709 806 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 807 /*
a8ad4709 808 if(mcArray){
809 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
810 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 811 pdg2 = TMath::Abs(mcp2->GetPdgCode());
812 }
813 */
a8ad4709 814
815 fPidKn = Int_t(probN[2]*100);
816
817 if(tofMatch2){
818 if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
819
820 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
821
bc629e38 822 if(fIsMC){
823 Float_t mismAdd = addMismatchForMC;
824 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
825
826 if(gRandom->Rndm() < mismAdd){
827 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
828 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
829 channel = channel % 8736;
830 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
831
832 // generate random time
833 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
fc9b31a7 834 Double_t times[AliPID::kSPECIESC];
835 KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
bc629e38 836 nSigmaTOF2 = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[4], AliPID::kProton);
837 }
838 }
839
840
841
a8ad4709 842 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
843
844 if(nSigmaTOF2 < 2) fPidKn += 256;
845 else if(nSigmaTOF2 < 3) fPidKn += 512;
846 }
847
848 px = KpTrack->Px() + KnTrack->Px();
849 py = KpTrack->Py() + KnTrack->Py();
850 pz = KpTrack->Pz() + KnTrack->Pz();
851 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 0.938e-01*0.938e-01);
852 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
853
854 ksV.SetPxPyPzE(px,py,pz,E);
855 fMassV0 = fMassLambda[iks];
856
857 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
858
859
860 fPtLambdaC = ksV.Pt();
861 fEtaLambda = ksV.Eta();
862 fPhiLambdaC = ksV.Phi();
863
864 if(tofMatch2){
865 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
866 } else {
867 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
868 }
869
870 // use sigmaTOF instead of sigmaComb
e34b28fe 871 if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
a8ad4709 872
873 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
874 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
875 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
876 else nSigmaComb2 = 6;
877
878 Bool_t isTrue = kFALSE;
879
880 if(mcArray){
881 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
882 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
883
884 if(labelP > -1 && labelN > -1){
885 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
886 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
887
888 Int_t mP = partP->GetMother();
889 Int_t mN = partN->GetMother();
890 if(mP == mN && mP > -1){
891 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
892 Int_t pdgM = partM->GetPdgCode();
893 if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
894 }
895 }
896
897 }
898
899 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
900 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
901
f8f53ae0 902 if(gRandom->Rndm() < 0.5){
903 deltaphi1 += TMath::Pi();
904 deltaphi2 += TMath::Pi();
905 }
906
a8ad4709 907 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
908 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
909 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
910 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
911
912 if(fPtKn > 4.299) fPtKn = 4.299;
913
2942f542 914 Float_t xTOfill[] = {static_cast<Float_t>(fPtLambdaC),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn,(fPidKp%128)*0.01),static_cast<Float_t>((fPidKn%128)*0.01),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
915 Float_t xTOfill2[] = {static_cast<Float_t>(fPtLambdaC),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKn),static_cast<Float_t>(fPtKp),static_cast<Float_t>((fPidKn%128)*0.01),static_cast<Float_t>((fPidKp%128)*0.01),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(fPsi)};
a8ad4709 916
917
918 Int_t ipt = 0;
919 while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
920 ipt++;
921 }
922 ipt--;
923 if(ipt < 0) ipt = 0; // just to be sure
924
925 if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
bc629e38 926 if(fSpeciesRef != 4){
927 xTOfill[4] = probP[fSpeciesRef];
928 xTOfill2[5] = probP[fSpeciesRef];
929
930 xTOfill[9] = nSigmaCombRef;
931 xTOfill2[10] = nSigmaCombRef;
932
933 }
934
a8ad4709 935 if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
936 else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
520899fb 937
938
939
940 if(fPIDuserCut){
2942f542 941 Float_t xUser[] = {static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(isTrue),0,static_cast<Float_t>(deltaphi1),static_cast<Float_t>(fPsi)};
520899fb 942
943 if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled
944 xUser[3] = 1;
945 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
946 xUser[3] = 2;
947 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
948 xUser[3] = 3;
949 }
950
951 if(isLambda) fContUser->Fill(0,fMassV0,fCentrality,xUser);
952 else fContUser2->Fill(0,fMassV0,fCentrality,xUser);
953 }
a8ad4709 954 }
955
956
957 }
958 } // end analysi Lambda
959}
960
961//_____________________________________________________________________________
962Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
963{
964
965 Float_t zvtx = -999;
966
967 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
968 if (!vtxAOD)
969 return zvtx;
970 if(vtxAOD->GetNContributors()>0)
971 zvtx = vtxAOD->GetZ();
972
973 return zvtx;
974}
975//_____________________________________________________________________________
976void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
977{
978 // Terminate loop
979 Printf("Terminate()");
980}
981//=======================================================================
982void AliAnalysisTaskLambdaBayes::SelectLambda(){
983 fNLambda=0;
984 fNpPos=0;
985 fNpNeg=0;
986
987 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
988 AliAODv0 *myV0;
989 Double_t dMASS=0.0;
990 for (Int_t i=0; i!=nV0s; ++i) {
991 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
992 if(!myV0) continue;
993 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
994 Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
995 if(pass) {
996 if(pass==1) dMASS = myV0->MassLambda();
997 else dMASS = myV0->MassAntiLambda();
998
f4c4ac9e 999 Float_t massKs = myV0->MassK0Short();
1000
1001 if(TMath::Abs(dMASS-1.115)/0.005 < 8 && TMath::Abs(massKs - 0.497)/0.005 > 8){
a8ad4709 1002 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1003 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1004 if(iT->Charge()<0){
1005 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1006 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1007 }
1008 fPhiLambda[fNLambda] = myV0->Phi();
1009 fPtLambda[fNLambda] = myV0->Pt();
1010 fIpP[fNLambda] = FindDaugheterIndex(iT);
1011 fIpN[fNLambda] = FindDaugheterIndex(jT);
1012
1013 if(pass==2){
1014 Int_t itemp = fIpP[fNLambda];
1015 fIpP[fNLambda] = fIpN[fNLambda];
1016 fIpN[fNLambda] = itemp;
1017 }
1018
1019 fMassLambda[fNLambda] = dMASS;
1020 if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
1021 fNLambda++;
1022 }
1023
1024 }
1025 }
1026 }
1027}
1028
1029//=======================================================================
1030Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1031{
1032 if (myV0->GetOnFlyStatus() ) return 0;
1033
1034 //the following is needed in order to evualuate track-quality
1035 AliAODTrack *iT, *jT;
1036 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1037 Double_t pos[3],cov[6];
1038 vV0s->GetXYZ(pos);
1039 vV0s->GetCovarianceMatrix(cov);
1040 const AliESDVertex vESD(pos,cov,100.,100);
1041
1042 // TESTING CHARGE
1043 int iPos, iNeg;
1044 iT=(AliAODTrack*) myV0->GetDaughter(0);
1045 if(iT->Charge()>0) {
1046 iPos = 0; iNeg = 1;
1047 } else {
1048 iPos = 1; iNeg = 0;
1049 }
1050 // END OF TEST
1051
1052 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1053
1054 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1055
1056 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1057 if(!trkFlag) return 0;
1058 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1059 if(!trkFlag2) return 0;
1060
1061 Double_t pvertex[3];
1062 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1063 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1064 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1065
1066 Double_t dDL=myV0->DecayLengthV0( pvertex );
1067 if(dDL < 0.5 || dDL > 25) return 0;
1068
1069 Double_t dDCA=myV0->DcaV0Daughters();
1070 if(dDCA >0.5) return 0;
1071
1072 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1073 if(dCTP < -1) return 0;
1074
1075// AliESDtrack ieT( iT );
1076// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1077// if(TMath::Abs(dD0P) < 0]) return 0;
1078
1079// AliESDtrack jeT( jT );
1080// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1081// if(TMath::Abs(dD0M) < 0) return 0;
1082
1083// Double_t dD0D0=dD0P*dD0M;
1084// if(dD0D0>0) return 0;
1085
1086// Double_t dETA=myV0->Eta();
1087// if(dETA <-0.8) return 0;
1088// if(dETA >0.8) return 0;
1089
1090// Double_t dQT=myV0->PtArmV0();
1091// if(specie==0) if(dQT<???) return 0;
1092
1093 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1094 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1095
1096 if(specie==1 && dALPHA<0) return 2; // antilambda
1097 return 1; // Ks or lambda
1098}
1099//-------------------------------------------------------------------------
1100Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
1101 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1102
1103 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1104 AliAODTrack* track = fOutputAOD->GetTrack(i);
1105 if(track == trk) return i;
1106 }
1107
1108 printf("Daughter for %p not found\n",trk);
1109 return -1;
1110}
1111//-------------------------------------------------------------------------
1112Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
1113 if(!fIsMC) return 1; // used only on MC
1114
9c668341 1115 if( fTypeCol==2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 1116 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1117
1118 if(!(channel%20)) return 0; // 5% additional loss in MC
1119 }
1120
1121 return 1;
1122}