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