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[AliPID::kSPECIESC];
 716             KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 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     nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
 743
 744     if(tofMatch1){
 745       nSigmaComb = nSigmaTOF;
 746       nSigmaCombRef = nSigmaTOFRef;
 747     }
 748
 749     if(nSigmaComb < 2) nSigmaComb = 2;
 750     else if(nSigmaComb < 3) nSigmaComb = 3;
 751     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 752     else nSigmaComb = 6;
 753
 754     if(nSigmaCombRef < 2) nSigmaCombRef = 2;
 755     else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
 756     else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
 757     else nSigmaCombRef = 6;
 758
 759     Int_t iks=-1;
 760     for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
 761       if(i == fIpiP[k]) iks = k;
 762     }
 763
 764     if(fPtKp > 4.299) fPtKp = 4.299;
 765
 766     if(iks > -1 && fIpiN[iks] > -1){
 767       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 768       Int_t j = fIpiN[iks];
 769       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 770
 771       if (!KnTrack){
 772         continue;
 773       }
 774
 775       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 776
 777       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 778       fPidKn=0;
 779
 780       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 781       Double_t oldpN[10];
 782       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 783
 784       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 785       nSigmaTPC2Ref = PIDResponse->NumberOfSigmasTPC(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 786
 787       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 788
 789       nSigmaComb2=5;
 790       nSigmaTOF2=5;
 791       nSigmaComb2Ref=5;
 792       nSigmaTOF2Ref=5;
 793
 794       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 795       /*
 796       if(mcArray){
 797         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 798         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 799         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 800       }
 801       */
 802
 803       fPidKn = Int_t(probN[2]*100);
 804
 805       if(tofMatch2){
 806         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 807           // remove this tof hit
 808           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 809           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 810           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 811           fPidKn = Int_t(probN[4]*100);
 812           tofMatch2=0;
 813         }
 814         else{
 815           if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 816
 817           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 818           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 819           nSigmaTOF2Ref = PIDResponse->NumberOfSigmasTOF(KnTrack,(AliPID::EParticleType) fSpeciesRef);
 820           nSigmaTOF2Ref = TMath::Abs(nSigmaTOF2Ref);
 821
 822           if(fIsMC){
 823             Float_t mismAdd = addMismatchForMC;
 824             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 825
 826             if(gRandom->Rndm() < mismAdd){
 827               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 828               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 829               channel = channel % 8736;
 830               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 831
 832               // generate random time
 833               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 834               Double_t times[AliPID::kSPECIESC];
 835               KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 836               nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[2], AliPID::kPion);
 837               nSigmaTOF2Ref = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[fSpeciesRef], fSpeciesRef);
 838             }
 839           }
 840
 841           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 842
 843           if(nSigmaTOF2 < 2) fPidKn += 256;
 844           else if(nSigmaTOF2 < 3) fPidKn += 512;
 845         }
 846       }
 847
 848       px = KpTrack->Px() + KnTrack->Px();
 849       py = KpTrack->Py() + KnTrack->Py();
 850       pz = KpTrack->Pz() + KnTrack->Pz();
 851       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
 852       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
 853
 854       ksV.SetPxPyPzE(px,py,pz,E);
 855       fMassV0 = fMassKs[iks];
 856
 857       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 858
 859
 860       fPtKs = ksV.Pt();
 861       fEtaKs = ksV.Eta();
 862       fPhiKs = ksV.Phi();
 863
 864       if(tofMatch2){
 865         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 866         nSigmaComb2Ref = TMath::Sqrt(nSigmaTOF2Ref*nSigmaTOF2Ref+ nSigmaTPC2Ref*nSigmaTPC2Ref);
 867         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 868           fCombsignal->Fill(nSigmaComb2);
 869         }
 870       } else {
 871         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 872         nSigmaComb2Ref = TMath::Abs(nSigmaTPC2Ref);
 873       }
 874
 875       // use sigmaTOF instead of sigmaComb
 876       if(tofMatch2){
 877         nSigmaComb2 = nSigmaTOF2;
 878         nSigmaComb2Ref = nSigmaTOF2Ref;
 879       }
 880
 881       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 882       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 883       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 884       else nSigmaComb2 = 6;
 885       if(nSigmaComb2Ref < 2) nSigmaComb2Ref = 2;
 886       else if(nSigmaComb2Ref < 3) nSigmaComb2Ref = 3;
 887       else if(nSigmaComb2Ref < 5) nSigmaComb2Ref = 4.99;
 888       else nSigmaComb2Ref = 6;
 889
 890       Bool_t isTrue = kFALSE;
 891
 892       if(mcArray){
 893         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 894         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 895
 896         if(labelP > -1 && labelN > -1){
 897           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 898           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 899
 900           Int_t mP = partP->GetMother();
 901           Int_t mN = partN->GetMother();
 902           if(mP == mN && mP > -1){
 903             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 904             Int_t pdgM = partM->GetPdgCode();
 905             if(pdgM == 310) isTrue = kTRUE;
 906           }
 907         }
 908
 909       }
 910
 911       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 912       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 913
 914       if(gRandom->Rndm() < 0.5){
 915         deltaphi1 += TMath::Pi();
 916         deltaphi2 += TMath::Pi();
 917       }
 918
 919       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 920       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 921       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 922       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 923
 924       if(fPtKn > 4.299) fPtKn = 4.299;
 925
 926       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)};
 927       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)};
 928
 929       Int_t ipt = 0;
 930       while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
 931         ipt++;
 932       }
 933       ipt--;
 934       if(ipt < 0) ipt = 0; // just to be sure
 935
 936       if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 937         if(fSpeciesRef != 2){
 938           xTOfill[4] = probP[fSpeciesRef];
 939           xTOfill2[5] = probN[fSpeciesRef];
 940
 941           xTOfill[9] = nSigmaCombRef;
 942           xTOfill2[10] = nSigmaComb2Ref;
 943
 944         }
 945
 946         fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 947         xTOfill[1] = KnTrack->Eta();
 948         fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
 949
 950         if(fPIDuserCut){
 951           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)};
 952           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)};
 953
 954           if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled for positive
 955             xUser[3] = 1;
 956           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 957             xUser[3] = 2;
 958           } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 959             xUser[3] = 3;
 960           }
 961           if(fPIDuserCut->IsSelected(KnTrack,AliPID::kPion)){ // to be filled for negative
 962             xUser2[3] = 1;
 963           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kKaon)){
 964             xUser2[3] = 2;
 965           } else if(fPIDuserCut->IsSelected(KnTrack,AliPID::kProton)){
 966             xUser2[3] = 3;
 967           }
 968           fContUser->Fill(0,fMassV0,fCentrality,xUser);
 969           fContUser2->Fill(0,fMassV0,fCentrality,xUser2);
 970         }
 971
 972       }
 973
 974
 975     }
 976   } // end analysi K0s
 977 }
 978
 979 //_____________________________________________________________________________
 980 Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
 981 {
 982
 983   Float_t zvtx = -999;
 984
 985   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 986   if (!vtxAOD)
 987     return zvtx;
 988   if(vtxAOD->GetNContributors()>0)
 989     zvtx = vtxAOD->GetZ();
 990
 991   return zvtx;
 992 }
 993 //_____________________________________________________________________________
 994 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
 995 {
 996   // Terminate loop
 997   Printf("Terminate()");
 998 }
 999 //=======================================================================
