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