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