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