6224755eb523bb3aa4cf4278ab11c4f58d81ab74
[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(!fTypeCol){
398         v0Centr=100./(fOutputAOD->GetNumberOfTracks()/12.+1);
399         trkCentr=v0Centr;
400       }
401
402       if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
403         fCentrality = v0Centr;
404         Analyze(fOutputAOD); // Do analysis!!!!
405
406       }
407     }
408     
409 }
410
411 //________________________________________________________________________
412 void AliAnalysisTaskPhiBayes::Analyze(AliAODEvent* aodEvent)
413 {
414   Int_t ntrack = aodEvent->GetNumberOfTracks();
415
416   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
417   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
418   fPidKp=0,fPidKn=0;
419   fMassV0=-1;
420   
421   TLorentzVector phiV;
422   
423   Double_t px,py,pz,E;
424
425   Float_t invMmin = 0.985;
426   Float_t invMmax = 1.045;
427   
428   Int_t icentr = 8;
429   if(fCentrality < 0) icentr = 8;
430   else if(fCentrality < 10) icentr = 0;
431   else if(fCentrality < 20) icentr = 1;
432   else if(fCentrality < 30) icentr = 2;
433   else if(fCentrality < 40) icentr = 3;
434   else if(fCentrality < 50) icentr = 4;
435   else if(fCentrality < 60) icentr = 5;
436   else if(fCentrality < 70) icentr = 6;
437   else if(fCentrality < 80) icentr = 7;
438
439   Float_t addMismatchForMC = 0.005;
440   if(fCentrality < 50) addMismatchForMC += 0.005;
441   if(fCentrality < 20) addMismatchForMC += 0.02;
442
443   if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
444
445   fPsi = 0;
446   /* Compute TPC EP */
447   Double_t Qx2 = 0, Qy2 = 0;
448   Double_t Qx3 = 0, Qy3 = 0;
449   for(Int_t iT = 0; iT < ntrack; iT++) {
450     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
451     
452     if (!aodTrack){
453       continue;
454     }
455     
456     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
457     
458     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag) 
459       continue;
460     
461     Double_t b[2] = {-99., -99.};
462     Double_t bCov[3] = {-99., -99., -99.};
463     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
464       continue;
465     
466     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
467       continue;
468     
469     Qx2 += TMath::Cos(2*aodTrack->Phi()); 
470     Qy2 += TMath::Sin(2*aodTrack->Phi());
471     Qx3 += TMath::Cos(3*aodTrack->Phi()); 
472     Qy3 += TMath::Sin(3*aodTrack->Phi());
473     
474   }
475
476   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
477
478   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
479   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
480   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
481
482 //   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
483 //   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
484 //   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
485 //   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
486
487   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
488
489   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
490   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
491   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
492
493   
494   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
495   TClonesArray *mcArray = NULL;
496   if (mcHeader)
497     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
498
499   Int_t nmc = 0;
500   if(mcArray)
501     nmc = mcArray->GetEntries();
502
503   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
504     AliAODTrack* track = aodEvent->GetTrack(i);
505     
506     AliAODMCParticle *mcp = NULL;
507     Int_t pdg = 0;
508     
509     if (!track){
510       continue;
511     }
512     
513     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
514     
515     Int_t label = -1;
516     if(mcArray){
517       label = track->GetLabel();
518       if(label != -1 && label < nmc){
519         label = TMath::Abs(label);
520         mcp = (AliAODMCParticle*)mcArray->At(label);
521         pdg = TMath::Abs(mcp->GetPdgCode());
522       }
523       else
524         label = -1;
525     }
526     else{
527       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
528     }
529     
530     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
531       hTracking[0]->Fill(track->Pt(),fCentrality);
532       if(pdg == 211)
533         hTracking[1]->Fill(track->Pt(),fCentrality);
534       else if(pdg == 321)
535         hTracking[2]->Fill(track->Pt(),fCentrality);
536       else if(pdg == 2212)
537         hTracking[3]->Fill(track->Pt(),fCentrality);
538       else if(! mcp){ // Fill matching histo with the prob
539         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
540         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
541         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
542       }
543     }
544     
545     if(!tofMatch) continue;
546     
547     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
548       hMatching[0]->Fill(track->Pt(),fCentrality);
549       if(pdg == 211)
550         hMatching[1]->Fill(track->Pt(),fCentrality);
551       else if(pdg == 321)
552         hMatching[2]->Fill(track->Pt(),fCentrality);
553       else if(pdg == 2212)
554         hMatching[3]->Fill(track->Pt(),fCentrality);
555       else if(! mcp){ // Fill matching histo with the prob
556         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
557         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
558         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
559       }
560     }
561   }
562   
563 //   Int_t pdg1 = -1;
564 //   Int_t pdg2 = -1;
565
566
567   // start analysis phi
568   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
569     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
570         
571     if (!KpTrack){
572       continue;
573     }
574     
575     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
576
577     nSigmaTOF=5;
578     nSigmaComb=5;
579     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
580     fPidKp=0;
581
582     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
583
584     Double_t oldpP[10];
585     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
586
587     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
588
589     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
590
591     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
592     /*
593     if(mcArray){
594       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
595       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
596       pdg1 = TMath::Abs(mcp1->GetPdgCode());
597     }
598     */
599
600     fPidKp = Int_t(probP[3]*100);
601
602     if(tofMatch1){
603       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
604         // remove this tof hit
605         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
606         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
607         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
608         fPidKp = Int_t(probP[4]*100);
609         tofMatch1=0;
610       }
611       else{
612         if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
613         
614         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
615         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
616         if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
617         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
618         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
619         if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
620         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
621         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
622         if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
623         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
624         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
625         if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
626         
627         if(fIsMC){
628           Float_t mismAdd = addMismatchForMC;
629           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
630           
631           if(gRandom->Rndm() < mismAdd){
632             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
633             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
634             channel = channel % 8736;
635             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
636             
637             // generate random time
638             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
639             Double_t times[10];
640             KpTrack->GetIntegratedTimes(times);
641             nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
642           }
643         }
644
645         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
646         nSigmaTOF = TMath::Abs(nSigmaTOF);
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 }