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 = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(iT));
 526     if(!aodTrack) AliFatal("Not a standard AOD");
 527
 528     if (!aodTrack){
 529       continue;
 530     }
 531
 532     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 533
 534     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 535       continue;
 536
 537     Double_t b[2] = {-99., -99.};
 538     Double_t bCov[3] = {-99., -99., -99.};
 539     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 540       continue;
 541
 542     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 543       continue;
 544
 545     Qx2 += TMath::Cos(2*aodTrack->Phi());
 546     Qy2 += TMath::Sin(2*aodTrack->Phi());
 547     Qx3 += TMath::Cos(3*aodTrack->Phi());
 548     Qy3 += TMath::Sin(3*aodTrack->Phi());
 549
 550   }
 551
 552   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 553
 554   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 555   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 556   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 557
 558   PIDResponse->SetTOFResponse(aodEvent,AliPIDResponse::kTOF_T0);
 559
 560   PIDResponse->GetTOFResponse().SetTOFtailAllPara(-23,1.1);
 561
 562   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 563
 564   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 565   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 566   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 567   Double_t nSigmaTPCRef,nSigmaTOFRef=6,nSigmaTPC2Ref,nSigmaTOF2Ref=6,nSigmaCombRef,nSigmaComb2Ref;
 568
 569
 570   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 571   TClonesArray *mcArray = NULL;
 572   if (mcHeader)
 573     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 574
 575   Int_t nmc = 0;
 576   if(mcArray)
 577     nmc = mcArray->GetEntries();
 578
 579   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 580     AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
 581     if(!track) AliFatal("Not a standard AOD");
 582
 583     AliAODMCParticle *mcp = NULL;
 584     Int_t pdg = 0;
 585
 586     if (!track){
 587       continue;
 588     }
 589
 590     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 591
 592     Int_t label = -1;
 593     if(mcArray){
 594       label = track->GetLabel();
 595       if(label != -1 && label < nmc){
 596         label = TMath::Abs(label);
 597         mcp = (AliAODMCParticle*)mcArray->At(label);
 598         pdg = TMath::Abs(mcp->GetPdgCode());
 599       }
 600       else
 601         label = -1;
 602     }
 603     else{
 604       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 605     }
 606
 607     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 608       hTracking[0]->Fill(track->Pt(),fCentrality);
 609       if(pdg == 211)
 610         hTracking[1]->Fill(track->Pt(),fCentrality);
 611       else if(pdg == 321)
 612         hTracking[2]->Fill(track->Pt(),fCentrality);
 613       else if(pdg == 2212)
 614         hTracking[3]->Fill(track->Pt(),fCentrality);
 615       else if(! mcp){ // Fill matching histo with the prob
 616         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 617         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 618         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 619       }
 620     }
 621
 622     if(!tofMatch) continue;
 623
 624     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 625       hMatching[0]->Fill(track->Pt(),fCentrality);
 626       if(pdg == 211)
 627         hMatching[1]->Fill(track->Pt(),fCentrality);
 628       else if(pdg == 321)
 629         hMatching[2]->Fill(track->Pt(),fCentrality);
 630       else if(pdg == 2212)
 631         hMatching[3]->Fill(track->Pt(),fCentrality);
 632       else if(! mcp){ // Fill matching histo with the prob
 633         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 634         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 635         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 636       }
 637     }
 638   }
 639
 640 //   Int_t pdg1 = -1;
 641 //   Int_t pdg2 = -1;
 642
 643
 644   // start analysis K0s
 645   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 646     AliAODTrack* KpTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
 647     if(!KpTrack) AliFatal("Not a standard AOD");
 648
 649     if (!KpTrack){
 650       continue;
 651     }
 652
 653     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 654
 655     nSigmaComb=5;
 656     nSigmaCombRef=5;
 657     nSigmaTOF = 5;
 658     nSigmaTOFRef = 5;
 659     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 660     fPidKp=0;
 661
 662     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 663
 664     Double_t oldpP[10];
 665     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 666
 667     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
 668     nSigmaTPCRef = PIDResponse->NumberOfSigmasTPC(KpTrack,(AliPID::EParticleType) fSpeciesRef);
 669
 670     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 671
 672     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 673
 674     /*
 675     if(mcArray){
 676       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 677       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 678       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 679     }
 680     */
 681
 682     fPidKp = Int_t(probP[2]*100);
 683
 684     if(tofMatch1){
 685       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 686         // remove this tof hit
 687         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 688         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 689         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 690         fPidKp = Int_t(probP[4]*100);
 691         tofMatch1=0;
 692       }
 693       else{
 694         if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
 695
 696         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 697         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 698         if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
 699         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 700         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 701         if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
 702         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 703         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 704         if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
 705         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 706         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 707         if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
 708         nSigmaTOFRef = PIDResponse->NumberOfSigmasTOF(KpTrack,(AliPID::EParticleType) fSpeciesRef);
 709
 710         if(fIsMC){
 711           Float_t mismAdd = addMismatchForMC;
 712           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 713
 714           if(gRandom->Rndm() < mismAdd){
 715             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 716             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 717             channel = channel % 8736;
 718             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 719
 720             // generate random time
 721             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 722             Double_t times[AliPID::kSPECIESC];
 723             KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 724             nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[2], AliPID::kPion);
 725             nSigmaTOFRef = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[fSpeciesRef], fSpeciesRef);
 726           }
 727         }
 728
 729         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 730         nSigmaTOF = TMath::Abs(nSigmaTOF);
 731
 732         if(nSigmaTOF < 2) fPidKp += 256;
 733         else if(nSigmaTOF < 3) fPidKp += 512;
 734       }
 735     }
 736
 737     if(tofMatch1){
 738       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 739       nSigmaCombRef = TMath::Sqrt(nSigmaTOFRef*nSigmaTOFRef + nSigmaTPCRef*nSigmaTPCRef);
 740       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 741         fCombsignal->Fill(nSigmaComb);
 742       }
 743     } else {
 744       nSigmaComb = TMath::Abs(nSigmaTPC);
 745       nSigmaCombRef = TMath::Abs(nSigmaTPCRef);
 746     }
 747
 748     // use sigmaTOF instead of sigmaComb
 749     nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
 750
 751     if(tofMatch1){
 752       nSigmaComb = nSigmaTOF;
 753       nSigmaCombRef = nSigmaTOFRef;
 754     }
 755
 756     if(nSigmaComb < 2) nSigmaComb = 2;
 757     else if(nSigmaComb < 3) nSigmaComb = 3;
 758     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 759     else nSigmaComb = 6;
 760
 761     if(nSigmaCombRef < 2) nSigmaCombRef = 2;
 762     else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
 763     else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
 764     else nSigmaCombRef = 6;
 765
 766     Int_t iks=-1;
 767     for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
 768       if(i == fIpiP[k]) iks = k;
 769     }
 770
 771     if(fPtKp > 4.299) fPtKp = 4.299;
 772
 773     if(iks > -1 && fIpiN[iks] > -1){
 774       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 775       Int_t j = fIpiN[iks];
 776       AliAODTrack* KnTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(j));
 777       if(!KnTrack) AliFatal("Not a standard AOD");
 778
 779       if (!KnTrack){
 780         continue;
 781       }
 782
 783       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 784
 785       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 786       fPidKn=0;
 787
 788       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 789       Double_t oldpN[10];
 790       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 791
 792       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 793       nSigmaTPC2Ref = PIDResponse->NumberOfSigmasTPC(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 794
 795       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 796
 797       nSigmaComb2=5;
 798       nSigmaTOF2=5;
 799       nSigmaComb2Ref=5;
 800       nSigmaTOF2Ref=5;
 801
 802       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 803       /*
 804       if(mcArray){
 805         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 806         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 807         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 808       }
 809       */
 810
 811       fPidKn = Int_t(probN[2]*100);
 812
 813       if(tofMatch2){
 814         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 815           // remove this tof hit
 816           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 817           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 818           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 819           fPidKn = Int_t(probN[4]*100);
 820           tofMatch2=0;
 821         }
 822         else{
 823           if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 824
 825           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 826           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 827           nSigmaTOF2Ref = PIDResponse->NumberOfSigmasTOF(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 828           nSigmaTOF2Ref = TMath::Abs(nSigmaTOF2Ref);
 829
 830           if(fIsMC){
 831             Float_t mismAdd = addMismatchForMC;
 832             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 833
 834             if(gRandom->Rndm() < mismAdd){
 835               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 836               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 837               channel = channel % 8736;
 838               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 839
 840               // generate random time
 841               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 842               Double_t times[AliPID::kSPECIESC];
 843               KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 844               nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[2], AliPID::kPion);
 845               nSigmaTOF2Ref = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[fSpeciesRef], fSpeciesRef);
 846             }
 847           }
 848
 849           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 850
 851           if(nSigmaTOF2 < 2) fPidKn += 256;
 852           else if(nSigmaTOF2 < 3) fPidKn += 512;
 853         }
 854       }
 855
 856       px = KpTrack->Px() + KnTrack->Px();
 857       py = KpTrack->Py() + KnTrack->Py();
 858       pz = KpTrack->Pz() + KnTrack->Pz();
 859       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
 860       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
 861
 862       ksV.SetPxPyPzE(px,py,pz,E);
 863       fMassV0 = fMassKs[iks];
 864
 865       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 866
 867
 868       fPtKs = ksV.Pt();
 869       fEtaKs = ksV.Eta();
 870       fPhiKs = ksV.Phi();
 871
 872       if(tofMatch2){
 873         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 874         nSigmaComb2Ref = TMath::Sqrt(nSigmaTOF2Ref*nSigmaTOF2Ref+ nSigmaTPC2Ref*nSigmaTPC2Ref);
 875         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 876           fCombsignal->Fill(nSigmaComb2);
 877         }
 878       } else {
 879         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 880         nSigmaComb2Ref = TMath::Abs(nSigmaTPC2Ref);
 881       }
 882
 883       // use sigmaTOF instead of sigmaComb
 884       if(tofMatch2){
 885         nSigmaComb2 = nSigmaTOF2;
 886         nSigmaComb2Ref = nSigmaTOF2Ref;
 887       }
 888
 889       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 890       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 891       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 892       else nSigmaComb2 = 6;
 893       if(nSigmaComb2Ref < 2) nSigmaComb2Ref = 2;
 894       else if(nSigmaComb2Ref < 3) nSigmaComb2Ref = 3;
 895       else if(nSigmaComb2Ref < 5) nSigmaComb2Ref = 4.99;
 896       else nSigmaComb2Ref = 6;
 897
 898       Bool_t isTrue = kFALSE;
 899
 900       if(mcArray){
 901         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 902         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 903
 904         if(labelP > -1 && labelN > -1){
 905           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 906           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 907
 908           Int_t mP = partP->GetMother();
 909           Int_t mN = partN->GetMother();
 910           if(mP == mN && mP > -1){
 911             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 912             Int_t pdgM = partM->GetPdgCode();
 913             if(pdgM == 310) isTrue = kTRUE;
 914           }
 915         }
 916
 917       }
 918
 919       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 920       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 921
 922       if(gRandom->Rndm() < 0.5){
 923         deltaphi1 += TMath::Pi();
 924         deltaphi2 += TMath::Pi();
 925       }
 926
 927       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 928       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 929       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 930       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 931
 932       if(fPtKn > 4.299) fPtKn = 4.299;
 933
 934       Float_t xTOfill[] = {static_cast<Float_t>(fPtKs),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn),static_cast<Float_t>(probP[2]),static_cast<Float_t>(probN[2]),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
 935       Float_t xTOfill2[] = {static_cast<Float_t>(fPtKs),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn),static_cast<Float_t>(probP[2]),static_cast<Float_t>(probN[2]),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
 936
 937       Int_t ipt = 0;
 938       while(ipt < nPtBin && fPtKsMin[ipt] < fPtKs){
 939         ipt++;
 940       }
 941       ipt--;
 942       if(ipt < 0) ipt = 0; // just to be sure
 943
 944       if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 945         if(fSpeciesRef != 2){
 946           xTOfill[4] = probP[fSpeciesRef];
 947           xTOfill2[5] = probN[fSpeciesRef];
 948
 949           xTOfill[9] = nSigmaCombRef;
 950           xTOfill2[10] = nSigmaComb2Ref;
 951
 952         }
 953
 954         fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 955         xTOfill[1] = KnTrack->Eta();
 956         fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
 957
 958         if(fPIDuserCut){
 959           Float_t xUser[] = {static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(isTrue),0,static_cast<Float_t>(deltaphi1),static_cast<Float_t>(fPsi)};
 960           Float_t xUser2[] = {static_cast<Float_t>(KnTrack->Eta()),static_cast<Float_t>(fPtKn),static_cast<Float_t>(isTrue),0,static_cast<Float_t>(deltaphi2),static_cast<Float_t>(fPsi)};
 961
 962           if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
 963             xUser[3] = 1;
 964           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 965             xUser[3] = 2;
 966           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 967             xUser[3] = 3;
 968           }
 969           if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
 970             xUser2[3] = 1;
 971           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
 972             xUser2[3] = 2;
 973           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
 974             xUser2[3] = 3;
 975           }
 976           fContUser->Fill(0,fMassV0,fCentrality,xUser);
 977           fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
 978         }
 979
 980       }
 981
 982
 983     }
 984   } // end analysi K0s
 985 }
 986
 987 //_____________________________________________________________________________
 988 Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
 989 {
 990
 991   Float_t zvtx = -999;
 992
 993   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 994   if (!vtxAOD)
 995     return zvtx;
 996   if(vtxAOD->GetNContributors()>0)
 997     zvtx = vtxAOD->GetZ();
 998
 999   return zvtx;
1000 }
1001 //_____________________________________________________________________________
1002 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
1003 {
1004   // Terminate loop
1005   Printf("Terminate()");
1006 }
1007 //=======================================================================
1008 void AliAnalysisTaskK0sBayes::SelectK0s(){
1009   fNK0s=0;
1010   fNpiPos=0;
1011   fNpiNeg=0;
1012
1013   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1014   AliAODv0 *myV0;
1015   Double_t dMASS=0.0;
1016   for (Int_t i=0; i!=nV0s; ++i) {
1017     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1018     if(!myV0) continue;
1019     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
1020     Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
1021     if(pass) {
1022       dMASS = myV0->MassK0Short();
1023       Float_t massLambda = myV0->MassLambda();
1024       Float_t massAntiLambda = myV0->MassAntiLambda();
1025
1026       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){
1027         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1028         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1029         if(iT->Charge()<0){
1030           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1031           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1032         }
1033         fPhiK0s[fNK0s] = myV0->Phi();
1034         fPtK0s[fNK0s] = myV0->Pt();
1035         fIpiP[fNK0s] = FindDaugheterIndex(iT);
1036         fIpiN[fNK0s] = FindDaugheterIndex(jT);
1037         fMassKs[fNK0s] = dMASS;
1038         if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
1039           fNK0s++;
1040       }
1041     }
1042   }
1043
1044   /* My V0 code
1045   // fill pion stacks
1046   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
1047   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1048     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fOutputAOD->GetTrack(iT));
1049     if(!aodTrack) AliFatal("Not a standard AOD");
1050
1051     if (!aodTrack){
1052       continue;
1053     }
1054
1055     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
1056
1057     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1058       continue;
1059     }
1060
1061     Double_t b[2] = {-99., -99.};
1062     Double_t bCov[3] = {-99., -99., -99.};
1063     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1064       continue;
1065
1066     if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
1067
1068
1069     Int_t charge = aodTrack->Charge();
1070     if(charge > 0){
1071       fIPiPos[fNpiPos] = iT;
1072       fNpiPos++;
1073     }
1074     else{
1075       fIPiNeg[fNpiNeg] = iT;
1076       fNpiNeg++;
1077     }
1078   }
1079
1080   for(Int_t i=0;i < fNpiPos;i++){
1081     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
1082     AliESDtrack pipE(pip);
1083
1084     for(Int_t j=0;j < fNpiNeg;j++){
1085       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
1086       AliESDtrack pinE(pin);
1087
1088       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
1089
1090       Double_t pPos[3];
1091       Double_t pNeg[3];
1092       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
1093       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
1094
1095       Float_t length = (xp+xn)*0.5;
1096
1097       Float_t pxs = pPos[0] + pNeg[0];
1098       Float_t pys = pPos[1] + pNeg[1];
1099       Float_t pzs = pPos[2] + pNeg[2];
1100       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);
1101
1102       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
1103       Float_t phi = TMath::ATan2(pys,pxs);
1104       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
1105
1106       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
1107
1108       if(mindist < 0.4&& length > 0.7 && length < 25){
1109         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);
1110         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);
1111
1112         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
1113         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
1114
1115         if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
1116           fPhiK0s[fNK0s] = phi;
1117           fPtK0s[fNK0s] = pt;
1118           fIpiP[fNK0s] =fIPiPos[i] ;
1119           fIpiN[fNK0s] = fIPiNeg[j];
1120           fMassKs[fNK0s] = mass;
1121           fNK0s++;
1122         }
1123       }
1124     }
1125   }
1126   */
1127 }
1128
1129 //=======================================================================
1130 Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1131 {
1132   if (myV0->GetOnFlyStatus() ) return 0;
1133
1134   //the following is needed in order to evualuate track-quality
1135   AliAODTrack *iT, *jT;
1136   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1137   Double_t pos[3],cov[6];
1138   vV0s->GetXYZ(pos);
1139   vV0s->GetCovarianceMatrix(cov);
1140   const AliESDVertex vESD(pos,cov,100.,100);
1141
1142   // TESTING CHARGE
1143   int iPos, iNeg;
1144   iT=(AliAODTrack*) myV0->GetDaughter(0);
1145   if(iT->Charge()>0) {
1146     iPos = 0; iNeg = 1;
1147   } else {
1148     iPos = 1; iNeg = 0;
1149   }
1150   // END OF TEST
1151
1152   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1153
1154   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1155
1156   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1157   if(!trkFlag) return 0;
1158   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1159   if(!trkFlag2) return 0;
1160
1161   Double_t pvertex[3];
1162   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1163   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1164   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1165
1166   Double_t dDL=myV0->DecayLengthV0( pvertex );
1167   if(dDL  < 0.5 || dDL > 25) return 0;
1168
1169   Double_t dDCA=myV0->DcaV0Daughters();
1170   if(dDCA >0.5) return 0;
1171
1172   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1173   if(dCTP < -1) return 0;
1174
1175 //   AliESDtrack ieT( iT );
1176 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1177 //   if(TMath::Abs(dD0P) < 0]) return 0;
1178
1179 //   AliESDtrack jeT( jT );
1180 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1181 //   if(TMath::Abs(dD0M) < 0) return 0;
1182
1183 //   Double_t dD0D0=dD0P*dD0M;
1184 //   if(dD0D0>0) return 0;
1185
1186 //   Double_t dETA=myV0->Eta();
1187 //   if(dETA <-0.8) return 0;
1188 //   if(dETA >0.8) return 0;
1189
1190 //   Double_t dQT=myV0->PtArmV0();
1191 //   if(specie==0) if(dQT<???) return 0;
1192
1193   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1194   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1195
1196   if(specie==1 && dALPHA<0) return 2; // antilambda
1197   return 1; // K0s or lambda
1198 }
1199 //-------------------------------------------------------------------------
1200 Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1201   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1202
1203   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1204     AliAODTrack* track = dynamic_cast<AliAODTrack*>(fOutputAOD->GetTrack(i));
1205     if(!track) AliFatal("Not a standard AOD");
1206     if(track == trk) return i;
1207   }
1208
1209   printf("Daughter for %p not found\n",trk);
1210   return -1;
1211 }
1212 //-------------------------------------------------------------------------
1213 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1214   if(!fIsMC) return 1; // used only on MC
1215
1216   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1217     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1218
1219     if(!(channel%20)) return 0; // 5% additional loss in MC
1220   }
1221
1222   return 1;
1223 }