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