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