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   PIDResponse->SetTOFResponse(aodEvent,AliPIDResponse::kTOF_T0);
 561   PIDResponse->GetTOFResponse().SetTOFtailAllPara(-23,1.1);
 562
 563 //   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
 564 //   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
 565 //   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
 566 //   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
 567
 568   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 569
 570
 571   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 572   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 573   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 574   Double_t nSigmaTPCRef,nSigmaTOFRef=6,nSigmaCombRef;
 575
 576   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 577   TClonesArray *mcArray = NULL;
 578   if (mcHeader)
 579     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 580
 581   Int_t nmc = 0;
 582   if(mcArray)
 583     nmc = mcArray->GetEntries();
 584
 585   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 586     AliAODTrack* track = aodEvent->GetTrack(i);
 587
 588     AliAODMCParticle *mcp = NULL;
 589     Int_t pdg = 0;
 590
 591     if (!track){
 592       continue;
 593     }
 594
 595     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 596
 597     Int_t label = -1;
 598     if(mcArray){
 599       label = track->GetLabel();
 600       if(label != -1 && label < nmc){
 601         label = TMath::Abs(label);
 602         mcp = (AliAODMCParticle*)mcArray->At(label);
 603         pdg = TMath::Abs(mcp->GetPdgCode());
 604       }
 605       else
 606         label = -1;
 607     }
 608     else{
 609       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 610     }
 611
 612     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 613       hTracking[0]->Fill(track->Pt(),fCentrality);
 614       if(pdg == 211)
 615         hTracking[1]->Fill(track->Pt(),fCentrality);
 616       else if(pdg == 321)
 617         hTracking[2]->Fill(track->Pt(),fCentrality);
 618       else if(pdg == 2212)
 619         hTracking[3]->Fill(track->Pt(),fCentrality);
 620       else if(! mcp){ // Fill matching histo with the prob
 621         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 622         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 623         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 624       }
 625     }
 626
 627     if(!tofMatch) continue;
 628
 629     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 630       hMatching[0]->Fill(track->Pt(),fCentrality);
 631       if(pdg == 211)
 632         hMatching[1]->Fill(track->Pt(),fCentrality);
 633       else if(pdg == 321)
 634         hMatching[2]->Fill(track->Pt(),fCentrality);
 635       else if(pdg == 2212)
 636         hMatching[3]->Fill(track->Pt(),fCentrality);
 637       else if(! mcp){ // Fill matching histo with the prob
 638         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 639         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 640         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 641       }
 642     }
 643   }
 644
 645 //   Int_t pdg1 = -1;
 646 //   Int_t pdg2 = -1;
 647
 648
 649   // start analysis Lambda
 650   for(Int_t i=0;i < ntrack;i++){ // loop on proton candidate tracks
 651     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 652
 653     if (!KpTrack){
 654       continue;
 655     }
 656
 657     if(!(KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 658
 659     Bool_t isLambda = (KpTrack->Charge() > 0);
 660
 661     nSigmaComb=5;
 662     nSigmaCombRef=5;
 663     nSigmaTOF=5;
 664     nSigmaTOFRef=5;
 665     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 666     fPidKp=0;
 667
 668     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 669
 670     Double_t oldpP[10];
 671     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 672
 673
 674     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
 675     nSigmaTPCRef = PIDResponse->NumberOfSigmasTPC(KpTrack,(AliPID::EParticleType) fSpeciesRef);
 676
 677     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 678
 679     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 680     /*
 681     if(mcArray){
 682       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 683       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 684       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 685     }
 686     */
 687
 688     fPidKp = Int_t(probP[4]*100);
 689
 690     if(tofMatch1){
 691       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 692         // remove this tof hit
 693         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 694         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 695         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 696         fPidKp = Int_t(probP[4]*100);
 697         tofMatch1=0;
 698       }
 699       else{
 700         if(probP[4] > probP[3] && probP[4] > probP[2] && probP[4] > probP[0]) fPidKp += 128; // max prob
 701
 702         if(isLambda){
 703           nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 704           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 705           if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
 706           nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 707           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 708           if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
 709           nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 710           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 711           if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
 712         }
 713         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 714         nSigmaTOFRef = PIDResponse->NumberOfSigmasTOF(KpTrack,(AliPID::EParticleType) fSpeciesRef);
 715         if(isLambda){
 716           if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 717           if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
 718         }
 719
 720         if(fIsMC){
 721           Float_t mismAdd = addMismatchForMC;
 722           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 723
 724           if(gRandom->Rndm() < mismAdd){
 725             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 726             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 727             channel = channel % 8736;
 728             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 729
 730             // generate random time
 731             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 732             Double_t times[AliPID::kSPECIESC];
 733             KpTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 734             nSigmaTOF = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[4], AliPID::kProton);
 735             nSigmaTOFRef = TMath::Abs(timeRandom - times[fSpeciesRef])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[fSpeciesRef], (AliPID::EParticleType) fSpeciesRef);
 736           }
 737         }
 738
 739         if(fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 740         nSigmaTOF = TMath::Abs(nSigmaTOF);
 741
 742         if(nSigmaTOF < 2) fPidKp += 256;
 743         else if(nSigmaTOF < 3) fPidKp += 512;
 744       }
 745     }
 746
 747     if(tofMatch1){
 748       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 749       nSigmaCombRef = TMath::Sqrt(nSigmaTOFRef*nSigmaTOFRef + nSigmaTPCRef*nSigmaTPCRef);
 750       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2){
 751         fCombsignal->Fill(nSigmaComb);
 752       }
 753     } else {
 754       nSigmaComb = TMath::Abs(nSigmaTPC);
 755       nSigmaCombRef = TMath::Abs(nSigmaTPCRef);
 756     }
 757
 758     // use sigmaTOF instead of sigmaComb
 759     nSigmaTOFRef = TMath::Abs(nSigmaTOFRef);
 760
 761     if(tofMatch1){
 762       nSigmaComb = nSigmaTOF;
 763       nSigmaCombRef = nSigmaTOFRef;
 764     }
 765
 766     if(nSigmaComb < 2) nSigmaComb = 2;
 767     else if(nSigmaComb < 3) nSigmaComb = 3;
 768     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 769     else nSigmaComb = 6;
 770
 771     if(nSigmaCombRef < 2) nSigmaCombRef = 2;
 772     else if(nSigmaCombRef < 3) nSigmaCombRef = 3;
 773     else if(nSigmaCombRef < 5) nSigmaCombRef = 4.99;
 774     else nSigmaCombRef = 6;
 775
 776     Int_t iks=-1;
 777     for(Int_t k=0;k < fNLambda;k++){ // find the Lambda which contains the positive track
 778       if(i == fIpP[k]) iks = k;
 779     }
 780
 781     if(fPtKp > 4.299) fPtKp = 4.299;
 782
 783     if(iks > -1 && fIpN[iks] > -1){
 784       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 785       Int_t j = fIpN[iks];
 786       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 787
 788       if (!KnTrack){
 789         continue;
 790       }
 791
 792       nSigmaComb2 = 5;
 793       nSigmaTOF2 = 5;
 794
 795       if(!(KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 796
 797       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 798       fPidKn=0;
 799
 800       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 801       Double_t oldpN[10];
 802       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 803
 804       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 805
 806       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 807
 808       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 809       /*
 810       if(mcArray){
 811         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 812         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 813         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 814       }
 815       */
 816
 817       fPidKn = Int_t(probN[2]*100);
 818
 819       if(tofMatch2){
 820         if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 821
 822         nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 823
 824         if(fIsMC){
 825           Float_t mismAdd = addMismatchForMC;
 826           if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 827
 828           if(gRandom->Rndm() < mismAdd){
 829             Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 830             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 831             channel = channel % 8736;
 832             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 833
 834             // generate random time
 835             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 836             Double_t times[AliPID::kSPECIESC];
 837             KnTrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
 838             nSigmaTOF2 = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[4], AliPID::kProton);
 839           }
 840         }
 841
 842
 843
 844         nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 845
 846         if(nSigmaTOF2 < 2) fPidKn += 256;
 847         else if(nSigmaTOF2 < 3) fPidKn += 512;
 848       }
 849
 850       px = KpTrack->Px() + KnTrack->Px();
 851       py = KpTrack->Py() + KnTrack->Py();
 852       pz = KpTrack->Pz() + KnTrack->Pz();
 853       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 0.938e-01*0.938e-01);
 854       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
 855
 856       ksV.SetPxPyPzE(px,py,pz,E);
 857       fMassV0 = fMassLambda[iks];
 858
 859       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 860
 861
 862       fPtLambdaC = ksV.Pt();
 863       fEtaLambda = ksV.Eta();
 864       fPhiLambdaC = ksV.Phi();
 865
 866       if(tofMatch2){
 867         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 868       } else {
 869         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 870       }
 871
 872       // use sigmaTOF instead of sigmaComb
 873       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 874
 875       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 876       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 877       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 878       else nSigmaComb2 = 6;
 879
 880       Bool_t isTrue = kFALSE;
 881
 882       if(mcArray){
 883         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 884         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 885
 886         if(labelP > -1 && labelN > -1){
 887           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 888           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 889
 890           Int_t mP = partP->GetMother();
 891           Int_t mN = partN->GetMother();
 892           if(mP == mN && mP > -1){
 893             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 894             Int_t pdgM = partM->GetPdgCode();
 895             if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
 896           }
 897         }
 898
 899       }
 900
 901       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 902       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 903
 904       if(gRandom->Rndm() < 0.5){
 905         deltaphi1 += TMath::Pi();
 906         deltaphi2 += TMath::Pi();
 907       }
 908
 909       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 910       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 911       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 912       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 913
 914       if(fPtKn > 4.299) fPtKn = 4.299;
 915
 916       Float_t xTOfill[] =  {static_cast<Float_t>(fPtLambdaC),static_cast<Float_t>(KpTrack->Eta()),static_cast<Float_t>(fPtKp),static_cast<Float_t>(fPtKn),static_cast<Float_t>((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)};
 917       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)};
 918
 919
 920       Int_t ipt = 0;
 921       while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
 922         ipt++;
 923       }
 924       ipt--;
 925       if(ipt < 0) ipt = 0; // just to be sure
 926
 927       if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 928         if(fSpeciesRef != 4){
 929           xTOfill[4] = probP[fSpeciesRef];
 930           xTOfill2[5] = probP[fSpeciesRef];
 931
 932           xTOfill[9] = nSigmaCombRef;
 933           xTOfill2[10] = nSigmaCombRef;
 934
 935         }
 936
 937         if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 938         else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
 939
 940
 941
 942         if(fPIDuserCut){
 943           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)};
 944
 945          if(fPIDuserCut->IsSelected(KpTrack,AliPID::kPion)){ // to be filled
 946            xUser[3] = 1;
 947          } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kKaon)){
 948            xUser[3] = 2;
 949          } else if(fPIDuserCut->IsSelected(KpTrack,AliPID::kProton)){
 950            xUser[3] = 3;
 951          }
 952
 953          if(isLambda) fContUser->Fill(0,fMassV0,fCentrality,xUser);
 954          else fContUser2->Fill(0,fMassV0,fCentrality,xUser);
 955         }
 956       }
 957
 958
 959     }
 960   } // end analysi Lambda
 961 }
 962
 963 //_____________________________________________________________________________
 964 Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
 965 {
 966
 967   Float_t zvtx = -999;
 968
 969   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 970   if (!vtxAOD)
 971     return zvtx;
 972   if(vtxAOD->GetNContributors()>0)
 973     zvtx = vtxAOD->GetZ();
 974
 975   return zvtx;
 976 }
 977 //_____________________________________________________________________________
 978 void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
 979 {
 980   // Terminate loop
 981   Printf("Terminate()");
 982 }
 983 //=======================================================================
 984 void AliAnalysisTaskLambdaBayes::SelectLambda(){
 985   fNLambda=0;
 986   fNpPos=0;
 987   fNpNeg=0;
 988
 989   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
 990   AliAODv0 *myV0;
 991   Double_t dMASS=0.0;
 992   for (Int_t i=0; i!=nV0s; ++i) {
 993     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
 994     if(!myV0) continue;
 995     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
 996     Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
 997     if(pass) {
 998       if(pass==1) dMASS = myV0->MassLambda();
 999       else dMASS = myV0->MassAntiLambda();
1000
1001       Float_t massKs = myV0->MassK0Short();
1002
1003       if(TMath::Abs(dMASS-1.115)/0.005 < 8 && TMath::Abs(massKs - 0.497)/0.005 > 8){
1004         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
1005         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
1006         if(iT->Charge()<0){
1007           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
1008           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
1009         }
1010         fPhiLambda[fNLambda] = myV0->Phi();
1011         fPtLambda[fNLambda] = myV0->Pt();
1012         fIpP[fNLambda] = FindDaugheterIndex(iT);
1013         fIpN[fNLambda] = FindDaugheterIndex(jT);
1014
1015         if(pass==2){
1016           Int_t itemp = fIpP[fNLambda];
1017           fIpP[fNLambda] = fIpN[fNLambda];
1018           fIpN[fNLambda] = itemp;
1019         }
1020
1021         fMassLambda[fNLambda] = dMASS;
1022         if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
1023           fNLambda++;
1024         }
1025
1026       }
1027     }
1028   }
1029 }
1030
1031 //=======================================================================
1032 Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1033 {
1034   if (myV0->GetOnFlyStatus() ) return 0;
1035
1036   //the following is needed in order to evualuate track-quality
1037   AliAODTrack *iT, *jT;
1038   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1039   Double_t pos[3],cov[6];
1040   vV0s->GetXYZ(pos);
1041   vV0s->GetCovarianceMatrix(cov);
1042   const AliESDVertex vESD(pos,cov,100.,100);
1043
1044   // TESTING CHARGE
1045   int iPos, iNeg;
1046   iT=(AliAODTrack*) myV0->GetDaughter(0);
1047   if(iT->Charge()>0) {
1048     iPos = 0; iNeg = 1;
1049   } else {
1050     iPos = 1; iNeg = 0;
1051   }
1052   // END OF TEST
1053
1054   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1055
1056   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1057
1058   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1059   if(!trkFlag) return 0;
1060   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1061   if(!trkFlag2) return 0;
1062
1063   Double_t pvertex[3];
1064   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1065   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1066   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1067
1068   Double_t dDL=myV0->DecayLengthV0( pvertex );
1069   if(dDL  < 0.5 || dDL > 25) return 0;
1070
1071   Double_t dDCA=myV0->DcaV0Daughters();
1072   if(dDCA >0.5) return 0;
1073
1074   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1075   if(dCTP < -1) return 0;
1076
1077 //   AliESDtrack ieT( iT );
1078 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1079 //   if(TMath::Abs(dD0P) < 0]) return 0;
1080
1081 //   AliESDtrack jeT( jT );
1082 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1083 //   if(TMath::Abs(dD0M) < 0) return 0;
1084
1085 //   Double_t dD0D0=dD0P*dD0M;
1086 //   if(dD0D0>0) return 0;
1087
1088 //   Double_t dETA=myV0->Eta();
1089 //   if(dETA <-0.8) return 0;
1090 //   if(dETA >0.8) return 0;
1091
1092 //   Double_t dQT=myV0->PtArmV0();
1093 //   if(specie==0) if(dQT<???) return 0;
1094
1095   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1096   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1097
1098   if(specie==1 && dALPHA<0) return 2; // antilambda
1099   return 1; // Ks or lambda
1100 }
1101 //-------------------------------------------------------------------------
1102 Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
1103   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1104
1105   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1106     AliAODTrack* track = fOutputAOD->GetTrack(i);
1107     if(track == trk) return i;
1108   }
1109
1110   printf("Daughter for %p not found\n",trk);
1111   return -1;
1112 }
1113 //-------------------------------------------------------------------------
1114 Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
1115   if(!fIsMC) return 1; // used only on MC
1116
1117   if( fTypeCol==2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1118     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1119
1120     if(!(channel%20)) return 0; // 5% additional loss in MC
1121   }
1122
1123   return 1;
1124 }