PWGPP/pid/AliAnalysisTaskK0sBayes.cxx

   1 #include "AliAnalysisTaskK0sBayes.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 AliAnalysisTaskK0sBayes::fPtKsMin[AliAnalysisTaskK0sBayes::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 AliAnalysisTaskK0sBayes::fPtKsMax[AliAnalysisTaskK0sBayes::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(AliAnalysisTaskK0sBayes)
  39 //_____________________________________________________________________________
  40 AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes():
  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   fPtKs(0.),
  55   fPhiKs(0.),
  56   fEtaKs(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   fNK0s(0),
  75   fNpiPos(0),
  76   fNpiNeg(0),
  77   fHmismTOF(0),
  78   fHchannelTOFdistr(0),
  79   fTypeCol(2),
  80   fPIDuserCut(NULL)
  81 {
  82   // Default constructor (should not be used)
  83   fList->SetName("contKsBayes1");
  84   fList2->SetName("contKsBayes2");
  85   fList3->SetName("contKsBayes3");
  86
  87   fList->SetOwner(kTRUE);
  88   fList2->SetOwner(kTRUE);
  89   fList3->SetOwner(kTRUE);
  90
  91   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
  92   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
  93
  94   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
  95   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
  96
  97   for(Int_t i=0;i < nCentrBin;i++){
  98     fElTOF[i] = NULL;
  99     fElTPC[i] = NULL;
 100     fPiTOF[i] = NULL;
 101     fPiTPC[i] = NULL;
 102     fKaTOF[i] = NULL;
 103     fKaTPC[i] = NULL;
 104     fPrTOF[i] = NULL;
 105     fPrTPC[i] = NULL;
 106   }
 107   for(Int_t i=0;i < 4;i++){
 108     hMatching[i] = NULL;
 109     hTracking[i] = NULL;
 110   }
 111   for(Int_t i=0;i < 1000;i++){
 112     fPhiK0s[i] = 0.0;
 113     fPtK0s[i] = 0.0;
 114     fIPiPos[i] = 0;
 115     fIPiNeg[i] = 0;
 116     fIpiP[i] = 0;
 117     fIpiN[i] = 0;
 118     fMassKs[i] = 0.0;
 119   }
 120 }
 121
 122 //______________________________________________________________________________
 123 AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes(const char *name):
 124   AliAnalysisTaskSE(name),
 125   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
 126   fEtaCut(0.8),   // cut on |eta| < fEtaCut
 127   fMinPt(0.15),   // cut on pt > fMinPt
 128   fIsMC(kFALSE),
 129   fQAsw(kFALSE),
 130   fNcluster(70),
 131   fFilterBit(1),
 132   fList(new TList()),
 133   fList2(new TList()),
 134   fList3(new TList()),
 135   fCentrality(-1),
 136   fPsi(0),
 137   fPtKs(0.),
 138   fPhiKs(0.),
 139   fEtaKs(0.),
 140   fMassV0(0.),
 141   fPtKp(0.),
 142   fPhiKp(0.),
 143   fEtaKp(0.),
 144   fPtKn(0.),
 145   fPhiKn(0.),
 146   fEtaKn(0.),
 147   fPidKp(0),
 148   fPidKn(0),
 149   fTOFTPCsignal(0),
 150   fCombsignal(0),
 151   fCutsDaughter(NULL),
 152   fPIDCombined(NULL),
 153   fContPid(NULL),
 154   fContPid2(NULL),
 155   fContUser(NULL),
 156   fContUser2(NULL),
 157   fNK0s(0),
 158   fNpiPos(0),
 159   fNpiNeg(0),
 160   fHmismTOF(0),
 161   fHchannelTOFdistr(0),
 162   fTypeCol(2),
 163   fPIDuserCut(NULL)
 164 {
 165
 166   DefineOutput(1, TList::Class());
 167   DefineOutput(2, TList::Class());
 168   DefineOutput(3, TList::Class());
 169   DefineOutput(4, TList::Class());
 170
 171   // Output slot #1 writes into a TTree
 172   fList->SetName("contKsBayes1");
 173   fList2->SetName("contKsBayes2");
 174   fList3->SetName("contKsBayes3");
 175
 176   fList->SetOwner(kTRUE);
 177   fList2->SetOwner(kTRUE);
 178   fList3->SetOwner(kTRUE);
 179
 180   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
 181   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
 182
 183   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
 184   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
 185
 186   for(Int_t i=0;i < nCentrBin;i++){
 187     fElTOF[i] = NULL;
 188     fElTPC[i] = NULL;
 189     fPiTOF[i] = NULL;
 190     fPiTPC[i] = NULL;
 191     fKaTOF[i] = NULL;
 192     fKaTPC[i] = NULL;
 193     fPrTOF[i] = NULL;
 194     fPrTPC[i] = NULL;
 195   }
 196   for(Int_t i=0;i < 4;i++){
 197     hMatching[i] = NULL;
 198     hTracking[i] = NULL;
 199   }
 200   for(Int_t i=0;i < 1000;i++){
 201     fPhiK0s[i] = 0.0;
 202     fPtK0s[i] = 0.0;
 203     fIPiPos[i] = 0;
 204     fIPiNeg[i] = 0;
 205     fIpiP[i] = 0;
 206     fIpiN[i] = 0;
 207     fMassKs[i] = 0.0;
 208   }
 209 }
 210 //_____________________________________________________________________________
 211 AliAnalysisTaskK0sBayes::~AliAnalysisTaskK0sBayes()
 212 {
 213
 214 }
 215
 216 //______________________________________________________________________________
 217 void AliAnalysisTaskK0sBayes::UserCreateOutputObjects()
 218 {
 219
 220   fPIDCombined=new AliPIDCombined;
 221   fPIDCombined->SetDefaultTPCPriors();
 222   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
 223
 224   Float_t invMmin = 0.497-0.005*8;
 225   Float_t invMmax = 0.497+0.005*8;
 226
 227   const Int_t nBinPid = 14;
 228
 229   Int_t binPid[nBinPid] = {1/*ptKs*/,8/*EtaPiP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
 230
 231   Int_t binPid2[nBinPid] = {1/*ptKs*/,8/*EtaPiN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
 232
 233   fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
 234   fContPid->SetTitleX("M_{K^{0}_{s}}");
 235   fContPid->SetTitleY("centrality (%)");
 236   fContPid->SetVarName(0,"p_{T}^{K^{0}_{s}}");
 237   fContPid->SetVarRange(0,0,10);
 238   fContPid->SetVarName(1,"#eta^{K^{0}_{s}}");
 239   fContPid->SetVarRange(1,-0.8,0.8);
 240   fContPid->SetVarName(2,"p_{T}^{Kp}");
 241   fContPid->SetVarRange(2,0.3,4.3);
 242   fContPid->SetVarName(3,"p_{T}^{Kn}");
 243   fContPid->SetVarRange(3,0.3,4.3);
 244   fContPid->SetVarName(4,"BayesProb^{Kp}");
 245   fContPid->SetVarRange(4,0,1.);
 246   fContPid->SetVarName(5,"BayesProb^{Kn}");
 247   fContPid->SetVarRange(5,0,1.);
 248   fContPid->SetVarName(6,"isTOFmatch^{Kp}");
 249   fContPid->SetVarRange(6,-0.5,1.5);
 250   fContPid->SetVarName(7,"isTOFmatch^{Kn}");
 251   fContPid->SetVarRange(7,-0.5,1.5);
 252   fContPid->SetVarName(8,"isKsTrue^{Kn}");
 253   fContPid->SetVarRange(8,-0.5,1.5);
 254   fContPid->SetVarName(9,"N#sigma^{Kp}");
 255   fContPid->SetVarRange(9,1.25,6.25);
 256   fContPid->SetVarName(10,"N#sigma^{Kn}");
 257   fContPid->SetVarRange(10,1.25,6.25);
 258   fContPid->SetVarName(11,"#Delta#phi^{Kp}");
 259   fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
 260   fContPid->SetVarName(12,"#Delta#phi^{Kn}");
 261   fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
 262   fContPid->SetVarName(13,"#Psi");
 263   fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
 264
 265   fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
 266   fContPid2->SetTitleX("M_{K^{0}_{s}}");
 267   fContPid2->SetTitleY("centrality (%)");
 268   fContPid2->SetVarName(0,"p_{T}^{K^{0}_{s}}");
 269   fContPid2->SetVarRange(0,0,10);
 270   fContPid2->SetVarName(1,"#eta^{K^{0}_{s}}");
 271   fContPid2->SetVarRange(1,-0.8,0.8);
 272   fContPid2->SetVarName(2,"p_{T}^{Kp}");
 273   fContPid2->SetVarRange(2,0.3,4.3);
 274   fContPid2->SetVarName(3,"p_{T}^{Kn}");
 275   fContPid2->SetVarRange(3,0.3,4.3);
 276   fContPid2->SetVarName(4,"BayesProb^{Kp}");
 277   fContPid2->SetVarRange(4,0,1.);
 278   fContPid2->SetVarName(5,"BayesProb^{Kn}");
 279   fContPid2->SetVarRange(5,0,1.);
 280   fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
 281   fContPid2->SetVarRange(6,-0.5,1.5);
 282   fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
 283   fContPid2->SetVarRange(7,-0.5,1.5);
 284   fContPid2->SetVarName(8,"isKsTrue^{Kn}");
 285   fContPid2->SetVarRange(8,-0.5,1.5);
 286   fContPid2->SetVarName(9,"N#sigma^{Kp}");
 287   fContPid2->SetVarRange(9,1.25,6.25);
 288   fContPid2->SetVarName(10,"N#sigma^{Kn}");
 289   fContPid2->SetVarRange(10,1.25,6.25);
 290   fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
 291   fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
 292   fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
 293   fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
 294   fContPid2->SetVarName(13,"#Psi");
 295   fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
 296
 297   const Int_t nDETsignal = 100; // mass
 298   Double_t binDETsignal[nDETsignal+1];
 299   for(Int_t i=0;i<nDETsignal+1;i++){
 300     binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
 301   }
 302   const Int_t nDETsignal2 = 10; // centrality
 303   Double_t binDETsignal2[nDETsignal2+1];
 304   for(Int_t i=0;i<nDETsignal2+1;i++){
 305     binDETsignal2[i] = i*100./ nDETsignal2;
 306   }
 307   fContPid->AddSpecies("K0s",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 308   fContPid2->AddSpecies("K0s2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 309
 310   fList->Add(fContPid);
 311   fList->Add(fContPid2);
 312
 313   const Int_t nBinUser = 6;
 314   Int_t binUser[nBinUser] = {8/*Eta*/,20/*pt*/,2/*istrue*/,4/*whatSelection*/,1/*DeltaPhi*/,1/*Psi*/};
 315   fContUser = new AliPIDperfContainer("contUserPID",nBinUser,binUser);
 316   fContUser->SetTitleX("M_{K^{0}_{s}}");
 317   fContUser->SetTitleY("centrality (%)");
 318   fContUser->SetVarName(0,"#eta^{K^{0}_{s}}");
 319   fContUser->SetVarRange(0,-0.8,0.8);
 320   fContUser->SetVarName(1,"p_{T}");
 321   fContUser->SetVarRange(1,0.3,4.3);
 322   fContUser->SetVarName(2,"isKsTrue^{Kn}");
 323   fContUser->SetVarRange(2,-0.5,1.5);
 324   fContUser->SetVarName(3,"whatSelected"); // 0=no, 1=pi, 2=K, 3=p
 325   fContUser->SetVarRange(3,-0.5,3.5);
 326   fContUser->SetVarName(4,"#Delta#phi");
 327   fContUser->SetVarRange(4,-TMath::Pi(),TMath::Pi());
 328   fContUser->SetVarName(5,"#Psi");
 329   fContUser->SetVarRange(5,-TMath::Pi()/2,TMath::Pi()/2);
 330
 331   fContUser2 = new AliPIDperfContainer("contUserPID2",nBinUser,binUser);
 332   fContUser2->SetTitleX("M_{K^{0}_{s}}");
 333   fContUser2->SetTitleY("centrality (%)");
 334   fContUser2->SetVarName(0,"#eta^{K^{0}_{s}}");
 335   fContUser2->SetVarRange(0,-0.8,0.8);
 336   fContUser2->SetVarName(1,"p_{T}");
 337   fContUser2->SetVarRange(1,0.3,4.3);
 338   fContUser2->SetVarName(2,"isKsTrue^{Kn}");
 339   fContUser2->SetVarRange(2,-0.5,1.5);
 340   fContUser2->SetVarName(3,"whatSelected");
 341   fContUser2->SetVarRange(3,-0.5,3.5);
 342   fContUser2->SetVarName(4,"#Delta#phi");
 343   fContUser2->SetVarRange(4,-TMath::Pi(),TMath::Pi());
 344   fContUser2->SetVarName(5,"#Psi");
 345   fContUser2->SetVarRange(5,-TMath::Pi()/2,TMath::Pi()/2);
 346
 347   fContUser->AddSpecies("K0s",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 348   fContUser2->AddSpecies("K0s2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 349
 350   fList->Add(fContUser);
 351   fList->Add(fContUser2);
 352
 353   hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 354   hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 355   hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 356   hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 357
 358   hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 359   hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 360   hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 361   hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 362
 363   fList2->Add(hMatching[0]);
 364   fList2->Add(hMatching[1]);
 365   fList2->Add(hMatching[2]);
 366   fList2->Add(hMatching[3]);
 367   fList2->Add(hTracking[0]);
 368   fList2->Add(hTracking[1]);
 369   fList2->Add(hTracking[2]);
 370   fList2->Add(hTracking[3]);
 371
 372   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);
 373   fList2->Add(fTOFTPCsignal);
 374   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);
 375   fList2->Add(fCombsignal);
 376
 377   // QA plots
 378   char name[100],title[400];
 379   for(Int_t i=0;i < nCentrBin;i++){
 380     snprintf(name,100,"hPiTPCQA_%i",i);
 381     snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 382     fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 383     fList3->Add(fPiTPC[i]);
 384
 385     snprintf(name,100,"hKaTPCQA_%i",i);
 386     snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 387     fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 388     fList3->Add(fKaTPC[i]);
 389
 390     snprintf(name,100,"hPrTPCQA_%i",i);
 391     snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 392     fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 393     fList3->Add(fPrTPC[i]);
 394
 395     snprintf(name,100,"hElTPCQA_%i",i);
 396     snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 397     fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 398     fList3->Add(fElTPC[i]);
 399
 400     snprintf(name,100,"hPiTOFQA_%i",i);
 401     snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 402     fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 403     fList3->Add(fPiTOF[i]);
 404
 405     snprintf(name,100,"hKaTOFQA_%i",i);
 406     snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 407     fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 408     fList3->Add(fKaTOF[i]);
 409
 410     snprintf(name,100,"hPrTOFQA_%i",i);
 411     snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 412     fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 413     fList3->Add(fPrTOF[i]);
 414
 415     snprintf(name,100,"hElTOFQA_%i",i);
 416     snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 417     fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 418     fList3->Add(fElTOF[i]);
 419
 420   }
 421
 422   // Post output data.
 423   PostData(1, fList);
 424   PostData(2, fList2);
 425   PostData(3, fList3);
 426 }
 427
 428 //______________________________________________________________________________
 429 void AliAnalysisTaskK0sBayes::UserExec(Option_t *)
 430 {
 431     // Main loop
 432     // Called for each event
 433
 434     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 435     if(!fOutputAOD){
 436         Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 437         this->Dump();
 438         return;
 439     }
 440
 441     Float_t zvtx = GetVertex(fOutputAOD);
 442
 443
 444
 445     //Get the MC object
 446     if(fIsMC){
 447       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 448       if (!mcHeader) {
 449         AliError("Could not find MC Header in AOD");
 450         return;
 451       }
 452     }
 453
 454     if (TMath::Abs(zvtx) < fVtxCut) {
 455
 456       SelectK0s();
 457
 458       //Centrality
 459       Float_t v0Centr  = -10.;
 460       Float_t trkCentr  = -10.;
 461       AliCentrality *centrality = fOutputAOD->GetCentrality();
 462       if (centrality){
 463         trkCentr  = centrality->GetCentralityPercentile("V0M");
 464         v0Centr = centrality->GetCentralityPercentile("TRK");
 465       }
 466
 467       if(!fTypeCol){
 468         v0Centr=100./(fOutputAOD->GetNumberOfTracks()/12.+1);
 469         trkCentr=v0Centr;
 470       }
 471
 472       if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
 473         fCentrality = v0Centr;
 474         Analyze(fOutputAOD); // Do analysis!!!!
 475
 476       }
 477     }
 478
 479 }
 480
 481 //________________________________________________________________________
 482 void AliAnalysisTaskK0sBayes::Analyze(AliAODEvent* aodEvent)
 483 {
 484
 485   Int_t ntrack = aodEvent->GetNumberOfTracks();
 486
 487   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
 488   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
 489   fPidKp=0,fPidKn=0;
 490   fMassV0=-1;
 491
 492   TLorentzVector ksV;
 493
 494   Double_t px,py,pz,E;
 495
 496   Float_t invMmin = 0.497-0.005*8;
 497   Float_t invMmax = 0.497+0.005*8;
 498
 499   Int_t icentr = 8;
 500   if(fCentrality < 0) icentr = 8;
 501   else if(fCentrality < 10) icentr = 0;
 502   else if(fCentrality < 20) icentr = 1;
 503   else if(fCentrality < 30) icentr = 2;
 504   else if(fCentrality < 40) icentr = 3;
 505   else if(fCentrality < 50) icentr = 4;
 506   else if(fCentrality < 60) icentr = 5;
 507   else if(fCentrality < 70) icentr = 6;
 508   else if(fCentrality < 80) icentr = 7;
 509
 510   Float_t addMismatchForMC = 0.005;
 511   if(fCentrality < 50) addMismatchForMC += 0.005;
 512   if(fCentrality < 20) addMismatchForMC += 0.02;
 513
 514   if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
 515
 516   fPsi = 0;
 517   /* Compute TPC EP */
 518   Double_t Qx2 = 0, Qy2 = 0;
 519   Double_t Qx3 = 0, Qy3 = 0;
 520   for(Int_t iT = 0; iT < ntrack; iT++) {
 521     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
 522
 523     if (!aodTrack){
 524       continue;
 525     }
 526
 527     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 528
 529     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 530       continue;
 531
 532     Double_t b[2] = {-99., -99.};
 533     Double_t bCov[3] = {-99., -99., -99.};
 534     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 535       continue;
 536
 537     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 538       continue;
 539
 540     Qx2 += TMath::Cos(2*aodTrack->Phi());
 541     Qy2 += TMath::Sin(2*aodTrack->Phi());
 542     Qx3 += TMath::Cos(3*aodTrack->Phi());
 543     Qy3 += TMath::Sin(3*aodTrack->Phi());
 544
 545   }
 546
 547   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 548
 549   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 550   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 551   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 552
 553   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 554
 555   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 556   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 557   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 558
 559
 560   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 561   TClonesArray *mcArray = NULL;
 562   if (mcHeader)
 563     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 564
 565   Int_t nmc = 0;
 566   if(mcArray)
 567     nmc = mcArray->GetEntries();
 568
 569   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 570     AliAODTrack* track = aodEvent->GetTrack(i);
 571
 572     AliAODMCParticle *mcp = NULL;
 573     Int_t pdg = 0;
 574
 575     if (!track){
 576       continue;
 577     }
 578
 579     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 580
 581     Int_t label = -1;
 582     if(mcArray){
 583       label = track->GetLabel();
 584       if(label != -1 && label < nmc){
 585         label = TMath::Abs(label);
 586         mcp = (AliAODMCParticle*)mcArray->At(label);
 587         pdg = TMath::Abs(mcp->GetPdgCode());
 588       }
 589       else
 590         label = -1;
 591     }
 592     else{
 593       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 594     }
 595
 596     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 597       hTracking[0]->Fill(track->Pt(),fCentrality);
 598       if(pdg == 211)
 599         hTracking[1]->Fill(track->Pt(),fCentrality);
 600       else if(pdg == 321)
 601         hTracking[2]->Fill(track->Pt(),fCentrality);
 602       else if(pdg == 2212)
 603         hTracking[3]->Fill(track->Pt(),fCentrality);
 604       else if(! mcp){ // Fill matching histo with the prob
 605         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 606         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 607         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 608       }
 609     }
 610
 611     if(!tofMatch) continue;
 612
 613     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 614       hMatching[0]->Fill(track->Pt(),fCentrality);
 615       if(pdg == 211)
 616         hMatching[1]->Fill(track->Pt(),fCentrality);
 617       else if(pdg == 321)
 618         hMatching[2]->Fill(track->Pt(),fCentrality);
 619       else if(pdg == 2212)
 620         hMatching[3]->Fill(track->Pt(),fCentrality);
 621       else if(! mcp){ // Fill matching histo with the prob
 622         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 623         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 624         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 625       }
 626     }
 627   }
 628
 629 //   Int_t pdg1 = -1;
 630 //   Int_t pdg2 = -1;
 631
 632
 633   // start analysis K0s
 634   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 635     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 636
 637     if (!KpTrack){
 638       continue;
 639     }
 640
 641     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 642
 643     nSigmaComb=5;
 644     nSigmaTOF = 5;
 645     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 646     fPidKp=0;
 647
 648     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 649
 650     Double_t oldpP[10];
 651     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 652
 653     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
 654
 655     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 656
 657     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 658
 659     /*
 660     if(mcArray){
 661       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 662       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 663       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 664     }
 665     */
 666
 667     fPidKp = Int_t(probP[2]*100);
 668
 669     if(tofMatch1){
 670       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 671         // remove this tof hit
 672         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 673         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 674         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 675         fPidKp = Int_t(probP[4]*100);
 676         tofMatch1=0;
 677       }
 678       else{
 679         if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
 680
 681         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 682         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 683         if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
 684         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 685         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 686         if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
 687         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 688         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 689         if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
 690         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 691         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 692         if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
 693
 694         if(fIsMC){
 695           Float_t mismAdd = addMismatchForMC;
 696           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 697
 698           if(gRandom->Rndm() < mismAdd){
 699             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 700             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 701             channel = channel % 8736;
 702             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 703
 704             // generate random time
 705             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 706             Double_t times[10];
 707             KpTrack->GetIntegratedTimes(times);
 708             nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kPion);
 709           }
 710         }
 711
 712         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 713         nSigmaTOF = TMath::Abs(nSigmaTOF);
 714
 715         if(nSigmaTOF < 2) fPidKp += 256;
 716         else if(nSigmaTOF < 3) fPidKp += 512;
 717       }
 718     }
 719
 720     if(tofMatch1){
 721       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 722       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 723         fCombsignal->Fill(nSigmaComb);
 724       }
 725     } else {
 726       nSigmaComb = TMath::Abs(nSigmaTPC);
 727     }
 728
 729     // use sigmaTOF instead of sigmaComb
 730     if(tofMatch1) nSigmaComb = nSigmaTOF;
 731
 732     if(nSigmaComb < 2) nSigmaComb = 2;
 733     else if(nSigmaComb < 3) nSigmaComb = 3;
 734     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 735     else nSigmaComb = 6;
 736
 737     Int_t iks=-1;
 738     for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
 739       if(i == fIpiP[k]) iks = k;
 740     }
 741
 742     if(fPtKp > 4.299) fPtKp = 4.299;
 743
 744     if(iks > -1 && fIpiN[iks] > -1){
 745       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 746       Int_t j = fIpiN[iks];
 747       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 748
 749       if (!KnTrack){
 750         continue;
 751       }
 752
 753       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 754
 755       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 756       fPidKn=0;
 757
 758       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 759       Double_t oldpN[10];
 760       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 761
 762       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 763
 764       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 765
 766       nSigmaComb2=5;
 767       nSigmaTOF2=5;
 768
 769       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 770       /*
 771       if(mcArray){
 772         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 773         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 774         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 775       }
 776       */
 777
 778       fPidKn = Int_t(probN[2]*100);
 779
 780       if(tofMatch2){
 781         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 782           // remove this tof hit
 783           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 784           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 785           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 786           fPidKn = Int_t(probN[4]*100);
 787           tofMatch2=0;
 788         }
 789         else{
 790           if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 791
 792           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 793
 794           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 795
 796           if(fIsMC){
 797             Float_t mismAdd = addMismatchForMC;
 798             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 799
 800             if(gRandom->Rndm() < mismAdd){
 801               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 802               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 803               channel = channel % 8736;
 804               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 805
 806               // generate random time
 807               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 808               Double_t times[10];
 809               KnTrack->GetIntegratedTimes(times);
 810               nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
 811             }
 812           }
 813
 814           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 815
 816           if(nSigmaTOF2 < 2) fPidKn += 256;
 817           else if(nSigmaTOF2 < 3) fPidKn += 512;
 818         }
 819       }
 820
 821       px = KpTrack->Px() + KnTrack->Px();
 822       py = KpTrack->Py() + KnTrack->Py();
 823       pz = KpTrack->Pz() + KnTrack->Pz();
 824       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
 825       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
 826
 827       ksV.SetPxPyPzE(px,py,pz,E);
 828       fMassV0 = fMassKs[iks];
 829
 830       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 831
 832
 833       fPtKs = ksV.Pt();
 834       fEtaKs = ksV.Eta();
 835       fPhiKs = ksV.Phi();
 836
 837       if(tofMatch2){
 838         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 839         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 840           fCombsignal->Fill(nSigmaComb2);
 841         }
 842       } else {
 843         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 844       }
 845
 846       // use sigmaTOF instead of sigmaComb
 847       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 848
 849       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 850       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 851       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 852       else nSigmaComb2 = 6;
 853
 854       Bool_t isTrue = kFALSE;
 855
 856       if(mcArray){
 857         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 858         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 859
 860         if(labelP > -1 && labelN > -1){
 861           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 862           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 863
 864           Int_t mP = partP->GetMother();
 865           Int_t mN = partN->GetMother();
 866           if(mP == mN && mP > -1){
 867             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 868             Int_t pdgM = partM->GetPdgCode();
 869             if(pdgM == 310) isTrue = kTRUE;
 870           }
 871         }
 872
 873       }
 874
 875       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 876       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 877
 878       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 879       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 880       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 881       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 882
 883       if(fPtKn > 4.299) fPtKn = 4.299;
 884
 885       Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 886
 887
 888       Int_t ipt = 0;
 889       while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
 890         ipt++;
 891       }
 892       ipt--;
 893       if(ipt < 0) ipt = 0; // just to be sure
 894
 895       if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 896         fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 897         xTOfill[1] = KnTrack->Eta();
 898         fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
 899
 900         if(fPIDuserCut){
 901           Float_t xUser[] = {KpTrack->Eta(),fPtKp,isTrue,0,deltaphi1,fPsi};
 902           Float_t xUser2[] = {KnTrack->Eta(),fPtKn,isTrue,0,deltaphi2,fPsi};
 903
 904           if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
 905             xUser[3] = 1;
 906           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 907             xUser[3] = 2;
 908           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 909             xUser[3] = 3;
 910           }
 911           if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
 912             xUser2[3] = 1;
 913           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
 914             xUser2[3] = 2;
 915           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
 916             xUser2[3] = 3;
 917           }
 918           fContUser->Fill(0,fMassV0,fCentrality,xUser);
 919           fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
 920         }
 921
 922       }
 923
 924
 925     }
 926   } // end analysi K0s
 927 }
 928
 929 //_____________________________________________________________________________
 930 Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
 931 {
 932
 933   Float_t zvtx = -999;
 934
 935   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 936   if (!vtxAOD)
 937     return zvtx;
 938   if(vtxAOD->GetNContributors()>0)
 939     zvtx = vtxAOD->GetZ();
 940
 941   return zvtx;
 942 }
 943 //_____________________________________________________________________________
 944 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
 945 {
 946   // Terminate loop
 947   Printf("Terminate()");
 948 }
 949 //=======================================================================
 950 void AliAnalysisTaskK0sBayes::SelectK0s(){
 951   fNK0s=0;
 952   fNpiPos=0;
 953   fNpiNeg=0;
 954
 955   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
 956   AliAODv0 *myV0;
 957   Double_t dMASS=0.0;
 958   for (Int_t i=0; i!=nV0s; ++i) {
 959     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
 960     if(!myV0) continue;
 961     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
 962     Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
 963     if(pass) {
 964       dMASS = myV0->MassK0Short();
 965       Float_t massLambda = myV0->MassLambda();
 966       Float_t massAntiLambda = myV0->MassAntiLambda();
 967
 968       if(TMath::Abs(dMASS-0.497)/0.005 < 8 && TMath::Abs(massLambda-1.115)/0.005 > 8 && TMath::Abs(massAntiLambda-1.115)/0.005 > 8){
 969         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
 970         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
 971         if(iT->Charge()<0){
 972           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
 973           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
 974         }
 975         fPhiK0s[fNK0s] = myV0->Phi();
 976         fPtK0s[fNK0s] = myV0->Pt();
 977         fIpiP[fNK0s] = FindDaugheterIndex(iT);
 978         fIpiN[fNK0s] = FindDaugheterIndex(jT);
 979         fMassKs[fNK0s] = dMASS;
 980         if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
 981           fNK0s++;
 982       }
 983     }
 984   }
 985
 986   /* My V0 code
 987   // fill pion stacks
 988   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
 989   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
 990     AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
 991
 992     if (!aodTrack){
 993       continue;
 994     }
 995
 996     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 997
 998     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
 999       continue;
1000     }
1001
1002     Double_t b[2] = {-99., -99.};
1003     Double_t bCov[3] = {-99., -99., -99.};
1004     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1005       continue;
1006
1007     if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
1008
1009
1010     Int_t charge = aodTrack->Charge();
1011     if(charge > 0){
1012       fIPiPos[fNpiPos] = iT;
1013       fNpiPos++;
1014     }
1015     else{
1016       fIPiNeg[fNpiNeg] = iT;
1017       fNpiNeg++;
1018     }
1019   }
1020
1021   for(Int_t i=0;i < fNpiPos;i++){
1022     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
1023     AliESDtrack pipE(pip);
1024
1025     for(Int_t j=0;j < fNpiNeg;j++){
1026       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
1027       AliESDtrack pinE(pin);
1028
1029       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
1030
1031       Double_t pPos[3];
1032       Double_t pNeg[3];
1033       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
1034       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
1035
1036       Float_t length = (xp+xn)*0.5;
1037
1038       Float_t pxs = pPos[0] + pNeg[0];
1039       Float_t pys = pPos[1] + pNeg[1];
1040       Float_t pzs = pPos[2] + pNeg[2];
1041       Float_t es = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
1042
1043       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
1044       Float_t phi = TMath::ATan2(pys,pxs);
1045       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
1046
1047       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
1048
1049       if(mindist < 0.4&& length > 0.7 && length < 25){
1050         Float_t esL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.938*0.938) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
1051         Float_t esAL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.938*0.938);
1052
1053         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
1054         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
1055
1056         if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
1057           fPhiK0s[fNK0s] = phi;
1058           fPtK0s[fNK0s] = pt;
1059           fIpiP[fNK0s] =fIPiPos[i] ;
1060           fIpiN[fNK0s] = fIPiNeg[j];
1061           fMassKs[fNK0s] = mass;
1062           fNK0s++;
1063         }
1064       }
1065     }
1066   }
1067   */
1068 }
1069
1070 //=======================================================================
1071 Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1072 {
1073   if (myV0->GetOnFlyStatus() ) return 0;
1074
1075   //the following is needed in order to evualuate track-quality
1076   AliAODTrack *iT, *jT;
1077   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1078   Double_t pos[3],cov[6];
1079   vV0s->GetXYZ(pos);
1080   vV0s->GetCovarianceMatrix(cov);
1081   const AliESDVertex vESD(pos,cov,100.,100);
1082
1083   // TESTING CHARGE
1084   int iPos, iNeg;
1085   iT=(AliAODTrack*) myV0->GetDaughter(0);
1086   if(iT->Charge()>0) {
1087     iPos = 0; iNeg = 1;
1088   } else {
1089     iPos = 1; iNeg = 0;
1090   }
1091   // END OF TEST
1092
1093   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1094
1095   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1096
1097   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1098   if(!trkFlag) return 0;
1099   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1100   if(!trkFlag2) return 0;
1101
1102   Double_t pvertex[3];
1103   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1104   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1105   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1106
1107   Double_t dDL=myV0->DecayLengthV0( pvertex );
1108   if(dDL  < 0.5 || dDL > 25) return 0;
1109
1110   Double_t dDCA=myV0->DcaV0Daughters();
1111   if(dDCA >0.5) return 0;
1112
1113   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1114   if(dCTP < -1) return 0;
1115
1116 //   AliESDtrack ieT( iT );
1117 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1118 //   if(TMath::Abs(dD0P) < 0]) return 0;
1119
1120 //   AliESDtrack jeT( jT );
1121 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1122 //   if(TMath::Abs(dD0M) < 0) return 0;
1123
1124 //   Double_t dD0D0=dD0P*dD0M;
1125 //   if(dD0D0>0) return 0;
1126
1127 //   Double_t dETA=myV0->Eta();
1128 //   if(dETA <-0.8) return 0;
1129 //   if(dETA >0.8) return 0;
1130
1131 //   Double_t dQT=myV0->PtArmV0();
1132 //   if(specie==0) if(dQT<???) return 0;
1133
1134   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1135   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1136
1137   if(specie==1 && dALPHA<0) return 2; // antilambda
1138   return 1; // K0s or lambda
1139 }
1140 //-------------------------------------------------------------------------
1141 Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1142   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1143
1144   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1145     AliAODTrack* track = fOutputAOD->GetTrack(i);
1146     if(track == trk) return i;
1147   }
1148
1149   printf("Daughter for %p not found\n",trk);
1150   return -1;
1151 }
1152 //-------------------------------------------------------------------------
1153 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1154   if(!fIsMC) return 1; // used only on MC
1155
1156   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1157     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1158
1159     if(!(channel%20)) return 0; // 5% additional loss in MC
1160   }
1161
1162   return 1;
1163 }