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