]> 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
Force use of T0-TOF for V0 analyses: PID perf
[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
55022e85 557 PIDResponse->SetTOFResponse(aodEvent,AliPIDResponse::kTOF_T0);
558
a8ad4709 559 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
560
f4c4ac9e 561 Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
562 Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
a8ad4709 563 Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
9dd528f5 564 Double_t nSigmaTPCRef,nSigmaTOFRef=6,nSigmaTPC2Ref,nSigmaTOF2Ref=6,nSigmaCombRef,nSigmaComb2Ref;
a8ad4709 565
566
567 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
568 TClonesArray *mcArray = NULL;
569 if (mcHeader)
570 mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
571
572 Int_t nmc = 0;
573 if(mcArray)
574 nmc = mcArray->GetEntries();
575
576 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
577 AliAODTrack* track = aodEvent->GetTrack(i);
578
579 AliAODMCParticle *mcp = NULL;
580 Int_t pdg = 0;
581
582 if (!track){
583 continue;
584 }
585
586 Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
587
588 Int_t label = -1;
589 if(mcArray){
590 label = track->GetLabel();
591 if(label != -1 && label < nmc){
592 label = TMath::Abs(label);
593 mcp = (AliAODMCParticle*)mcArray->At(label);
594 pdg = TMath::Abs(mcp->GetPdgCode());
595 }
596 else
597 label = -1;
598 }
599 else{
600 /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
601 }
602
603 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
604 hTracking[0]->Fill(track->Pt(),fCentrality);
605 if(pdg == 211)
606 hTracking[1]->Fill(track->Pt(),fCentrality);
607 else if(pdg == 321)
608 hTracking[2]->Fill(track->Pt(),fCentrality);
609 else if(pdg == 2212)
610 hTracking[3]->Fill(track->Pt(),fCentrality);
611 else if(! mcp){ // Fill matching histo with the prob
612 hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
613 hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
614 hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
615 }
616 }
617
618 if(!tofMatch) continue;
619
620 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
621 hMatching[0]->Fill(track->Pt(),fCentrality);
622 if(pdg == 211)
623 hMatching[1]->Fill(track->Pt(),fCentrality);
624 else if(pdg == 321)
625 hMatching[2]->Fill(track->Pt(),fCentrality);
626 else if(pdg == 2212)
627 hMatching[3]->Fill(track->Pt(),fCentrality);
628 else if(! mcp){ // Fill matching histo with the prob
629 hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
630 hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
631 hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
632 }
633 }
634 }
635
007df888 636// Int_t pdg1 = -1;
637// Int_t pdg2 = -1;
a8ad4709 638
639
640 // start analysis K0s
641 for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
642 AliAODTrack* KpTrack = aodEvent->GetTrack(i);
643
644 if (!KpTrack){
645 continue;
646 }
647
648 if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3 && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
649
650 nSigmaComb=5;
9dd528f5 651 nSigmaCombRef=5;
e34b28fe 652 nSigmaTOF = 5;
9dd528f5 653 nSigmaTOFRef = 5;
a8ad4709 654 fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
655 fPidKp=0;
656
657 UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
658
659 Double_t oldpP[10];
660 fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
661
a8ad4709 662 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
9dd528f5 663 nSigmaTPCRef = PIDResponse->NumberOfSigmasTPC(KpTrack,(AliPID::EParticleType) fSpeciesRef);
a8ad4709 664
665 if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
666
667 Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
668
9c668341 669 /*
a8ad4709 670 if(mcArray){
671 Int_t labelK = TMath::Abs(KpTrack->GetLabel());
672 AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 673 pdg1 = TMath::Abs(mcp1->GetPdgCode());
a8ad4709 674 }
9c668341 675 */
a8ad4709 676
677 fPidKp = Int_t(probP[2]*100);
678
679 if(tofMatch1){
680 if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
681 // remove this tof hit
682 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
683 detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
684 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
685 fPidKp = Int_t(probP[4]*100);
686 tofMatch1=0;
687 }
688 else{
689 if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
690
691 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
cd71f9f4 692 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
693 if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
a8ad4709 694 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
cd71f9f4 695 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
696 if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
a8ad4709 697 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
cd71f9f4 698 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
699 if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
a8ad4709 700 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
cd71f9f4 701 if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
702 if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
9dd528f5 703 nSigmaTOFRef = PIDResponse->NumberOfSigmasTOF(KpTrack,(AliPID::EParticleType) fSpeciesRef);
704
a8ad4709 705 if(fIsMC){
706 Float_t mismAdd = addMismatchForMC;
707 if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
708
709 if(gRandom->Rndm() < mismAdd){
710 Float_t etaAbs = TMath::Abs(KpTrack->Eta());
711 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
712 channel = channel % 8736;
713 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
714
715 // generate random time
a16898fb 716 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
fc9b31a7 717 Double_t times[AliPID::kSPECIESC];
718 KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
9dd528f5 719 nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[2], AliPID::kPion);
720 nSigmaTOFRef = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[fSpeciesRef], fSpeciesRef);
a8ad4709 721 }
722 }
723
724 if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
63d1d548 725 nSigmaTOF = TMath::Abs(nSigmaTOF);
a8ad4709 726
727 if(nSigmaTOF < 2) fPidKp += 256;
728 else if(nSigmaTOF < 3) fPidKp += 512;
729 }
730 }
731
732 if(tofMatch1){
733 nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
9dd528f5 734 nSigmaCombRef = TMath::Sqrt(nSigmaTOFRef*nSigmaTOFRef + nSigmaTPCRef*nSigmaTPCRef);
a8ad4709 735 if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
736 fCombsignal->Fill(nSigmaComb);
737 }
738 } else {
739 nSigmaComb = TMath::Abs(nSigmaTPC);
9dd528f5 740 nSigmaCombRef = TMath::Abs(nSigmaTPCRef);
a8ad4709 741 }
742
743 // use sigmaTOF instead of sigmaComb
d82096df 744 nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
745
9dd528f5 746 if(tofMatch1){
747 nSigmaComb = nSigmaTOF;
748 nSigmaCombRef = nSigmaTOFRef;
749 }
a8ad4709 750
751 if(nSigmaComb < 2) nSigmaComb = 2;
752 else if(nSigmaComb < 3) nSigmaComb = 3;
753 else if(nSigmaComb < 5) nSigmaComb = 4.99;
754 else nSigmaComb = 6;
755
9dd528f5 756 if(nSigmaCombRef < 2) nSigmaCombRef = 2;
757 else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
758 else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
759 else nSigmaCombRef = 6;
760
a8ad4709 761 Int_t iks=-1;
762 for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
763 if(i == fIpiP[k]) iks = k;
764 }
765
766 if(fPtKp > 4.299) fPtKp = 4.299;
767
768 if(iks > -1 && fIpiN[iks] > -1){
769 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
770 Int_t j = fIpiN[iks];
771 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
772
773 if (!KnTrack){
774 continue;
775 }
776
777 if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
778
779 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
780 fPidKn=0;
781
782 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
783 Double_t oldpN[10];
784 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
785
786 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
9dd528f5 787 nSigmaTPC2Ref = PIDResponse->NumberOfSigmasTPC(KnTrack,(AliPID::EParticleType) fSpeciesRef);
a8ad4709 788
789 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
790
791 nSigmaComb2=5;
e34b28fe 792 nSigmaTOF2=5;
9dd528f5 793 nSigmaComb2Ref=5;
794 nSigmaTOF2Ref=5;
a8ad4709 795
796 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 797 /*
a8ad4709 798 if(mcArray){
799 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
800 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 801 pdg2 = TMath::Abs(mcp2->GetPdgCode());
802 }
803 */
a8ad4709 804
805 fPidKn = Int_t(probN[2]*100);
806
807 if(tofMatch2){
808 if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
809 // remove this tof hit
810 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
811 detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
812 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
813 fPidKn = Int_t(probN[4]*100);
814 tofMatch2=0;
815 }
816 else{
817 if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
818
9dd528f5 819 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
a8ad4709 820 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
9dd528f5 821 nSigmaTOF2Ref = PIDResponse->NumberOfSigmasTOF(KnTrack,(AliPID::EParticleType) fSpeciesRef);
822 nSigmaTOF2Ref = TMath::Abs(nSigmaTOF2Ref);
a8ad4709 823
824 if(fIsMC){
825 Float_t mismAdd = addMismatchForMC;
826 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
827
828 if(gRandom->Rndm() < mismAdd){
829 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
830 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
831 channel = channel % 8736;
832 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
833
834 // generate random time
835 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
fc9b31a7 836 Double_t times[AliPID::kSPECIESC];
837 KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
9dd528f5 838 nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[2], AliPID::kPion);
839 nSigmaTOF2Ref = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[fSpeciesRef], fSpeciesRef);
a8ad4709 840 }
841 }
842
843 if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
844
845 if(nSigmaTOF2 < 2) fPidKn += 256;
846 else if(nSigmaTOF2 < 3) fPidKn += 512;
847 }
848 }
849
850 px = KpTrack->Px() + KnTrack->Px();
851 py = KpTrack->Py() + KnTrack->Py();
852 pz = KpTrack->Pz() + KnTrack->Pz();
853 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
854 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
855
856 ksV.SetPxPyPzE(px,py,pz,E);
857 fMassV0 = fMassKs[iks];
858
859 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
860
861
862 fPtKs = ksV.Pt();
863 fEtaKs = ksV.Eta();
864 fPhiKs = ksV.Phi();
865
866 if(tofMatch2){
867 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
9dd528f5 868 nSigmaComb2Ref = TMath::Sqrt(nSigmaTOF2Ref*nSigmaTOF2Ref+ nSigmaTPC2Ref*nSigmaTPC2Ref);
a8ad4709 869 if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
870 fCombsignal->Fill(nSigmaComb2);
871 }
872 } else {
873 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
9dd528f5 874 nSigmaComb2Ref = TMath::Abs(nSigmaTPC2Ref);
a8ad4709 875 }
876
877 // use sigmaTOF instead of sigmaComb
9dd528f5 878 if(tofMatch2){
879 nSigmaComb2 = nSigmaTOF2;
880 nSigmaComb2Ref = nSigmaTOF2Ref;
881 }
a8ad4709 882
883 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
884 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
885 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
886 else nSigmaComb2 = 6;
9dd528f5 887 if(nSigmaComb2Ref < 2) nSigmaComb2Ref = 2;
888 else if(nSigmaComb2Ref < 3) nSigmaComb2Ref = 3;
889 else if(nSigmaComb2Ref < 5) nSigmaComb2Ref = 4.99;
890 else nSigmaComb2Ref = 6;
a8ad4709 891
892 Bool_t isTrue = kFALSE;
893
894 if(mcArray){
895 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
896 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
897
898 if(labelP > -1 && labelN > -1){
899 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
900 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
901
902 Int_t mP = partP->GetMother();
903 Int_t mN = partN->GetMother();
904 if(mP == mN && mP > -1){
905 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
906 Int_t pdgM = partM->GetPdgCode();
907 if(pdgM == 310) isTrue = kTRUE;
908 }
909 }
910
911 }
912
913 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
914 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
915
f8f53ae0 916 if(gRandom->Rndm() < 0.5){
917 deltaphi1 += TMath::Pi();
918 deltaphi2 += TMath::Pi();
919 }
920
a8ad4709 921 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
922 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
923 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
924 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
925
926 if(fPtKn > 4.299) fPtKn = 4.299;
927
2942f542 928 Float_t xTOfill[] = {static_cast<Float_t>(fPtKs),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn),static_cast<Float_t>(probP[2]),static_cast<Float_t>(probN[2]),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
929 Float_t xTOfill2[] = {static_cast<Float_t>(fPtKs),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn),static_cast<Float_t>(probP[2]),static_cast<Float_t>(probN[2]),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
a8ad4709 930
a8ad4709 931 Int_t ipt = 0;
932 while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
933 ipt++;
934 }
935 ipt--;
936 if(ipt < 0) ipt = 0; // just to be sure
937
938 if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
9dd528f5 939 if(fSpeciesRef != 2){
940 xTOfill[4] = probP[fSpeciesRef];
941 xTOfill2[5] = probN[fSpeciesRef];
942
943 xTOfill[9] = nSigmaCombRef;
944 xTOfill2[10] = nSigmaComb2Ref;
945
946 }
947
a8ad4709 948 fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
949 xTOfill[1] = KnTrack->Eta();
9dd528f5 950 fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
520899fb 951
952 if(fPIDuserCut){
2942f542 953 Float_t xUser[] = {static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(isTrue),0,static_cast<Float_t>(deltaphi1),static_cast<Float_t>(fPsi)};
954 Float_t xUser2[] = {static_cast<Float_t>(KnTrack->Eta()),static_cast<Float_t>(fPtKn),static_cast<Float_t>(isTrue),0,static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
520899fb 955
956 if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
957 xUser[3] = 1;
958 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
959 xUser[3] = 2;
960 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
961 xUser[3] = 3;
962 }
963 if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
964 xUser2[3] = 1;
965 } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
966 xUser2[3] = 2;
967 } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
968 xUser2[3] = 3;
969 }
970 fContUser->Fill(0,fMassV0,fCentrality,xUser);
971 fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
972 }
973
a8ad4709 974 }
975
976
977 }
978 } // end analysi K0s
979}
980
981//_____________________________________________________________________________
982Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
983{
984
985 Float_t zvtx = -999;
986
987 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
988 if (!vtxAOD)
989 return zvtx;
990 if(vtxAOD->GetNContributors()>0)
991 zvtx = vtxAOD->GetZ();
992
993 return zvtx;
994}
995//_____________________________________________________________________________
996void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
997{
998 // Terminate loop
999 Printf("Terminate()");
1000}
1001//=======================================================================
1002void AliAnalysisTaskK0sBayes::SelectK0s(){
1003 fNK0s=0;
1004 fNpiPos=0;
1005 fNpiNeg=0;
1006
1007 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1008 AliAODv0 *myV0;
1009 Double_t dMASS=0.0;
1010 for (Int_t i=0; i!=nV0s; ++i) {
1011 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1012 if(!myV0) continue;
1013 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
1014 Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
1015 if(pass) {
1016 dMASS = myV0->MassK0Short();
f4c4ac9e 1017 Float_t massLambda = myV0->MassLambda();
1018 Float_t massAntiLambda = myV0->MassAntiLambda();
1019
1020 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 1021 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1022 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1023 if(iT->Charge()<0){
1024 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1025 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1026 }
1027 fPhiK0s[fNK0s] = myV0->Phi();
1028 fPtK0s[fNK0s] = myV0->Pt();
1029 fIpiP[fNK0s] = FindDaugheterIndex(iT);
1030 fIpiN[fNK0s] = FindDaugheterIndex(jT);
1031 fMassKs[fNK0s] = dMASS;
1032 if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
1033 fNK0s++;
1034 }
1035 }
1036 }
1037
1038 /* My V0 code
1039 // fill pion stacks
1040 Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
1041 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1042 AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
1043
1044 if (!aodTrack){
1045 continue;
1046 }
1047
1048 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
1049
1050 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1051 continue;
1052 }
1053
1054 Double_t b[2] = {-99., -99.};
1055 Double_t bCov[3] = {-99., -99., -99.};
1056 if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1057 continue;
1058
1059 if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
1060
1061
1062 Int_t charge = aodTrack->Charge();
1063 if(charge > 0){
1064 fIPiPos[fNpiPos] = iT;
1065 fNpiPos++;
1066 }
1067 else{
1068 fIPiNeg[fNpiNeg] = iT;
1069 fNpiNeg++;
1070 }
1071 }
1072
1073 for(Int_t i=0;i < fNpiPos;i++){
1074 AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
1075 AliESDtrack pipE(pip);
1076
1077 for(Int_t j=0;j < fNpiNeg;j++){
1078 AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
1079 AliESDtrack pinE(pin);
1080
1081 Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
1082
1083 Double_t pPos[3];
1084 Double_t pNeg[3];
1085 pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
1086 pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
1087
1088 Float_t length = (xp+xn)*0.5;
1089
1090 Float_t pxs = pPos[0] + pNeg[0];
1091 Float_t pys = pPos[1] + pNeg[1];
1092 Float_t pzs = pPos[2] + pNeg[2];
1093 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);
1094
1095 Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
1096 Float_t phi = TMath::ATan2(pys,pxs);
1097 Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
1098
1099 // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
1100
1101 if(mindist < 0.4&& length > 0.7 && length < 25){
1102 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);
1103 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);
1104
1105 Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
1106 Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
1107
1108 if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
1109 fPhiK0s[fNK0s] = phi;
1110 fPtK0s[fNK0s] = pt;
1111 fIpiP[fNK0s] =fIPiPos[i] ;
1112 fIpiN[fNK0s] = fIPiNeg[j];
1113 fMassKs[fNK0s] = mass;
1114 fNK0s++;
1115 }
1116 }
1117 }
1118 }
1119 */
1120}
1121
1122//=======================================================================
1123Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1124{
1125 if (myV0->GetOnFlyStatus() ) return 0;
1126
1127 //the following is needed in order to evualuate track-quality
1128 AliAODTrack *iT, *jT;
1129 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1130 Double_t pos[3],cov[6];
1131 vV0s->GetXYZ(pos);
1132 vV0s->GetCovarianceMatrix(cov);
1133 const AliESDVertex vESD(pos,cov,100.,100);
1134
1135 // TESTING CHARGE
1136 int iPos, iNeg;
1137 iT=(AliAODTrack*) myV0->GetDaughter(0);
1138 if(iT->Charge()>0) {
1139 iPos = 0; iNeg = 1;
1140 } else {
1141 iPos = 1; iNeg = 0;
1142 }
1143 // END OF TEST
1144
1145 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1146
1147 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1148
1149 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1150 if(!trkFlag) return 0;
1151 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1152 if(!trkFlag2) return 0;
1153
1154 Double_t pvertex[3];
1155 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1156 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1157 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1158
1159 Double_t dDL=myV0->DecayLengthV0( pvertex );
1160 if(dDL < 0.5 || dDL > 25) return 0;
1161
1162 Double_t dDCA=myV0->DcaV0Daughters();
1163 if(dDCA >0.5) return 0;
1164
1165 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1166 if(dCTP < -1) return 0;
1167
1168// AliESDtrack ieT( iT );
1169// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1170// if(TMath::Abs(dD0P) < 0]) return 0;
1171
1172// AliESDtrack jeT( jT );
1173// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1174// if(TMath::Abs(dD0M) < 0) return 0;
1175
1176// Double_t dD0D0=dD0P*dD0M;
1177// if(dD0D0>0) return 0;
1178
1179// Double_t dETA=myV0->Eta();
1180// if(dETA <-0.8) return 0;
1181// if(dETA >0.8) return 0;
1182
1183// Double_t dQT=myV0->PtArmV0();
1184// if(specie==0) if(dQT<???) return 0;
1185
1186 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1187 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1188
1189 if(specie==1 && dALPHA<0) return 2; // antilambda
1190 return 1; // K0s or lambda
1191}
1192//-------------------------------------------------------------------------
1193Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1194 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1195
1196 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1197 AliAODTrack* track = fOutputAOD->GetTrack(i);
1198 if(track == trk) return i;
1199 }
1200
1201 printf("Daughter for %p not found\n",trk);
1202 return -1;
1203}
1204//-------------------------------------------------------------------------
1205Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1206 if(!fIsMC) return 1; // used only on MC
1207
9c668341 1208 if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 1209 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1210
1211 if(!(channel%20)) return 0; // 5% additional loss in MC
1212 }
1213
1214 return 1;
1215}