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