improved the macro to extract PID perf (automaticly through a new macro)
[u/mrichter/AliRoot.git] / PWGPP / pid / AliAnalysisTaskPhiBayes.cxx
1 #include "AliAnalysisTaskPhiBayes.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 AliAnalysisTaskPhiBayes::fPtPhiMin[AliAnalysisTaskPhiBayes::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 AliAnalysisTaskPhiBayes::fPtPhiMax[AliAnalysisTaskPhiBayes::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(AliAnalysisTaskPhiBayes)
39 //_____________________________________________________________________________
40 AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes():
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   fPtPhi(0.),
55   fPhiPhi(0.),
56   fEtaPhi(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   fHmismTOF(0),
73   fHchannelTOFdistr(0),
74   fTypeCol(2)
75 {
76   // Default constructor (should not be used)
77   fList->SetName("contPhiBayes1");
78   fList2->SetName("contPhiBayes2");
79   fList3->SetName("contPhiBayes3");
80
81   fList->SetOwner(kTRUE); 
82   fList2->SetOwner(kTRUE); 
83   fList3->SetOwner(kTRUE); 
84
85   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
86   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
87
88   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
89   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
90
91   for(Int_t i=0;i < nCentrBin;i++){
92     fElTOF[i] = NULL; 
93     fElTPC[i] = NULL; 
94     fPiTOF[i] = NULL; 
95     fPiTPC[i] = NULL; 
96     fKaTOF[i] = NULL; 
97     fKaTPC[i] = NULL; 
98     fPrTOF[i] = NULL; 
99     fPrTPC[i] = NULL; 
100   }
101   for(Int_t i=0;i < 4;i++){
102     hMatching[i] = NULL;
103     hTracking[i] = NULL;
104   }
105 }
106
107 //______________________________________________________________________________
108 AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes(const char *name):
109   AliAnalysisTaskSE(name),
110   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
111   fEtaCut(0.8),   // cut on |eta| < fEtaCut
112   fMinPt(0.15),   // cut on pt > fMinPt
113   fIsMC(kFALSE),
114   fQAsw(kFALSE),
115   fNcluster(70),
116   fFilterBit(1),
117   fList(new TList()),
118   fList2(new TList()),
119   fList3(new TList()),
120   fCentrality(-1),
121   fPsi(0),
122   fPtPhi(0.),
123   fPhiPhi(0.),
124   fEtaPhi(0.),
125   fMassV0(0.),
126   fPtKp(0.),
127   fPhiKp(0.),
128   fEtaKp(0.),
129   fPtKn(0.),
130   fPhiKn(0.),
131   fEtaKn(0.),
132   fPidKp(0),
133   fPidKn(0),
134   fTOFTPCsignal(0),
135   fCombsignal(0),
136   fCutsDaughter(NULL),
137   fPIDCombined(NULL),
138   fContPid(NULL),
139   fContPid2(NULL),
140   fHmismTOF(0),
141   fHchannelTOFdistr(0),
142   fTypeCol(2)
143 {
144
145   DefineOutput(1, TList::Class());
146   DefineOutput(2, TList::Class());
147   DefineOutput(3, TList::Class());
148   DefineOutput(4, TList::Class());
149
150   // Output slot #1 writes into a TTree
151   fList->SetName("contPhiBayes1");
152   fList2->SetName("contPhiBayes2");
153   fList3->SetName("contPhiBayes3");
154
155   fList->SetOwner(kTRUE); 
156   fList2->SetOwner(kTRUE); 
157   fList3->SetOwner(kTRUE); 
158
159   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
160   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
161
162   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
163   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
164
165   for(Int_t i=0;i < nCentrBin;i++){
166     fElTOF[i] = NULL; 
167     fElTPC[i] = NULL; 
168     fPiTOF[i] = NULL; 
169     fPiTPC[i] = NULL; 
170     fKaTOF[i] = NULL; 
171     fKaTPC[i] = NULL; 
172     fPrTOF[i] = NULL; 
173     fPrTPC[i] = NULL; 
174   }
175   for(Int_t i=0;i < 4;i++){
176     hMatching[i] = NULL;
177     hTracking[i] = NULL;
178   }
179 }
180 //_____________________________________________________________________________
181 AliAnalysisTaskPhiBayes::~AliAnalysisTaskPhiBayes()
182 {
183
184 }
185
186 //______________________________________________________________________________
187 void AliAnalysisTaskPhiBayes::UserCreateOutputObjects()
188 {
189
190   fPIDCombined=new AliPIDCombined;
191   fPIDCombined->SetDefaultTPCPriors();
192   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
193
194   Float_t invMmin = 0.985;
195   Float_t invMmax = 1.045;
196
197   const Int_t nBinPid = 14;
198
199   Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaKP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
200
201   Int_t binPid2[nBinPid] = {1/*ptPhi*/,8/*EtaKN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
202
203   fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
204   fContPid->SetTitleX("M_{#phi}");
205   fContPid->SetTitleY("centrality (%)");
206   fContPid->SetVarName(0,"p_{T}^#phi}");
207   fContPid->SetVarRange(0,0,10);
208   fContPid->SetVarName(1,"#eta^{#phi}");
209   fContPid->SetVarRange(1,-0.8,0.8);
210   fContPid->SetVarName(2,"p_{T}^{Kp}");
211   fContPid->SetVarRange(2,0.3,4.3);
212   fContPid->SetVarName(3,"p_{T}^{Kn}");
213   fContPid->SetVarRange(3,0.3,4.3);
214   fContPid->SetVarName(4,"BayesProb^{Kp}");
215   fContPid->SetVarRange(4,0,1.);
216   fContPid->SetVarName(5,"BayesProb^{Kn}");
217   fContPid->SetVarRange(5,0,1.);
218   fContPid->SetVarName(6,"isTOFmatch^{Kp}");
219   fContPid->SetVarRange(6,-0.5,1.5);
220   fContPid->SetVarName(7,"isTOFmatch^{Kn}");
221   fContPid->SetVarRange(7,-0.5,1.5);
222   fContPid->SetVarName(8,"isPhiTrue^{Kn}");
223   fContPid->SetVarRange(8,-0.5,1.5);
224   fContPid->SetVarName(9,"N#sigma^{Kp}");
225   fContPid->SetVarRange(9,1.25,6.25);
226   fContPid->SetVarName(10,"N#sigma^{Kn}");
227   fContPid->SetVarRange(10,1.25,6.25);
228   fContPid->SetVarName(11,"#Delta#phi^{Kp}");
229   fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
230   fContPid->SetVarName(12,"#Delta#phi^{Kn}");
231   fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
232   fContPid->SetVarName(13,"#Psi");
233   fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
234
235   fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
236   fContPid2->SetTitleX("M_{#phi}");
237   fContPid2->SetTitleY("centrality (%)");
238   fContPid2->SetVarName(0,"p_{T}^{#phi}");
239   fContPid2->SetVarRange(0,0,10);
240   fContPid2->SetVarName(1,"#eta^{#phi}");
241   fContPid2->SetVarRange(1,-0.8,0.8);
242   fContPid2->SetVarName(2,"p_{T}^{Kp}");
243   fContPid2->SetVarRange(2,0.3,4.3);
244   fContPid2->SetVarName(3,"p_{T}^{Kn}");
245   fContPid2->SetVarRange(3,0.3,4.3);
246   fContPid2->SetVarName(4,"BayesProb^{Kp}");
247   fContPid2->SetVarRange(4,0,1.);
248   fContPid2->SetVarName(5,"BayesProb^{Kn}");
249   fContPid2->SetVarRange(5,0,1.);
250   fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
251   fContPid2->SetVarRange(6,-0.5,1.5);
252   fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
253   fContPid2->SetVarRange(7,-0.5,1.5);
254   fContPid2->SetVarName(8,"isPhiTrue^{Kn}");
255   fContPid2->SetVarRange(8,-0.5,1.5);
256   fContPid2->SetVarName(9,"N#sigma^{Kp}");
257   fContPid2->SetVarRange(9,1.25,6.25);
258   fContPid2->SetVarName(10,"N#sigma^{Kn}");
259   fContPid2->SetVarRange(10,1.25,6.25);
260   fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
261   fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
262   fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
263   fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
264   fContPid2->SetVarName(13,"#Psi");
265   fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
266
267
268
269   const Int_t nDETsignal = 100; // mass
270   Double_t binDETsignal[nDETsignal+1];
271   for(Int_t i=0;i<nDETsignal+1;i++){
272     binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
273   }
274   const Int_t nDETsignal2 = 10; // centrality
275   Double_t binDETsignal2[nDETsignal2+1];
276   for(Int_t i=0;i<nDETsignal2+1;i++){
277     binDETsignal2[i] = i*100./ nDETsignal2;
278   }
279   fContPid->AddSpecies("phi",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
280   fContPid2->AddSpecies("phi2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
281
282   fList->Add(fContPid);
283   fList->Add(fContPid2);
284
285   hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
286   hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
287   hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
288   hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
289
290   hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
291   hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
292   hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
293   hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
294
295   fList2->Add(hMatching[0]);
296   fList2->Add(hMatching[1]);
297   fList2->Add(hMatching[2]);
298   fList2->Add(hMatching[3]);
299   fList2->Add(hTracking[0]);
300   fList2->Add(hTracking[1]);
301   fList2->Add(hTracking[2]);
302   fList2->Add(hTracking[3]);
303
304   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);
305   fList2->Add(fTOFTPCsignal);
306   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);
307   fList2->Add(fCombsignal);
308
309   // QA plots
310   char name[100],title[400];
311   for(Int_t i=0;i < nCentrBin;i++){
312     snprintf(name,100,"hPiTPCQA_%i",i);
313     snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
314     fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
315     fList3->Add(fPiTPC[i]);
316
317     snprintf(name,100,"hKaTPCQA_%i",i);
318     snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
319     fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
320     fList3->Add(fKaTPC[i]);
321
322     snprintf(name,100,"hPrTPCQA_%i",i);
323     snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
324     fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
325     fList3->Add(fPrTPC[i]);
326
327     snprintf(name,100,"hElTPCQA_%i",i);
328     snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
329     fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
330     fList3->Add(fElTPC[i]);
331
332     snprintf(name,100,"hPiTOFQA_%i",i);
333     snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
334     fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
335     fList3->Add(fPiTOF[i]);
336
337     snprintf(name,100,"hKaTOFQA_%i",i);
338     snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
339     fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
340     fList3->Add(fKaTOF[i]);
341
342     snprintf(name,100,"hPrTOFQA_%i",i);
343     snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
344     fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
345     fList3->Add(fPrTOF[i]);
346
347     snprintf(name,100,"hElTOFQA_%i",i);
348     snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
349     fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
350     fList3->Add(fElTOF[i]);
351
352   }
353
354   // Post output data.
355   PostData(1, fList);
356   PostData(2, fList2);
357   PostData(3, fList3);
358 }
359
360 //______________________________________________________________________________
361 void AliAnalysisTaskPhiBayes::UserExec(Option_t *) 
362 {
363     // Main loop
364     // Called for each event
365     
366     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
367     if(!fOutputAOD){
368         Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
369         this->Dump();
370         return;
371     }
372     
373     Float_t zvtx = GetVertex(fOutputAOD);
374
375
376
377     //Get the MC object
378     if(fIsMC){
379       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
380       if (!mcHeader) {
381         AliError("Could not find MC Header in AOD");
382         return;
383       }
384     }
385
386     if (TMath::Abs(zvtx) < fVtxCut) {
387
388       //Centrality
389       Float_t v0Centr  = -10.;
390       Float_t trkCentr  = -10.;
391       AliCentrality *centrality = fOutputAOD->GetCentrality();
392       if (centrality){
393         trkCentr  = centrality->GetCentralityPercentile("V0M");
394         v0Centr = centrality->GetCentralityPercentile("TRK"); 
395       }
396
397       if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
398         fCentrality = v0Centr;
399         Analyze(fOutputAOD); // Do analysis!!!!
400
401       }
402     }
403     
404 }
405
406 //________________________________________________________________________
407 void AliAnalysisTaskPhiBayes::Analyze(AliAODEvent* aodEvent)
408 {
409   Int_t ntrack = aodEvent->GetNumberOfTracks();
410
411   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
412   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
413   fPidKp=0,fPidKn=0;
414   fMassV0=-1;
415   
416   TLorentzVector phiV;
417   
418   Double_t px,py,pz,E;
419
420   Float_t invMmin = 0.985;
421   Float_t invMmax = 1.045;
422   
423   Int_t icentr = 8;
424   if(fCentrality < 0) icentr = 8;
425   else if(fCentrality < 10) icentr = 0;
426   else if(fCentrality < 20) icentr = 1;
427   else if(fCentrality < 30) icentr = 2;
428   else if(fCentrality < 40) icentr = 3;
429   else if(fCentrality < 50) icentr = 4;
430   else if(fCentrality < 60) icentr = 5;
431   else if(fCentrality < 70) icentr = 6;
432   else if(fCentrality < 80) icentr = 7;
433
434   Float_t addMismatchForMC = 0.005;
435   if(fCentrality < 50) addMismatchForMC += 0.01;
436   if(fCentrality < 20) addMismatchForMC += 0.02;
437
438   if(fTypeCol == 0) addMismatchForMC = 0.005;
439   else if(fTypeCol == 1) addMismatchForMC *= 0.5;
440
441   fPsi = 0;
442   /* Compute TPC EP */
443   Double_t Qx2 = 0, Qy2 = 0;
444   Double_t Qx3 = 0, Qy3 = 0;
445   for(Int_t iT = 0; iT < ntrack; iT++) {
446     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
447     
448     if (!aodTrack){
449       continue;
450     }
451     
452     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
453     
454     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag) 
455       continue;
456     
457     Double_t b[2] = {-99., -99.};
458     Double_t bCov[3] = {-99., -99., -99.};
459     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
460       continue;
461     
462     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
463       continue;
464     
465     Qx2 += TMath::Cos(2*aodTrack->Phi()); 
466     Qy2 += TMath::Sin(2*aodTrack->Phi());
467     Qx3 += TMath::Cos(3*aodTrack->Phi()); 
468     Qy3 += TMath::Sin(3*aodTrack->Phi());
469     
470   }
471
472   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
473
474   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
475   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
476   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
477
478   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
479   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
480   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
481   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
482
483   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
484
485   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
486   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
487   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
488
489   
490   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
491   TClonesArray *mcArray = NULL;
492   if (mcHeader)
493     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
494
495   Int_t nmc = 0;
496   if(mcArray)
497     nmc = mcArray->GetEntries();
498
499   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
500     AliAODTrack* track = aodEvent->GetTrack(i);
501     
502     AliAODMCParticle *mcp = NULL;
503     Int_t pdg = 0;
504     
505     if (!track){
506       continue;
507     }
508     
509     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
510     
511     Int_t label = -1;
512     if(mcArray){
513       label = track->GetLabel();
514       if(label != -1 && label < nmc){
515         label = TMath::Abs(label);
516         mcp = (AliAODMCParticle*)mcArray->At(label);
517         pdg = TMath::Abs(mcp->GetPdgCode());
518       }
519       else
520         label = -1;
521     }
522     else{
523       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
524     }
525     
526     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
527       hTracking[0]->Fill(track->Pt(),fCentrality);
528       if(pdg == 211)
529         hTracking[1]->Fill(track->Pt(),fCentrality);
530       else if(pdg == 321)
531         hTracking[2]->Fill(track->Pt(),fCentrality);
532       else if(pdg == 2212)
533         hTracking[3]->Fill(track->Pt(),fCentrality);
534       else if(! mcp){ // Fill matching histo with the prob
535         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
536         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
537         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
538       }
539     }
540     
541     if(!tofMatch) continue;
542     
543     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
544       hMatching[0]->Fill(track->Pt(),fCentrality);
545       if(pdg == 211)
546         hMatching[1]->Fill(track->Pt(),fCentrality);
547       else if(pdg == 321)
548         hMatching[2]->Fill(track->Pt(),fCentrality);
549       else if(pdg == 2212)
550         hMatching[3]->Fill(track->Pt(),fCentrality);
551       else if(! mcp){ // Fill matching histo with the prob
552         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
553         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
554         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
555       }
556     }
557   }
558   
559 //   Int_t pdg1 = -1;
560 //   Int_t pdg2 = -1;
561
562
563   // start analysis phi
564   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
565     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
566         
567     if (!KpTrack){
568       continue;
569     }
570     
571     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
572
573     nSigmaTOF=5;
574     nSigmaComb=5;
575     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
576     fPidKp=0;
577
578     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
579
580     Double_t oldpP[10];
581     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
582
583     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
584     fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
585     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
586     fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
587     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
588     fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
589     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
590     fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
591
592     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
593
594     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
595     /*
596     if(mcArray){
597       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
598       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
599       pdg1 = TMath::Abs(mcp1->GetPdgCode());
600     }
601     */
602
603     fPidKp = Int_t(probP[3]*100);
604
605     if(tofMatch1){
606       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
607         // remove this tof hit
608         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
609         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
610         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
611         fPidKp = Int_t(probP[4]*100);
612         tofMatch1=0;
613       }
614       else{
615         if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
616         
617         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
618         fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
619         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
620         fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
621         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
622         fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
623         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
624         fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
625                 
626         nSigmaTOF = TMath::Abs(nSigmaTOF);
627         
628         if(fIsMC){
629           Float_t mismAdd = addMismatchForMC;
630           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
631           
632           if(gRandom->Rndm() < mismAdd){
633             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
634             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
635             channel = channel % 8736;
636             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
637             
638             // generate random time
639             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
640             Double_t times[10];
641             KpTrack->GetIntegratedTimes(times);
642             nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
643           }
644         }
645
646         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
647
648         if(nSigmaTOF < 2) fPidKp += 256;
649         else if(nSigmaTOF < 3) fPidKp += 512;
650       }
651     }
652     
653     if(tofMatch1){
654       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
655       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
656         fCombsignal->Fill(nSigmaComb);
657       }
658     } else {
659       nSigmaComb = TMath::Abs(nSigmaTPC);
660     }
661
662     // use sigmaTOF instead of sigmaComb
663     if(tofMatch1) nSigmaComb = nSigmaTOF;
664
665     if(nSigmaComb < 2) nSigmaComb = 2;
666     else if(nSigmaComb < 3) nSigmaComb = 3;
667     else if(nSigmaComb < 5) nSigmaComb = 4.99;
668     else nSigmaComb = 6;
669
670     if(fPtKp > 4.299) fPtKp = 4.299;
671
672     for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
673       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
674       
675       if (!KnTrack){
676         continue;
677       }
678
679       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
680
681       px = KpTrack->Px() + KnTrack->Px();
682       py = KpTrack->Py() + KnTrack->Py();
683       pz = KpTrack->Pz() + KnTrack->Pz();
684       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
685       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
686
687       phiV.SetPxPyPzE(px,py,pz,E);
688       fMassV0 = phiV.M();
689       
690       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
691
692       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
693       fPidKn=0;
694
695       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
696       Double_t oldpN[10];
697       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
698
699       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
700       
701       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
702       
703       nSigmaTOF2=5;
704       nSigmaComb2=5;
705
706       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
707       /*
708       if(mcArray){
709         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
710         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
711         pdg2 = TMath::Abs(mcp2->GetPdgCode());
712       }
713       */
714
715       fPidKn = Int_t(probN[3]*100);
716
717       if(tofMatch2){
718         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
719           // remove this tof hit
720           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
721           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
722           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
723           fPidKn = Int_t(probN[4]*100);
724           tofMatch2=0;
725         }
726         else{
727           if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
728           
729           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
730                   
731           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
732           
733           if(fIsMC){
734             Float_t mismAdd = addMismatchForMC;
735             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
736             
737             if(gRandom->Rndm() < mismAdd){
738               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
739               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
740               channel = channel % 8736;
741               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
742               
743               // generate random time
744               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
745               Double_t times[10];
746               KnTrack->GetIntegratedTimes(times);
747               nSigmaTOF2 = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
748             }
749           }
750
751           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
752
753           if(nSigmaTOF2 < 2) fPidKn += 256;
754           else if(nSigmaTOF2 < 3) fPidKn += 512;
755         }
756       }
757
758       fPtPhi = phiV.Pt();
759       fEtaPhi = phiV.Eta();
760       fPhiPhi = phiV.Phi();
761
762       if(tofMatch2){
763         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
764         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
765           fCombsignal->Fill(nSigmaComb2);
766         }
767       } else {
768         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
769       }
770
771       // use sigmaTOF instead of sigmaComb
772       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
773
774       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
775       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
776       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
777       else nSigmaComb2 = 6;  
778
779       Bool_t isTrue = kFALSE;
780
781       if(mcArray){
782         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
783         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
784
785         if(labelP > -1 && labelN > -1){
786           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
787           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
788
789           Int_t mP = partP->GetMother();
790           Int_t mN = partN->GetMother();
791           if(mP == mN && mP > -1){
792             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
793             Int_t pdgM = partM->GetPdgCode();
794             if(pdgM == 333) isTrue = kTRUE;
795           }
796         }
797
798       }
799
800       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
801       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
802
803       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
804       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
805       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
806       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
807
808       if(fPtKn > 4.299) fPtKn = 4.299;
809
810       Float_t xTOfill[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
811       
812
813       Int_t ipt = 0;
814       while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
815         ipt++;
816       }
817       ipt--;
818       if(ipt < 0) ipt = 0; // just to be sure
819
820       if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
821         if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
822         xTOfill[1] = KnTrack->Eta();
823         if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);// use tagging on positive track
824
825       }
826
827
828     }
829   } // end analysi phi
830
831 }
832
833 //_____________________________________________________________________________
834 Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
835 {
836
837   Float_t zvtx = -999;
838
839   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
840   if (!vtxAOD)
841     return zvtx;
842   if(vtxAOD->GetNContributors()>0)
843     zvtx = vtxAOD->GetZ();
844   
845   return zvtx;
846 }
847 //_____________________________________________________________________________
848 void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
849
850   // Terminate loop
851   Printf("Terminate()");
852 }
853 //-------------------------------------------------------------------------
854 Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
855   if(!fIsMC) return 1; // used only on MC
856
857   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
858     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs); 
859   
860     if(!(channel%20)) return 0; // 5% additional loss in MC
861   }
862
863   return 1;
864 }