Added the possibility to the user to manage own PID
[u/mrichter/AliRoot.git] / PWGPP / pid / AliAnalysisTaskLambdaBayes.cxx
1 #include "AliAnalysisTaskLambdaBayes.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 AliAnalysisTaskLambdaBayes::fPtLambdaMin[AliAnalysisTaskLambdaBayes::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 AliAnalysisTaskLambdaBayes::fPtLambdaMax[AliAnalysisTaskLambdaBayes::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(AliAnalysisTaskLambdaBayes)
39 //_____________________________________________________________________________
40 AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes():
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   fPtLambdaC(0.),
55   fPhiLambdaC(0.),
56   fEtaLambda(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   fNLambda(0),
75   fNpPos(0),
76   fNpNeg(0),
77   fHmismTOF(0),
78   fHchannelTOFdistr(0),
79   fTypeCol(2),
80   fPIDuserCut(NULL)
81 {
82   // Default constructor (should not be used)
83   fList->SetName("contLambdaBayes1");
84   fList2->SetName("contLambdaBayes2");
85   fList3->SetName("contLambdaBayes3");
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
112   for(Int_t i=0;i < 1000;i++){
113     fPhiLambda[i] = 0.0;
114     fPtLambda[i] = 0.0;
115     fIPPos[i] = 0;
116     fIPNeg[i] = 0;
117     fIpP[i] = 0;
118     fIpN[i] = 0;
119     fMassLambda[i] = 0.0;
120   }
121 }
122
123 //______________________________________________________________________________
124 AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes(const char *name):
125   AliAnalysisTaskSE(name),
126   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
127   fEtaCut(0.8),   // cut on |eta| < fEtaCut
128   fMinPt(0.15),   // cut on pt > fMinPt
129   fIsMC(kFALSE),
130   fQAsw(kFALSE),
131   fNcluster(70),
132   fFilterBit(1),
133   fList(new TList()),
134   fList2(new TList()),
135   fList3(new TList()),
136   fCentrality(-1),
137   fPsi(0),
138   fPtLambdaC(0.),
139   fPhiLambdaC(0.),
140   fEtaLambda(0.),
141   fMassV0(0.),
142   fPtKp(0.),
143   fPhiKp(0.),
144   fEtaKp(0.),
145   fPtKn(0.),
146   fPhiKn(0.),
147   fEtaKn(0.),
148   fPidKp(0),
149   fPidKn(0),
150   fTOFTPCsignal(0),
151   fCombsignal(0),
152   fCutsDaughter(NULL),
153   fPIDCombined(NULL),
154   fContPid(NULL),
155   fContPid2(NULL),
156   fContUser(NULL),
157   fContUser2(NULL),
158   fNLambda(0),
159   fNpPos(0),
160   fNpNeg(0),
161   fHmismTOF(0),
162   fHchannelTOFdistr(0),
163   fTypeCol(2),
164   fPIDuserCut(NULL)
165 {
166
167   DefineOutput(1, TList::Class());
168   DefineOutput(2, TList::Class());
169   DefineOutput(3, TList::Class());
170   DefineOutput(4, TList::Class());
171
172   // Output slot #1 writes into a TTree
173   fList->SetName("contLambdaBayes1");
174   fList2->SetName("contLambdaBayes2");
175   fList3->SetName("contLambdaBayes3");
176
177   fList->SetOwner(kTRUE); 
178   fList2->SetOwner(kTRUE); 
179   fList3->SetOwner(kTRUE); 
180
181   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
182   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
183
184   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
185   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
186
187   for(Int_t i=0;i < nCentrBin;i++){
188     fElTOF[i] = NULL; 
189     fElTPC[i] = NULL; 
190     fPiTOF[i] = NULL; 
191     fPiTPC[i] = NULL; 
192     fKaTOF[i] = NULL; 
193     fKaTPC[i] = NULL; 
194     fPrTOF[i] = NULL; 
195     fPrTPC[i] = NULL; 
196   }
197   for(Int_t i=0;i < 4;i++){
198     hMatching[i] = NULL;
199     hTracking[i] = NULL;
200   }
201
202   for(Int_t i=0;i < 1000;i++){
203     fPhiLambda[i] = 0.0;
204     fPtLambda[i] = 0.0;
205     fIPPos[i] = 0;
206     fIPNeg[i] = 0;
207     fIpP[i] = 0;
208     fIpN[i] = 0;
209     fMassLambda[i] = 0.0;
210   }
211 }
212 //_____________________________________________________________________________
213 AliAnalysisTaskLambdaBayes::~AliAnalysisTaskLambdaBayes()
214 {
215
216 }
217
218 //______________________________________________________________________________
219 void AliAnalysisTaskLambdaBayes::UserCreateOutputObjects()
220 {
221
222   fPIDCombined=new AliPIDCombined;
223   fPIDCombined->SetDefaultTPCPriors();
224   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
225
226   Float_t invMmin = 1.115-0.005*8;
227   Float_t invMmax = 1.115+0.005*8;
228
229   const Int_t nBinPid = 14;
230
231   Int_t binPid[nBinPid] = {1/*ptLambda*/,8/*EtaPP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
232
233   Int_t binPid2[nBinPid] = {1/*ptLambda*/,8/*EtaPN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
234
235   fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
236   fContPid->SetTitleX("M_{#Lambda}");
237   fContPid->SetTitleY("centrality (%)");
238   fContPid->SetVarName(0,"p_{T}^{#Lambda}");
239   fContPid->SetVarRange(0,0,10);
240   fContPid->SetVarName(1,"#eta^{#Lambda}");
241   fContPid->SetVarRange(1,-0.8,0.8);
242   fContPid->SetVarName(2,"p_{T}^{Kp}");
243   fContPid->SetVarRange(2,0.3,4.3);
244   fContPid->SetVarName(3,"p_{T}^{Kn}");
245   fContPid->SetVarRange(3,0.3,4.3);
246   fContPid->SetVarName(4,"BayesProb^{Kp}");
247   fContPid->SetVarRange(4,0,1.);
248   fContPid->SetVarName(5,"BayesProb^{Kn}");
249   fContPid->SetVarRange(5,0,1.);
250   fContPid->SetVarName(6,"isTOFmatch^{Kp}");
251   fContPid->SetVarRange(6,-0.5,1.5);
252   fContPid->SetVarName(7,"isTOFmatch^{Kn}");
253   fContPid->SetVarRange(7,-0.5,1.5);
254   fContPid->SetVarName(8,"isLambdaTrue^{Kn}");
255   fContPid->SetVarRange(8,-0.5,1.5);
256   fContPid->SetVarName(9,"N#sigma^{Kp}");
257   fContPid->SetVarRange(9,1.25,6.25);
258   fContPid->SetVarName(10,"N#sigma^{Kn}");
259   fContPid->SetVarRange(10,1.25,6.25);
260   fContPid->SetVarName(11,"#Delta#phi^{Kp}");
261   fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
262   fContPid->SetVarName(12,"#Delta#phi^{Kn}");
263   fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
264   fContPid->SetVarName(13,"#Psi");
265   fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
266
267   fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
268   fContPid2->SetTitleX("M_{#Lambda}");
269   fContPid2->SetTitleY("centrality (%)");
270   fContPid2->SetVarName(0,"p_{T}^{#Lambda}");
271   fContPid2->SetVarRange(0,0,10);
272   fContPid2->SetVarName(1,"#eta^{#Lambda}");
273   fContPid2->SetVarRange(1,-0.8,0.8);
274   fContPid2->SetVarName(2,"p_{T}^{Kp}");
275   fContPid2->SetVarRange(2,0.3,4.3);
276   fContPid2->SetVarName(3,"p_{T}^{Kn}");
277   fContPid2->SetVarRange(3,0.3,4.3);
278   fContPid2->SetVarName(4,"BayesProb^{Kp}");
279   fContPid2->SetVarRange(4,0,1.);
280   fContPid2->SetVarName(5,"BayesProb^{Kn}");
281   fContPid2->SetVarRange(5,0,1.);
282   fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
283   fContPid2->SetVarRange(6,-0.5,1.5);
284   fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
285   fContPid2->SetVarRange(7,-0.5,1.5);
286   fContPid2->SetVarName(8,"isLambdaTrue^{Kn}");
287   fContPid2->SetVarRange(8,-0.5,1.5);
288   fContPid2->SetVarName(9,"N#sigma^{Kp}");
289   fContPid2->SetVarRange(9,1.25,6.25);
290   fContPid2->SetVarName(10,"N#sigma^{Kn}");
291   fContPid2->SetVarRange(10,1.25,6.25);
292   fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
293   fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
294   fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
295   fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
296   fContPid2->SetVarName(13,"#Psi");
297   fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
298
299
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("Lambda",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
312   fContPid2->AddSpecies("Lambda2",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_{#Lambda}");
321   fContUser->SetTitleY("centrality (%)");
322   fContUser->SetVarName(0,"#eta");
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,"isLambdaTrue");
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_{#Lambda}");
337   fContUser2->SetTitleY("centrality (%)");
338   fContUser2->SetVarName(0,"#eta");
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,"isLambdaTrue");
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("Lambda",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
352   fContUser2->AddSpecies("Lambda2",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 protons (1.2 < p_{T} < 1.3 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 protons (1.2 < p_{T} < 1.3 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 AliAnalysisTaskLambdaBayes::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       SelectLambda();
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 AliAnalysisTaskLambdaBayes::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 = 1.115-0.005*8;
501   Float_t invMmax = 1.115+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 //   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
558 //   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
559 //   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
560 //   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
561
562   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
563
564
565   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
566   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
567   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
568
569   
570   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
571   TClonesArray *mcArray = NULL;
572   if (mcHeader)
573     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
574
575   Int_t nmc = 0;
576   if(mcArray)
577     nmc = mcArray->GetEntries();
578
579   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
580     AliAODTrack* track = aodEvent->GetTrack(i);
581     
582     AliAODMCParticle *mcp = NULL;
583     Int_t pdg = 0;
584     
585     if (!track){
586       continue;
587     }
588     
589     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
590     
591     Int_t label = -1;
592     if(mcArray){
593       label = track->GetLabel();
594       if(label != -1 && label < nmc){
595         label = TMath::Abs(label);
596         mcp = (AliAODMCParticle*)mcArray->At(label);
597         pdg = TMath::Abs(mcp->GetPdgCode());
598       }
599       else
600         label = -1;
601     }
602     else{
603       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
604     }
605     
606     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
607       hTracking[0]->Fill(track->Pt(),fCentrality);
608       if(pdg == 211)
609         hTracking[1]->Fill(track->Pt(),fCentrality);
610       else if(pdg == 321)
611         hTracking[2]->Fill(track->Pt(),fCentrality);
612       else if(pdg == 2212)
613         hTracking[3]->Fill(track->Pt(),fCentrality);
614       else if(! mcp){ // Fill matching histo with the prob
615         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
616         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
617         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
618       }
619     }
620     
621     if(!tofMatch) continue;
622     
623     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
624       hMatching[0]->Fill(track->Pt(),fCentrality);
625       if(pdg == 211)
626         hMatching[1]->Fill(track->Pt(),fCentrality);
627       else if(pdg == 321)
628         hMatching[2]->Fill(track->Pt(),fCentrality);
629       else if(pdg == 2212)
630         hMatching[3]->Fill(track->Pt(),fCentrality);
631       else if(! mcp){ // Fill matching histo with the prob
632         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
633         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
634         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
635       }
636     }
637   }
638   
639 //   Int_t pdg1 = -1;
640 //   Int_t pdg2 = -1;
641
642
643   // start analysis Lambda
644   for(Int_t i=0;i < ntrack;i++){ // loop on proton candidate tracks
645     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
646         
647     if (!KpTrack){
648       continue;
649     }
650     
651     if(!(KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
652
653     Bool_t isLambda = (KpTrack->Charge() > 0);
654
655     nSigmaComb=5;
656     nSigmaTOF=5;
657     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
658     fPidKp=0;
659
660     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
661
662     Double_t oldpP[10];
663     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
664
665     
666     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
667
668     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
669
670     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
671     /*
672     if(mcArray){
673       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
674       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
675       pdg1 = TMath::Abs(mcp1->GetPdgCode());
676     }
677     */
678
679     fPidKp = Int_t(probP[4]*100);
680
681     if(tofMatch1){
682       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
683         // remove this tof hit
684         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
685         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
686         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
687         fPidKp = Int_t(probP[4]*100);
688         tofMatch1=0;
689       }
690       else{
691         if(probP[4] > probP[3] && probP[4] > probP[2] && probP[4] > probP[0]) fPidKp += 128; // max prob
692         
693         if(isLambda){
694           nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
695           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
696           if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
697           nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
698           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
699           if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
700           nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
701           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
702           if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
703         }
704         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
705         if(isLambda){
706           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
707           if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
708         }
709         
710         if(fIsMC){
711           Float_t mismAdd = addMismatchForMC;
712           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
713           
714           if(gRandom->Rndm() < mismAdd){
715             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
716             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
717             channel = channel % 8736;
718             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
719             
720             // generate random time
721             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
722             Double_t times[10];
723             KpTrack->GetIntegratedTimes(times);
724             nSigmaTOF = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kProton);
725           }
726         }
727
728         if(fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
729         nSigmaTOF = TMath::Abs(nSigmaTOF);
730
731         if(nSigmaTOF < 2) fPidKp += 256;
732         else if(nSigmaTOF < 3) fPidKp += 512;
733       }
734     }
735     
736     if(tofMatch1){
737       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
738       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2){
739         fCombsignal->Fill(nSigmaComb);
740       }
741     } else {
742       nSigmaComb = TMath::Abs(nSigmaTPC);
743     }
744
745     // use sigmaTOF instead of sigmaComb
746     if(tofMatch1) nSigmaComb = nSigmaTOF;
747
748     if(nSigmaComb < 2) nSigmaComb = 2;
749     else if(nSigmaComb < 3) nSigmaComb = 3;
750     else if(nSigmaComb < 5) nSigmaComb = 4.99;
751     else nSigmaComb = 6;
752
753     Int_t iks=-1;
754     for(Int_t k=0;k < fNLambda;k++){ // find the Lambda which contains the positive track
755       if(i == fIpP[k]) iks = k;
756     }
757
758     if(fPtKp > 4.299) fPtKp = 4.299;
759
760     if(iks > -1 && fIpN[iks] > -1){
761       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
762       Int_t j = fIpN[iks];
763       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
764       
765       if (!KnTrack){
766         continue;
767       }
768
769       nSigmaComb2 = 5;
770       nSigmaTOF2 = 5;
771
772       if(!(KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
773
774       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
775       fPidKn=0;
776
777       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
778       Double_t oldpN[10];
779       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
780
781       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
782       
783       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
784       
785       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
786       /*
787       if(mcArray){
788         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
789         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
790         pdg2 = TMath::Abs(mcp2->GetPdgCode());
791       }
792       */
793
794       fPidKn = Int_t(probN[2]*100);
795
796       if(tofMatch2){
797         if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
798
799         nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
800
801         nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
802
803         if(nSigmaTOF2 < 2) fPidKn += 256;
804         else if(nSigmaTOF2 < 3) fPidKn += 512;
805       }
806
807       px = KpTrack->Px() + KnTrack->Px();
808       py = KpTrack->Py() + KnTrack->Py();
809       pz = KpTrack->Pz() + KnTrack->Pz();
810       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 0.938e-01*0.938e-01);
811       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
812
813       ksV.SetPxPyPzE(px,py,pz,E);
814       fMassV0 = fMassLambda[iks];
815       
816       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
817
818
819       fPtLambdaC = ksV.Pt();
820       fEtaLambda = ksV.Eta();
821       fPhiLambdaC = ksV.Phi();
822
823       if(tofMatch2){
824         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
825       } else {
826         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
827       }
828
829       // use sigmaTOF instead of sigmaComb
830       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
831
832       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
833       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
834       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
835       else nSigmaComb2 = 6;  
836
837       Bool_t isTrue = kFALSE;
838
839       if(mcArray){
840         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
841         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
842
843         if(labelP > -1 && labelN > -1){
844           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
845           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
846
847           Int_t mP = partP->GetMother();
848           Int_t mN = partN->GetMother();
849           if(mP == mN && mP > -1){
850             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
851             Int_t pdgM = partM->GetPdgCode();
852             if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
853           }
854         }
855
856       }
857
858       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
859       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
860
861       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
862       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
863       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
864       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
865
866       if(fPtKn > 4.299) fPtKn = 4.299;
867
868       Float_t xTOfill[] = {fPtLambdaC,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
869       Float_t xTOfill2[] = {fPtLambdaC,KpTrack->Eta(),fPtKn,fPtKp,(fPidKn%128)*0.01,(fPidKp%128)*0.01,tofMatch2,tofMatch1,isTrue,nSigmaComb2,nSigmaComb,deltaphi2,deltaphi1,fPsi};
870       
871
872       Int_t ipt = 0;
873       while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
874         ipt++;
875       }
876       ipt--;
877       if(ipt < 0) ipt = 0; // just to be sure
878
879       if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
880         if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
881         else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
882
883
884
885         if(fPIDuserCut){
886          Float_t xUser[] = {KpTrack->Eta(),fPtKp,isTrue,0,deltaphi1,fPsi};
887          
888          if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled
889            xUser[3] = 1;
890          } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
891            xUser[3] = 2;
892          } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
893            xUser[3] = 3;
894          }
895          
896          if(isLambda) fContUser->Fill(0,fMassV0,fCentrality,xUser);
897          else fContUser2->Fill(0,fMassV0,fCentrality,xUser);
898         }
899       }
900
901
902     }
903   } // end analysi Lambda
904 }
905
906 //_____________________________________________________________________________
907 Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
908 {
909
910   Float_t zvtx = -999;
911
912   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
913   if (!vtxAOD)
914     return zvtx;
915   if(vtxAOD->GetNContributors()>0)
916     zvtx = vtxAOD->GetZ();
917   
918   return zvtx;
919 }
920 //_____________________________________________________________________________
921 void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
922
923   // Terminate loop
924   Printf("Terminate()");
925 }
926 //=======================================================================
927 void AliAnalysisTaskLambdaBayes::SelectLambda(){
928   fNLambda=0;
929   fNpPos=0;
930   fNpNeg=0;
931
932   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
933   AliAODv0 *myV0;
934   Double_t dMASS=0.0;
935   for (Int_t i=0; i!=nV0s; ++i) {
936     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
937     if(!myV0) continue;
938     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
939     Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
940     if(pass) {
941       if(pass==1) dMASS = myV0->MassLambda();
942       else dMASS = myV0->MassAntiLambda();
943
944       Float_t massKs = myV0->MassK0Short();
945
946       if(TMath::Abs(dMASS-1.115)/0.005 < 8 && TMath::Abs(massKs - 0.497)/0.005 > 8){
947         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
948         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
949         if(iT->Charge()<0){
950           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
951           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
952         }
953         fPhiLambda[fNLambda] = myV0->Phi();
954         fPtLambda[fNLambda] = myV0->Pt();
955         fIpP[fNLambda] = FindDaugheterIndex(iT);
956         fIpN[fNLambda] = FindDaugheterIndex(jT);
957
958         if(pass==2){
959           Int_t itemp = fIpP[fNLambda];
960           fIpP[fNLambda] = fIpN[fNLambda];
961           fIpN[fNLambda] = itemp;
962         }
963
964         fMassLambda[fNLambda] = dMASS;
965         if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
966           fNLambda++;
967         }
968
969       }
970     }
971   }
972 }
973
974 //=======================================================================
975 Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
976
977   if (myV0->GetOnFlyStatus() ) return 0;
978   
979   //the following is needed in order to evualuate track-quality
980   AliAODTrack *iT, *jT;
981   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
982   Double_t pos[3],cov[6];
983   vV0s->GetXYZ(pos);
984   vV0s->GetCovarianceMatrix(cov);
985   const AliESDVertex vESD(pos,cov,100.,100);
986   
987   // TESTING CHARGE
988   int iPos, iNeg;
989   iT=(AliAODTrack*) myV0->GetDaughter(0);
990   if(iT->Charge()>0) {
991     iPos = 0; iNeg = 1;
992   } else {
993     iPos = 1; iNeg = 0;
994   }
995   // END OF TEST
996
997   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
998
999   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1000
1001   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1002   if(!trkFlag) return 0;
1003   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1004   if(!trkFlag2) return 0;
1005
1006   Double_t pvertex[3];
1007   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1008   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1009   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1010
1011   Double_t dDL=myV0->DecayLengthV0( pvertex );
1012   if(dDL  < 0.5 || dDL > 25) return 0;
1013
1014   Double_t dDCA=myV0->DcaV0Daughters();
1015   if(dDCA >0.5) return 0;
1016
1017   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1018   if(dCTP < -1) return 0;
1019
1020 //   AliESDtrack ieT( iT );
1021 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1022 //   if(TMath::Abs(dD0P) < 0]) return 0;
1023
1024 //   AliESDtrack jeT( jT );
1025 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1026 //   if(TMath::Abs(dD0M) < 0) return 0;
1027
1028 //   Double_t dD0D0=dD0P*dD0M;
1029 //   if(dD0D0>0) return 0;
1030
1031 //   Double_t dETA=myV0->Eta();
1032 //   if(dETA <-0.8) return 0;
1033 //   if(dETA >0.8) return 0;
1034
1035 //   Double_t dQT=myV0->PtArmV0();
1036 //   if(specie==0) if(dQT<???) return 0;
1037
1038   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1039   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1040
1041   if(specie==1 && dALPHA<0) return 2; // antilambda
1042   return 1; // Ks or lambda
1043 }
1044 //-------------------------------------------------------------------------
1045 Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
1046   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1047
1048   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1049     AliAODTrack* track = fOutputAOD->GetTrack(i);
1050     if(track == trk) return i;
1051   }
1052   
1053   printf("Daughter for %p not found\n",trk);
1054   return -1;
1055 }
1056 //-------------------------------------------------------------------------
1057 Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
1058   if(!fIsMC) return 1; // used only on MC
1059
1060   if( fTypeCol==2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1061     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs); 
1062   
1063     if(!(channel%20)) return 0; // 5% additional loss in MC
1064   }
1065
1066   return 1;
1067 }