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 = aodEvent->GetTrack(iT);
 501
 502     if (!aodTrack){
 503       continue;
 504     }
 505
 506     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 507
 508     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 509       continue;
 510
 511     Double_t b[2] = {-99., -99.};
 512     Double_t bCov[3] = {-99., -99., -99.};
 513     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 514       continue;
 515
 516     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 517       continue;
 518
 519     Qx2 += TMath::Cos(2*aodTrack->Phi());
 520     Qy2 += TMath::Sin(2*aodTrack->Phi());
 521     Qx3 += TMath::Cos(3*aodTrack->Phi());
 522     Qy3 += TMath::Sin(3*aodTrack->Phi());
 523
 524   }
 525
 526   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 527
 528   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 529   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 530   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 531   PIDResponse->SetTOFResponse(aodEvent,AliPIDResponse::kTOF_T0);
 532   PIDResponse->GetTOFResponse().SetTOFtailAllPara(-3,1.1);
 533
 534 //   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
 535 //   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
 536 //   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
 537 //   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
 538
 539   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 540
 541   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 542   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 543   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 544   Double_t nSigmaTPCRef,nSigmaTOFRef=6,nSigmaTPC2Ref,nSigmaTOF2Ref=6,nSigmaCombRef,nSigmaComb2Ref;
 545
 546   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 547   TClonesArray *mcArray = NULL;
 548   if (mcHeader)
 549     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 550
 551   Int_t nmc = 0;
 552   if(mcArray)
 553     nmc = mcArray->GetEntries();
 554
 555   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 556     AliAODTrack* track = aodEvent->GetTrack(i);
 557
 558     AliAODMCParticle *mcp = NULL;
 559     Int_t pdg = 0;
 560
 561     if (!track){
 562       continue;
 563     }
 564
 565     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 566
 567     Int_t label = -1;
 568     if(mcArray){
 569       label = track->GetLabel();
 570       if(label != -1 && label < nmc){
 571         label = TMath::Abs(label);
 572         mcp = (AliAODMCParticle*)mcArray->At(label);
 573         pdg = TMath::Abs(mcp->GetPdgCode());
 574       }
 575       else
 576         label = -1;
 577     }
 578     else{
 579       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 580     }
 581
 582     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 583       hTracking[0]->Fill(track->Pt(),fCentrality);
 584       if(pdg == 211)
 585         hTracking[1]->Fill(track->Pt(),fCentrality);
 586       else if(pdg == 321)
 587         hTracking[2]->Fill(track->Pt(),fCentrality);
 588       else if(pdg == 2212)
 589         hTracking[3]->Fill(track->Pt(),fCentrality);
 590       else if(! mcp){ // Fill matching histo with the prob
 591         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 592         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 593         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 594       }
 595     }
 596
 597     if(!tofMatch) continue;
 598
 599     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 600       hMatching[0]->Fill(track->Pt(),fCentrality);
 601       if(pdg == 211)
 602         hMatching[1]->Fill(track->Pt(),fCentrality);
 603       else if(pdg == 321)
 604         hMatching[2]->Fill(track->Pt(),fCentrality);
 605       else if(pdg == 2212)
 606         hMatching[3]->Fill(track->Pt(),fCentrality);
 607       else if(! mcp){ // Fill matching histo with the prob
 608         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 609         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 610         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 611       }
 612     }
 613   }
 614
 615 //   Int_t pdg1 = -1;
 616 //   Int_t pdg2 = -1;
 617
 618
 619   // start analysis phi
 620   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 621     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 622
 623     if (!KpTrack){
 624       continue;
 625     }
 626
 627     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 628
 629     nSigmaTOF=5;
 630     nSigmaTOFRef=5;
 631     nSigmaComb=5;
 632     nSigmaCombRef=5;
 633
 634     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 635     fPidKp=0;
 636
 637     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 638
 639     Double_t oldpP[10];
 640     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 641
 642     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
 643     nSigmaTPCRef = PIDResponse->NumberOfSigmasTPC(KpTrack,(AliPID::EParticleType) fSpeciesRef);
 644
 645     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 646
 647     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 648     /*
 649     if(mcArray){
 650       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 651       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 652       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 653     }
 654     */
 655
 656     fPidKp = Int_t(probP[3]*100);
 657
 658     if(tofMatch1){
 659       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 660         // remove this tof hit
 661         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 662         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 663         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 664         fPidKp = Int_t(probP[4]*100);
 665         tofMatch1=0;
 666       }
 667       else{
 668         if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
 669
 670         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 671         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 672         if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
 673         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 674         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 675         if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
 676         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 677         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 678         if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
 679         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 680         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 681         if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
 682
 683         nSigmaTOFRef = PIDResponse->NumberOfSigmasTOF(KpTrack,(AliPID::EParticleType) fSpeciesRef);
 684
 685         if(fIsMC){
 686           Float_t mismAdd = addMismatchForMC;
 687           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 688
 689           if(gRandom->Rndm() < mismAdd){
 690             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 691             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 692             channel = channel % 8736;
 693             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 694
 695             // generate random time
 696             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 697             Double_t times[AliPID::kSPECIESC];
 698             KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 699             nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
 700             nSigmaTOFRef = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[fSpeciesRef], (AliPID::EParticleType) fSpeciesRef);
 701           }
 702         }
 703
 704         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 705         nSigmaTOF = TMath::Abs(nSigmaTOF);
 706
 707         if(nSigmaTOF < 2) fPidKp += 256;
 708         else if(nSigmaTOF < 3) fPidKp += 512;
 709       }
 710     }
 711
 712     if(tofMatch1){
 713       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 714       nSigmaCombRef = TMath::Sqrt(nSigmaTOFRef*nSigmaTOFRef + nSigmaTPCRef*nSigmaTPCRef);
 715       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 716         fCombsignal->Fill(nSigmaComb);
 717       }
 718     } else {
 719       nSigmaComb = TMath::Abs(nSigmaTPC);
 720       nSigmaCombRef = TMath::Abs(nSigmaTPCRef);
 721     }
 722
 723     // use sigmaTOF instead of sigmaComb
 724     nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
 725
 726     if(tofMatch1){
 727       nSigmaComb = nSigmaTOF;
 728       nSigmaCombRef = nSigmaTOFRef;
 729     }
 730
 731     if(nSigmaComb < 2) nSigmaComb = 2;
 732     else if(nSigmaComb < 3) nSigmaComb = 3;
 733     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 734     else nSigmaComb = 6;
 735
 736     if(nSigmaCombRef < 2) nSigmaCombRef = 2;
 737     else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
 738     else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
 739     else nSigmaCombRef = 6;
 740
 741     if(fPtKp > 4.299) fPtKp = 4.299;
 742
 743     for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 744       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 745
 746       if (!KnTrack){
 747         continue;
 748       }
 749
 750       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 751
 752       px = KpTrack->Px() + KnTrack->Px();
 753       py = KpTrack->Py() + KnTrack->Py();
 754       pz = KpTrack->Pz() + KnTrack->Pz();
 755       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
 756       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
 757
 758       phiV.SetPxPyPzE(px,py,pz,E);
 759       fMassV0 = phiV.M();
 760
 761       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 762
 763       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 764       fPidKn=0;
 765
 766       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 767       Double_t oldpN[10];
 768       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 769
 770       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
 771       nSigmaTPC2Ref = PIDResponse->NumberOfSigmasTPC(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 772
 773       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 774
 775       nSigmaTOF2=5;
 776       nSigmaTOF2Ref=5;
 777       nSigmaComb2=5;
 778       nSigmaComb2Ref=5;
 779
 780       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 781       /*
 782       if(mcArray){
 783         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 784         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 785         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 786       }
 787       */
 788
 789       fPidKn = Int_t(probN[3]*100);
 790
 791       if(tofMatch2){
 792         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 793           // remove this tof hit
 794           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 795           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 796           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 797           fPidKn = Int_t(probN[4]*100);
 798           tofMatch2=0;
 799         }
 800         else{
 801           if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
 802
 803           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
 804           nSigmaTOF2Ref = PIDResponse->NumberOfSigmasTOF(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 805
 806           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 807           nSigmaTOF2Ref = TMath::Abs(nSigmaTOF2Ref);
 808
 809           if(fIsMC){
 810             Float_t mismAdd = addMismatchForMC;
 811             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 812
 813             if(gRandom->Rndm() < mismAdd){
 814               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 815               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 816               channel = channel % 8736;
 817               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 818
 819               // generate random time
 820               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 821               Double_t times[AliPID::kSPECIESC];
 822               KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 823               nSigmaTOF2 = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
 824               nSigmaTOF2Ref = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[fSpeciesRef], (AliPID::EParticleType) fSpeciesRef);
 825             }
 826           }
 827
 828           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 829
 830           if(nSigmaTOF2 < 2) fPidKn += 256;
 831           else if(nSigmaTOF2 < 3) fPidKn += 512;
 832         }
 833       }
 834
 835       fPtPhi = phiV.Pt();
 836       fEtaPhi = phiV.Eta();
 837       fPhiPhi = phiV.Phi();
 838
 839       if(tofMatch2){
 840         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 841         nSigmaComb2Ref = TMath::Sqrt(nSigmaTOF2Ref*nSigmaTOF2Ref+ nSigmaTPC2Ref*nSigmaTPC2Ref);
 842         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 843           fCombsignal->Fill(nSigmaComb2);
 844         }
 845       } else {
 846         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 847         nSigmaComb2Ref = TMath::Abs(nSigmaTPC2Ref);
 848       }
 849
 850       // use sigmaTOF instead of sigmaComb
 851       if(tofMatch2){
 852         nSigmaComb2 = nSigmaTOF2;
 853         nSigmaComb2Ref = nSigmaTOF2Ref;
 854       }
 855
 856       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 857       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 858       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 859       else nSigmaComb2 = 6;
 860
 861       if(nSigmaComb2Ref < 2) nSigmaComb2Ref = 2;
 862       else if(nSigmaComb2Ref < 3) nSigmaComb2Ref = 3;
 863       else if(nSigmaComb2Ref < 5) nSigmaComb2Ref = 4.99;
 864       else nSigmaComb2Ref = 6;
 865
 866       Bool_t isTrue = kFALSE;
 867
 868       if(mcArray){
 869         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 870         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 871
 872         if(labelP > -1 && labelN > -1){
 873           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 874           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 875
 876           Int_t mP = partP->GetMother();
 877           Int_t mN = partN->GetMother();
 878           if(mP == mN && mP > -1){
 879             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 880             Int_t pdgM = partM->GetPdgCode();
 881             if(pdgM == 333) isTrue = kTRUE;
 882           }
 883         }
 884
 885       }
 886
 887       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 888       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 889
 890       if(gRandom->Rndm() < 0.5){
 891         deltaphi1 += TMath::Pi();
 892         deltaphi2 += TMath::Pi();
 893       }
 894
 895       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 896       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 897       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 898       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 899
 900       if(fPtKn > 4.299) fPtKn = 4.299;
 901
 902       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)};
 903       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)};
 904
 905       Int_t ipt = 0;
 906       while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
 907         ipt++;
 908       }
 909       ipt--;
 910       if(ipt < 0) ipt = 0; // just to be sure
 911
 912       if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 913         if(fSpeciesRef != 3){
 914           xTOfill[4] = probP[fSpeciesRef];
 915           xTOfill2[5] = probN[fSpeciesRef];
 916
 917           xTOfill[9] = nSigmaCombRef;
 918           xTOfill2[10] = nSigmaComb2Ref;
 919         }
 920
 921
 922         if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
 923         xTOfill[1] = KnTrack->Eta();
 924         if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);// use tagging on positive track
 925
 926         if(fPIDuserCut){
 927           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)};
 928           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)};
 929
 930           if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
 931             xUser[3] = 1;
 932           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 933             xUser[3] = 2;
 934           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 935             xUser[3] = 3;
 936           }
 937           if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
 938             xUser2[3] = 1;
 939           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
 940             xUser2[3] = 2;
 941           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
 942             xUser2[3] = 3;
 943           }
 944           if((fPidKn%128) > 80) fContUser->Fill(0,fMassV0,fCentrality,xUser);
 945           if((fPidKp%128) > 80) fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
 946         }
 947
 948       }
 949     }
 950   } // end analysi phi
 951
 952 }
 953
 954 //_____________________________________________________________________________
 955 Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
 956 {
 957
 958   Float_t zvtx = -999;
 959
 960   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 961   if (!vtxAOD)
 962     return zvtx;
 963   if(vtxAOD->GetNContributors()>0)
 964     zvtx = vtxAOD->GetZ();
 965
 966   return zvtx;
 967 }
 968 //_____________________________________________________________________________
 969 void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
 970 {
 971   // Terminate loop
 972   Printf("Terminate()");
 973 }
 974 //-------------------------------------------------------------------------
 975 Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
 976   if(!fIsMC) return 1; // used only on MC
 977
 978   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
 979     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 980
 981     if(!(channel%20)) return 0; // 5% additional loss in MC
 982   }
 983
 984   return 1;
 985 }