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