]> 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
RelVal: return exitcode of Makeflow to the caller
[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;
fc9b31a7 715 Double_t times[AliPID::kSPECIESC];
716 KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
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
d82096df 742 nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
743
9dd528f5 744 if(tofMatch1){
745 nSigmaComb = nSigmaTOF;
746 nSigmaCombRef = nSigmaTOFRef;
747 }
a8ad4709 748
749 if(nSigmaComb < 2) nSigmaComb = 2;
750 else if(nSigmaComb < 3) nSigmaComb = 3;
751 else if(nSigmaComb < 5) nSigmaComb = 4.99;
752 else nSigmaComb = 6;
753
9dd528f5 754 if(nSigmaCombRef < 2) nSigmaCombRef = 2;
755 else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
756 else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
757 else nSigmaCombRef = 6;
758
a8ad4709 759 Int_t iks=-1;
760 for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
761 if(i == fIpiP[k]) iks = k;
762 }
763
764 if(fPtKp > 4.299) fPtKp = 4.299;
765
766 if(iks > -1 && fIpiN[iks] > -1){
767 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
768 Int_t j = fIpiN[iks];
769 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
770
771 if (!KnTrack){
772 continue;
773 }
774
775 if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
776
777 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
778 fPidKn=0;
779
780 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
781 Double_t oldpN[10];
782 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
783
784 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
9dd528f5 785 nSigmaTPC2Ref = PIDResponse->NumberOfSigmasTPC(KnTrack,(AliPID::EParticleType) fSpeciesRef);
a8ad4709 786
787 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
788
789 nSigmaComb2=5;
e34b28fe 790 nSigmaTOF2=5;
9dd528f5 791 nSigmaComb2Ref=5;
792 nSigmaTOF2Ref=5;
a8ad4709 793
794 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 795 /*
a8ad4709 796 if(mcArray){
797 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
798 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 799 pdg2 = TMath::Abs(mcp2->GetPdgCode());
800 }
801 */
a8ad4709 802
803 fPidKn = Int_t(probN[2]*100);
804
805 if(tofMatch2){
806 if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
807 // remove this tof hit
808 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
809 detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
810 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
811 fPidKn = Int_t(probN[4]*100);
812 tofMatch2=0;
813 }
814 else{
815 if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
816
9dd528f5 817 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
a8ad4709 818 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
9dd528f5 819 nSigmaTOF2Ref = PIDResponse->NumberOfSigmasTOF(KnTrack,(AliPID::EParticleType) fSpeciesRef);
820 nSigmaTOF2Ref = TMath::Abs(nSigmaTOF2Ref);
a8ad4709 821
822 if(fIsMC){
823 Float_t mismAdd = addMismatchForMC;
824 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
825
826 if(gRandom->Rndm() < mismAdd){
827 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
828 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
829 channel = channel % 8736;
830 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
831
832 // generate random time
833 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
fc9b31a7 834 Double_t times[AliPID::kSPECIESC];
835 KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
9dd528f5 836 nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[2], AliPID::kPion);
837 nSigmaTOF2Ref = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[fSpeciesRef], fSpeciesRef);
a8ad4709 838 }
839 }
840
841 if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
842
843 if(nSigmaTOF2 < 2) fPidKn += 256;
844 else if(nSigmaTOF2 < 3) fPidKn += 512;
845 }
846 }
847
848 px = KpTrack->Px() + KnTrack->Px();
849 py = KpTrack->Py() + KnTrack->Py();
850 pz = KpTrack->Pz() + KnTrack->Pz();
851 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
852 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
853
854 ksV.SetPxPyPzE(px,py,pz,E);
855 fMassV0 = fMassKs[iks];
856
857 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
858
859
860 fPtKs = ksV.Pt();
861 fEtaKs = ksV.Eta();
862 fPhiKs = ksV.Phi();
863
864 if(tofMatch2){
865 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
9dd528f5 866 nSigmaComb2Ref = TMath::Sqrt(nSigmaTOF2Ref*nSigmaTOF2Ref+ nSigmaTPC2Ref*nSigmaTPC2Ref);
a8ad4709 867 if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
868 fCombsignal->Fill(nSigmaComb2);
869 }
870 } else {
871 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
9dd528f5 872 nSigmaComb2Ref = TMath::Abs(nSigmaTPC2Ref);
a8ad4709 873 }
874
875 // use sigmaTOF instead of sigmaComb
9dd528f5 876 if(tofMatch2){
877 nSigmaComb2 = nSigmaTOF2;
878 nSigmaComb2Ref = nSigmaTOF2Ref;
879 }
a8ad4709 880
881 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
882 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
883 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
884 else nSigmaComb2 = 6;
9dd528f5 885 if(nSigmaComb2Ref < 2) nSigmaComb2Ref = 2;
886 else if(nSigmaComb2Ref < 3) nSigmaComb2Ref = 3;
887 else if(nSigmaComb2Ref < 5) nSigmaComb2Ref = 4.99;
888 else nSigmaComb2Ref = 6;
a8ad4709 889
890 Bool_t isTrue = kFALSE;
891
892 if(mcArray){
893 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
894 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
895
896 if(labelP > -1 && labelN > -1){
897 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
898 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
899
900 Int_t mP = partP->GetMother();
901 Int_t mN = partN->GetMother();
902 if(mP == mN && mP > -1){
903 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
904 Int_t pdgM = partM->GetPdgCode();
905 if(pdgM == 310) isTrue = kTRUE;
906 }
907 }
908
909 }
910
911 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
912 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
913
f8f53ae0 914 if(gRandom->Rndm() < 0.5){
915 deltaphi1 += TMath::Pi();
916 deltaphi2 += TMath::Pi();
917 }
918
a8ad4709 919 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
920 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
921 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
922 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
923
924 if(fPtKn > 4.299) fPtKn = 4.299;
925
2942f542 926 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)};
927 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 928
a8ad4709 929 Int_t ipt = 0;
930 while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
931 ipt++;
932 }
933 ipt--;
934 if(ipt < 0) ipt = 0; // just to be sure
935
936 if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
9dd528f5 937 if(fSpeciesRef != 2){
938 xTOfill[4] = probP[fSpeciesRef];
939 xTOfill2[5] = probN[fSpeciesRef];
940
941 xTOfill[9] = nSigmaCombRef;
942 xTOfill2[10] = nSigmaComb2Ref;
943
944 }
945
a8ad4709 946 fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
947 xTOfill[1] = KnTrack->Eta();
9dd528f5 948 fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
520899fb 949
950 if(fPIDuserCut){
2942f542 951 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)};
952 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 953
954 if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
955 xUser[3] = 1;
956 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
957 xUser[3] = 2;
958 } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
959 xUser[3] = 3;
960 }
961 if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
962 xUser2[3] = 1;
963 } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
964 xUser2[3] = 2;
965 } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
966 xUser2[3] = 3;
967 }
968 fContUser->Fill(0,fMassV0,fCentrality,xUser);
969 fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
970 }
971
a8ad4709 972 }
973
974
975 }
976 } // end analysi K0s
977}
978
979//_____________________________________________________________________________
980Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
981{
982
983 Float_t zvtx = -999;
984
985 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
986 if (!vtxAOD)
987 return zvtx;
988 if(vtxAOD->GetNContributors()>0)
989 zvtx = vtxAOD->GetZ();
990
991 return zvtx;
992}
993//_____________________________________________________________________________
994void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
995{
996 // Terminate loop
997 Printf("Terminate()");
998}
999//=======================================================================
1000void AliAnalysisTaskK0sBayes::SelectK0s(){
1001 fNK0s=0;
1002 fNpiPos=0;
1003 fNpiNeg=0;
1004
1005 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1006 AliAODv0 *myV0;
1007 Double_t dMASS=0.0;
1008 for (Int_t i=0; i!=nV0s; ++i) {
1009 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1010 if(!myV0) continue;
1011 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
1012 Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
1013 if(pass) {
1014 dMASS = myV0->MassK0Short();
f4c4ac9e 1015 Float_t massLambda = myV0->MassLambda();
1016 Float_t massAntiLambda = myV0->MassAntiLambda();
1017
1018 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 1019 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1020 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1021 if(iT->Charge()<0){
1022 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1023 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1024 }
1025 fPhiK0s[fNK0s] = myV0->Phi();
1026 fPtK0s[fNK0s] = myV0->Pt();
1027 fIpiP[fNK0s] = FindDaugheterIndex(iT);
1028 fIpiN[fNK0s] = FindDaugheterIndex(jT);
1029 fMassKs[fNK0s] = dMASS;
1030 if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
1031 fNK0s++;
1032 }
1033 }
1034 }
1035
1036 /* My V0 code
1037 // fill pion stacks
1038 Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
1039 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1040 AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
1041
1042 if (!aodTrack){
1043 continue;
1044 }
1045
1046 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
1047
1048 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1049 continue;
1050 }
1051
1052 Double_t b[2] = {-99., -99.};
1053 Double_t bCov[3] = {-99., -99., -99.};
1054 if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1055 continue;
1056
1057 if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
1058
1059
1060 Int_t charge = aodTrack->Charge();
1061 if(charge > 0){
1062 fIPiPos[fNpiPos] = iT;
1063 fNpiPos++;
1064 }
1065 else{
1066 fIPiNeg[fNpiNeg] = iT;
1067 fNpiNeg++;
1068 }
1069 }
1070
1071 for(Int_t i=0;i < fNpiPos;i++){
1072 AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
1073 AliESDtrack pipE(pip);
1074
1075 for(Int_t j=0;j < fNpiNeg;j++){
1076 AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
1077 AliESDtrack pinE(pin);
1078
1079 Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
1080
1081 Double_t pPos[3];
1082 Double_t pNeg[3];
1083 pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
1084 pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
1085
1086 Float_t length = (xp+xn)*0.5;
1087
1088 Float_t pxs = pPos[0] + pNeg[0];
1089 Float_t pys = pPos[1] + pNeg[1];
1090 Float_t pzs = pPos[2] + pNeg[2];
1091 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);
1092
1093 Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
1094 Float_t phi = TMath::ATan2(pys,pxs);
1095 Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
1096
1097 // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
1098
1099 if(mindist < 0.4&& length > 0.7 && length < 25){
1100 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);
1101 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);
1102
1103 Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
1104 Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
1105
1106 if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
1107 fPhiK0s[fNK0s] = phi;
1108 fPtK0s[fNK0s] = pt;
1109 fIpiP[fNK0s] =fIPiPos[i] ;
1110 fIpiN[fNK0s] = fIPiNeg[j];
1111 fMassKs[fNK0s] = mass;
1112 fNK0s++;
1113 }
1114 }
1115 }
1116 }
1117 */
1118}
1119
1120//=======================================================================
1121Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1122{
1123 if (myV0->GetOnFlyStatus() ) return 0;
1124
1125 //the following is needed in order to evualuate track-quality
1126 AliAODTrack *iT, *jT;
1127 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1128 Double_t pos[3],cov[6];
1129 vV0s->GetXYZ(pos);
1130 vV0s->GetCovarianceMatrix(cov);
1131 const AliESDVertex vESD(pos,cov,100.,100);
1132
1133 // TESTING CHARGE
1134 int iPos, iNeg;
1135 iT=(AliAODTrack*) myV0->GetDaughter(0);
1136 if(iT->Charge()>0) {
1137 iPos = 0; iNeg = 1;
1138 } else {
1139 iPos = 1; iNeg = 0;
1140 }
1141 // END OF TEST
1142
1143 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1144
1145 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1146
1147 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1148 if(!trkFlag) return 0;
1149 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1150 if(!trkFlag2) return 0;
1151
1152 Double_t pvertex[3];
1153 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1154 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1155 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1156
1157 Double_t dDL=myV0->DecayLengthV0( pvertex );
1158 if(dDL < 0.5 || dDL > 25) return 0;
1159
1160 Double_t dDCA=myV0->DcaV0Daughters();
1161 if(dDCA >0.5) return 0;
1162
1163 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1164 if(dCTP < -1) return 0;
1165
1166// AliESDtrack ieT( iT );
1167// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1168// if(TMath::Abs(dD0P) < 0]) return 0;
1169
1170// AliESDtrack jeT( jT );
1171// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1172// if(TMath::Abs(dD0M) < 0) return 0;
1173
1174// Double_t dD0D0=dD0P*dD0M;
1175// if(dD0D0>0) return 0;
1176
1177// Double_t dETA=myV0->Eta();
1178// if(dETA <-0.8) return 0;
1179// if(dETA >0.8) return 0;
1180
1181// Double_t dQT=myV0->PtArmV0();
1182// if(specie==0) if(dQT<???) return 0;
1183
1184 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1185 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1186
1187 if(specie==1 && dALPHA<0) return 2; // antilambda
1188 return 1; // K0s or lambda
1189}
1190//-------------------------------------------------------------------------
1191Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1192 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1193
1194 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1195 AliAODTrack* track = fOutputAOD->GetTrack(i);
1196 if(track == trk) return i;
1197 }
1198
1199 printf("Daughter for %p not found\n",trk);
1200 return -1;
1201}
1202//-------------------------------------------------------------------------
1203Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1204 if(!fIsMC) return 1; // used only on MC
1205
9c668341 1206 if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 1207 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1208
1209 if(!(channel%20)) return 0; // 5% additional loss in MC
1210 }
1211
1212 return 1;
1213}