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     if(tofMatch1){
 723       nSigmaComb = nSigmaTOF;
 724       nSigmaCombRef = nSigmaTOFRef;
 725     }
 726
 727     if(nSigmaComb < 2) nSigmaComb = 2;
 728     else if(nSigmaComb < 3) nSigmaComb = 3;
 729     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 730     else nSigmaComb = 6;
 731
 732     if(nSigmaCombRef < 2) nSigmaCombRef = 2;
 733     else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
 734     else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
 735     else nSigmaCombRef = 6;
 736
 737     if(fPtKp > 4.299) fPtKp = 4.299;
 738
 739     for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 740       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 741
 742       if (!KnTrack){
 743         continue;
 744       }
 745
 746       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 747
 748       px = KpTrack->Px() + KnTrack->Px();
 749       py = KpTrack->Py() + KnTrack->Py();
 750       pz = KpTrack->Pz() + KnTrack->Pz();
 751       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
 752       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
 753
 754       phiV.SetPxPyPzE(px,py,pz,E);
 755       fMassV0 = phiV.M();
 756
 757       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 758
 759       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 760       fPidKn=0;
 761
 762       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 763       Double_t oldpN[10];
 764       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 765
 766       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
 767       nSigmaTPC2Ref = PIDResponse->NumberOfSigmasTPC(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 768
 769       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 770
 771       nSigmaTOF2=5;
 772       nSigmaTOF2Ref=5;
 773       nSigmaComb2=5;
 774       nSigmaComb2Ref=5;
 775
 776       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 777       /*
 778       if(mcArray){
 779         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 780         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 781         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 782       }
 783       */
 784
 785       fPidKn = Int_t(probN[3]*100);
 786
 787       if(tofMatch2){
 788         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 789           // remove this tof hit
 790           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 791           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 792           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 793           fPidKn = Int_t(probN[4]*100);
 794           tofMatch2=0;
 795         }
 796         else{
 797           if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
 798
 799           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
 800           nSigmaTOF2Ref = PIDResponse->NumberOfSigmasTOF(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 801
 802           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 803           nSigmaTOF2Ref = TMath::Abs(nSigmaTOF2Ref);
 804
 805           if(fIsMC){
 806             Float_t mismAdd = addMismatchForMC;
 807             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 808
 809             if(gRandom->Rndm() < mismAdd){
 810               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 811               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 812               channel = channel % 8736;
 813               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 814
 815               // generate random time
 816               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 817               Double_t times[AliPID::kSPECIESC];
 818               KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 819               nSigmaTOF2 = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
 820               nSigmaTOF2Ref = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[fSpeciesRef], (AliPID::EParticleType) fSpeciesRef);
 821             }
 822           }
 823
 824           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 825
 826           if(nSigmaTOF2 < 2) fPidKn += 256;
 827           else if(nSigmaTOF2 < 3) fPidKn += 512;
 828         }
 829       }
 830
 831       fPtPhi = phiV.Pt();
 832       fEtaPhi = phiV.Eta();
 833       fPhiPhi = phiV.Phi();
 834
 835       if(tofMatch2){
 836         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 837         nSigmaComb2Ref = TMath::Sqrt(nSigmaTOF2Ref*nSigmaTOF2Ref+ nSigmaTPC2Ref*nSigmaTPC2Ref);
 838         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 839           fCombsignal->Fill(nSigmaComb2);
 840         }
 841       } else {
 842         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 843         nSigmaComb2Ref = TMath::Abs(nSigmaTPC2Ref);
 844       }
 845
 846       // use sigmaTOF instead of sigmaComb
 847       if(tofMatch2){
 848         nSigmaComb2 = nSigmaTOF2;
 849         nSigmaComb2Ref = nSigmaTOF2Ref;
 850       }
 851
 852       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 853       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 854       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 855       else nSigmaComb2 = 6;
 856
 857       if(nSigmaComb2Ref < 2) nSigmaComb2Ref = 2;
 858       else if(nSigmaComb2Ref < 3) nSigmaComb2Ref = 3;
 859       else if(nSigmaComb2Ref < 5) nSigmaComb2Ref = 4.99;
 860       else nSigmaComb2Ref = 6;
 861
 862       Bool_t isTrue = kFALSE;
 863
 864       if(mcArray){
 865         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 866         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 867
 868         if(labelP > -1 && labelN > -1){
 869           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 870           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 871
 872           Int_t mP = partP->GetMother();
 873           Int_t mN = partN->GetMother();
 874           if(mP == mN && mP > -1){
 875             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 876             Int_t pdgM = partM->GetPdgCode();
 877             if(pdgM == 333) isTrue = kTRUE;
 878           }
 879         }
 880
 881       }
 882
 883       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 884       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 885
 886       if(gRandom->Rndm() < 0.5){
 887         deltaphi1 += TMath::Pi();
 888         deltaphi2 += TMath::Pi();
 889       }
 890
 891       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 892       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 893       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 894       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 895
 896       if(fPtKn > 4.299) fPtKn = 4.299;
 897
 898       Float_t xTOfill[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 899       Float_t xTOfill2[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 900
 901       Int_t ipt = 0;
 902       while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
 903         ipt++;
 904       }
 905       ipt--;
 906       if(ipt < 0) ipt = 0; // just to be sure
 907
 908       if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 909         if(fSpeciesRef != 3){
 910           xTOfill[4] = probP[fSpeciesRef];
 911           xTOfill2[5] = probN[fSpeciesRef];
 912
 913           xTOfill[9] = nSigmaCombRef;
 914           xTOfill2[10] = nSigmaComb2Ref;
 915         }
 916
 917
 918         if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
 919         xTOfill[1] = KnTrack->Eta();
 920         if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);// use tagging on positive track
 921
 922         if(fPIDuserCut){
 923           Float_t xUser[] = {KpTrack->Eta(),fPtKp,isTrue,0,deltaphi1,fPsi};
 924           Float_t xUser2[] = {KnTrack->Eta(),fPtKn,isTrue,0,deltaphi2,fPsi};
 925
 926           if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
 927             xUser[3] = 1;
 928           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 929             xUser[3] = 2;
 930           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 931             xUser[3] = 3;
 932           }
 933           if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
 934             xUser2[3] = 1;
 935           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
 936             xUser2[3] = 2;
 937           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
 938             xUser2[3] = 3;
 939           }
 940           if((fPidKn%128) > 80) fContUser->Fill(0,fMassV0,fCentrality,xUser);
 941           if((fPidKp%128) > 80) fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
 942         }
 943
 944       }
 945     }
 946   } // end analysi phi
 947
 948 }
 949
 950 //_____________________________________________________________________________
 951 Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
 952 {
 953
 954   Float_t zvtx = -999;
 955
 956   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 957   if (!vtxAOD)
 958     return zvtx;
 959   if(vtxAOD->GetNContributors()>0)
 960     zvtx = vtxAOD->GetZ();
 961
 962   return zvtx;
 963 }
 964 //_____________________________________________________________________________
 965 void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
 966 {
 967   // Terminate loop
 968   Printf("Terminate()");
 969 }
 970 //-------------------------------------------------------------------------
 971 Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
 972   if(!fIsMC) return 1; // used only on MC
 973
 974   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
 975     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 976
 977     if(!(channel%20)) return 0; // 5% additional loss in MC
 978   }
 979
 980   return 1;
 981 }