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