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