1000 void AliAnalysisTaskK0sBayes::SelectK0s(){
1001   fNK0s=0;
1002   fNpiPos=0;
1003   fNpiNeg=0;
1004
1005   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1006   AliAODv0 *myV0;
1007   Double_t dMASS=0.0;
1008   for (Int_t i=0; i!=nV0s; ++i) {
1009     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1010     if(!myV0) continue;
1011     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
1012     Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
1013     if(pass) {
1014       dMASS = myV0->MassK0Short();
1015       Float_t massLambda = myV0->MassLambda();
1016       Float_t massAntiLambda = myV0->MassAntiLambda();
1017
1018       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){
1019         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1020         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1021         if(iT->Charge()<0){
1022           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1023           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1024         }
1025         fPhiK0s[fNK0s] = myV0->Phi();
1026         fPtK0s[fNK0s] = myV0->Pt();
1027         fIpiP[fNK0s] = FindDaugheterIndex(iT);
1028         fIpiN[fNK0s] = FindDaugheterIndex(jT);
1029         fMassKs[fNK0s] = dMASS;
1030         if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
1031           fNK0s++;
1032       }
1033     }
1034   }
1035
1036   /* My V0 code
1037   // fill pion stacks
1038   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
1039   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1040     AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
1041
1042     if (!aodTrack){
1043       continue;
1044     }
1045
1046     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
1047
1048     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1049       continue;
1050     }
1051
1052     Double_t b[2] = {-99., -99.};
1053     Double_t bCov[3] = {-99., -99., -99.};
1054     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
1055       continue;
1056
1057     if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
1058
1059
1060     Int_t charge = aodTrack->Charge();
1061     if(charge > 0){
1062       fIPiPos[fNpiPos] = iT;
1063       fNpiPos++;
1064     }
1065     else{
1066       fIPiNeg[fNpiNeg] = iT;
1067       fNpiNeg++;
1068     }
1069   }
1070
1071   for(Int_t i=0;i < fNpiPos;i++){
1072     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
1073     AliESDtrack pipE(pip);
1074
1075     for(Int_t j=0;j < fNpiNeg;j++){
1076       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
1077       AliESDtrack pinE(pin);
1078
1079       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
1080
1081       Double_t pPos[3];
1082       Double_t pNeg[3];
1083       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
1084       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
1085
1086       Float_t length = (xp+xn)*0.5;
1087
1088       Float_t pxs = pPos[0] + pNeg[0];
1089       Float_t pys = pPos[1] + pNeg[1];
1090       Float_t pzs = pPos[2] + pNeg[2];
1091       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);
1092
1093       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
1094       Float_t phi = TMath::ATan2(pys,pxs);
1095       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
1096
1097       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
1098
1099       if(mindist < 0.4&& length > 0.7 && length < 25){
1100         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);
1101         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);
1102
1103         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
1104         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
1105
1106         if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
1107           fPhiK0s[fNK0s] = phi;
1108           fPtK0s[fNK0s] = pt;
1109           fIpiP[fNK0s] =fIPiPos[i] ;
1110           fIpiN[fNK0s] = fIPiNeg[j];
1111           fMassKs[fNK0s] = mass;
1112           fNK0s++;
1113         }
1114       }
1115     }
1116   }
1117   */
1118 }
1119
1120 //=======================================================================
1121 Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1122 {
1123   if (myV0->GetOnFlyStatus() ) return 0;
1124
1125   //the following is needed in order to evualuate track-quality
1126   AliAODTrack *iT, *jT;
1127   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1128   Double_t pos[3],cov[6];
1129   vV0s->GetXYZ(pos);
1130   vV0s->GetCovarianceMatrix(cov);
1131   const AliESDVertex vESD(pos,cov,100.,100);
1132
1133   // TESTING CHARGE
1134   int iPos, iNeg;
1135   iT=(AliAODTrack*) myV0->GetDaughter(0);
1136   if(iT->Charge()>0) {
1137     iPos = 0; iNeg = 1;
1138   } else {
1139     iPos = 1; iNeg = 0;
1140   }
1141   // END OF TEST
1142
1143   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1144
1145   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1146
1147   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1148   if(!trkFlag) return 0;
1149   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1150   if(!trkFlag2) return 0;
1151
1152   Double_t pvertex[3];
1153   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1154   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1155   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1156
1157   Double_t dDL=myV0->DecayLengthV0( pvertex );
1158   if(dDL  < 0.5 || dDL > 25) return 0;
1159
1160   Double_t dDCA=myV0->DcaV0Daughters();
1161   if(dDCA >0.5) return 0;
1162
1163   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1164   if(dCTP < -1) return 0;
1165
1166 //   AliESDtrack ieT( iT );
1167 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1168 //   if(TMath::Abs(dD0P) < 0]) return 0;
1169
1170 //   AliESDtrack jeT( jT );
1171 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1172 //   if(TMath::Abs(dD0M) < 0) return 0;
1173
1174 //   Double_t dD0D0=dD0P*dD0M;
1175 //   if(dD0D0>0) return 0;
1176
1177 //   Double_t dETA=myV0->Eta();
1178 //   if(dETA <-0.8) return 0;
1179 //   if(dETA >0.8) return 0;
1180
1181 //   Double_t dQT=myV0->PtArmV0();
1182 //   if(specie==0) if(dQT<???) return 0;
1183
1184   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1185   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1186
1187   if(specie==1 && dALPHA<0) return 2; // antilambda
1188   return 1; // K0s or lambda
1189 }
1190 //-------------------------------------------------------------------------
1191 Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1192   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1193
1194   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1195     AliAODTrack* track = fOutputAOD->GetTrack(i);
1196     if(track == trk) return i;
1197   }
1198
1199   printf("Daughter for %p not found\n",trk);
1200   return -1;
1201 }
1202 //-------------------------------------------------------------------------
1203 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1204   if(!fIsMC) return 1; // used only on MC
1205
1206   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1207     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1208
1209     if(!(channel%20)) return 0; // 5% additional loss in MC
1210   }
1211
1212   return 1;
1213 }