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       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 881       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 882       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 883       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 884
 885       if(fPtKn > 4.299) fPtKn = 4.299;
 886
 887       Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 888
 889
 890       Int_t ipt = 0;
 891       while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
 892         ipt++;
 893       }
 894       ipt--;
 895       if(ipt < 0) ipt = 0; // just to be sure
 896
 897       if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 898         fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 899         xTOfill[1] = KnTrack->Eta();
 900         fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
 901
 902         if(fPIDuserCut){
 903           Float_t xUser[] = {KpTrack->Eta(),fPtKp,isTrue,0,deltaphi1,fPsi};
 904           Float_t xUser2[] = {KnTrack->Eta(),fPtKn,isTrue,0,deltaphi2,fPsi};
 905
 906           if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
 907             xUser[3] = 1;
 908           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 909             xUser[3] = 2;
 910           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 911             xUser[3] = 3;
 912           }
 913           if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
 914             xUser2[3] = 1;
 915           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
 916             xUser2[3] = 2;
 917           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
 918             xUser2[3] = 3;
 919           }
 920           fContUser->Fill(0,fMassV0,fCentrality,xUser);
 921           fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
 922         }
 923
 924       }
 925
 926
 927     }
 928   } // end analysi K0s
 929 }
 930
 931 //_____________________________________________________________________________
 932 Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
 933 {
 934
 935   Float_t zvtx = -999;
 936
 937   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 938   if (!vtxAOD)
 939     return zvtx;
 940   if(vtxAOD->GetNContributors()>0)
 941     zvtx = vtxAOD->GetZ();
 942
 943   return zvtx;
 944 }
 945 //_____________________________________________________________________________
 946 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
 947 {
 948   // Terminate loop
 949   Printf("Terminate()");
 950 }
 951 //=======================================================================
 952 void AliAnalysisTaskK0sBayes::SelectK0s(){
 953   fNK0s=0;
 954   fNpiPos=0;
 955   fNpiNeg=0;
 956
 957   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
 958   AliAODv0 *myV0;
 959   Double_t dMASS=0.0;
 960   for (Int_t i=0; i!=nV0s; ++i) {
 961     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
 962     if(!myV0) continue;
 963     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
 964     Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
 965     if(pass) {
 966       dMASS = myV0->MassK0Short();
 967       Float_t massLambda = myV0->MassLambda();
 968       Float_t massAntiLambda = myV0->MassAntiLambda();
 969
 970       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){
 971         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
 972         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
 973         if(iT->Charge()<0){
 974           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
 975           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
 976         }
 977         fPhiK0s[fNK0s] = myV0->Phi();
 978         fPtK0s[fNK0s] = myV0->Pt();
 979         fIpiP[fNK0s] = FindDaugheterIndex(iT);
 980         fIpiN[fNK0s] = FindDaugheterIndex(jT);
 981         fMassKs[fNK0s] = dMASS;
 982         if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
 983           fNK0s++;
 984       }
 985     }
 986   }
 987
 988   /* My V0 code
 989   // fill pion stacks
 990   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
 991   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
 992     AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
 993
 994     if (!aodTrack){
 995       continue;
 996     }
 997
 998     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 999
1000     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1001       continue;
1002     }
1003
1004     Double_t b[2] = {-99., -99.};
1005     Double_t bCov[3] = {-99., -99., -99.};
1006     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1007       continue;
1008
1009     if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
1010
1011
1012     Int_t charge = aodTrack->Charge();
1013     if(charge > 0){
1014       fIPiPos[fNpiPos] = iT;
1015       fNpiPos++;
1016     }
1017     else{
1018       fIPiNeg[fNpiNeg] = iT;
1019       fNpiNeg++;
1020     }
1021   }
1022
1023   for(Int_t i=0;i < fNpiPos;i++){
1024     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
1025     AliESDtrack pipE(pip);
1026
1027     for(Int_t j=0;j < fNpiNeg;j++){
1028       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
1029       AliESDtrack pinE(pin);
1030
1031       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
1032
1033       Double_t pPos[3];
1034       Double_t pNeg[3];
1035       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
1036       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
1037
1038       Float_t length = (xp+xn)*0.5;
1039
1040       Float_t pxs = pPos[0] + pNeg[0];
1041       Float_t pys = pPos[1] + pNeg[1];
1042       Float_t pzs = pPos[2] + pNeg[2];
1043       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);
1044
1045       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
1046       Float_t phi = TMath::ATan2(pys,pxs);
1047       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
1048
1049       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
1050
1051       if(mindist < 0.4&& length > 0.7 && length < 25){
1052         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);
1053         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);
1054
1055         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
1056         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
1057
1058         if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
1059           fPhiK0s[fNK0s] = phi;
1060           fPtK0s[fNK0s] = pt;
1061           fIpiP[fNK0s] =fIPiPos[i] ;
1062           fIpiN[fNK0s] = fIPiNeg[j];
1063           fMassKs[fNK0s] = mass;
1064           fNK0s++;
1065         }
1066       }
1067     }
1068   }
1069   */
1070 }
1071
1072 //=======================================================================
1073 Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1074 {
1075   if (myV0->GetOnFlyStatus() ) return 0;
1076
1077   //the following is needed in order to evualuate track-quality
1078   AliAODTrack *iT, *jT;
1079   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1080   Double_t pos[3],cov[6];
1081   vV0s->GetXYZ(pos);
1082   vV0s->GetCovarianceMatrix(cov);
1083   const AliESDVertex vESD(pos,cov,100.,100);
1084
1085   // TESTING CHARGE
1086   int iPos, iNeg;
1087   iT=(AliAODTrack*) myV0->GetDaughter(0);
1088   if(iT->Charge()>0) {
1089     iPos = 0; iNeg = 1;
1090   } else {
1091     iPos = 1; iNeg = 0;
1092   }
1093   // END OF TEST
1094
1095   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1096
1097   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1098
1099   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1100   if(!trkFlag) return 0;
1101   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1102   if(!trkFlag2) return 0;
1103
1104   Double_t pvertex[3];
1105   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1106   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1107   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1108
1109   Double_t dDL=myV0->DecayLengthV0( pvertex );
1110   if(dDL  < 0.5 || dDL > 25) return 0;
1111
1112   Double_t dDCA=myV0->DcaV0Daughters();
1113   if(dDCA >0.5) return 0;
1114
1115   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1116   if(dCTP < -1) return 0;
1117
1118 //   AliESDtrack ieT( iT );
1119 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1120 //   if(TMath::Abs(dD0P) < 0]) return 0;
1121
1122 //   AliESDtrack jeT( jT );
1123 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1124 //   if(TMath::Abs(dD0M) < 0) return 0;
1125
1126 //   Double_t dD0D0=dD0P*dD0M;
1127 //   if(dD0D0>0) return 0;
1128
1129 //   Double_t dETA=myV0->Eta();
1130 //   if(dETA <-0.8) return 0;
1131 //   if(dETA >0.8) return 0;
1132
1133 //   Double_t dQT=myV0->PtArmV0();
1134 //   if(specie==0) if(dQT<???) return 0;
1135
1136   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1137   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1138
1139   if(specie==1 && dALPHA<0) return 2; // antilambda
1140   return 1; // K0s or lambda
1141 }
1142 //-------------------------------------------------------------------------
1143 Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1144   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1145
1146   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1147     AliAODTrack* track = fOutputAOD->GetTrack(i);
1148     if(track == trk) return i;
1149   }
1150
1151   printf("Daughter for %p not found\n",trk);
1152   return -1;
1153 }
1154 //-------------------------------------------------------------------------
1155 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1156   if(!fIsMC) return 1; // used only on MC
1157
1158   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1159     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1160
1161     if(!(channel%20)) return 0; // 5% additional loss in MC
1162   }
1163
1164   return 1;
1165 }