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