PWGPP/pid/AliAnalysisTaskLambdaBayes.cxx

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