Added the possibility to the user to manage own PID
[u/mrichter/AliRoot.git] / PWGPP / pid / AliAnalysisTaskK0sBayes.cxx
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
31 Float_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
32 Float_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
38 ClassImp(AliAnalysisTaskK0sBayes)
39 //_____________________________________________________________________________
40 AliAnalysisTaskK0sBayes::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),
72   fContUser(NULL),
73   fContUser2(NULL),
74   fNK0s(0),
75   fNpiPos(0),
76   fNpiNeg(0),
77   fHmismTOF(0),
78   fHchannelTOFdistr(0),
79   fTypeCol(2),
80   fPIDuserCut(NULL)
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"); 
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   }
120 }
121
122 //______________________________________________________________________________
123 AliAnalysisTaskK0sBayes::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),
155   fContUser(NULL),
156   fContUser2(NULL),
157   fNK0s(0),
158   fNpiPos(0),
159   fNpiNeg(0),
160   fHmismTOF(0),
161   fHchannelTOFdistr(0),
162   fTypeCol(2),
163   fPIDuserCut(NULL)
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"); 
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   }
209 }
210 //_____________________________________________________________________________
211 AliAnalysisTaskK0sBayes::~AliAnalysisTaskK0sBayes()
212 {
213
214 }
215
216 //______________________________________________________________________________
217 void 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
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
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
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 //______________________________________________________________________________
429 void 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
467       if(!fTypeCol){
468         v0Centr=100./(fOutputAOD->GetNumberOfTracks()/12.+1);
469         trkCentr=v0Centr;
470       }
471
472       if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
473         fCentrality = v0Centr;
474         Analyze(fOutputAOD); // Do analysis!!!!
475
476       }
477     }
478     
479 }
480
481 //________________________________________________________________________
482 void 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;
511   if(fCentrality < 50) addMismatchForMC += 0.005;
512   if(fCentrality < 20) addMismatchForMC += 0.02;
513
514   if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
515
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
553   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
554
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.};
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   
629 //   Int_t pdg1 = -1;
630 //   Int_t pdg2 = -1;
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;
644     nSigmaTOF = 5;
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
653     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
654
655     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
656
657     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
658
659     /*
660     if(mcArray){
661       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
662       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
663       pdg1 = TMath::Abs(mcp1->GetPdgCode());
664     }
665     */
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);
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));
684         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
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));
687         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
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));
690         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
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));
693                         
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
705             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
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);
713         nSigmaTOF = TMath::Abs(nSigmaTOF);
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
730     if(tofMatch1) nSigmaComb = nSigmaTOF;
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;
767       nSigmaTOF2=5;
768
769       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
770       /*
771       if(mcArray){
772         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
773         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
774         pdg2 = TMath::Abs(mcp2->GetPdgCode());
775       }
776       */
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);
810               nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
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
847       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
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);
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
922       }
923
924
925     }
926   } // end analysi K0s
927 }
928
929 //_____________________________________________________________________________
930 Float_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 //_____________________________________________________________________________
944 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
945
946   // Terminate loop
947   Printf("Terminate()");
948 }
949 //=======================================================================
950 void 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();
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){
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 //=======================================================================
1071 Int_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 //-------------------------------------------------------------------------
1141 Int_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 //-------------------------------------------------------------------------
1153 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1154   if(!fIsMC) return 1; // used only on MC
1155
1156   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
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 }