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