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