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