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[AliPID::kSPECIESC];
 731             KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 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     nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
 758
 759     if(tofMatch1){
 760       nSigmaComb = nSigmaTOF;
 761       nSigmaCombRef = nSigmaTOFRef;
 762     }
 763
 764     if(nSigmaComb < 2) nSigmaComb = 2;
 765     else if(nSigmaComb < 3) nSigmaComb = 3;
 766     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 767     else nSigmaComb = 6;
 768
 769     if(nSigmaCombRef < 2) nSigmaCombRef = 2;
 770     else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
 771     else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
 772     else nSigmaCombRef = 6;
 773
 774     Int_t iks=-1;
 775     for(Int_t k=0;k < fNLambda;k++){ // find the Lambda which contains the positive track
 776       if(i == fIpP[k]) iks = k;
 777     }
 778
 779     if(fPtKp > 4.299) fPtKp = 4.299;
 780
 781     if(iks > -1 && fIpN[iks] > -1){
 782       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 783       Int_t j = fIpN[iks];
 784       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 785
 786       if (!KnTrack){
 787         continue;
 788       }
 789
 790       nSigmaComb2 = 5;
 791       nSigmaTOF2 = 5;
 792
 793       if(!(KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 794
 795       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 796       fPidKn=0;
 797
 798       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 799       Double_t oldpN[10];
 800       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 801
 802       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 803
 804       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 805
 806       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 807       /*
 808       if(mcArray){
 809         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 810         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 811         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 812       }
 813       */
 814
 815       fPidKn = Int_t(probN[2]*100);
 816
 817       if(tofMatch2){
 818         if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 819
 820         nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 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+01;
 834             Double_t times[AliPID::kSPECIESC];
 835             KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 836             nSigmaTOF2 = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[4], AliPID::kProton);
 837           }
 838         }
 839
 840
 841
 842         nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 843
 844         if(nSigmaTOF2 < 2) fPidKn += 256;
 845         else if(nSigmaTOF2 < 3) fPidKn += 512;
 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() + 0.938e-01*0.938e-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 = fMassLambda[iks];
 856
 857       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 858
 859
 860       fPtLambdaC = ksV.Pt();
 861       fEtaLambda = ksV.Eta();
 862       fPhiLambdaC = ksV.Phi();
 863
 864       if(tofMatch2){
 865         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 866       } else {
 867         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 868       }
 869
 870       // use sigmaTOF instead of sigmaComb
 871       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 872
 873       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 874       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 875       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 876       else nSigmaComb2 = 6;
 877
 878       Bool_t isTrue = kFALSE;
 879
 880       if(mcArray){
 881         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 882         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 883
 884         if(labelP > -1 && labelN > -1){
 885           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 886           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 887
 888           Int_t mP = partP->GetMother();
 889           Int_t mN = partN->GetMother();
 890           if(mP == mN && mP > -1){
 891             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 892             Int_t pdgM = partM->GetPdgCode();
 893             if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
 894           }
 895         }
 896
 897       }
 898
 899       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 900       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 901
 902       if(gRandom->Rndm() < 0.5){
 903         deltaphi1 += TMath::Pi();
 904         deltaphi2 += TMath::Pi();
 905       }
 906
 907       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 908       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 909       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 910       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 911
 912       if(fPtKn > 4.299) fPtKn = 4.299;
 913
 914       Float_t xTOfill[] = {static_cast<Float_t>(fPtLambdaC),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn,(fPidKp%128)*0.01),static_cast<Float_t>((fPidKn%128)*0.01),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)};
 915       Float_t xTOfill2[] = {static_cast<Float_t>(fPtLambdaC),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKn),static_cast<Float_t>(fPtKp),static_cast<Float_t>((fPidKn%128)*0.01),static_cast<Float_t>((fPidKp%128)*0.01),static_cast<Float_t>(tofMatch2),static_cast<Float_t>(tofMatch1),static_cast<Float_t>(isTrue),static_cast<Float_t>(nSigmaComb2),static_cast<Float_t>(nSigmaComb),static_cast<Float_t>(deltaphi2),static_cast<Float_t>(deltaphi1),static_cast<Float_t>(fPsi)};
 916
 917
 918       Int_t ipt = 0;
 919       while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
 920         ipt++;
 921       }
 922       ipt--;
 923       if(ipt < 0) ipt = 0; // just to be sure
 924
 925       if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 926         if(fSpeciesRef != 4){
 927           xTOfill[4] = probP[fSpeciesRef];
 928           xTOfill2[5] = probP[fSpeciesRef];
 929
 930           xTOfill[9] = nSigmaCombRef;
 931           xTOfill2[10] = nSigmaCombRef;
 932
 933         }
 934
 935         if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 936         else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
 937
 938
 939
 940         if(fPIDuserCut){
 941           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)};
 942
 943          if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled
 944            xUser[3] = 1;
 945          } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 946            xUser[3] = 2;
 947          } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 948            xUser[3] = 3;
 949          }
 950
 951          if(isLambda) fContUser->Fill(0,fMassV0,fCentrality,xUser);
 952          else fContUser2->Fill(0,fMassV0,fCentrality,xUser);
 953         }
 954       }
 955
 956
 957     }
 958   } // end analysi Lambda
 959 }
 960
 961 //_____________________________________________________________________________
 962 Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
 963 {
 964
 965   Float_t zvtx = -999;
 966
 967   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 968   if (!vtxAOD)
 969     return zvtx;
 970   if(vtxAOD->GetNContributors()>0)
 971     zvtx = vtxAOD->GetZ();
 972
 973   return zvtx;
 974 }
 975 //_____________________________________________________________________________
 976 void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
 977 {
 978   // Terminate loop
 979   Printf("Terminate()");
 980 }
 981 //=======================================================================
 982 void AliAnalysisTaskLambdaBayes::SelectLambda(){
 983   fNLambda=0;
 984   fNpPos=0;
 985   fNpNeg=0;
 986
 987   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
 988   AliAODv0 *myV0;
 989   Double_t dMASS=0.0;
 990   for (Int_t i=0; i!=nV0s; ++i) {
 991     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
 992     if(!myV0) continue;
 993     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
 994     Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
 995     if(pass) {
 996       if(pass==1) dMASS = myV0->MassLambda();
 997       else dMASS = myV0->MassAntiLambda();
 998
 999       Float_t massKs = myV0->MassK0Short();
1000
1001       if(TMath::Abs(dMASS-1.115)/0.005 < 8 && TMath::Abs(massKs - 0.497)/0.005 > 8){
1002         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1003         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1004         if(iT->Charge()<0){
1005           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1006           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1007         }
1008         fPhiLambda[fNLambda] = myV0->Phi();
1009         fPtLambda[fNLambda] = myV0->Pt();
1010         fIpP[fNLambda] = FindDaugheterIndex(iT);
1011         fIpN[fNLambda] = FindDaugheterIndex(jT);
1012
1013         if(pass==2){
1014           Int_t itemp = fIpP[fNLambda];
1015           fIpP[fNLambda] = fIpN[fNLambda];
1016           fIpN[fNLambda] = itemp;
1017         }
1018
1019         fMassLambda[fNLambda] = dMASS;
1020         if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
1021           fNLambda++;
1022         }
1023
1024       }
1025     }
1026   }
1027 }
1028
1029 //=======================================================================
1030 Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1031 {
1032   if (myV0->GetOnFlyStatus() ) return 0;
1033
1034   //the following is needed in order to evualuate track-quality
1035   AliAODTrack *iT, *jT;
1036   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1037   Double_t pos[3],cov[6];
1038   vV0s->GetXYZ(pos);
1039   vV0s->GetCovarianceMatrix(cov);
1040   const AliESDVertex vESD(pos,cov,100.,100);
1041
1042   // TESTING CHARGE
1043   int iPos, iNeg;
1044   iT=(AliAODTrack*) myV0->GetDaughter(0);
1045   if(iT->Charge()>0) {
1046     iPos = 0; iNeg = 1;
1047   } else {
1048     iPos = 1; iNeg = 0;
1049   }
1050   // END OF TEST
1051
1052   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1053
1054   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1055
1056   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1057   if(!trkFlag) return 0;
1058   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1059   if(!trkFlag2) return 0;
1060
1061   Double_t pvertex[3];
1062   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1063   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1064   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1065
1066   Double_t dDL=myV0->DecayLengthV0( pvertex );
1067   if(dDL  < 0.5 || dDL > 25) return 0;
1068
1069   Double_t dDCA=myV0->DcaV0Daughters();
1070   if(dDCA >0.5) return 0;
1071
1072   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1073   if(dCTP < -1) return 0;
1074
1075 //   AliESDtrack ieT( iT );
1076 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1077 //   if(TMath::Abs(dD0P) < 0]) return 0;
1078
1079 //   AliESDtrack jeT( jT );
1080 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1081 //   if(TMath::Abs(dD0M) < 0) return 0;
1082
1083 //   Double_t dD0D0=dD0P*dD0M;
1084 //   if(dD0D0>0) return 0;
1085
1086 //   Double_t dETA=myV0->Eta();
1087 //   if(dETA <-0.8) return 0;
1088 //   if(dETA >0.8) return 0;
1089
1090 //   Double_t dQT=myV0->PtArmV0();
1091 //   if(specie==0) if(dQT<???) return 0;
1092
1093   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1094   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1095
1096   if(specie==1 && dALPHA<0) return 2; // antilambda
1097   return 1; // Ks or lambda
1098 }
1099 //-------------------------------------------------------------------------
1100 Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
1101   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1102
1103   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1104     AliAODTrack* track = fOutputAOD->GetTrack(i);
1105     if(track == trk) return i;
1106   }
1107
1108   printf("Daughter for %p not found\n",trk);
1109   return -1;
1110 }
1111 //-------------------------------------------------------------------------
1112 Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
1113   if(!fIsMC) return 1; // used only on MC
1114
1115   if( fTypeCol==2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1116     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1117
1118     if(!(channel%20)) return 0; // 5% additional loss in MC
1119   }
1120
1121   return 1;
1122 }