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