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