add pid dir with some task to run PID performance
authorfnoferin <fnoferin@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Jul 2013 14:36:13 +0000 (14:36 +0000)
committerfnoferin <fnoferin@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Jul 2013 14:36:13 +0000 (14:36 +0000)
14 files changed:
PWGPP/pid/AddTaskK0sBayes.C [new file with mode: 0644]
PWGPP/pid/AddTaskLambdaBayes.C [new file with mode: 0644]
PWGPP/pid/AddTaskPhiBayes.C [new file with mode: 0644]
PWGPP/pid/AliAnalysisTaskK0sBayes.cxx [new file with mode: 0644]
PWGPP/pid/AliAnalysisTaskK0sBayes.h [new file with mode: 0644]
PWGPP/pid/AliAnalysisTaskLambdaBayes.cxx [new file with mode: 0644]
PWGPP/pid/AliAnalysisTaskLambdaBayes.h [new file with mode: 0644]
PWGPP/pid/AliAnalysisTaskPhiBayes.cxx [new file with mode: 0644]
PWGPP/pid/AliAnalysisTaskPhiBayes.h [new file with mode: 0644]
PWGPP/pid/AliPIDperfContainer.cxx [new file with mode: 0644]
PWGPP/pid/AliPIDperfContainer.h [new file with mode: 0644]
PWGPP/pid/doeffKa.C [new file with mode: 0644]
PWGPP/pid/doeffPi.C [new file with mode: 0644]
PWGPP/pid/doeffPr.C [new file with mode: 0644]

diff --git a/PWGPP/pid/AddTaskK0sBayes.C b/PWGPP/pid/AddTaskK0sBayes.C
new file mode 100644 (file)
index 0000000..b88d8cc
--- /dev/null
@@ -0,0 +1,43 @@
+AliAnalysisTask *AddTaskK0sBayes(Bool_t ismc=kFALSE,Bool_t qa=kTRUE,Int_t filterbit=4){
+
+  //get the current analysis manager
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("No manager found in AddTaskVZERO. Why?");
+    return 0;
+  }
+  // currently don't accept AOD input
+  if (!mgr->GetInputEventHandler()->InheritsFrom(AliAODInputHandler::Class())) {
+    Error("AddTaskK0sBayes","This task works only with AOD input!");
+    return 0;
+  }
+
+  //========= Add tender to the ANALYSIS manager and set default storage =====
+  char mytaskName[100];
+  snprintf(mytaskName,100,"AliAnalysisTaskK0sBayes.cxx"); 
+
+  AliAnalysisTaskK0sBayes *task = new AliAnalysisTaskK0sBayes(mytaskName);
+  if(ismc) task->SetMC();
+  if(qa) task->SetQA();
+  task->SetEtaCut(0.8);
+  task->SetFilterBit(filterbit);
+
+  mgr->AddTask(task);
+
+  //Attach input to my tasks
+  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+
+  // Attach output to my tasks
+  AliAnalysisDataContainer *cOutputL= mgr->CreateContainer("contK0sBayes1",TList::Class(), AliAnalysisManager::kOutputContainer, AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 1, cOutputL);
+
+  AliAnalysisDataContainer *cOutputL2= mgr->CreateContainer("contK0sBayes2",TList::Class(), AliAnalysisManager::kOutputContainer, AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 2, cOutputL2);
+
+  AliAnalysisDataContainer *cOutputL3= mgr->CreateContainer("contK0sBayes3",TList::Class(), AliAnalysisManager::kOutputContainer, AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 3, cOutputL3);
+
+  return task;
+}
+
diff --git a/PWGPP/pid/AddTaskLambdaBayes.C b/PWGPP/pid/AddTaskLambdaBayes.C
new file mode 100644 (file)
index 0000000..c1dd5c6
--- /dev/null
@@ -0,0 +1,46 @@
+AliAnalysisTask *AddTaskLambdaBayes(Bool_t ismc=kFALSE,Bool_t qa=kTRUE,Int_t filterbit=4){
+
+  //get the current analysis manager
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("No manager found in AddTaskVZERO. Why?");
+    return 0;
+  }
+  // currently don't accept AOD input
+  if (!mgr->GetInputEventHandler()->InheritsFrom(AliAODInputHandler::Class())) {
+    Error("AddTaskLambdaBayes","This task works only with AOD input!");
+    return 0;
+  }
+
+  //========= Add tender to the ANALYSIS manager and set default storage =====
+  char mytaskName[100];
+  snprintf(mytaskName,100,"AliAnalysisTaskLambdaBayes.cxx"); 
+
+  AliAnalysisTaskLambdaBayes *task = new AliAnalysisTaskLambdaBayes(mytaskName);
+  if(ismc) task->SetMC();
+  if(qa) task->SetQA();
+  task->SetEtaCut(0.8);
+  task->SetFilterBit(filterbit);
+
+  mgr->AddTask(task);
+
+  //Attach input to my tasks
+  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+
+  // Attach output to my tasks
+  AliAnalysisDataContainer *cOutputL= mgr->CreateContainer("contLambdaBayes1",TList::Class(), AliAnalysisManager::kOutputContainer, 
+AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 1, cOutputL);
+
+  AliAnalysisDataContainer *cOutputL2= mgr->CreateContainer("contLambdaBayes2",TList::Class(), AliAnalysisManager::kOutputContainer, 
+AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 2, cOutputL2);
+
+  AliAnalysisDataContainer *cOutputL3= mgr->CreateContainer("contLambdaBayes3",TList::Class(), AliAnalysisManager::kOutputContainer, 
+AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 3, cOutputL3);
+
+  return task;
+}
+
diff --git a/PWGPP/pid/AddTaskPhiBayes.C b/PWGPP/pid/AddTaskPhiBayes.C
new file mode 100644 (file)
index 0000000..606b608
--- /dev/null
@@ -0,0 +1,43 @@
+AliAnalysisTask *AddTaskPhiBayes(Bool_t ismc=kFALSE,Bool_t qa=kTRUE,Int_t filterbit=16){
+
+  //get the current analysis manager
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("No manager found in AddTaskVZERO. Why?");
+    return 0;
+  }
+  // currently don't accept AOD input
+  if (!mgr->GetInputEventHandler()->InheritsFrom(AliAODInputHandler::Class())) {
+    Error("AddTaskPhiBayes","This task works only with AOD input!");
+    return 0;
+  }
+
+  //========= Add tender to the ANALYSIS manager and set default storage =====
+  char mytaskName[100];
+  snprintf(mytaskName,100,"AliAnalysisTaskPhiBayes.cxx"); 
+
+  AliAnalysisTaskPhiBayes *task = new AliAnalysisTaskPhiBayes(mytaskName);
+  if(ismc) task->SetMC();
+  if(qa) task->SetQA();
+  task->SetEtaCut(0.8);
+  task->SetFilterBit(filterbit);
+
+  mgr->AddTask(task);
+
+  //Attach input to my tasks
+  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+
+  // Attach output to my tasks
+  AliAnalysisDataContainer *cOutputL= mgr->CreateContainer("contPhiBayes1",TList::Class(), AliAnalysisManager::kOutputContainer, AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 1, cOutputL);
+
+  AliAnalysisDataContainer *cOutputL2= mgr->CreateContainer("contPhiBayes2",TList::Class(), AliAnalysisManager::kOutputContainer, AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 2, cOutputL2);
+
+  AliAnalysisDataContainer *cOutputL3= mgr->CreateContainer("contPhiBayes3",TList::Class(), AliAnalysisManager::kOutputContainer, AliAnalysisManager::GetCommonFileName());
+  mgr->ConnectOutput(task, 3, cOutputL3);
+
+  return task;
+}
+
diff --git a/PWGPP/pid/AliAnalysisTaskK0sBayes.cxx b/PWGPP/pid/AliAnalysisTaskK0sBayes.cxx
new file mode 100644 (file)
index 0000000..4ef22b2
--- /dev/null
@@ -0,0 +1,1043 @@
+#include "AliAnalysisTaskK0sBayes.h"
+
+// ROOT includes
+#include <TMath.h>
+
+// AliRoot includes
+#include "AliInputEventHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliCentrality.h"
+#include "AliVHeader.h"
+#include "AliAODVZERO.h"
+#include "TFile.h"
+#include "AliOADBContainer.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliMCEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "TChain.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDVertex.h"
+#include "AliEventplane.h"
+#include "AliAnalysisManager.h"
+#include "TRandom.h"
+
+
+Float_t AliAnalysisTaskK0sBayes::fPtKsMin[AliAnalysisTaskK0sBayes::nPtBin] = {0.,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.};// ptmin bin
+Float_t AliAnalysisTaskK0sBayes::fPtKsMax[AliAnalysisTaskK0sBayes::nPtBin] = {0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.,10.};// ptmax bin
+
+// STL includes
+//#include <iostream>
+//using namespace std;
+
+ClassImp(AliAnalysisTaskK0sBayes)
+//_____________________________________________________________________________
+AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes():
+  AliAnalysisTaskSE(),
+  fVtxCut(10.0),  // cut on |vertex| < fVtxCut
+  fEtaCut(0.8),   // cut on |eta| < fEtaCut
+  fMinPt(0.15),   // cut on pt > fMinPt
+  fIsMC(kFALSE),
+  fQAsw(kFALSE),
+  fNcluster(70),
+  fFilterBit(1),
+  fList(new TList()),
+  fList2(new TList()),
+  fList3(new TList()),
+  fCentrality(-1),
+  fPsi(0),
+  fPtKs(0.),
+  fPhiKs(0.),
+  fEtaKs(0.),
+  fMassV0(0.),
+  fPtKp(0.),
+  fPhiKp(0.),
+  fEtaKp(0.),
+  fPtKn(0.),
+  fPhiKn(0.),
+  fEtaKn(0.),
+  fPidKp(0),
+  fPidKn(0),
+  fTOFTPCsignal(0),
+  fCombsignal(0),
+  fCutsDaughter(NULL),
+  fPIDCombined(NULL),
+  fContPid(NULL),
+  fContPid2(NULL),
+  fNK0s(0),
+  fNpiPos(0),
+  fNpiNeg(0),
+  fHmismTOF(0),
+  fHchannelTOFdistr(0)
+{
+  // Default constructor (should not be used)
+  fList->SetName("contKsBayes1");
+  fList2->SetName("contKsBayes2");
+  fList3->SetName("contKsBayes3");
+
+  fList->SetOwner(kTRUE); 
+  fList2->SetOwner(kTRUE); 
+  fList3->SetOwner(kTRUE); 
+
+  TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
+  fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
+
+  TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+  fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
+}
+
+//______________________________________________________________________________
+AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes(const char *name):
+  AliAnalysisTaskSE(name),
+  fVtxCut(10.0),  // cut on |vertex| < fVtxCut
+  fEtaCut(0.8),   // cut on |eta| < fEtaCut
+  fMinPt(0.15),   // cut on pt > fMinPt
+  fIsMC(kFALSE),
+  fQAsw(kFALSE),
+  fNcluster(70),
+  fFilterBit(1),
+  fList(new TList()),
+  fList2(new TList()),
+  fList3(new TList()),
+  fCentrality(-1),
+  fPsi(0),
+  fPtKs(0.),
+  fPhiKs(0.),
+  fEtaKs(0.),
+  fMassV0(0.),
+  fPtKp(0.),
+  fPhiKp(0.),
+  fEtaKp(0.),
+  fPtKn(0.),
+  fPhiKn(0.),
+  fEtaKn(0.),
+  fPidKp(0),
+  fPidKn(0),
+  fTOFTPCsignal(0),
+  fCombsignal(0),
+  fCutsDaughter(NULL),
+  fPIDCombined(NULL),
+  fContPid(NULL),
+  fContPid2(NULL),
+  fNK0s(0),
+  fNpiPos(0),
+  fNpiNeg(0),
+  fHmismTOF(0),
+  fHchannelTOFdistr(0)
+{
+
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
+  DefineOutput(3, TList::Class());
+  DefineOutput(4, TList::Class());
+
+  // Output slot #1 writes into a TTree
+  fList->SetName("contKsBayes1");
+  fList2->SetName("contKsBayes2");
+  fList3->SetName("contKsBayes3");
+
+  fList->SetOwner(kTRUE); 
+  fList2->SetOwner(kTRUE); 
+  fList3->SetOwner(kTRUE); 
+
+  TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
+  fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
+
+  TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+  fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
+}
+//_____________________________________________________________________________
+AliAnalysisTaskK0sBayes::~AliAnalysisTaskK0sBayes()
+{
+
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskK0sBayes::UserCreateOutputObjects()
+{
+
+  fPIDCombined=new AliPIDCombined;
+  fPIDCombined->SetDefaultTPCPriors();
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+
+  Float_t invMmin = 0.497-0.005*8;
+  Float_t invMmax = 0.497+0.005*8;
+
+  const Int_t nBinPid = 14;
+
+  Int_t binPid[nBinPid] = {1/*ptKs*/,8/*EtaPiP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+
+  Int_t binPid2[nBinPid] = {1/*ptKs*/,8/*EtaPiN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+
+  fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
+  fContPid->SetTitleX("M_{K^{0}_{s}}");
+  fContPid->SetTitleY("centrality (%)");
+  fContPid->SetVarName(0,"p_{T}^{K^{0}_{s}}");
+  fContPid->SetVarRange(0,0,10);
+  fContPid->SetVarName(1,"#eta^{K^{0}_{s}}");
+  fContPid->SetVarRange(1,-0.8,0.8);
+  fContPid->SetVarName(2,"p_{T}^{Kp}");
+  fContPid->SetVarRange(2,0.3,4.3);
+  fContPid->SetVarName(3,"p_{T}^{Kn}");
+  fContPid->SetVarRange(3,0.3,4.3);
+  fContPid->SetVarName(4,"BayesProb^{Kp}");
+  fContPid->SetVarRange(4,0,1.);
+  fContPid->SetVarName(5,"BayesProb^{Kn}");
+  fContPid->SetVarRange(5,0,1.);
+  fContPid->SetVarName(6,"isTOFmatch^{Kp}");
+  fContPid->SetVarRange(6,-0.5,1.5);
+  fContPid->SetVarName(7,"isTOFmatch^{Kn}");
+  fContPid->SetVarRange(7,-0.5,1.5);
+  fContPid->SetVarName(8,"isKsTrue^{Kn}");
+  fContPid->SetVarRange(8,-0.5,1.5);
+  fContPid->SetVarName(9,"N#sigma^{Kp}");
+  fContPid->SetVarRange(9,1.25,6.25);
+  fContPid->SetVarName(10,"N#sigma^{Kn}");
+  fContPid->SetVarRange(10,1.25,6.25);
+  fContPid->SetVarName(11,"#Delta#phi^{Kp}");
+  fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
+  fContPid->SetVarName(12,"#Delta#phi^{Kn}");
+  fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
+  fContPid->SetVarName(13,"#Psi");
+  fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
+
+  fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
+  fContPid2->SetTitleX("M_{K^{0}_{s}}");
+  fContPid2->SetTitleY("centrality (%)");
+  fContPid2->SetVarName(0,"p_{T}^{K^{0}_{s}}");
+  fContPid2->SetVarRange(0,0,10);
+  fContPid2->SetVarName(1,"#eta^{K^{0}_{s}}");
+  fContPid2->SetVarRange(1,-0.8,0.8);
+  fContPid2->SetVarName(2,"p_{T}^{Kp}");
+  fContPid2->SetVarRange(2,0.3,4.3);
+  fContPid2->SetVarName(3,"p_{T}^{Kn}");
+  fContPid2->SetVarRange(3,0.3,4.3);
+  fContPid2->SetVarName(4,"BayesProb^{Kp}");
+  fContPid2->SetVarRange(4,0,1.);
+  fContPid2->SetVarName(5,"BayesProb^{Kn}");
+  fContPid2->SetVarRange(5,0,1.);
+  fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
+  fContPid2->SetVarRange(6,-0.5,1.5);
+  fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
+  fContPid2->SetVarRange(7,-0.5,1.5);
+  fContPid2->SetVarName(8,"isKsTrue^{Kn}");
+  fContPid2->SetVarRange(8,-0.5,1.5);
+  fContPid2->SetVarName(9,"N#sigma^{Kp}");
+  fContPid2->SetVarRange(9,1.25,6.25);
+  fContPid2->SetVarName(10,"N#sigma^{Kn}");
+  fContPid2->SetVarRange(10,1.25,6.25);
+  fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
+  fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
+  fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
+  fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
+  fContPid2->SetVarName(13,"#Psi");
+  fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
+
+
+
+  const Int_t nDETsignal = 100; // mass
+  Double_t binDETsignal[nDETsignal+1];
+  for(Int_t i=0;i<nDETsignal+1;i++){
+    binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
+  }
+  const Int_t nDETsignal2 = 10; // centrality
+  Double_t binDETsignal2[nDETsignal2+1];
+  for(Int_t i=0;i<nDETsignal2+1;i++){
+    binDETsignal2[i] = i*100./ nDETsignal2;
+  }
+  fContPid->AddSpecies("K0s",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
+  fContPid2->AddSpecies("K0s2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
+
+  fList->Add(fContPid);
+  fList->Add(fContPid2);
+
+  hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+
+  hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+
+  fList2->Add(hMatching[0]);
+  fList2->Add(hMatching[1]);
+  fList2->Add(hMatching[2]);
+  fList2->Add(hMatching[3]);
+  fList2->Add(hTracking[0]);
+  fList2->Add(hTracking[1]);
+  fList2->Add(hTracking[2]);
+  fList2->Add(hTracking[3]);
+
+  fTOFTPCsignal = new TH2F("hTOFTPCsignal","TOF-TPC signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma_{TOF};N#sigma_{TPC}",100,-5,5,100,-5,5);
+  fList2->Add(fTOFTPCsignal);
+  fCombsignal = new TH1F("hCombsignal","Combined signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma",100,0,5);
+  fList2->Add(fCombsignal);
+
+  // QA plots
+  char name[100],title[400];
+  for(Int_t i=0;i < nCentrBin;i++){
+    snprintf(name,100,"hPiTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPiTPC[i]);
+
+    snprintf(name,100,"hKaTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fKaTPC[i]);
+
+    snprintf(name,100,"hPrTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPrTPC[i]);
+
+    snprintf(name,100,"hElTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fElTPC[i]);
+
+    snprintf(name,100,"hPiTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPiTOF[i]);
+
+    snprintf(name,100,"hKaTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fKaTOF[i]);
+
+    snprintf(name,100,"hPrTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPrTOF[i]);
+
+    snprintf(name,100,"hElTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fElTOF[i]);
+
+  }
+
+  // Post output data.
+  PostData(1, fList);
+  PostData(2, fList2);
+  PostData(3, fList3);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskK0sBayes::UserExec(Option_t *) 
+{
+    // Main loop
+    // Called for each event
+    
+    fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    if(!fOutputAOD){
+       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+       this->Dump();
+       return;
+    }
+    
+    Float_t zvtx = GetVertex(fOutputAOD);
+
+
+
+    //Get the MC object
+    if(fIsMC){
+      AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+      if (!mcHeader) {
+       AliError("Could not find MC Header in AOD");
+       return;
+      }
+    }
+
+    if (TMath::Abs(zvtx) < fVtxCut) {
+
+      SelectK0s();
+
+      //Centrality
+      Float_t v0Centr  = -10.;
+      Float_t trkCentr  = -10.;
+      AliCentrality *centrality = fOutputAOD->GetCentrality();
+      if (centrality){
+       trkCentr  = centrality->GetCentralityPercentile("V0M");
+       v0Centr = centrality->GetCentralityPercentile("TRK"); 
+      }
+
+      if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr>0){ // consistency cut on centrality selection
+        fCentrality = v0Centr;
+       Analyze(fOutputAOD); // Do analysis!!!!
+
+      }
+    }
+    
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskK0sBayes::Analyze(AliAODEvent* aodEvent)
+{
+
+  Int_t ntrack = aodEvent->GetNumberOfTracks();
+
+  fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
+  fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
+  fPidKp=0,fPidKn=0;
+  fMassV0=-1;
+  
+  TLorentzVector ksV;
+  
+  Double_t px,py,pz,E;
+
+  Float_t invMmin = 0.497-0.005*8;
+  Float_t invMmax = 0.497+0.005*8;
+  
+  Int_t icentr = 8;
+  if(fCentrality < 0) icentr = 8;
+  else if(fCentrality < 10) icentr = 0;
+  else if(fCentrality < 20) icentr = 1;
+  else if(fCentrality < 30) icentr = 2;
+  else if(fCentrality < 40) icentr = 3;
+  else if(fCentrality < 50) icentr = 4;
+  else if(fCentrality < 60) icentr = 5;
+  else if(fCentrality < 70) icentr = 6;
+  else if(fCentrality < 80) icentr = 7;
+
+  Float_t addMismatchForMC = 0.005;
+  if(fCentrality < 50) addMismatchForMC += 0.01;
+  if(fCentrality < 20) addMismatchForMC += 0.02;
+
+  fPsi = 0;
+  /* Compute TPC EP */
+  Double_t Qx2 = 0, Qy2 = 0;
+  Double_t Qx3 = 0, Qy3 = 0;
+  for(Int_t iT = 0; iT < ntrack; iT++) {
+    AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
+    
+    if (!aodTrack){
+      continue;
+    }
+    
+    Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
+    
+    if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag) 
+      continue;
+    
+    Double_t b[2] = {-99., -99.};
+    Double_t bCov[3] = {-99., -99., -99.};
+    if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
+      continue;
+    
+    if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
+      continue;
+    
+    Qx2 += TMath::Cos(2*aodTrack->Phi()); 
+    Qy2 += TMath::Sin(2*aodTrack->Phi());
+    Qx3 += TMath::Cos(3*aodTrack->Phi()); 
+    Qy3 += TMath::Sin(3*aodTrack->Phi());
+    
+  }
+
+  fPsi = TMath::ATan2(Qy2, Qx2)/2.;
+
+  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+  AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
+
+  PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
+  PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
+  PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
+  PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
+
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+
+  Double_t probP[10],probN[10];
+  Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
+
+  
+  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+  TClonesArray *mcArray = NULL;
+  if (mcHeader)
+    mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+
+  Int_t nmc = 0;
+  if(mcArray)
+    nmc = mcArray->GetEntries();
+
+  for(Int_t i=0;i < ntrack;i++){ // loop on tracks
+    AliAODTrack* track = aodEvent->GetTrack(i);
+    
+    AliAODMCParticle *mcp = NULL;
+    Int_t pdg = 0;
+    
+    if (!track){
+      continue;
+    }
+    
+    Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
+    
+    Int_t label = -1;
+    if(mcArray){
+      label = track->GetLabel();
+      if(label != -1 && label < nmc){
+       label = TMath::Abs(label);
+       mcp = (AliAODMCParticle*)mcArray->At(label);
+       pdg = TMath::Abs(mcp->GetPdgCode());
+      }
+      else
+       label = -1;
+    }
+    else{
+      /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
+    }
+    
+    if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
+      hTracking[0]->Fill(track->Pt(),fCentrality);
+      if(pdg == 211)
+       hTracking[1]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 321)
+       hTracking[2]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 2212)
+       hTracking[3]->Fill(track->Pt(),fCentrality);
+      else if(! mcp){ // Fill matching histo with the prob
+       hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
+       hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
+       hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
+      }
+    }
+    
+    if(!tofMatch) continue;
+    
+    if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
+      hMatching[0]->Fill(track->Pt(),fCentrality);
+      if(pdg == 211)
+       hMatching[1]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 321)
+       hMatching[2]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 2212)
+       hMatching[3]->Fill(track->Pt(),fCentrality);
+      else if(! mcp){ // Fill matching histo with the prob
+       hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
+       hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
+       hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
+      }
+    }
+  }
+  
+  Int_t pdg1 = -1;
+  Int_t pdg2 = -1;
+
+
+  // start analysis K0s
+  for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
+    AliAODTrack* KpTrack = aodEvent->GetTrack(i);
+        
+    if (!KpTrack){
+      continue;
+    }
+    
+    if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
+
+    nSigmaComb=5;
+    fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
+    fPidKp=0;
+
+    UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
+
+    Double_t oldpP[10];
+    fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
+
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
+    fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
+    fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
+    fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
+    fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+
+    if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
+
+    Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
+
+    if(mcArray){
+      Int_t labelK = TMath::Abs(KpTrack->GetLabel());
+      AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
+      pdg1 = TMath::Abs(mcp1->GetPdgCode());
+    }
+
+    fPidKp = Int_t(probP[2]*100);
+
+    if(tofMatch1){
+      if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
+       // remove this tof hit
+       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+       detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
+       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+       fPidKp = Int_t(probP[4]*100);
+       tofMatch1=0;
+      }
+      else{
+       if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
+       
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
+       fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
+       fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
+       fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
+       fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+               
+       nSigmaTOF = TMath::Abs(nSigmaTOF);
+       
+       if(fIsMC){
+         Float_t mismAdd = addMismatchForMC;
+         if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
+         
+         if(gRandom->Rndm() < mismAdd){
+           Float_t etaAbs = TMath::Abs(KpTrack->Eta());
+           Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
+           channel = channel % 8736;
+           Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
+           
+           // generate random time
+           Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
+           Double_t times[10];
+           KpTrack->GetIntegratedTimes(times);
+           nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kPion);
+         }
+       }
+
+       if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
+
+       if(nSigmaTOF < 2) fPidKp += 256;
+       else if(nSigmaTOF < 3) fPidKp += 512;
+      }
+    }
+    
+    if(tofMatch1){
+      nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
+      if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
+       fCombsignal->Fill(nSigmaComb);
+      }
+    } else {
+      nSigmaComb = TMath::Abs(nSigmaTPC);
+    }
+
+    // use sigmaTOF instead of sigmaComb
+    //nSigmaComb = nSigmaTOF;
+
+    if(nSigmaComb < 2) nSigmaComb = 2;
+    else if(nSigmaComb < 3) nSigmaComb = 3;
+    else if(nSigmaComb < 5) nSigmaComb = 4.99;
+    else nSigmaComb = 6;
+
+    Int_t iks=-1;
+    for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
+      if(i == fIpiP[k]) iks = k;
+    }
+
+    if(fPtKp > 4.299) fPtKp = 4.299;
+
+    if(iks > -1 && fIpiN[iks] > -1){
+      //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
+      Int_t j = fIpiN[iks];
+      AliAODTrack* KnTrack = aodEvent->GetTrack(j);
+      
+      if (!KnTrack){
+       continue;
+      }
+
+      if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
+
+      fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
+      fPidKn=0;
+
+      UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
+      Double_t oldpN[10];
+      fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
+
+      nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
+      
+      if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
+      
+      nSigmaComb2=5;
+
+      Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
+
+      if(mcArray){
+       Int_t labelK = TMath::Abs(KnTrack->GetLabel());
+       AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
+       pdg2 = TMath::Abs(mcp2->GetPdgCode());
+     }
+
+      fPidKn = Int_t(probN[2]*100);
+
+      if(tofMatch2){
+       if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
+         // remove this tof hit
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+         detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+         fPidKn = Int_t(probN[4]*100);
+         tofMatch2=0;
+       }
+       else{
+         if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
+         
+         nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
+                 
+         nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
+         
+         if(fIsMC){
+           Float_t mismAdd = addMismatchForMC;
+           if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
+           
+           if(gRandom->Rndm() < mismAdd){
+             Float_t etaAbs = TMath::Abs(KnTrack->Eta());
+             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
+             channel = channel % 8736;
+             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
+             
+             // generate random time
+             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
+             Double_t times[10];
+             KnTrack->GetIntegratedTimes(times);
+             nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
+           }
+         }
+
+         if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
+
+         if(nSigmaTOF2 < 2) fPidKn += 256;
+         else if(nSigmaTOF2 < 3) fPidKn += 512;
+       }
+      }
+
+      px = KpTrack->Px() + KnTrack->Px();
+      py = KpTrack->Py() + KnTrack->Py();
+      pz = KpTrack->Pz() + KnTrack->Pz();
+      E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
+      E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
+
+      ksV.SetPxPyPzE(px,py,pz,E);
+      fMassV0 = fMassKs[iks];
+      
+      if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
+
+
+      fPtKs = ksV.Pt();
+      fEtaKs = ksV.Eta();
+      fPhiKs = ksV.Phi();
+
+      if(tofMatch2){
+       nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
+       if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
+         fCombsignal->Fill(nSigmaComb2);
+       }
+      } else {
+       nSigmaComb2 = TMath::Abs(nSigmaTPC2);
+      }
+
+      // use sigmaTOF instead of sigmaComb
+      //nSigmaComb2 = nSigmaTOF2;
+
+      if(nSigmaComb2 < 2) nSigmaComb2 = 2;
+      else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
+      else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
+      else nSigmaComb2 = 6;  
+
+      Bool_t isTrue = kFALSE;
+
+      if(mcArray){
+       Int_t labelP = TMath::Abs(KpTrack->GetLabel());
+       Int_t labelN = TMath::Abs(KnTrack->GetLabel());
+
+       if(labelP > -1 && labelN > -1){
+         AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
+         AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
+
+         Int_t mP = partP->GetMother();
+         Int_t mN = partN->GetMother();
+         if(mP == mN && mP > -1){
+           AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
+           Int_t pdgM = partM->GetPdgCode();
+           if(pdgM == 310) isTrue = kTRUE;
+         }
+       }
+
+      }
+
+      Double_t deltaphi1 = KpTrack->Phi() - fPsi;
+      Double_t deltaphi2 = KnTrack->Phi() - fPsi;
+
+      while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
+      while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
+      while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
+      while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
+
+      if(fPtKn > 4.299) fPtKn = 4.299;
+
+      Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
+      
+
+      Int_t ipt = 0;
+      while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
+       ipt++;
+      }
+      ipt--;
+      if(ipt < 0) ipt = 0; // just to be sure
+
+      if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
+       fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
+        xTOfill[1] = KnTrack->Eta();
+       fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
+      }
+
+
+    }
+  } // end analysi K0s
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
+{
+
+  Float_t zvtx = -999;
+
+  const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
+  if (!vtxAOD)
+    return zvtx;
+  if(vtxAOD->GetNContributors()>0)
+    zvtx = vtxAOD->GetZ();
+  
+  return zvtx;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
+{ 
+  // Terminate loop
+  Printf("Terminate()");
+}
+//=======================================================================
+void AliAnalysisTaskK0sBayes::SelectK0s(){
+  fNK0s=0;
+  fNpiPos=0;
+  fNpiNeg=0;
+
+  Int_t nV0s = fOutputAOD->GetNumberOfV0s();
+  AliAODv0 *myV0;
+  Double_t dMASS=0.0;
+  for (Int_t i=0; i!=nV0s; ++i) {
+    myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
+    if(!myV0) continue;
+    if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
+    Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
+    if(pass) {
+      dMASS = myV0->MassK0Short();
+      if(TMath::Abs(dMASS-0.497)/0.005 < 8){
+       AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
+       AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
+       if(iT->Charge()<0){
+         iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
+         jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
+       }
+       fPhiK0s[fNK0s] = myV0->Phi();
+       fPtK0s[fNK0s] = myV0->Pt();
+       fIpiP[fNK0s] = FindDaugheterIndex(iT);
+       fIpiN[fNK0s] = FindDaugheterIndex(jT);
+       fMassKs[fNK0s] = dMASS;
+       if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
+         fNK0s++;
+      }
+    }
+  }
+
+  /* My V0 code
+  // fill pion stacks
+  Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
+  for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
+    AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
+    
+    if (!aodTrack){
+      continue;
+    }
+    
+    Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
+
+    if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
+      continue;
+    }
+
+    Double_t b[2] = {-99., -99.};
+    Double_t bCov[3] = {-99., -99., -99.};
+    if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
+      continue;
+    
+    if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
+
+
+    Int_t charge = aodTrack->Charge();
+    if(charge > 0){
+      fIPiPos[fNpiPos] = iT;
+      fNpiPos++;
+    }
+    else{
+      fIPiNeg[fNpiNeg] = iT;
+      fNpiNeg++;
+    }     
+  }
+  
+  for(Int_t i=0;i < fNpiPos;i++){
+    AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
+    AliESDtrack pipE(pip);
+
+    for(Int_t j=0;j < fNpiNeg;j++){
+      AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
+      AliESDtrack pinE(pin);
+
+      Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
+
+      Double_t pPos[3];
+      Double_t pNeg[3];
+      pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
+      pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
+
+      Float_t length = (xp+xn)*0.5;
+
+      Float_t pxs = pPos[0] + pNeg[0];
+      Float_t pys = pPos[1] + pNeg[1];
+      Float_t pzs = pPos[2] + pNeg[2];
+      Float_t es = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
+
+      Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
+      Float_t phi = TMath::ATan2(pys,pxs);
+      Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
+      
+      //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
+
+      if(mindist < 0.4&& length > 0.7 && length < 25){
+        Float_t esL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.938*0.938) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
+        Float_t esAL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.938*0.938);
+
+        Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
+        Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
+
+        if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
+          fPhiK0s[fNK0s] = phi;
+          fPtK0s[fNK0s] = pt;
+         fIpiP[fNK0s] =fIPiPos[i] ;
+         fIpiN[fNK0s] = fIPiNeg[j];
+          fMassKs[fNK0s] = mass;
+          fNK0s++;
+        }
+      }
+    }
+  }
+  */
+}
+
+//=======================================================================
+Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
+{ 
+  if (myV0->GetOnFlyStatus() ) return 0;
+  
+  //the following is needed in order to evualuate track-quality
+  AliAODTrack *iT, *jT;
+  AliAODVertex *vV0s = myV0->GetSecondaryVtx();
+  Double_t pos[3],cov[6];
+  vV0s->GetXYZ(pos);
+  vV0s->GetCovarianceMatrix(cov);
+  const AliESDVertex vESD(pos,cov,100.,100);
+  
+  // TESTING CHARGE
+  int iPos, iNeg;
+  iT=(AliAODTrack*) myV0->GetDaughter(0);
+  if(iT->Charge()>0) {
+    iPos = 0; iNeg = 1;
+  } else {
+    iPos = 1; iNeg = 0;
+  }
+  // END OF TEST
+
+  iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
+
+  jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
+
+  Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
+  if(!trkFlag) return 0;
+  Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
+  if(!trkFlag2) return 0;
+
+  Double_t pvertex[3];
+  pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
+  pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
+  pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
+
+  Double_t dDL=myV0->DecayLengthV0( pvertex );
+  if(dDL  < 0.5 || dDL > 25) return 0;
+
+  Double_t dDCA=myV0->DcaV0Daughters();
+  if(dDCA >0.5) return 0;
+
+  Double_t dCTP=myV0->CosPointingAngle( pvertex );
+  if(dCTP < -1) return 0;
+
+//   AliESDtrack ieT( iT );
+//   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
+//   if(TMath::Abs(dD0P) < 0]) return 0;
+
+//   AliESDtrack jeT( jT );
+//   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
+//   if(TMath::Abs(dD0M) < 0) return 0;
+
+//   Double_t dD0D0=dD0P*dD0M;
+//   if(dD0D0>0) return 0;
+
+//   Double_t dETA=myV0->Eta();
+//   if(dETA <-0.8) return 0;
+//   if(dETA >0.8) return 0;
+
+//   Double_t dQT=myV0->PtArmV0();
+//   if(specie==0) if(dQT<???) return 0;
+
+  Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
+  if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
+
+  if(specie==1 && dALPHA<0) return 2; // antilambda
+  return 1; // K0s or lambda
+}
+//-------------------------------------------------------------------------
+Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
+  Int_t ntrack = fOutputAOD->GetNumberOfTracks();
+
+  for(Int_t i=0;i < ntrack;i++){ // loop on tracks
+    AliAODTrack* track = fOutputAOD->GetTrack(i);
+    if(track == trk) return i;
+  }
+  
+  printf("Daughter for %p not found\n",trk);
+  return -1;
+}
+//-------------------------------------------------------------------------
+Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
+  if(!fIsMC) return 1; // used only on MC
+
+  Int_t run = fOutputAOD->GetRunNumber();
+  if( (run>=136851&&run<=139846) || (run>=165772 && run<=170718) ){ // LHC10h or LHC11h because of TOF matching window at 3 cm
+    Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs); 
+  
+    if(!(channel%20)) return 0; // 5% additional loss in MC
+  }
+
+  return 1;
+}
diff --git a/PWGPP/pid/AliAnalysisTaskK0sBayes.h b/PWGPP/pid/AliAnalysisTaskK0sBayes.h
new file mode 100644 (file)
index 0000000..cb01dc6
--- /dev/null
@@ -0,0 +1,141 @@
+#ifndef ALIANALYSISTASKK0SBAYES_H
+#define ALIANALYSISTASKK0SBayes_H
+
+// ROOT includes
+#include <TObject.h>
+#include <TClonesArray.h>
+#include <TList.h>
+#include <TProfile.h>
+#include "AliPIDCombined.h"
+#include "TF1.h"
+
+// AliRoot includes
+#include <AliAnalysisTaskSE.h>
+#include <AliAODEvent.h>
+#include "AliPIDperfContainer.h"
+
+class TH1F;
+class TH2F;
+class AliESDtrackCuts;
+
+class AliAnalysisTaskK0sBayes : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskK0sBayes();
+  AliAnalysisTaskK0sBayes(const char *name);
+
+  virtual ~AliAnalysisTaskK0sBayes();
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *); 
+  
+  Double_t GetVtxCut() { return fVtxCut; }   
+  Double_t GetEtaCut() { return fEtaCut; }     
+  Double_t GetMinPt() { return fMinPt; }   
+
+  virtual void  SetVtxCut(Double_t vtxCut){fVtxCut = vtxCut;}
+  virtual void  SetEtaCut(Double_t etaCut){fEtaCut = etaCut;}
+  virtual void  SetMinPt(Double_t value) {fMinPt = value;}   
+
+  virtual void SetMC(Bool_t flag = kTRUE){fIsMC = flag;};
+  virtual void SetQA(Bool_t flag = kTRUE){fQAsw = flag;};
+
+  void SetTPCclusterN(Int_t ncl){fNcluster=ncl;};
+
+  static const Int_t nPtBin = 13;  //! # pt ks bins
+
+  void SetFilterBit(Int_t fb){fFilterBit=fb;};
+  Int_t GetFilterBit() const {return fFilterBit;};
+
+ private:
+  AliAnalysisTaskK0sBayes(const AliAnalysisTaskK0sBayes &old); 
+  AliAnalysisTaskK0sBayes& operator=(const AliAnalysisTaskK0sBayes &source); 
+
+  Int_t PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie);
+  Int_t FindDaugheterIndex(AliAODTrack *trk);
+
+  virtual Float_t GetVertex(AliAODEvent* aod) const;
+  virtual void Analyze(AliAODEvent* aodEvent); 
+  virtual void SelectK0s();
+  virtual Int_t IsChannelValid(Float_t etaAbs);
+
+  Double_t     fVtxCut;             // Vtx cut on z position in cm
+  Double_t     fEtaCut;             // Eta cut used to select particles
+  Double_t     fMinPt;              // Min pt - for histogram limits
+
+  Bool_t fIsMC; // if MC
+  Bool_t fQAsw;   // if QA
+
+  static const Int_t nCentrBin = 9; //! # cenrality bins
+
+  //
+  // Cuts and options
+  //
+
+  Int_t fNcluster;           // Numer of TPC cluster required
+  Int_t fFilterBit;          // filter bit to be used
+  TList *fList;              //! List for output objects
+  TList *fList2;             //! List for output objects
+  TList *fList3;             //! List for output objects
+
+  //
+  Float_t fCentrality;  //! current centrality for the tree
+  Float_t fPsi;         //! current Psi from TPC
+
+  Float_t fPtKs;   //! variable to fill the tree
+  Float_t fPhiKs;  //! variable to fill the tree
+  Float_t fEtaKs;  //! variable to fill the tree
+  Float_t fMassV0; //! variable to fill the tree
+  Float_t fPtKp;  //! variable to fill the tree
+  Float_t fPhiKp; //! variable to fill the tree
+  Float_t fEtaKp; //! variable to fill the tree
+  Float_t fPtKn;  //! variable to fill the tree
+  Float_t fPhiKn; //! variable to fill the tree
+  Float_t fEtaKn; //! variable to fill the tree
+  Int_t fPidKp;   //! variable to fill the tree
+  Int_t fPidKn;   //! variable to fill the tree
+
+  TH2F *hMatching[4]; //! matching (all, pi, k ,p)
+  TH2F *hTracking[4]; //! tracking (all, pi, k ,p)
+  TH2F *fTOFTPCsignal;   //! histo with tof signal
+  TH1F *fCombsignal;   //! histo with tof signal
+
+  static Float_t fPtKsMin[nPtBin];// ptmin bin
+  static Float_t fPtKsMax[nPtBin];// ptmax bin
+
+  // QA plots
+  TH2F *fPiTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fKaTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fPrTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fElTPC[nCentrBin];//! TPC dE/dx plot
+
+  TH2F *fPiTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fKaTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fPrTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fElTOF[nCentrBin];//! TPC dE/dx plot
+
+  AliESDtrackCuts *fCutsDaughter;
+
+  AliPIDCombined *fPIDCombined;  //! PID combined object
+  AliPIDperfContainer* fContPid;      //! results for positive
+  AliPIDperfContainer* fContPid2;     //! results for negative
+
+  Int_t fNK0s; //! number of K0s in my private selection
+  Float_t fPhiK0s[1000]; //! phi of K0s in my private selection
+  Float_t fPtK0s[1000];//! pt of K0s in my private selection
+  Int_t fNpiPos; //! number of positive pions for K0s selection
+  Int_t fNpiNeg; //! number of negative pions for K0s selection
+  Int_t fIPiPos[1000]; //! position in the AOD stack of positive pions candidates
+  Int_t fIPiNeg[1000]; //! position in the AOD stack of negative pions candidates
+  Int_t fIpiP[1000]; //! position in the AOD stack of positive pions for K0s
+  Int_t fIpiN[1000]; //! position in the AOD stack of negative pions for K0s
+  Float_t fMassKs[1000]; //! K0s mass
+
+  TH1F *fHmismTOF; //! TOF mismatch distribution
+  TH1D *fHchannelTOFdistr; //! TOF channel distance w.r.t. IP
+
+  ClassDef(AliAnalysisTaskK0sBayes, 1);    //Analysis task for bayesian (K0s)
+};
+
+#endif
diff --git a/PWGPP/pid/AliAnalysisTaskLambdaBayes.cxx b/PWGPP/pid/AliAnalysisTaskLambdaBayes.cxx
new file mode 100644 (file)
index 0000000..45124c3
--- /dev/null
@@ -0,0 +1,938 @@
+#include "AliAnalysisTaskLambdaBayes.h"
+
+// ROOT includes
+#include <TMath.h>
+
+// AliRoot includes
+#include "AliInputEventHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliCentrality.h"
+#include "AliVHeader.h"
+#include "AliAODVZERO.h"
+#include "TFile.h"
+#include "AliOADBContainer.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliMCEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "TChain.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDVertex.h"
+#include "AliEventplane.h"
+#include "AliAnalysisManager.h"
+#include "TRandom.h"
+
+
+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
+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
+
+// STL includes
+//#include <iostream>
+//using namespace std;
+
+ClassImp(AliAnalysisTaskLambdaBayes)
+//_____________________________________________________________________________
+AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes():
+  AliAnalysisTaskSE(),
+  fVtxCut(10.0),  // cut on |vertex| < fVtxCut
+  fEtaCut(0.8),   // cut on |eta| < fEtaCut
+  fMinPt(0.15),   // cut on pt > fMinPt
+  fIsMC(kFALSE),
+  fQAsw(kFALSE),
+  fNcluster(70),
+  fFilterBit(1),
+  fList(new TList()),
+  fList2(new TList()),
+  fList3(new TList()),
+  fCentrality(-1),
+  fPsi(0),
+  fPtLambdaC(0.),
+  fPhiLambdaC(0.),
+  fEtaLambda(0.),
+  fMassV0(0.),
+  fPtKp(0.),
+  fPhiKp(0.),
+  fEtaKp(0.),
+  fPtKn(0.),
+  fPhiKn(0.),
+  fEtaKn(0.),
+  fPidKp(0),
+  fPidKn(0),
+  fTOFTPCsignal(0),
+  fCombsignal(0),
+  fCutsDaughter(NULL),
+  fPIDCombined(NULL),
+  fContPid(NULL),
+  fContPid2(NULL),
+  fNLambda(0),
+  fNpPos(0),
+  fNpNeg(0),
+  fHmismTOF(0),
+  fHchannelTOFdistr(0)
+{
+  // Default constructor (should not be used)
+  fList->SetName("contLambdaBayes1");
+  fList2->SetName("contLambdaBayes2");
+  fList3->SetName("contLambdaBayes3");
+
+  fList->SetOwner(kTRUE); 
+  fList2->SetOwner(kTRUE); 
+  fList3->SetOwner(kTRUE); 
+
+  TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
+  fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
+
+  TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+  fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
+}
+
+//______________________________________________________________________________
+AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes(const char *name):
+  AliAnalysisTaskSE(name),
+  fVtxCut(10.0),  // cut on |vertex| < fVtxCut
+  fEtaCut(0.8),   // cut on |eta| < fEtaCut
+  fMinPt(0.15),   // cut on pt > fMinPt
+  fIsMC(kFALSE),
+  fQAsw(kFALSE),
+  fNcluster(70),
+  fFilterBit(1),
+  fList(new TList()),
+  fList2(new TList()),
+  fList3(new TList()),
+  fCentrality(-1),
+  fPsi(0),
+  fPtLambdaC(0.),
+  fPhiLambdaC(0.),
+  fEtaLambda(0.),
+  fMassV0(0.),
+  fPtKp(0.),
+  fPhiKp(0.),
+  fEtaKp(0.),
+  fPtKn(0.),
+  fPhiKn(0.),
+  fEtaKn(0.),
+  fPidKp(0),
+  fPidKn(0),
+  fTOFTPCsignal(0),
+  fCombsignal(0),
+  fCutsDaughter(NULL),
+  fPIDCombined(NULL),
+  fContPid(NULL),
+  fContPid2(NULL),
+  fNLambda(0),
+  fNpPos(0),
+  fNpNeg(0),
+  fHmismTOF(0),
+  fHchannelTOFdistr(0)
+{
+
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
+  DefineOutput(3, TList::Class());
+  DefineOutput(4, TList::Class());
+
+  // Output slot #1 writes into a TTree
+  fList->SetName("contLambdaBayes1");
+  fList2->SetName("contLambdaBayes2");
+  fList3->SetName("contLambdaBayes3");
+
+  fList->SetOwner(kTRUE); 
+  fList2->SetOwner(kTRUE); 
+  fList3->SetOwner(kTRUE); 
+
+  TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
+  fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
+
+  TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+  fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
+}
+//_____________________________________________________________________________
+AliAnalysisTaskLambdaBayes::~AliAnalysisTaskLambdaBayes()
+{
+
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskLambdaBayes::UserCreateOutputObjects()
+{
+
+  fPIDCombined=new AliPIDCombined;
+  fPIDCombined->SetDefaultTPCPriors();
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+
+  Float_t invMmin = 1.115-0.005*8;
+  Float_t invMmax = 1.115+0.005*8;
+
+  const Int_t nBinPid = 14;
+
+  Int_t binPid[nBinPid] = {1/*ptLambda*/,8/*EtaPP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+
+  Int_t binPid2[nBinPid] = {1/*ptLambda*/,8/*EtaPN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+
+  fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
+  fContPid->SetTitleX("M_{#Lambda}");
+  fContPid->SetTitleY("centrality (%)");
+  fContPid->SetVarName(0,"p_{T}^{#Lambda}");
+  fContPid->SetVarRange(0,0,10);
+  fContPid->SetVarName(1,"#eta^{#Lambda}");
+  fContPid->SetVarRange(1,-0.8,0.8);
+  fContPid->SetVarName(2,"p_{T}^{Kp}");
+  fContPid->SetVarRange(2,0.3,4.3);
+  fContPid->SetVarName(3,"p_{T}^{Kn}");
+  fContPid->SetVarRange(3,0.3,4.3);
+  fContPid->SetVarName(4,"BayesProb^{Kp}");
+  fContPid->SetVarRange(4,0,1.);
+  fContPid->SetVarName(5,"BayesProb^{Kn}");
+  fContPid->SetVarRange(5,0,1.);
+  fContPid->SetVarName(6,"isTOFmatch^{Kp}");
+  fContPid->SetVarRange(6,-0.5,1.5);
+  fContPid->SetVarName(7,"isTOFmatch^{Kn}");
+  fContPid->SetVarRange(7,-0.5,1.5);
+  fContPid->SetVarName(8,"isLambdaTrue^{Kn}");
+  fContPid->SetVarRange(8,-0.5,1.5);
+  fContPid->SetVarName(9,"N#sigma^{Kp}");
+  fContPid->SetVarRange(9,1.25,6.25);
+  fContPid->SetVarName(10,"N#sigma^{Kn}");
+  fContPid->SetVarRange(10,1.25,6.25);
+  fContPid->SetVarName(11,"#Delta#phi^{Kp}");
+  fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
+  fContPid->SetVarName(12,"#Delta#phi^{Kn}");
+  fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
+  fContPid->SetVarName(13,"#Psi");
+  fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
+
+  fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
+  fContPid2->SetTitleX("M_{#Lambda}");
+  fContPid2->SetTitleY("centrality (%)");
+  fContPid2->SetVarName(0,"p_{T}^{#Lambda}");
+  fContPid2->SetVarRange(0,0,10);
+  fContPid2->SetVarName(1,"#eta^{#Lambda}");
+  fContPid2->SetVarRange(1,-0.8,0.8);
+  fContPid2->SetVarName(2,"p_{T}^{Kp}");
+  fContPid2->SetVarRange(2,0.3,4.3);
+  fContPid2->SetVarName(3,"p_{T}^{Kn}");
+  fContPid2->SetVarRange(3,0.3,4.3);
+  fContPid2->SetVarName(4,"BayesProb^{Kp}");
+  fContPid2->SetVarRange(4,0,1.);
+  fContPid2->SetVarName(5,"BayesProb^{Kn}");
+  fContPid2->SetVarRange(5,0,1.);
+  fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
+  fContPid2->SetVarRange(6,-0.5,1.5);
+  fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
+  fContPid2->SetVarRange(7,-0.5,1.5);
+  fContPid2->SetVarName(8,"isLambdaTrue^{Kn}");
+  fContPid2->SetVarRange(8,-0.5,1.5);
+  fContPid2->SetVarName(9,"N#sigma^{Kp}");
+  fContPid2->SetVarRange(9,1.25,6.25);
+  fContPid2->SetVarName(10,"N#sigma^{Kn}");
+  fContPid2->SetVarRange(10,1.25,6.25);
+  fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
+  fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
+  fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
+  fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
+  fContPid2->SetVarName(13,"#Psi");
+  fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
+
+
+
+  const Int_t nDETsignal = 100; // mass
+  Double_t binDETsignal[nDETsignal+1];
+  for(Int_t i=0;i<nDETsignal+1;i++){
+    binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
+  }
+  const Int_t nDETsignal2 = 10; // centrality
+  Double_t binDETsignal2[nDETsignal2+1];
+  for(Int_t i=0;i<nDETsignal2+1;i++){
+    binDETsignal2[i] = i*100./ nDETsignal2;
+  }
+  fContPid->AddSpecies("Lambda",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
+  fContPid2->AddSpecies("Lambda2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
+
+  fList->Add(fContPid);
+  fList->Add(fContPid2);
+
+  hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+
+  hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+
+  fList2->Add(hMatching[0]);
+  fList2->Add(hMatching[1]);
+  fList2->Add(hMatching[2]);
+  fList2->Add(hMatching[3]);
+  fList2->Add(hTracking[0]);
+  fList2->Add(hTracking[1]);
+  fList2->Add(hTracking[2]);
+  fList2->Add(hTracking[3]);
+
+  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);
+  fList2->Add(fTOFTPCsignal);
+  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);
+  fList2->Add(fCombsignal);
+
+  // QA plots
+  char name[100],title[400];
+  for(Int_t i=0;i < nCentrBin;i++){
+    snprintf(name,100,"hPiTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPiTPC[i]);
+
+    snprintf(name,100,"hKaTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fKaTPC[i]);
+
+    snprintf(name,100,"hPrTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPrTPC[i]);
+
+    snprintf(name,100,"hElTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fElTPC[i]);
+
+    snprintf(name,100,"hPiTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPiTOF[i]);
+
+    snprintf(name,100,"hKaTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fKaTOF[i]);
+
+    snprintf(name,100,"hPrTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPrTOF[i]);
+
+    snprintf(name,100,"hElTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fElTOF[i]);
+
+  }
+
+  // Post output data.
+  PostData(1, fList);
+  PostData(2, fList2);
+  PostData(3, fList3);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskLambdaBayes::UserExec(Option_t *) 
+{
+    // Main loop
+    // Called for each event
+    
+    fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    if(!fOutputAOD){
+       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+       this->Dump();
+       return;
+    }
+    
+    Float_t zvtx = GetVertex(fOutputAOD);
+
+
+
+    //Get the MC object
+    if(fIsMC){
+      AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+      if (!mcHeader) {
+       AliError("Could not find MC Header in AOD");
+       return;
+      }
+    }
+
+    if (TMath::Abs(zvtx) < fVtxCut) {
+
+      SelectLambda();
+
+      //Centrality
+      Float_t v0Centr  = -10.;
+      Float_t trkCentr  = -10.;
+      AliCentrality *centrality = fOutputAOD->GetCentrality();
+      if (centrality){
+       trkCentr  = centrality->GetCentralityPercentile("V0M");
+       v0Centr = centrality->GetCentralityPercentile("TRK"); 
+      }
+
+      if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr>0){ // consistency cut on centrality selection
+        fCentrality = v0Centr;
+       Analyze(fOutputAOD); // Do analysis!!!!
+
+      }
+    }
+    
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskLambdaBayes::Analyze(AliAODEvent* aodEvent)
+{
+
+  Int_t ntrack = aodEvent->GetNumberOfTracks();
+
+  fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
+  fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
+  fPidKp=0,fPidKn=0;
+  fMassV0=-1;
+  
+  TLorentzVector ksV;
+  
+  Double_t px,py,pz,E;
+
+  Float_t invMmin = 1.115-0.005*8;
+  Float_t invMmax = 1.115+0.005*8;
+  
+  Int_t icentr = 8;
+  if(fCentrality < 0) icentr = 8;
+  else if(fCentrality < 10) icentr = 0;
+  else if(fCentrality < 20) icentr = 1;
+  else if(fCentrality < 30) icentr = 2;
+  else if(fCentrality < 40) icentr = 3;
+  else if(fCentrality < 50) icentr = 4;
+  else if(fCentrality < 60) icentr = 5;
+  else if(fCentrality < 70) icentr = 6;
+  else if(fCentrality < 80) icentr = 7;
+
+  Float_t addMismatchForMC = 0.005;
+  if(fCentrality < 50) addMismatchForMC += 0.01;
+  if(fCentrality < 20) addMismatchForMC += 0.02;
+
+  fPsi = 0;
+  /* Compute TPC EP */
+  Double_t Qx2 = 0, Qy2 = 0;
+  Double_t Qx3 = 0, Qy3 = 0;
+  for(Int_t iT = 0; iT < ntrack; iT++) {
+    AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
+    
+    if (!aodTrack){
+      continue;
+    }
+    
+    Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
+    
+    if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag) 
+      continue;
+    
+    Double_t b[2] = {-99., -99.};
+    Double_t bCov[3] = {-99., -99., -99.};
+    if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
+      continue;
+    
+    if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
+      continue;
+    
+    Qx2 += TMath::Cos(2*aodTrack->Phi()); 
+    Qy2 += TMath::Sin(2*aodTrack->Phi());
+    Qx3 += TMath::Cos(3*aodTrack->Phi()); 
+    Qy3 += TMath::Sin(3*aodTrack->Phi());
+    
+  }
+
+  fPsi = TMath::ATan2(Qy2, Qx2)/2.;
+
+  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+  AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
+
+  PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
+  PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
+  PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
+  PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
+
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+
+  Double_t probP[10],probN[10];
+  Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
+
+  
+  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+  TClonesArray *mcArray = NULL;
+  if (mcHeader)
+    mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+
+  Int_t nmc = 0;
+  if(mcArray)
+    nmc = mcArray->GetEntries();
+
+  for(Int_t i=0;i < ntrack;i++){ // loop on tracks
+    AliAODTrack* track = aodEvent->GetTrack(i);
+    
+    AliAODMCParticle *mcp = NULL;
+    Int_t pdg = 0;
+    
+    if (!track){
+      continue;
+    }
+    
+    Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
+    
+    Int_t label = -1;
+    if(mcArray){
+      label = track->GetLabel();
+      if(label != -1 && label < nmc){
+       label = TMath::Abs(label);
+       mcp = (AliAODMCParticle*)mcArray->At(label);
+       pdg = TMath::Abs(mcp->GetPdgCode());
+      }
+      else
+       label = -1;
+    }
+    else{
+      /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
+    }
+    
+    if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
+      hTracking[0]->Fill(track->Pt(),fCentrality);
+      if(pdg == 211)
+       hTracking[1]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 321)
+       hTracking[2]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 2212)
+       hTracking[3]->Fill(track->Pt(),fCentrality);
+      else if(! mcp){ // Fill matching histo with the prob
+       hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
+       hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
+       hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
+      }
+    }
+    
+    if(!tofMatch) continue;
+    
+    if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
+      hMatching[0]->Fill(track->Pt(),fCentrality);
+      if(pdg == 211)
+       hMatching[1]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 321)
+       hMatching[2]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 2212)
+       hMatching[3]->Fill(track->Pt(),fCentrality);
+      else if(! mcp){ // Fill matching histo with the prob
+       hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
+       hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
+       hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
+      }
+    }
+  }
+  
+  Int_t pdg1 = -1;
+  Int_t pdg2 = -1;
+
+
+  // start analysis Lambda
+  for(Int_t i=0;i < ntrack;i++){ // loop on proton candidate tracks
+    AliAODTrack* KpTrack = aodEvent->GetTrack(i);
+        
+    if (!KpTrack){
+      continue;
+    }
+    
+    if(!(KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
+
+    Bool_t isLambda = (KpTrack->Charge() > 0);
+
+    nSigmaComb=5;
+    fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
+    fPidKp=0;
+
+    UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
+
+    Double_t oldpP[10];
+    fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
+
+    
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
+    if(isLambda)fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
+    if(isLambda) fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
+    if(isLambda)fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
+    if(isLambda) fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+
+    if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
+
+    Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
+
+    if(mcArray){
+      Int_t labelK = TMath::Abs(KpTrack->GetLabel());
+      AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
+      pdg1 = TMath::Abs(mcp1->GetPdgCode());
+    }
+
+    fPidKp = Int_t(probP[4]*100);
+
+    if(tofMatch1){
+      if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
+       // remove this tof hit
+       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+       detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
+       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+       fPidKp = Int_t(probP[4]*100);
+       tofMatch1=0;
+      }
+      else{
+       if(probP[4] > probP[3] && probP[4] > probP[2] && probP[4] > probP[0]) fPidKp += 128; // max prob
+       
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
+       if(isLambda) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
+       if(isLambda) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
+       if(isLambda) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
+       if(isLambda) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       
+       nSigmaTOF = TMath::Abs(nSigmaTOF);
+       
+       if(fIsMC){
+         Float_t mismAdd = addMismatchForMC;
+         if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
+         
+         if(gRandom->Rndm() < mismAdd){
+           Float_t etaAbs = TMath::Abs(KpTrack->Eta());
+           Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
+           channel = channel % 8736;
+           Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
+           
+           // generate random time
+           Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
+           Double_t times[10];
+           KpTrack->GetIntegratedTimes(times);
+           nSigmaTOF = TMath::Abs(timeRandom - times[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kProton);
+         }
+       }
+
+       if(fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
+
+       if(nSigmaTOF < 2) fPidKp += 256;
+       else if(nSigmaTOF < 3) fPidKp += 512;
+      }
+    }
+    
+    if(tofMatch1){
+      nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
+      if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2){
+       fCombsignal->Fill(nSigmaComb);
+      }
+    } else {
+      nSigmaComb = TMath::Abs(nSigmaTPC);
+    }
+
+    // use sigmaTOF instead of sigmaComb
+    //nSigmaComb = nSigmaTOF;
+
+    if(nSigmaComb < 2) nSigmaComb = 2;
+    else if(nSigmaComb < 3) nSigmaComb = 3;
+    else if(nSigmaComb < 5) nSigmaComb = 4.99;
+    else nSigmaComb = 6;
+
+    Int_t iks=-1;
+    for(Int_t k=0;k < fNLambda;k++){ // find the Lambda which contains the positive track
+      if(i == fIpP[k]) iks = k;
+    }
+
+    if(fPtKp > 4.299) fPtKp = 4.299;
+
+    if(iks > -1 && fIpN[iks] > -1){
+      //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
+      Int_t j = fIpN[iks];
+      AliAODTrack* KnTrack = aodEvent->GetTrack(j);
+      
+      if (!KnTrack){
+       continue;
+      }
+
+      if(!(KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
+
+      fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
+      fPidKn=0;
+
+      UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
+      Double_t oldpN[10];
+      fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
+
+      nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
+      
+      if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
+      
+      nSigmaComb2=5;
+
+      Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
+
+      if(mcArray){
+       Int_t labelK = TMath::Abs(KnTrack->GetLabel());
+       AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
+       pdg2 = TMath::Abs(mcp2->GetPdgCode());
+     }
+
+      fPidKn = Int_t(probN[2]*100);
+
+      if(tofMatch2){
+       if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
+
+       nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
+
+       nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
+
+       if(nSigmaTOF2 < 2) fPidKn += 256;
+       else if(nSigmaTOF2 < 3) fPidKn += 512;
+      }
+
+      px = KpTrack->Px() + KnTrack->Px();
+      py = KpTrack->Py() + KnTrack->Py();
+      pz = KpTrack->Pz() + KnTrack->Pz();
+      E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 0.938e-01*0.938e-01);
+      E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
+
+      ksV.SetPxPyPzE(px,py,pz,E);
+      fMassV0 = fMassLambda[iks];
+      
+      if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
+
+
+      fPtLambdaC = ksV.Pt();
+      fEtaLambda = ksV.Eta();
+      fPhiLambdaC = ksV.Phi();
+
+      if(tofMatch2){
+       nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
+      } else {
+       nSigmaComb2 = TMath::Abs(nSigmaTPC2);
+      }
+
+      // use sigmaTOF instead of sigmaComb
+      //nSigmaComb2 = nSigmaTOF2;
+
+      if(nSigmaComb2 < 2) nSigmaComb2 = 2;
+      else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
+      else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
+      else nSigmaComb2 = 6;  
+
+      Bool_t isTrue = kFALSE;
+
+      if(mcArray){
+       Int_t labelP = TMath::Abs(KpTrack->GetLabel());
+       Int_t labelN = TMath::Abs(KnTrack->GetLabel());
+
+       if(labelP > -1 && labelN > -1){
+         AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
+         AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
+
+         Int_t mP = partP->GetMother();
+         Int_t mN = partN->GetMother();
+         if(mP == mN && mP > -1){
+           AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
+           Int_t pdgM = partM->GetPdgCode();
+           if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
+         }
+       }
+
+      }
+
+      Double_t deltaphi1 = KpTrack->Phi() - fPsi;
+      Double_t deltaphi2 = KnTrack->Phi() - fPsi;
+
+      while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
+      while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
+      while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
+      while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
+
+      if(fPtKn > 4.299) fPtKn = 4.299;
+
+      Float_t xTOfill[] = {fPtLambdaC,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
+      Float_t xTOfill2[] = {fPtLambdaC,KpTrack->Eta(),fPtKn,fPtKp,(fPidKn%128)*0.01,(fPidKp%128)*0.01,tofMatch2,tofMatch1,isTrue,nSigmaComb2,nSigmaComb,deltaphi2,deltaphi1,fPsi};
+      
+
+      Int_t ipt = 0;
+      while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
+       ipt++;
+      }
+      ipt--;
+      if(ipt < 0) ipt = 0; // just to be sure
+
+      if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
+       if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
+       else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
+      }
+
+
+    }
+  } // end analysi Lambda
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
+{
+
+  Float_t zvtx = -999;
+
+  const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
+  if (!vtxAOD)
+    return zvtx;
+  if(vtxAOD->GetNContributors()>0)
+    zvtx = vtxAOD->GetZ();
+  
+  return zvtx;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
+{ 
+  // Terminate loop
+  Printf("Terminate()");
+}
+//=======================================================================
+void AliAnalysisTaskLambdaBayes::SelectLambda(){
+  fNLambda=0;
+  fNpPos=0;
+  fNpNeg=0;
+
+  Int_t nV0s = fOutputAOD->GetNumberOfV0s();
+  AliAODv0 *myV0;
+  Double_t dMASS=0.0;
+  for (Int_t i=0; i!=nV0s; ++i) {
+    myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
+    if(!myV0) continue;
+    if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
+    Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
+    if(pass) {
+      if(pass==1) dMASS = myV0->MassLambda();
+      else dMASS = myV0->MassAntiLambda();
+
+      if(TMath::Abs(dMASS-1.115)/0.005 < 8){
+       AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
+       AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
+       if(iT->Charge()<0){
+         iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
+         jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
+       }
+       fPhiLambda[fNLambda] = myV0->Phi();
+       fPtLambda[fNLambda] = myV0->Pt();
+       fIpP[fNLambda] = FindDaugheterIndex(iT);
+       fIpN[fNLambda] = FindDaugheterIndex(jT);
+
+       if(pass==2){
+         Int_t itemp = fIpP[fNLambda];
+         fIpP[fNLambda] = fIpN[fNLambda];
+         fIpN[fNLambda] = itemp;
+       }
+
+       fMassLambda[fNLambda] = dMASS;
+       if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
+         fNLambda++;
+       }
+
+      }
+    }
+  }
+}
+
+//=======================================================================
+Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
+{ 
+  if (myV0->GetOnFlyStatus() ) return 0;
+  
+  //the following is needed in order to evualuate track-quality
+  AliAODTrack *iT, *jT;
+  AliAODVertex *vV0s = myV0->GetSecondaryVtx();
+  Double_t pos[3],cov[6];
+  vV0s->GetXYZ(pos);
+  vV0s->GetCovarianceMatrix(cov);
+  const AliESDVertex vESD(pos,cov,100.,100);
+  
+  // TESTING CHARGE
+  int iPos, iNeg;
+  iT=(AliAODTrack*) myV0->GetDaughter(0);
+  if(iT->Charge()>0) {
+    iPos = 0; iNeg = 1;
+  } else {
+    iPos = 1; iNeg = 0;
+  }
+  // END OF TEST
+
+  iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
+
+  jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
+
+  Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
+  if(!trkFlag) return 0;
+  Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
+  if(!trkFlag2) return 0;
+
+  Double_t pvertex[3];
+  pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
+  pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
+  pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
+
+  Double_t dDL=myV0->DecayLengthV0( pvertex );
+  if(dDL  < 0.5 || dDL > 25) return 0;
+
+  Double_t dDCA=myV0->DcaV0Daughters();
+  if(dDCA >0.5) return 0;
+
+  Double_t dCTP=myV0->CosPointingAngle( pvertex );
+  if(dCTP < -1) return 0;
+
+//   AliESDtrack ieT( iT );
+//   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
+//   if(TMath::Abs(dD0P) < 0]) return 0;
+
+//   AliESDtrack jeT( jT );
+//   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
+//   if(TMath::Abs(dD0M) < 0) return 0;
+
+//   Double_t dD0D0=dD0P*dD0M;
+//   if(dD0D0>0) return 0;
+
+//   Double_t dETA=myV0->Eta();
+//   if(dETA <-0.8) return 0;
+//   if(dETA >0.8) return 0;
+
+//   Double_t dQT=myV0->PtArmV0();
+//   if(specie==0) if(dQT<???) return 0;
+
+  Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
+  if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
+
+  if(specie==1 && dALPHA<0) return 2; // antilambda
+  return 1; // Ks or lambda
+}
+//-------------------------------------------------------------------------
+Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
+  Int_t ntrack = fOutputAOD->GetNumberOfTracks();
+
+  for(Int_t i=0;i < ntrack;i++){ // loop on tracks
+    AliAODTrack* track = fOutputAOD->GetTrack(i);
+    if(track == trk) return i;
+  }
+  
+  printf("Daughter for %p not found\n",trk);
+  return -1;
+}
+//-------------------------------------------------------------------------
+Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
+  if(!fIsMC) return 1; // used only on MC
+
+  Int_t run = fOutputAOD->GetRunNumber();
+  if( (run>=136851&&run<=139846) || (run>=165772 && run<=170718) ){ // LHC10h or LHC11h because of TOF matching window at 3 cm
+    Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs); 
+  
+    if(!(channel%20)) return 0; // 5% additional loss in MC
+  }
+
+  return 1;
+}
diff --git a/PWGPP/pid/AliAnalysisTaskLambdaBayes.h b/PWGPP/pid/AliAnalysisTaskLambdaBayes.h
new file mode 100644 (file)
index 0000000..5b17a96
--- /dev/null
@@ -0,0 +1,141 @@
+#ifndef ALIANALYSISTASKLAMBDABAYES_H
+#define ALIANALYSISTASKLAMBDABayes_H
+
+// ROOT includes
+#include <TObject.h>
+#include <TClonesArray.h>
+#include <TList.h>
+#include <TProfile.h>
+#include "AliPIDCombined.h"
+#include "TF1.h"
+
+// AliRoot includes
+#include <AliAnalysisTaskSE.h>
+#include <AliAODEvent.h>
+#include "AliPIDperfContainer.h"
+
+class TH1F;
+class TH2F;
+class AliESDtrackCuts;
+
+class AliAnalysisTaskLambdaBayes : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskLambdaBayes();
+  AliAnalysisTaskLambdaBayes(const char *name);
+
+  virtual ~AliAnalysisTaskLambdaBayes();
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *); 
+  
+  Double_t GetVtxCut() { return fVtxCut; }   
+  Double_t GetEtaCut() { return fEtaCut; }     
+  Double_t GetMinPt() { return fMinPt; }   
+
+  virtual void  SetVtxCut(Double_t vtxCut){fVtxCut = vtxCut;}
+  virtual void  SetEtaCut(Double_t etaCut){fEtaCut = etaCut;}
+  virtual void  SetMinPt(Double_t value) {fMinPt = value;}   
+
+  virtual void SetMC(Bool_t flag = kTRUE){fIsMC = flag;};
+  virtual void SetQA(Bool_t flag = kTRUE){fQAsw = flag;};
+
+  void SetTPCclusterN(Int_t ncl){fNcluster=ncl;};
+
+  static const Int_t nPtBin = 13;  //! # pt ks bins
+
+  void SetFilterBit(Int_t fb){fFilterBit=fb;};
+  Int_t GetFilterBit() const {return fFilterBit;};
+
+ private:
+  AliAnalysisTaskLambdaBayes(const AliAnalysisTaskLambdaBayes &old); 
+  AliAnalysisTaskLambdaBayes& operator=(const AliAnalysisTaskLambdaBayes &source); 
+
+  Int_t PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie);
+  Int_t FindDaugheterIndex(AliAODTrack *trk);
+
+  virtual Float_t GetVertex(AliAODEvent* aod) const;
+  virtual void Analyze(AliAODEvent* aodEvent); 
+  virtual Int_t IsChannelValid(Float_t etaAbs);
+  virtual void SelectLambda();
+
+  Double_t     fVtxCut;             // Vtx cut on z position in cm
+  Double_t     fEtaCut;             // Eta cut used to select particles
+  Double_t     fMinPt;              // Min pt - for histogram limits
+
+  Bool_t fIsMC; // if MC
+  Bool_t fQAsw;   // if QA
+
+  static const Int_t nCentrBin = 9; //! # cenrality bins
+
+  //
+  // Cuts and options
+  //
+
+  Int_t fNcluster;           // Numer of TPC cluster required
+  Int_t fFilterBit;          // filter bit to be used
+  TList *fList;              //! List for output objects
+  TList *fList2;             //! List for output objects
+  TList *fList3;             //! List for output objects
+
+  //
+  Float_t fCentrality;  //! current centrality for the tree
+  Float_t fPsi;         //! current Psi from TPC
+
+  Float_t fPtLambdaC;   //! variable to fill the tree
+  Float_t fPhiLambdaC;  //! variable to fill the tree
+  Float_t fEtaLambda;  //! variable to fill the tree
+  Float_t fMassV0; //! variable to fill the tree
+  Float_t fPtKp;  //! variable to fill the tree
+  Float_t fPhiKp; //! variable to fill the tree
+  Float_t fEtaKp; //! variable to fill the tree
+  Float_t fPtKn;  //! variable to fill the tree
+  Float_t fPhiKn; //! variable to fill the tree
+  Float_t fEtaKn; //! variable to fill the tree
+  Int_t fPidKp;   //! variable to fill the tree
+  Int_t fPidKn;   //! variable to fill the tree
+
+  TH2F *hMatching[4]; //! matching (all, pi, k ,p)
+  TH2F *hTracking[4]; //! tracking (all, pi, k ,p)
+  TH2F *fTOFTPCsignal;   //! histo with tof signal
+  TH1F *fCombsignal;   //! histo with tof signal
+
+  static Float_t fPtLambdaMin[nPtBin];// ptmin bin
+  static Float_t fPtLambdaMax[nPtBin];// ptmax bin
+
+  // QA plots
+  TH2F *fPiTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fKaTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fPrTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fElTPC[nCentrBin];//! TPC dE/dx plot
+
+  TH2F *fPiTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fKaTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fPrTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fElTOF[nCentrBin];//! TPC dE/dx plot
+
+  AliESDtrackCuts *fCutsDaughter;
+
+  AliPIDCombined *fPIDCombined;  //! PID combined object
+  AliPIDperfContainer* fContPid;      //! results for positive
+  AliPIDperfContainer* fContPid2;     //! results for negative
+
+  Int_t fNLambda; //! number of Lambda in my private selection
+  Float_t fPhiLambda[1000]; //! phi of Lambda in my private selection
+  Float_t fPtLambda[1000];//! pt of Lambda in my private selection
+  Int_t fNpPos; //! number of positive pions for Lambda selection
+  Int_t fNpNeg; //! number of negative pions for Lambda selection
+  Int_t fIPPos[1000]; //! position in the AOD stack of positive pions candidates
+  Int_t fIPNeg[1000]; //! position in the AOD stack of negative pions candidates
+  Int_t fIpP[1000]; //! position in the AOD stack of positive pions for Lambda
+  Int_t fIpN[1000]; //! position in the AOD stack of negative pions for Lambda
+  Float_t fMassLambda[1000]; //! Lambda mass
+  
+  TH1F *fHmismTOF; //! TOF mismatch distribution
+  TH1D *fHchannelTOFdistr; //! TOF channel distance w.r.t. IP
+
+  ClassDef(AliAnalysisTaskLambdaBayes, 1);    //Analysis task for Bayesian (Lambda)
+};
+
+#endif
diff --git a/PWGPP/pid/AliAnalysisTaskPhiBayes.cxx b/PWGPP/pid/AliAnalysisTaskPhiBayes.cxx
new file mode 100644 (file)
index 0000000..c18c12c
--- /dev/null
@@ -0,0 +1,825 @@
+#include "AliAnalysisTaskPhiBayes.h"
+
+// ROOT includes
+#include <TMath.h>
+
+// AliRoot includes
+#include "AliInputEventHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliCentrality.h"
+#include "AliVHeader.h"
+#include "AliAODVZERO.h"
+#include "TFile.h"
+#include "AliOADBContainer.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliMCEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "TChain.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDVertex.h"
+#include "AliEventplane.h"
+#include "AliAnalysisManager.h"
+#include "TRandom.h"
+
+
+Float_t AliAnalysisTaskPhiBayes::fPtPhiMin[AliAnalysisTaskPhiBayes::nPtBin] = {0.,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.};// ptmin bin
+Float_t AliAnalysisTaskPhiBayes::fPtPhiMax[AliAnalysisTaskPhiBayes::nPtBin] = {0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.,10.};// ptmax bin
+
+// STL includes
+//#include <iostream>
+//using namespace std;
+
+ClassImp(AliAnalysisTaskPhiBayes)
+//_____________________________________________________________________________
+AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes():
+  AliAnalysisTaskSE(),
+  fVtxCut(10.0),  // cut on |vertex| < fVtxCut
+  fEtaCut(0.8),   // cut on |eta| < fEtaCut
+  fMinPt(0.15),   // cut on pt > fMinPt
+  fIsMC(kFALSE),
+  fQAsw(kFALSE),
+  fNcluster(70),
+  fFilterBit(1),
+  fList(new TList()),
+  fList2(new TList()),
+  fList3(new TList()),
+  fCentrality(-1),
+  fPsi(0),
+  fPtPhi(0.),
+  fPhiPhi(0.),
+  fEtaPhi(0.),
+  fMassV0(0.),
+  fPtKp(0.),
+  fPhiKp(0.),
+  fEtaKp(0.),
+  fPtKn(0.),
+  fPhiKn(0.),
+  fEtaKn(0.),
+  fPidKp(0),
+  fPidKn(0),
+  fTOFTPCsignal(0),
+  fCombsignal(0),
+  fCutsDaughter(NULL),
+  fPIDCombined(NULL),
+  fContPid(NULL),
+  fContPid2(NULL),
+  fHmismTOF(0),
+  fHchannelTOFdistr(0)
+{
+  // Default constructor (should not be used)
+  fList->SetName("contPhiBayes1");
+  fList2->SetName("contPhiBayes2");
+  fList3->SetName("contPhiBayes3");
+
+  fList->SetOwner(kTRUE); 
+  fList2->SetOwner(kTRUE); 
+  fList3->SetOwner(kTRUE); 
+
+  TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
+  fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
+
+  TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+  fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
+}
+
+//______________________________________________________________________________
+AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes(const char *name):
+  AliAnalysisTaskSE(name),
+  fVtxCut(10.0),  // cut on |vertex| < fVtxCut
+  fEtaCut(0.8),   // cut on |eta| < fEtaCut
+  fMinPt(0.15),   // cut on pt > fMinPt
+  fIsMC(kFALSE),
+  fQAsw(kFALSE),
+  fNcluster(70),
+  fFilterBit(1),
+  fList(new TList()),
+  fList2(new TList()),
+  fList3(new TList()),
+  fCentrality(-1),
+  fPsi(0),
+  fPtPhi(0.),
+  fPhiPhi(0.),
+  fEtaPhi(0.),
+  fMassV0(0.),
+  fPtKp(0.),
+  fPhiKp(0.),
+  fEtaKp(0.),
+  fPtKn(0.),
+  fPhiKn(0.),
+  fEtaKn(0.),
+  fPidKp(0),
+  fPidKn(0),
+  fTOFTPCsignal(0),
+  fCombsignal(0),
+  fCutsDaughter(NULL),
+  fPIDCombined(NULL),
+  fContPid(NULL),
+  fContPid2(NULL),
+  fHmismTOF(0),
+  fHchannelTOFdistr(0)
+{
+
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
+  DefineOutput(3, TList::Class());
+  DefineOutput(4, TList::Class());
+
+  // Output slot #1 writes into a TTree
+  fList->SetName("contPhiBayes1");
+  fList2->SetName("contPhiBayes2");
+  fList3->SetName("contPhiBayes3");
+
+  fList->SetOwner(kTRUE); 
+  fList2->SetOwner(kTRUE); 
+  fList3->SetOwner(kTRUE); 
+
+  TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
+  fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
+
+  TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+  fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist"); 
+}
+//_____________________________________________________________________________
+AliAnalysisTaskPhiBayes::~AliAnalysisTaskPhiBayes()
+{
+
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskPhiBayes::UserCreateOutputObjects()
+{
+
+  fPIDCombined=new AliPIDCombined;
+  fPIDCombined->SetDefaultTPCPriors();
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+
+  Float_t invMmin = 0.985;
+  Float_t invMmax = 1.045;
+
+  const Int_t nBinPid = 14;
+
+  Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaKP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+
+  Int_t binPid2[nBinPid] = {1/*ptPhi*/,8/*EtaKN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+
+  fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
+  fContPid->SetTitleX("M_{#phi}");
+  fContPid->SetTitleY("centrality (%)");
+  fContPid->SetVarName(0,"p_{T}^#phi}");
+  fContPid->SetVarRange(0,0,10);
+  fContPid->SetVarName(1,"#eta^{#phi}");
+  fContPid->SetVarRange(1,-0.8,0.8);
+  fContPid->SetVarName(2,"p_{T}^{Kp}");
+  fContPid->SetVarRange(2,0.3,4.3);
+  fContPid->SetVarName(3,"p_{T}^{Kn}");
+  fContPid->SetVarRange(3,0.3,4.3);
+  fContPid->SetVarName(4,"BayesProb^{Kp}");
+  fContPid->SetVarRange(4,0,1.);
+  fContPid->SetVarName(5,"BayesProb^{Kn}");
+  fContPid->SetVarRange(5,0,1.);
+  fContPid->SetVarName(6,"isTOFmatch^{Kp}");
+  fContPid->SetVarRange(6,-0.5,1.5);
+  fContPid->SetVarName(7,"isTOFmatch^{Kn}");
+  fContPid->SetVarRange(7,-0.5,1.5);
+  fContPid->SetVarName(8,"isPhiTrue^{Kn}");
+  fContPid->SetVarRange(8,-0.5,1.5);
+  fContPid->SetVarName(9,"N#sigma^{Kp}");
+  fContPid->SetVarRange(9,1.25,6.25);
+  fContPid->SetVarName(10,"N#sigma^{Kn}");
+  fContPid->SetVarRange(10,1.25,6.25);
+  fContPid->SetVarName(11,"#Delta#phi^{Kp}");
+  fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
+  fContPid->SetVarName(12,"#Delta#phi^{Kn}");
+  fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
+  fContPid->SetVarName(13,"#Psi");
+  fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
+
+  fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
+  fContPid2->SetTitleX("M_{#phi}");
+  fContPid2->SetTitleY("centrality (%)");
+  fContPid2->SetVarName(0,"p_{T}^{#phi}");
+  fContPid2->SetVarRange(0,0,10);
+  fContPid2->SetVarName(1,"#eta^{#phi}");
+  fContPid2->SetVarRange(1,-0.8,0.8);
+  fContPid2->SetVarName(2,"p_{T}^{Kp}");
+  fContPid2->SetVarRange(2,0.3,4.3);
+  fContPid2->SetVarName(3,"p_{T}^{Kn}");
+  fContPid2->SetVarRange(3,0.3,4.3);
+  fContPid2->SetVarName(4,"BayesProb^{Kp}");
+  fContPid2->SetVarRange(4,0,1.);
+  fContPid2->SetVarName(5,"BayesProb^{Kn}");
+  fContPid2->SetVarRange(5,0,1.);
+  fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
+  fContPid2->SetVarRange(6,-0.5,1.5);
+  fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
+  fContPid2->SetVarRange(7,-0.5,1.5);
+  fContPid2->SetVarName(8,"isPhiTrue^{Kn}");
+  fContPid2->SetVarRange(8,-0.5,1.5);
+  fContPid2->SetVarName(9,"N#sigma^{Kp}");
+  fContPid2->SetVarRange(9,1.25,6.25);
+  fContPid2->SetVarName(10,"N#sigma^{Kn}");
+  fContPid2->SetVarRange(10,1.25,6.25);
+  fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
+  fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
+  fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
+  fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
+  fContPid2->SetVarName(13,"#Psi");
+  fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
+
+
+
+  const Int_t nDETsignal = 100; // mass
+  Double_t binDETsignal[nDETsignal+1];
+  for(Int_t i=0;i<nDETsignal+1;i++){
+    binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
+  }
+  const Int_t nDETsignal2 = 10; // centrality
+  Double_t binDETsignal2[nDETsignal2+1];
+  for(Int_t i=0;i<nDETsignal2+1;i++){
+    binDETsignal2[i] = i*100./ nDETsignal2;
+  }
+  fContPid->AddSpecies("phi",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
+  fContPid2->AddSpecies("phi2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
+
+  fList->Add(fContPid);
+  fList->Add(fContPid2);
+
+  hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+
+  hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+  hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
+
+  fList2->Add(hMatching[0]);
+  fList2->Add(hMatching[1]);
+  fList2->Add(hMatching[2]);
+  fList2->Add(hMatching[3]);
+  fList2->Add(hTracking[0]);
+  fList2->Add(hTracking[1]);
+  fList2->Add(hTracking[2]);
+  fList2->Add(hTracking[3]);
+
+  fTOFTPCsignal = new TH2F("hTOFTPCsignal","TOF-TPC signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma_{TOF};N#sigma_{TPC}",100,-5,5,100,-5,5);
+  fList2->Add(fTOFTPCsignal);
+  fCombsignal = new TH1F("hCombsignal","Combined signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma",100,0,5);
+  fList2->Add(fCombsignal);
+
+  // QA plots
+  char name[100],title[400];
+  for(Int_t i=0;i < nCentrBin;i++){
+    snprintf(name,100,"hPiTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPiTPC[i]);
+
+    snprintf(name,100,"hKaTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fKaTPC[i]);
+
+    snprintf(name,100,"hPrTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPrTPC[i]);
+
+    snprintf(name,100,"hElTPCQA_%i",i);
+    snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fElTPC[i]);
+
+    snprintf(name,100,"hPiTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPiTOF[i]);
+
+    snprintf(name,100,"hKaTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fKaTOF[i]);
+
+    snprintf(name,100,"hPrTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fPrTOF[i]);
+
+    snprintf(name,100,"hElTOFQA_%i",i);
+    snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
+    fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
+    fList3->Add(fElTOF[i]);
+
+  }
+
+  // Post output data.
+  PostData(1, fList);
+  PostData(2, fList2);
+  PostData(3, fList3);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskPhiBayes::UserExec(Option_t *) 
+{
+    // Main loop
+    // Called for each event
+    
+    fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    if(!fOutputAOD){
+       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+       this->Dump();
+       return;
+    }
+    
+    Float_t zvtx = GetVertex(fOutputAOD);
+
+
+
+    //Get the MC object
+    if(fIsMC){
+      AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+      if (!mcHeader) {
+       AliError("Could not find MC Header in AOD");
+       return;
+      }
+    }
+
+    if (TMath::Abs(zvtx) < fVtxCut) {
+
+      //Centrality
+      Float_t v0Centr  = -10.;
+      Float_t trkCentr  = -10.;
+      AliCentrality *centrality = fOutputAOD->GetCentrality();
+      if (centrality){
+       trkCentr  = centrality->GetCentralityPercentile("V0M");
+       v0Centr = centrality->GetCentralityPercentile("TRK"); 
+      }
+
+      if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr>0){ // consistency cut on centrality selection
+        fCentrality = v0Centr;
+       Analyze(fOutputAOD); // Do analysis!!!!
+
+      }
+    }
+    
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPhiBayes::Analyze(AliAODEvent* aodEvent)
+{
+  Int_t ntrack = aodEvent->GetNumberOfTracks();
+
+  fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
+  fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
+  fPidKp=0,fPidKn=0;
+  fMassV0=-1;
+  
+  TLorentzVector phiV;
+  
+  Double_t px,py,pz,E;
+
+  Float_t invMmin = 0.985;
+  Float_t invMmax = 1.045;
+  
+  Int_t icentr = 8;
+  if(fCentrality < 0) icentr = 8;
+  else if(fCentrality < 10) icentr = 0;
+  else if(fCentrality < 20) icentr = 1;
+  else if(fCentrality < 30) icentr = 2;
+  else if(fCentrality < 40) icentr = 3;
+  else if(fCentrality < 50) icentr = 4;
+  else if(fCentrality < 60) icentr = 5;
+  else if(fCentrality < 70) icentr = 6;
+  else if(fCentrality < 80) icentr = 7;
+
+  Float_t addMismatchForMC = 0.005;
+  if(fCentrality < 50) addMismatchForMC += 0.01;
+  if(fCentrality < 20) addMismatchForMC += 0.02;
+
+  fPsi = 0;
+  /* Compute TPC EP */
+  Double_t Qx2 = 0, Qy2 = 0;
+  Double_t Qx3 = 0, Qy3 = 0;
+  for(Int_t iT = 0; iT < ntrack; iT++) {
+    AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
+    
+    if (!aodTrack){
+      continue;
+    }
+    
+    Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
+    
+    if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag) 
+      continue;
+    
+    Double_t b[2] = {-99., -99.};
+    Double_t bCov[3] = {-99., -99., -99.};
+    if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
+      continue;
+    
+    if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
+      continue;
+    
+    Qx2 += TMath::Cos(2*aodTrack->Phi()); 
+    Qy2 += TMath::Sin(2*aodTrack->Phi());
+    Qx3 += TMath::Cos(3*aodTrack->Phi()); 
+    Qy3 += TMath::Sin(3*aodTrack->Phi());
+    
+  }
+
+  fPsi = TMath::ATan2(Qy2, Qx2)/2.;
+
+  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+  AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
+
+  PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
+  PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
+  PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
+  PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
+
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+
+  Double_t probP[10],probN[10];
+  Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
+
+  
+  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+  TClonesArray *mcArray = NULL;
+  if (mcHeader)
+    mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+
+  Int_t nmc = 0;
+  if(mcArray)
+    nmc = mcArray->GetEntries();
+
+  for(Int_t i=0;i < ntrack;i++){ // loop on tracks
+    AliAODTrack* track = aodEvent->GetTrack(i);
+    
+    AliAODMCParticle *mcp = NULL;
+    Int_t pdg = 0;
+    
+    if (!track){
+      continue;
+    }
+    
+    Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
+    
+    Int_t label = -1;
+    if(mcArray){
+      label = track->GetLabel();
+      if(label != -1 && label < nmc){
+       label = TMath::Abs(label);
+       mcp = (AliAODMCParticle*)mcArray->At(label);
+       pdg = TMath::Abs(mcp->GetPdgCode());
+      }
+      else
+       label = -1;
+    }
+    else{
+      /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
+    }
+    
+    if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
+      hTracking[0]->Fill(track->Pt(),fCentrality);
+      if(pdg == 211)
+       hTracking[1]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 321)
+       hTracking[2]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 2212)
+       hTracking[3]->Fill(track->Pt(),fCentrality);
+      else if(! mcp){ // Fill matching histo with the prob
+       hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
+       hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
+       hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
+      }
+    }
+    
+    if(!tofMatch) continue;
+    
+    if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
+      hMatching[0]->Fill(track->Pt(),fCentrality);
+      if(pdg == 211)
+       hMatching[1]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 321)
+       hMatching[2]->Fill(track->Pt(),fCentrality);
+      else if(pdg == 2212)
+       hMatching[3]->Fill(track->Pt(),fCentrality);
+      else if(! mcp){ // Fill matching histo with the prob
+       hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
+       hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
+       hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
+      }
+    }
+  }
+  
+  Int_t pdg1 = -1;
+  Int_t pdg2 = -1;
+
+
+  // start analysis phi
+  for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
+    AliAODTrack* KpTrack = aodEvent->GetTrack(i);
+        
+    if (!KpTrack){
+      continue;
+    }
+    
+    if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
+
+    nSigmaComb=5;
+    fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
+    fPidKp=0;
+
+    UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
+
+    Double_t oldpP[10];
+    fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
+
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
+    fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
+    fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
+    fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+    nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
+    fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
+
+    if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
+
+    Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
+
+    if(mcArray){
+      Int_t labelK = TMath::Abs(KpTrack->GetLabel());
+      AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
+      pdg1 = TMath::Abs(mcp1->GetPdgCode());
+    }
+
+    fPidKp = Int_t(probP[3]*100);
+
+    if(tofMatch1){
+      if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
+       // remove this tof hit
+       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+       detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
+       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+       fPidKp = Int_t(probP[4]*100);
+       tofMatch1=0;
+      }
+      else{
+       if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
+       
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
+       fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
+       fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
+       fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+       nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
+       fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
+               
+       nSigmaTOF = TMath::Abs(nSigmaTOF);
+       
+       if(fIsMC){
+         Float_t mismAdd = addMismatchForMC;
+         if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
+         
+         if(gRandom->Rndm() < mismAdd){
+           Float_t etaAbs = TMath::Abs(KpTrack->Eta());
+           Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
+           channel = channel % 8736;
+           Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
+           
+           // generate random time
+           Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
+           Double_t times[10];
+           KpTrack->GetIntegratedTimes(times);
+           nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
+         }
+       }
+
+       if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
+
+       if(nSigmaTOF < 2) fPidKp += 256;
+       else if(nSigmaTOF < 3) fPidKp += 512;
+      }
+    }
+    
+    if(tofMatch1){
+      nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
+      if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
+       fCombsignal->Fill(nSigmaComb);
+      }
+    } else {
+      nSigmaComb = TMath::Abs(nSigmaTPC);
+    }
+
+    // use sigmaTOF instead of sigmaComb
+    //nSigmaComb = nSigmaTOF;
+
+    if(nSigmaComb < 2) nSigmaComb = 2;
+    else if(nSigmaComb < 3) nSigmaComb = 3;
+    else if(nSigmaComb < 5) nSigmaComb = 4.99;
+    else nSigmaComb = 6;
+
+    if(fPtKp > 4.299) fPtKp = 4.299;
+
+    for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
+      AliAODTrack* KnTrack = aodEvent->GetTrack(j);
+      
+      if (!KnTrack){
+       continue;
+      }
+
+      if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
+
+      px = KpTrack->Px() + KnTrack->Px();
+      py = KpTrack->Py() + KnTrack->Py();
+      pz = KpTrack->Pz() + KnTrack->Pz();
+      E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
+      E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
+
+      phiV.SetPxPyPzE(px,py,pz,E);
+      fMassV0 = phiV.M();
+      
+      if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
+
+      fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
+      fPidKn=0;
+
+      UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
+      Double_t oldpN[10];
+      fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
+
+      nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
+      
+      if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
+      
+      nSigmaComb2=5;
+
+      Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
+
+      if(mcArray){
+       Int_t labelK = TMath::Abs(KnTrack->GetLabel());
+       AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
+       pdg2 = TMath::Abs(mcp2->GetPdgCode());
+     }
+
+      fPidKn = Int_t(probN[3]*100);
+
+      if(tofMatch2){
+       if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
+         // remove this tof hit
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+         detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+         fPidKn = Int_t(probN[4]*100);
+         tofMatch2=0;
+       }
+       else{
+         if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
+         
+         nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
+                 
+         nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
+         
+         if(fIsMC){
+           Float_t mismAdd = addMismatchForMC;
+           if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
+           
+           if(gRandom->Rndm() < mismAdd){
+             Float_t etaAbs = TMath::Abs(KnTrack->Eta());
+             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
+             channel = channel % 8736;
+             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
+             
+             // generate random time
+             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
+             Double_t times[10];
+             KnTrack->GetIntegratedTimes(times);
+             nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
+           }
+         }
+
+         if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
+
+         if(nSigmaTOF2 < 2) fPidKn += 256;
+         else if(nSigmaTOF2 < 3) fPidKn += 512;
+       }
+      }
+
+      fPtPhi = phiV.Pt();
+      fEtaPhi = phiV.Eta();
+      fPhiPhi = phiV.Phi();
+
+      if(tofMatch2){
+       nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
+       if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
+         fCombsignal->Fill(nSigmaComb2);
+       }
+      } else {
+       nSigmaComb2 = TMath::Abs(nSigmaTPC2);
+      }
+
+      // use sigmaTOF instead of sigmaComb
+      //nSigmaComb2 = nSigmaTOF2;
+
+      if(nSigmaComb2 < 2) nSigmaComb2 = 2;
+      else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
+      else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
+      else nSigmaComb2 = 6;  
+
+      Bool_t isTrue = kFALSE;
+
+      if(mcArray){
+       Int_t labelP = TMath::Abs(KpTrack->GetLabel());
+       Int_t labelN = TMath::Abs(KnTrack->GetLabel());
+
+       if(labelP > -1 && labelN > -1){
+         AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
+         AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
+
+         Int_t mP = partP->GetMother();
+         Int_t mN = partN->GetMother();
+         if(mP == mN && mP > -1){
+           AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
+           Int_t pdgM = partM->GetPdgCode();
+           if(pdgM == 333) isTrue = kTRUE;
+         }
+       }
+
+      }
+
+      Double_t deltaphi1 = KpTrack->Phi() - fPsi;
+      Double_t deltaphi2 = KnTrack->Phi() - fPsi;
+
+      while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
+      while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
+      while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
+      while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
+
+      if(fPtKn > 4.299) fPtKn = 4.299;
+
+      Float_t xTOfill[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
+      
+
+      Int_t ipt = 0;
+      while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
+       ipt++;
+      }
+      ipt--;
+      if(ipt < 0) ipt = 0; // just to be sure
+
+      if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
+       if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
+        xTOfill[1] = KnTrack->Eta();
+       if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);// use tagging on positive track
+
+      }
+
+
+    }
+  } // end analysi phi
+
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
+{
+
+  Float_t zvtx = -999;
+
+  const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
+  if (!vtxAOD)
+    return zvtx;
+  if(vtxAOD->GetNContributors()>0)
+    zvtx = vtxAOD->GetZ();
+  
+  return zvtx;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
+{ 
+  // Terminate loop
+  Printf("Terminate()");
+}
+//-------------------------------------------------------------------------
+Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
+  if(!fIsMC) return 1; // used only on MC
+
+  Int_t run = fOutputAOD->GetRunNumber();
+  if( (run>=136851&&run<=139846) || (run>=165772 && run<=170718) ){ // LHC10h or LHC11h because of TOF matching window at 3 cm
+    Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs); 
+  
+    if(!(channel%20)) return 0; // 5% additional loss in MC
+  }
+
+  return 1;
+}
diff --git a/PWGPP/pid/AliAnalysisTaskPhiBayes.h b/PWGPP/pid/AliAnalysisTaskPhiBayes.h
new file mode 100644 (file)
index 0000000..6fd3dc2
--- /dev/null
@@ -0,0 +1,126 @@
+#ifndef ALIANALYSISTASKPHIBAYES_H
+#define ALIANALYSISTASKPHIBayes_H
+
+// ROOT includes
+#include <TObject.h>
+#include <TClonesArray.h>
+#include <TList.h>
+#include <TProfile.h>
+#include "AliPIDCombined.h"
+#include "TF1.h"
+
+// AliRoot includes
+#include <AliAnalysisTaskSE.h>
+#include <AliAODEvent.h>
+#include "AliPIDperfContainer.h"
+
+class TH1F;
+class TH2F;
+class AliESDtrackCuts;
+
+class AliAnalysisTaskPhiBayes : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskPhiBayes();
+  AliAnalysisTaskPhiBayes(const char *name);
+
+  virtual ~AliAnalysisTaskPhiBayes();
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *); 
+  
+  Double_t GetVtxCut() { return fVtxCut; }   
+  Double_t GetEtaCut() { return fEtaCut; }     
+  Double_t GetMinPt() { return fMinPt; }   
+
+  virtual void  SetVtxCut(Double_t vtxCut){fVtxCut = vtxCut;}
+  virtual void  SetEtaCut(Double_t etaCut){fEtaCut = etaCut;}
+  virtual void  SetMinPt(Double_t value) {fMinPt = value;}   
+
+  virtual void SetMC(Bool_t flag = kTRUE){fIsMC = flag;};
+  virtual void SetQA(Bool_t flag = kTRUE){fQAsw = flag;};
+
+  void SetTPCclusterN(Int_t ncl){fNcluster=ncl;};
+
+  static const Int_t nPtBin = 13;  //! # pt phi bins
+
+  void SetFilterBit(Int_t fb){fFilterBit=fb;};
+  Int_t GetFilterBit() const {return fFilterBit;};
+
+ private:
+  AliAnalysisTaskPhiBayes(const AliAnalysisTaskPhiBayes &old); 
+  AliAnalysisTaskPhiBayes& operator=(const AliAnalysisTaskPhiBayes &source); 
+
+  virtual Float_t GetVertex(AliAODEvent* aod) const;
+  virtual void Analyze(AliAODEvent* aodEvent); 
+  virtual Int_t IsChannelValid(Float_t etaAbs);
+
+  Double_t     fVtxCut;             // Vtx cut on z position in cm
+  Double_t     fEtaCut;             // Eta cut used to select particles
+  Double_t     fMinPt;              // Min pt - for histogram limits
+
+  Bool_t fIsMC; // if MC
+  Bool_t fQAsw;   // if QA
+
+  static const Int_t nCentrBin = 9; //! # cenrality bins
+
+  //
+  // Cuts and options
+  //
+
+  Int_t fNcluster;           // Numer of TPC cluster required
+  Int_t fFilterBit;          // filter bit to be used
+  TList *fList;              //! List for output objects
+  TList *fList2;             //! List for output objects
+  TList *fList3;             //! List for output objects
+
+  //
+  Float_t fCentrality;  //! current centrality for the tree
+  Float_t fPsi;         //! current Psi from TPC
+
+  Float_t fPtPhi;   //! variable to fill the tree
+  Float_t fPhiPhi;  //! variable to fill the tree
+  Float_t fEtaPhi;  //! variable to fill the tree
+  Float_t fMassV0; //! variable to fill the tree
+  Float_t fPtKp;  //! variable to fill the tree
+  Float_t fPhiKp; //! variable to fill the tree
+  Float_t fEtaKp; //! variable to fill the tree
+  Float_t fPtKn;  //! variable to fill the tree
+  Float_t fPhiKn; //! variable to fill the tree
+  Float_t fEtaKn; //! variable to fill the tree
+  Int_t fPidKp;   //! variable to fill the tree
+  Int_t fPidKn;   //! variable to fill the tree
+
+  TH2F *hMatching[4]; //! matching (all, pi, k ,p)
+  TH2F *hTracking[4]; //! tracking (all, pi, k ,p)
+  TH2F *fTOFTPCsignal;   //! histo with tof signal
+  TH1F *fCombsignal;   //! histo with tof signal
+
+  static Float_t fPtPhiMin[nPtBin];// ptmin bin
+  static Float_t fPtPhiMax[nPtBin];// ptmax bin
+
+  // QA plots
+  TH2F *fPiTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fKaTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fPrTPC[nCentrBin];//! TPC dE/dx plot
+  TH2F *fElTPC[nCentrBin];//! TPC dE/dx plot
+
+  TH2F *fPiTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fKaTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fPrTOF[nCentrBin];//! TPC dE/dx plot
+  TH2F *fElTOF[nCentrBin];//! TPC dE/dx plot
+
+  AliESDtrackCuts *fCutsDaughter;
+
+  AliPIDCombined *fPIDCombined;  //! PID combined object
+  AliPIDperfContainer* fContPid;      //! results for positive
+  AliPIDperfContainer* fContPid2;     //! results for negative
+
+  TH1F *fHmismTOF;         //! TOF mismatch distribution
+  TH1D *fHchannelTOFdistr; //! TOF channel distance w.r.t. IP
+
+  ClassDef(AliAnalysisTaskPhiBayes, 1);    //Analysis task for bayesian (K0s)
+};
+
+#endif
diff --git a/PWGPP/pid/AliPIDperfContainer.cxx b/PWGPP/pid/AliPIDperfContainer.cxx
new file mode 100644 (file)
index 0000000..8549f2a
--- /dev/null
@@ -0,0 +1,316 @@
+#include<stdio.h>
+#include "AliPIDperfContainer.h"
+#include "TList.h"
+
+ClassImp(AliPIDperfContainer);
+
+AliPIDperfContainer::AliPIDperfContainer(const char *name,const Int_t nvar,const Int_t* binVar) :
+  TNamed(name,name),
+  fNbinVar(new TArrayI(nvar)),
+  fXmin(new TArrayF(nvar)),
+  fXmax(new TArrayF(nvar)),
+  fNameVar(new TClonesArray("TNamed")),
+  fQA(new TClonesArray("TH2F"))
+{
+
+  for(Int_t i=0;i < GetNvar();i++){
+    (*fNbinVar)[i] = binVar[i];
+  }
+  
+  for(Int_t i=0; i<GetNvar();i++){
+    new((*fNameVar)[i]) TNamed("","");
+  }
+
+  SetTitleX("");
+  SetTitleY("");
+}
+
+AliPIDperfContainer::AliPIDperfContainer() :
+  TNamed("qa","qa"),
+  fNbinVar(new TArrayI(0)),
+  fXmin(new TArrayF(0)),
+  fXmax(new TArrayF(0)),
+  fNameVar(new TClonesArray("TNamed")),
+  fQA(new TClonesArray("TH2F"))
+{
+  SetTitleX("");
+  SetTitleY("");  
+}
+
+AliPIDperfContainer::~AliPIDperfContainer(){
+  for(Int_t i=fNameVar->GetEntries();i>0;i--){
+    delete fNameVar->At(i-1);
+    fNameVar->RemoveAt(i-1);
+  }
+
+  for(Int_t i=fQA->GetEntries();i>0;i--){
+    delete fQA->At(i-1);
+    fQA->RemoveAt(i-1);
+  }
+
+  delete fNbinVar;
+  delete fXmin;
+  delete fXmax;
+}
+
+void AliPIDperfContainer::Reset(){
+  for(Int_t i=fNameVar->GetEntries();i>0;i--){
+    delete fNameVar->At(i-1);
+    fNameVar->RemoveAt(i-1);
+  }
+
+  for(Int_t i=fQA->GetEntries();i>0;i--){
+    delete fQA->At(i-1);
+    fQA->RemoveAt(i-1);
+  }
+
+  delete fNbinVar;
+  delete fXmin;
+  delete fXmax;
+
+  fNbinVar = new TArrayI(0);
+  fXmin = new TArrayF(0);
+  fXmax = new TArrayF(0);
+
+}
+
+AliPIDperfContainer::AliPIDperfContainer(const AliPIDperfContainer &old) :
+  TNamed(old),
+  fNbinVar(NULL),
+  fXmin(NULL),
+  fXmax(NULL),
+  fNameVar(new TClonesArray("TNamed")),
+  fQA(new TClonesArray("TH2F"))
+{
+
+  fNbinVar = new TArrayI(old.GetNvar());
+  fXmin = new TArrayF(old.GetNvar());
+  fXmax = new TArrayF(old.GetNvar());
+
+  for(Int_t i=0; i<old.GetNhistos();i++){
+    new((*fQA)[i]) TH2F(*((TH2F *) old.GetQA(i)));
+  }
+
+  for(Int_t i=0; i<old.GetNvar();i++){
+    new((*fNameVar)[i]) TNamed(old.GetVarName(i),old.GetVarName(i));
+  }
+
+  for(Int_t i=0;i < GetNvar();i++){
+    (*fNbinVar)[i] = (*old.fNbinVar)[i];
+    (*fXmin)[i] = (*old.fXmin)[i];
+    (*fXmax)[i] = (*old.fXmax)[i];
+   }
+
+}
+
+AliPIDperfContainer& AliPIDperfContainer::operator=(const AliPIDperfContainer &old){
+
+  if(this != &old){
+     printf("different\n");
+   }
+
+   for(Int_t i=0; i<old.GetNhistos();i++){
+     new((*fQA)[i]) TH2F(*((TH2F *) old.GetQA(i)));
+   }
+   
+   fNbinVar = new TArrayI(old.GetNvar());
+   fXmin = new TArrayF(old.GetNvar());
+   fXmax = new TArrayF(old.GetNvar());
+   
+   for(Int_t i=0;i < old.GetNvar();i++){
+     (*fNbinVar)[i] = (*old.fNbinVar)[i];
+     (*fXmin)[i] = (*old.fXmin)[i];
+     (*fXmax)[i] = (*old.fXmax)[i];
+   }
+   
+   fNameVar = new TClonesArray("TNamed");
+   for(Int_t i=0; i<old.GetNvar();i++){
+     new((*fNameVar)[i]) TNamed(old.GetVarName(i),old.GetVarName(i));
+   }
+
+  return *this;
+}
+  
+Int_t AliPIDperfContainer::GetNspecies() const{
+  Int_t n = fQA->GetEntries();
+
+  for(Int_t i=0;i < GetNvar();i++){
+    n /= (*fNbinVar)[i];
+  }
+
+  return n;
+}
+
+void AliPIDperfContainer::AddSpecies(const char *name,Int_t nXbin,const Double_t *xbin,Int_t nYbin,const Double_t *ybin){
+  
+  Bool_t kErr = kFALSE;
+  for(Int_t i=0;i < GetNvar();i++){ // check the var ranges are set properly    
+    if((*fNbinVar)[i] < 1 || (*fXmin)[i] >= (*fXmax)[i]){
+    printf("var ranges are not set properly for variable %i please chek it before to define the species\n",i);
+       kErr = kTRUE;
+    }
+  }
+  if(kErr){
+    printf("AddSpecies: NOTHING DONE\n");
+    return;
+  }
+
+  Int_t ncomb = 1;
+  for(Int_t i=0;i < GetNvar();i++){
+    ncomb *= (*fNbinVar)[i];
+  }
+
+  char nameHisto[200];
+  char title[300];
+  char title2[300];
+  for(Int_t i=0; i < ncomb;i++){
+    snprintf(nameHisto,200,"%s_%s_%i",GetName(),name,i);
+    snprintf(title,300,"%s",name);
+    Int_t ncombTemp = i;
+    for(Int_t j=0;j < GetNvar();j++){
+      Int_t ibin = ncombTemp%(*fNbinVar)[j];
+      snprintf(title2,300,"%s",title);
+      snprintf(title,300,"%s_%04.1f<%s<%04.1f",title2,(*fXmin)[j] + ((*fXmax)[j]-(*fXmin)[j])/(*fNbinVar)[j]*ibin,fNameVar->At(j)->GetName(),(*fXmin)[j] + ((*fXmax)[j]-(*fXmin)[j])/(*fNbinVar)[j]*(ibin+1));
+      ncombTemp /= (*fNbinVar)[j];
+    }
+
+    new((*fQA)[GetNhistos()]) TH2F(nameHisto,title,nXbin,xbin,nYbin,ybin);
+    ((TH2F *) GetQA(GetNhistos()-1))->GetXaxis()->SetTitle(fTitleX);
+    ((TH2F *) GetQA(GetNhistos()-1))->GetYaxis()->SetTitle(fTitleY);
+  }
+}
+
+Int_t AliPIDperfContainer::Add(const AliPIDperfContainer *oth){
+  if(GetNhistos() == oth->GetNhistos()){
+    for(Int_t i=0;i < GetNhistos();i++){
+      GetQA(i)->Add(oth->GetQA(i));
+    }
+    return 0;
+  }
+  else{
+    printf("ADD error: number of objects is different (%i != %i)\n",GetNhistos(),oth->GetNhistos());
+    return 1;
+  }
+}
+
+void AliPIDperfContainer::SetVarRange(Int_t ivar,Float_t xMin,Float_t xMax){
+  if(!GetNhistos()){
+    (*fXmin)[ivar]=xMin;
+    (*fXmax)[ivar]=xMax;
+  }
+  else{ // to avoid different range among the histos
+    printf("Ranges should be set before to define the species\nNOTHING DONE\n");
+  }
+
+}
+
+void AliPIDperfContainer::Fill(Int_t species,Float_t var1,Float_t var2,Float_t x[]){
+  Int_t ncomb = 1;
+  Int_t histo = 0;
+  for(Int_t i=0;i < GetNvar();i++){
+    Int_t ibin = GetBin(i,x[i]);
+    if(ibin < 0 || ibin >= (*fNbinVar)[i]){
+      printf("%i) %i not good w.r.t. %i (%f) (%f,%f)\n",i,ibin,(*fNbinVar)[i],x[i],(*fXmin)[i],(*fXmax)[i]);
+      return;
+    }
+    histo += ncomb * ibin;
+    ncomb *= (*fNbinVar)[i];
+  }
+  histo += species*ncomb;
+  DirectFill(histo,var1,var2);
+ };
+
+TH2F *AliPIDperfContainer::GetQA(Int_t species,Float_t x[]) const{
+  Int_t ncomb = 1;
+  Int_t histo = 0;
+  for(Int_t i=0;i < GetNvar();i++){
+    Int_t ibin = GetBin(i,x[i]);
+    if(ibin < 0 || ibin >= (*fNbinVar)[i]){
+      printf("%i) %i not good w.r.t. %i (%f) (%f,%f)\n",i,ibin,(*fNbinVar)[i],x[i],(*fXmin)[i],(*fXmax)[i]);
+      return NULL;
+    }
+    histo += ncomb * ibin;
+    ncomb *= (*fNbinVar)[i];
+  }
+  histo += species*ncomb;
+
+
+  return GetQA(histo);
+
+}
+
+TH2F *AliPIDperfContainer::GetQA(Int_t species,Float_t xMin[],Float_t xMax[]) const{
+  if(GetNvar()){
+    char title[300];
+    char title2[300];
+    Int_t ncomb = 1;
+    for(Int_t i=0;i < GetNvar();i++){
+      ncomb *= (*fNbinVar)[i];
+    }
+
+    TH2F *htemplate = GetQA(species*ncomb);
+    TH2F *temp = new TH2F(*htemplate);
+    temp->SetName("histo");
+    temp->Reset();
+    snprintf(title,300,"%i",species);
+    for(Int_t i=0;i < GetNvar();i++){
+      Int_t imin = GetBin(i,xMin[i]);
+      if(imin < 0) imin = 0;
+      else if(imin >= (*fNbinVar)[i]) imin = (*fNbinVar)[i]-1;
+      Int_t imax = GetBin(i,xMax[i]);
+      if(imax < imin) imax = imin;
+      else if(imax >= (*fNbinVar)[i]) imax = (*fNbinVar)[i]-1;
+      snprintf(title2,300,"%s",title);
+      snprintf(title,300,"%s_%04.1f<%s<%04.1f",title2,
+                (*fXmin)[i] + ((*fXmax)[i]-(*fXmin)[i])/(*fNbinVar)[i]*imin,
+                fNameVar->At(i)->GetName(),
+                (*fXmin)[i] + ((*fXmax)[i]-(*fXmin)[i])/(*fNbinVar)[i]*(imax+1));
+    }
+    temp->SetTitle(title);
+
+    for(Int_t i=species*ncomb;i < (species+1)*ncomb;i++){
+      Bool_t kGood = kTRUE;
+
+      Int_t ncombTemp = i;
+      for(Int_t j=0;j < GetNvar();j++){
+       Int_t imin = GetBin(j,xMin[j]);
+       if(imin < 0) imin = 0;
+       else if(imin >= (*fNbinVar)[j]) imin = (*fNbinVar)[j]-1;
+       Int_t imax = GetBin(j,xMax[j]);
+       if(imax < imin) imax = imin;
+       else if(imax >= (*fNbinVar)[j]) imax = (*fNbinVar)[j]-1;
+       
+       Int_t ibin = ncombTemp%(*fNbinVar)[j];
+       ncombTemp /= (*fNbinVar)[j];
+
+       if(ibin < imin || ibin > imax){
+         kGood = kFALSE;
+         j = GetNvar();
+       }
+      }
+
+      if(kGood) temp->Add(GetQA(i));
+    }
+    return temp;
+  }
+
+  return GetQA(species);
+}
+
+Long64_t AliPIDperfContainer::Merge(TCollection* list){
+  Long64_t res=0;
+  if (!list) return 0;
+  if (list->IsEmpty()) return 0;
+
+  TList *listObj = new TList();
+  listObj->AddAll(list);
+
+  for(Int_t i=0;i < listObj->GetEntries();i++){
+    AliPIDperfContainer *obj = (AliPIDperfContainer *) listObj->At(i);
+    Add(obj);
+    res++;
+  }
+  return res;
+}
diff --git a/PWGPP/pid/AliPIDperfContainer.h b/PWGPP/pid/AliPIDperfContainer.h
new file mode 100644 (file)
index 0000000..5efaa45
--- /dev/null
@@ -0,0 +1,77 @@
+#ifndef ALIPIDPERFCONTAINER_H
+#define ALIPIDPERFCONTAINER_H
+
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliPIDperfContainer.h 49869 2013-07-10 04:49:51Z fnoferin $ */
+
+/////////////////////////////////////////////////
+//                                             //
+//         output v2-VZERO Class               //
+//           noferini@bo.infn.it               //
+/////////////////////////////////////////////////
+
+#include "TH2F.h"
+#include "TClonesArray.h"
+#include "TArrayI.h"
+#include "TArrayF.h"
+
+class AliPIDperfContainer : public TNamed
+{
+ public:
+  AliPIDperfContainer(const char *name,const Int_t nvar,const Int_t* binVar);
+  AliPIDperfContainer();
+  ~AliPIDperfContainer();
+  AliPIDperfContainer(const AliPIDperfContainer &old);
+  AliPIDperfContainer& operator=(const AliPIDperfContainer &source);
+
+  Int_t GetNhistos() const {return fQA->GetEntries();};
+  Int_t GetNspecies() const;
+  TH2F *GetQA(Int_t histo) const {return ((TH2F *) fQA->At(histo));};
+  TH2F *GetQA(Int_t species,Float_t x[]) const;
+  TH2F *GetQA(Int_t species,Float_t xMin[],Float_t xMax[]) const;
+  void DirectFill(Int_t histo,Float_t var1,Float_t var2){GetQA(histo)->Fill(var1,var2);};
+  void Fill(Int_t species,Float_t var1,Float_t var2,Float_t x[]);
+
+  void AddSpecies(const char *name,Int_t nXbin,const Double_t *xbin,Int_t nYbin,const Double_t *ybin);
+
+  const char *GetSpeciesName(Int_t species){if(species<GetNspecies()) return GetQA(species*GetNhistos()/GetNspecies())->GetName();else return "";};
+
+  Int_t Add(const AliPIDperfContainer *oth);
+
+  Int_t GetNvar() const {return fNbinVar->GetSize();};
+  Int_t GetNbinVar(Int_t ivar) const {return (*fNbinVar)[ivar];};
+
+  void SetVarRange(Int_t ivar,Float_t xMin,Float_t xMax);
+  void SetVarName(Int_t ivar,const char *name){TNamed *atemp = (TNamed *) fNameVar->At(ivar); atemp->SetName(name);};
+
+  Float_t GetXmin(Int_t ivar) const {return (*fXmin)[ivar];};
+  Float_t GetXmax(Int_t ivar) const {return (*fXmax)[ivar];};
+  const char *GetVarName(Int_t ivar) const {TNamed *atemp = (TNamed *) fNameVar->At(ivar); return atemp->GetName();};
+
+  Int_t GetBin(Int_t ivar,Float_t x) const {return Int_t((x-(*fXmin)[ivar])/((*fXmax)[ivar]-(*fXmin)[ivar])*(*fNbinVar)[ivar]);};
+
+  Long64_t Merge(TCollection* list);
+
+  void Reset();
+
+  void SetTitleX(const char *title){snprintf(fTitleX,100,"%s",title);};
+  void SetTitleY(const char *title){snprintf(fTitleY,100,"%s",title);};
+
+ private:
+  TArrayI *fNbinVar;
+  TArrayF *fXmin,*fXmax;
+  TClonesArray *fNameVar;
+
+  TClonesArray *fQA;
+
+  char fTitleX[100]; // title for X-axis
+  char fTitleY[100]; // title for Y-axis
+
+  ClassDef(AliPIDperfContainer,1)  // qa vzero outuput object
+};
+#endif
+
+
diff --git a/PWGPP/pid/doeffKa.C b/PWGPP/pid/doeffKa.C
new file mode 100644 (file)
index 0000000..28ed25b
--- /dev/null
@@ -0,0 +1,525 @@
+TObject* fContPid1;
+TObject* fContPid2;
+const Int_t nBinPid = 14; // pt,eta, ptPip, ptPin, PPip, PPin, TOF3sigmaPip, TOF3sigmaPin, isPhiTrue, nsigmaPip, nsigmaPin
+// 0.985 < mass < 1.045 (60) and 0 < centrality < 100 (10)
+Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaPi*/,20/*pt+*/,20/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+Float_t xmin[nBinPid] = {1,-0.8,0.3,0.3,0,0,-0.5,-0.5,-0.5,0,0,-TMath::Pi(),-TMath::Pi(),-TMath::Pi()/2};
+Float_t xmax[nBinPid] = {5,0.8,4.3,4.3,1,1,1.5,1.5,1.5,7.5,7.5,TMath::Pi(),TMath::Pi(),TMath::Pi()/2};
+
+TF1 *fsign;
+TF1 *fall;
+TF1 *fback;
+
+Int_t ifunc=0;
+
+Float_t fitmin = 0.99;
+Float_t fitmax = 1.045;
+
+Int_t cmin = 1;
+Int_t cmax = 8;
+
+Float_t weightS = -1.;
+
+Int_t rebinsize = 4;
+
+Int_t parplotted = 2;
+
+Bool_t isMC = kFALSE; // don't change this (is set automatically)
+Bool_t selectTrue = kTRUE; // put it to true to remove background (only for MC)
+Bool_t keepTrue = kFALSE; // put it to false to fit only background (only for MC)
+
+Bool_t kGoodMatch = kFALSE; // to check good matching
+
+Bool_t kSigma2vs3 = kFALSE; // to check good matching
+
+Bool_t require5sigma = kFALSE; // don't touch this flag
+
+Bool_t bayesVsigma = kFALSE; // only to do checks
+
+Bool_t kTOFmatch = kFALSE; // for combined PID requires TOF matching
+
+Bool_t kOverAll = kFALSE;
+
+Bool_t kLoaded=kFALSE;
+LoadLib(){
+  weightS = -1.;
+
+  if(! kLoaded){
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libPhysics.so");
+    gSystem->Load("libTree.so");
+    gSystem->Load("libMinuit.so");
+    gSystem->Load("libSTEERBase.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libAOD.so");
+    gSystem->Load("libESD.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libANALYSISalice.so");
+    gSystem->Load("libCORRFW.so");
+    gSystem->Load("libNetx.so");
+    gSystem->Load("libPWGPPpid.so");
+
+    TFile *f = new TFile("AnalysisResults.root");
+    f->ls();
+    TList *l = (TList *) f->Get("contPhiBayes1");
+    l->ls();
+    fContPid1 = (AliPIDperfContainer *) l->At(0);
+    fContPid2 = (AliPIDperfContainer *) l->At(1);
+  }
+  kLoaded = kTRUE;
+
+  // check if MC
+  Float_t x[] = {xmin[0]+0.001,xmin[1]+0.001,xmin[2]+0.001,xmin[3]+0.001,xmin[4]+0.001,xmin[5]+0.001,xmin[6]+0.001,xmin[7]+0.001,1/*trueMC*/,xmin[9],xmin[10]};
+  Float_t x2[] = {xmax[0],xmax[1],xmax[2],xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],xmax[8],xmax[9],xmax[10]};
+
+  AliPIDperfContainer *tmp = fContPid1;
+  TH1D *h = tmp->GetQA(0, x, x2)->ProjectionX("checkMC");
+
+  if(h->GetEntries()) isMC = kTRUE;
+  else isMC=kFALSE;
+
+  if(!isMC){
+    selectTrue = kFALSE;
+    keepTrue = kTRUE;
+  }
+  else{
+    printf("MC truth found!!!!!!\nIt is MC!!!!!!");
+  }
+
+  fsign = new TF1("fsign","[0]*TMath::Voigt(x-[1],[3],[2])*(x>0.987)*(x > 1.005 && x < 1.035 || [4])",fitmin,fitmax);
+  fback = new TF1("fback","([0]*sqrt(x-0.987) + [1]*(x-0.987) + [2]*sqrt(x-0.987)*(x-0.987) +[3]*(x-0.987)*(x-0.987)+[4]*(x-0.987)*(x-0.987)*sqrt(x-0.987))*(x>0.987)",fitmin,fitmax);
+  fall = new TF1("fall","([0]*TMath::Voigt(x-[1],[3],[2])*(x > 1.005 && x < 1.035 || [9]) + [4]*sqrt(x-0.987) + [5]*(x-0.987) + [6]*sqrt(x-0.987)*(x-0.987) +[7]*(x-0.987)*(x-0.987)+[8]*(x-0.987)*(x-0.987)*sqrt(x-0.987))*(x>0.987)",0.987,1.05);
+
+  if(isMC){
+    fsign->SetParameter(4,0);
+    fall->FixParameter(9,0);
+  }
+  else{
+    fsign->SetParameter(4,1);
+    fall->FixParameter(9,1);
+  }
+
+  fsign->SetLineColor(2);
+  fback->SetLineColor(4);
+
+  if(kSigma2vs3){
+    kGoodMatch=kFALSE;
+    kOverAll = 0;
+  }
+
+  if(bayesVsigma){
+    kOverAll = 0;
+    kGoodMatch=kFALSE;
+    kSigma2vs3=kFALSE;
+    kTOFmatch=kTRUE;
+    weightS = -0.7;
+  }
+}
+
+doeffKa(Int_t pos=1,Float_t prob=0.1,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+  LoadLib();
+
+  Int_t nptbin = binPid[2];
+  Float_t minptbin = xmin[2];
+  Float_t maxptbin = xmax[2];
+
+  if(pos == 0){
+    nptbin = binPid[3];
+    minptbin = xmin[3];
+    maxptbin = xmax[3];
+  }
+
+  if(prob > 0.1999){
+    kGoodMatch = kFALSE;
+    kSigma2vs3 = kFALSE;
+    if(! kOverAll) require5sigma = kTRUE;
+    if(!isMC) weightS = -0.95;
+  }
+
+  TCanvas *c = new TCanvas();
+  c->Divide((nptbin+1)/2,2);
+  TH2F *hh.*hh2;
+  TH1D *h,*h2;
+  char name[100];
+  Float_t b[50][3];
+
+  Double_t xx[50],yy[50];
+  Double_t exx[50],eyy[50];
+
+  for(Int_t i=0;i < nptbin;i++){
+    c->cd(i+1)->SetLogy();
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    xx[i] = (ptmin+ptmax)/2;
+    exx[i] = (-ptmin+ptmax)/2;
+
+    Float_t pp=0.1;
+    if(prob < 0.2) pp = 0.;
+    if(pos) hh=GetHistoKap(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh=GetHistoKan(ptmin,ptmax,pp,0.0);
+    sprintf(name,"TOF matched: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
+    hh->SetTitle(name);
+    sprintf(name,"hNoPid%i",i);
+    
+    pp=prob;
+    if(prob < 0.2) pp = 0.1;
+    if(pos) hh2=GetHistoKap(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh2=GetHistoKan(ptmin,ptmax,pp,0.0);
+    AddHisto(hh,hh2,weightS);
+
+    h = hh->ProjectionX(name,cmin,cmax);
+    h->RebinX(rebinsize);
+    h->Draw("ERR");
+    h->SetMarkerStyle(24);
+    b[i][0]=-1;
+    Int_t ntrial = 0;
+    Float_t chi2 = 10000;
+    while(ntrial < 10 && (chi2 > 20 + 1000*selectTrue)){
+      fit(h,b[i],"WW","",xx[i]);
+      c1->Update();
+//       getchar();
+      fit(h,b[i],"","",xx[i]);
+      ntrial++;
+      chi2 = b[i][2];
+      printf("chi2 = %f\n",chi2);
+      c1->Update();
+//       getchar();
+      
+    }
+
+    yy[i] = fall->GetParameter(parplotted);
+    eyy[i] = fall->GetParError(parplotted);
+  }
+
+  TGraphErrors *gpar = new TGraphErrors(nptbin,xx,yy,exx,eyy);
+  c->cd(8);
+  gpar->Draw("AP");
+  gpar->SetMarkerStyle(20);
+
+  TCanvas *c2 = new TCanvas();
+  c2->Divide((nptbin+1)/2,2);
+  Float_t b2[50][3];
+
+  for(Int_t i=0;i < nptbin;i++){
+    c2->cd(i+1);
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    Float_t pp=prob;
+    if(prob < 0.2) pp = 0.1;
+    if(pos) hh=GetHistoKap(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh=GetHistoKan(ptmin,ptmax,pp,0.0);
+    sprintf(name,"P_{TOF} > 0.8: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
+    hh->SetTitle(name);
+    sprintf(name,"hPid60_%i",i);
+    h = hh->ProjectionX(name,cmin,cmax);
+    h->RebinX(rebinsize);
+    h->Draw("ERR");
+    h->SetMarkerStyle(24);
+    b2[i][0]=-1;
+    Int_t ntrial = 0;
+    Float_t chi2 = 10000;
+    while(ntrial < 40 && (chi2 > 20 + 1000*selectTrue)){
+      fit(h,b2[i],"WW","");
+      fit(h,b2[i],"","");
+      ntrial++;
+      chi2 = b2[i][2];
+      printf("chi2 = %f\n",chi2);
+    }
+    yy[i] = fall->GetParameter(parplotted);
+    eyy[i] = fall->GetParError(parplotted);
+
+  }
+
+  TGraphErrors *gpar2 = new TGraphErrors(nptbin,xx,yy,exx,eyy);
+  c2->cd(8);
+  gpar2->Draw("AP");
+  gpar2->SetMarkerStyle(20);
+  
+  Double_t xpt[50],expt[50],eff[50],efferr[50];
+  for(Int_t i=0;i<nptbin;i++){
+    printf("%f +/- %f -  %f +/- %f\n",b[i][0],b[i][1],b2[i][0],b2[i][1]);
+
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    xpt[i] = (ptmin+ptmax)/2;
+    expt[i] = (-ptmin+ptmax)/2;
+    eff[i] = b2[i][0]/(b[i][0]-b2[i][0]*weightS);
+
+    b[i][0] = b[i][0]-b2[i][0]*weightS;
+
+    efferr[i] = TMath::Sqrt(b[i][1]*b[i][1]/b[i][0]/b[i][0] + b2[i][1]*b2[i][1]/b2[i][0]/b2[i][0])*(b2[i][0]+b2[i][1])*(1+weightS*(b2[i][0]-b2[i][1])/b[i][0])/b[i][0];//*(1-eff[i]);//der2*der2*(b[i][1]*b[i][1] - b2[i][1]*b2[i][1]));
+
+    if(TMath::Abs(efferr[i]) > 1)efferr[i]=1;
+  }
+  new TCanvas();
+  TGraphErrors *geff = new TGraphErrors(nptbin,xpt,eff,expt,efferr);
+  geff->Draw("AP");
+
+  char flag[100];
+  sprintf(flag,"");
+
+  if(isMC){
+    if(selectTrue) sprintf(flag,"true");
+    else if(!keepTrue) sprintf(flag,"back");
+  }
+
+  char flag2[100];
+  sprintf(flag2,"");
+
+  char etarange[100];
+  sprintf(etarange,"_%.1f-%.1f_",etaminkp,etamaxkp);
+
+  if(kGoodMatch)
+    sprintf(flag2,"GM");
+
+  if(bayesVsigma)
+    sprintf(flag2,"BayesVsSigma");
+
+  if(kSigma2vs3)
+    sprintf(flag2,"Sigma2vs3");
+
+  if(kOverAll)
+    sprintf(flag2,"OverAll");
+
+  if(pos){
+    if(prob >=0.2) sprintf(name,"kaonPos%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
+    else sprintf(name,"kaonPos%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
+  }
+  else{
+    if(prob >=0.2) sprintf(name,"kaonNeg%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
+    else sprintf(name,"kaonNeg%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
+  }
+
+  TFile *fout = new TFile(name,"RECREATE");
+  geff->Write();
+  fout->Close();
+
+
+  TF1 *ff = new TF1("ff","[0] - [1]*TMath::Exp([2]*x)",0,3);
+  ff->SetParameter(0,0.67);
+  ff->SetParameter(1,1.14383e+00);
+  ff->SetParameter(2,-2.29910);
+  ff->SetLineColor(4);
+  ff->SetLineColor(2);
+  ff->Draw("SAME");
+
+  TF1 *ff2 = new TF1("ff2","[0] - [1]*TMath::Exp([2]*x)",0,3);
+  ff2->SetParameter(0,0.67);
+  ff2->SetParameter(1,9.23126e-01);
+  ff2->SetParameter(2,-1.851);
+  ff2->SetLineColor(4);
+  ff2->Draw("SAME");
+}
+
+TH2F *GetHistoKap(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkp=0,Float_t pMinkn=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+
+  Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,pt+0.001,xmin[3]+0.001,pMinkp+0.001,pMinkn+0.001,(pMinkp>0.09)+0.001,kTOFmatch+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
+  Float_t x2[] = {xmax[0],etamaxkp-0.001,ptM-0.001,xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
+
+  if(kGoodMatch){
+    x[6] = 1.0001;
+    if(pMinkp > 0)
+      x2[9] = 4.9;
+      
+  }
+
+  if(kTOFmatch){
+    x[6] = 1.0001;
+  }
+
+  if(kSigma2vs3){
+    x[6] = 1.0001;
+    x2[9] = 3;
+    if(pMinkp > 0)
+      x2[9] = 2;
+  }
+
+  if(bayesVsigma){
+    if(pMinkp > 0){
+      x[4] = 0.2001;
+      x2[9] = 5;
+    }
+    else{
+      x2[9] = 3;
+    }
+
+    
+  }
+
+  if(require5sigma) x2[9] = 4.9;
+
+  AliPIDperfContainer *tmp = fContPid1;
+
+  TH2F *h = tmp->GetQA(0, x, x2);
+
+  h->GetXaxis()->SetTitle("M_{#phi} (GeV/#it{c}^{2})");
+  h->GetYaxis()->SetTitle("centrality [%]");
+
+  return h;
+}
+
+TH2F *GetHistoKan(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkn=0,Float_t pMinkp=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+
+  Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,xmin[2]+0.001,pt+0.001,pMinkp+0.001,pMinkn+0.001,kTOFmatch+0.001,(pMinkn>0.09)+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
+  Float_t x2[] = {xmax[0],etamaxkp-0.001,xmax[2],ptM-0.001,xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
+
+  if(kGoodMatch){
+    x[7] = 1.0001;
+    if(pMinkn > 0)
+      x2[10] = 4.9;
+      
+  }
+
+  if(kTOFmatch){
+    x[7] = 1.0001;
+  }
+
+  if(kSigma2vs3){
+    x[7] = 1.0001;
+    x2[10] = 3;
+    if(pMinkn > 0)
+      x2[10] = 2;
+  }
+ if(bayesVsigma){
+    if(pMinkn > 0){
+      x[5] = 0.2001;
+      x2[10] = 5;
+    }
+    else{
+      x2[10] = 3;
+    }    
+  }
+
+  if(require5sigma) x2[10] = 4.9;
+
+  AliPIDperfContainer *tmp = fContPid2;
+
+  TH2F *h = tmp->GetQA(0, x, x2);
+
+  h->GetXaxis()->SetTitle("M_{#phi} (GeV/#it{c}^{2})");
+  h->GetYaxis()->SetTitle("centrality [%]");
+
+  return h;
+}
+
+
+fit(TH1D *h,Float_t *a=NULL,char *opt="",char *opt2="",Float_t pt=1.5){
+  if(h->GetEntries() < 1){
+    if(a){
+      a[0]=0.01;
+      a[1]=1;
+    }
+    return;
+  }
+
+
+ fall->SetParameter(0,100);
+ fall->SetParameter(0,1.01898 + 2.4e-04*pt);
+ fall->SetParameter(2,0.0044);
+ fall->SetParameter(3,0.0015);
+
+ fall->SetParLimits(0,-100,100000);
+ fall->SetParLimits(1,1.01898 + 2.4e-04*pt-1e-03,1.01898 + 2.4e-04*pt+1e-03);
+ fall->SetParLimits(2,0.0005,0.006);
+ fall->SetParLimits(3,0.001,0.0017);
+
+ fall->FixParameter(1,1.01884 + 2.9891e-04*pt);
+ fall->FixParameter(2,0.0044);
+ fall->FixParameter(3,7.57574e-04 + 3.85408e-04*pt);
+
+ fall->ReleaseParameter(4);
+ fall->ReleaseParameter(5);
+ fall->ReleaseParameter(6);
+ fall->ReleaseParameter(7);
+ fall->ReleaseParameter(8);
+
+
+ if(!kGoodMatch && !kSigma2vs3){
+   if(pt > 1.5){
+     fall->FixParameter(7,0);   
+     fall->FixParameter(8,0);   
+   }
+   if(pt > 1.7){
+     fall->FixParameter(6,0);   
+   }
+ }
+
+ if(selectTrue){
+   fall->FixParameter(4,0);
+   fall->FixParameter(5,0);
+   fall->FixParameter(6,0);
+   fall->FixParameter(7,0);
+   fall->FixParameter(8,0);
+ }
+
+ char name[100];
+ TF1 *ftmp=fall;
+
+ TF1 *ftmp2=new TF1(*fsign);
+ sprintf(name,"fsign%i",ifunc);
+ ftmp2->SetName(name);
+
+ TF1 *ftmp3=new TF1(*fback);
+ sprintf(name,"ftmp3%i",ifunc);
+ ftmp3->SetName(name);
+
+ ifunc++;
+
+ h->Fit(ftmp,opt,opt2,fitmin,fitmax);
+ h->Draw("ERR");
+
+ ftmp2->SetParameter(0,ftmp->GetParameter(0));
+ ftmp2->SetParameter(1,ftmp->GetParameter(1));
+ ftmp2->SetParameter(2,ftmp->GetParameter(2));
+ ftmp2->SetParameter(3,ftmp->GetParameter(3));
+ ftmp2->Draw("SAME");
+ ftmp3->SetParameter(0,ftmp->GetParameter(4));
+ ftmp3->SetParameter(1,ftmp->GetParameter(5));
+ ftmp3->SetParameter(2,ftmp->GetParameter(6));
+ ftmp3->SetParameter(3,ftmp->GetParameter(7));
+ ftmp3->SetParameter(4,ftmp->GetParameter(8));
+ ftmp3->Draw("SAME");
+
+ Float_t mean = ftmp->GetParameter(1);
+ Float_t sigma = 0.0044;//TMath::Abs(ftmp->GetParameter(2));
+
+ Float_t signI = ftmp2->Integral(mean-10*sigma,mean+10*sigma)/h->GetBinWidth(1);
+ Float_t backI = ftmp3->Integral(mean-3*sigma,mean+3*sigma)/h->GetBinWidth(1);
+
+ Float_t errI = TMath::Sqrt(ftmp->GetParError(0)*ftmp->GetParError(0)/(0.001+ftmp->GetParameter(0))/(0.001+ftmp->GetParameter(0)));
+
+ printf("signal(5 sigma) = %f +/- %f(fit) +/- %f(stat)\n",signI,errI*signI,TMath::Sqrt(signI));
+ printf("backgr(3sigma) = %f\n",backI);
+ printf("significance(3 sigma) = %f\n",signI/sqrt(signI+backI));
+
+ if(a){
+   a[0]=signI;
+   a[1]=signI*errI*signI*errI + signI;
+   a[1] = TMath::Sqrt(a[1]);
+   if(ftmp->GetNDF()) a[2] = ftmp->GetChisquare()/ftmp->GetNDF();
+
+
+   if(selectTrue){
+     a[0] = h->GetEntries();
+     a[1] = TMath::Sqrt(a[0]);
+   }
+ }
+}
+
+AddHisto(TH2F *h1,TH2F *h2,Float_t w){
+  Int_t nbinx = h1->GetNbinsX();
+  Int_t nbiny = h1->GetNbinsY();
+
+  for(Int_t i=1;i<=nbinx;i++){
+    for(Int_t j=1;j<=nbiny;j++){
+      Float_t val = h1->GetBinContent(i,j) + h2->GetBinContent(i,j)*w;
+      Float_t err = TMath::Min(TMath::Sqrt(val),val);
+      h1->SetBinContent(i,j,val);
+      h1->SetBinError(i,j,err);
+    }
+  }
+}
diff --git a/PWGPP/pid/doeffPi.C b/PWGPP/pid/doeffPi.C
new file mode 100644 (file)
index 0000000..0497c91
--- /dev/null
@@ -0,0 +1,492 @@
+TObject* fContPid1;
+TObject* fContPid2;
+const Int_t nBinPid = 14; // pt,eta, ptPip, ptPin, PPip, PPin, TOF3sigmaPip, TOF3sigmaPin, isPhiTrue, nsigmaPip, nsigmaPin
+// 0.985 < mass < 1.045 (60) and 0 < centrality < 100 (10)
+Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaPi*/,20/*pt+*/,20/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+Float_t xmin[nBinPid] = {1,-0.8,0.3,0.3,0,0,-0.5,-0.5,-0.5,0,0,-TMath::Pi(),-TMath::Pi(),-TMath::Pi()/2};
+Float_t xmax[nBinPid] = {5,0.8,4.3,4.3,1,1,1.5,1.5,1.5,7.5,7.5,TMath::Pi(),TMath::Pi(),TMath::Pi()/2};
+
+TF1 *fsign;
+TF1 *fall;
+TF1 *fback;
+
+Int_t ifunc=0;
+
+Float_t fitmin = 0.3;
+Float_t fitmax = 0.7;
+
+Int_t cmin = 1;
+Int_t cmax = 8;
+
+Float_t weightS = -1.;
+
+Int_t rebinsize = 4;
+
+Int_t parplotted = 2;
+
+Bool_t isMC = kFALSE; // don't change this (is set automatically)
+Bool_t selectTrue = kTRUE; // put it to true to remove background (only for MC)
+Bool_t keepTrue = kFALSE; // put it to false to fit only background (only for MC)
+
+Bool_t kGoodMatch = kFALSE; // to check good matching
+
+Bool_t kSigma2vs3 = kFALSE; // to check good matching
+
+Bool_t require5sigma = kFALSE; // don't touch this flag
+
+Bool_t bayesVsigma = kFALSE; // only to do checks
+
+Bool_t kTOFmatch = kFALSE; // for combined PID requires TOF matching
+
+Bool_t kOverAll = kFALSE;
+
+Bool_t kLoaded=kFALSE;
+LoadLib(){
+  weightS = -1.;
+
+  if(! kLoaded){
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libPhysics.so");
+    gSystem->Load("libTree.so");
+    gSystem->Load("libMinuit.so");
+    gSystem->Load("libSTEERBase.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libAOD.so");
+    gSystem->Load("libESD.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libANALYSISalice.so");
+    gSystem->Load("libCORRFW.so");
+    gSystem->Load("libNetx.so");
+    gSystem->Load("libPWGPPpid.so");
+
+    TFile *f = new TFile("AnalysisResults.root");
+    f->ls();
+    TList *l = (TList *) f->Get("contK0sBayes1");
+    l->ls();
+    fContPid1 = (AliPIDperfContainer *) l->At(0);
+    fContPid2 = (AliPIDperfContainer *) l->At(1);
+  }
+  kLoaded = kTRUE;
+
+  // check if MC
+  Float_t x[] = {xmin[0]+0.001,xmin[1]+0.001,xmin[2]+0.001,xmin[3]+0.001,xmin[4]+0.001,xmin[5]+0.001,xmin[6]+0.001,xmin[7]+0.001,1/*trueMC*/,xmin[9],xmin[10]};
+  Float_t x2[] = {xmax[0],xmax[1],xmax[2],xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],xmax[8],xmax[9],xmax[10]};
+
+  AliPIDperfContainer *tmp = fContPid1;
+  TH1D *h = tmp->GetQA(0, x, x2)->ProjectionX("checkMC");
+
+  if(h->GetEntries()) isMC = kTRUE;
+  else isMC=kFALSE;
+
+  if(!isMC){
+    selectTrue = kFALSE;
+    keepTrue = kTRUE;
+  }
+  else{
+    printf("MC truth found!!!!!!\nIt is MC!!!!!!");
+  }
+
+  fsign = new TF1("fsign","[0]*TMath::Voigt(x-[1],[3],[2])",fitmin,fitmax);
+  fback = new TF1("fback","pol1",fitmin,fitmax);
+  fall = new TF1("fall","[0]*TMath::Voigt(x-[1],[3],[2]) + pol1(4)",fitmin,fitmax);
+
+  fsign->SetLineColor(2);
+  fback->SetLineColor(4);
+
+  if(kSigma2vs3){
+    kGoodMatch=kFALSE;
+    kOverAll = 0;
+  }
+
+  if(bayesVsigma){
+    kOverAll = 0;
+    kGoodMatch=kFALSE;
+    kSigma2vs3=kFALSE;
+    kTOFmatch=kTRUE;
+    weightS = -0.7;
+  }
+}
+
+doeffPi(Int_t pos=1,Float_t prob=0.1,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+  LoadLib();
+
+  Int_t nptbin = binPid[2];
+  Float_t minptbin = xmin[2];
+  Float_t maxptbin = xmax[2];
+
+  if(pos == 0){
+    nptbin = binPid[3];
+    minptbin = xmin[3];
+    maxptbin = xmax[3];
+  }
+
+  if(prob > 0.1999){
+    kGoodMatch = kFALSE;
+    kSigma2vs3 = kFALSE;
+    if(! kOverAll) require5sigma = kTRUE;
+    if(!isMC) weightS = -0.95;
+  }
+
+  TCanvas *c = new TCanvas();
+  c->Divide((nptbin+1)/2,2);
+  TH2F *hh.*hh2;
+  TH1D *h,*h2;
+  char name[100];
+  Float_t b[50][3];
+
+  Double_t xx[50],yy[50];
+  Double_t exx[50],eyy[50];
+
+  for(Int_t i=0;i < nptbin;i++){
+    c->cd(i+1)->SetLogy();
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    xx[i] = (ptmin+ptmax)/2;
+    exx[i] = (-ptmin+ptmax)/2;
+
+    Float_t pp=0.1;
+    if(prob < 0.2) pp = 0.;
+    if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
+    sprintf(name,"TOF matched: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
+    hh->SetTitle(name);
+    sprintf(name,"hNoPid%i",i);
+    
+    pp=prob;
+    if(prob < 0.2) pp = 0.1;
+    if(pos) hh2=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh2=GetHistoPin(ptmin,ptmax,pp,0.0);
+    AddHisto(hh,hh2,weightS);
+
+    h = hh->ProjectionX(name,cmin,cmax);
+    h->RebinX(rebinsize);
+    h->Draw("ERR");
+    h->SetMarkerStyle(24);
+    b[i][0]=-1;
+    Int_t ntrial = 0;
+    Float_t chi2 = 10000;
+    while(ntrial < 10 && (chi2 > 20 + 1000*selectTrue)){
+      fit(h,b[i],"WW","",xx[i]);
+      c1->Update();
+//       getchar();
+      fit(h,b[i],"","",xx[i]);
+      ntrial++;
+      chi2 = b[i][2];
+      printf("chi2 = %f\n",chi2);
+      c1->Update();
+//       getchar();
+      
+    }
+
+    yy[i] = fall->GetParameter(parplotted);
+    eyy[i] = fall->GetParError(parplotted);
+  }
+
+  TGraphErrors *gpar = new TGraphErrors(nptbin,xx,yy,exx,eyy);
+  c->cd(8);
+  gpar->Draw("AP");
+  gpar->SetMarkerStyle(20);
+
+  TCanvas *c2 = new TCanvas();
+  c2->Divide((nptbin+1)/2,2);
+  Float_t b2[50][3];
+
+  for(Int_t i=0;i < nptbin;i++){
+    c2->cd(i+1);
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    Float_t pp=prob;
+    if(prob < 0.2) pp = 0.1;
+    if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
+    sprintf(name,"P_{TOF} > 0.8: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
+    hh->SetTitle(name);
+    sprintf(name,"hPid60_%i",i);
+    h = hh->ProjectionX(name,cmin,cmax);
+    h->RebinX(rebinsize);
+    h->Draw("ERR");
+    h->SetMarkerStyle(24);
+    b2[i][0]=-1;
+    Int_t ntrial = 0;
+    Float_t chi2 = 10000;
+    while(ntrial < 40 && (chi2 > 20 + 1000*selectTrue)){
+      fit(h,b2[i],"WW","");
+      fit(h,b2[i],"","");
+      ntrial++;
+      chi2 = b2[i][2];
+      printf("chi2 = %f\n",chi2);
+    }
+    yy[i] = fall->GetParameter(parplotted);
+    eyy[i] = fall->GetParError(parplotted);
+
+  }
+
+  TGraphErrors *gpar2 = new TGraphErrors(nptbin,xx,yy,exx,eyy);
+  c2->cd(8);
+  gpar2->Draw("AP");
+  gpar2->SetMarkerStyle(20);
+  
+  Double_t xpt[50],expt[50],eff[50],efferr[50];
+  for(Int_t i=0;i<nptbin;i++){
+    printf("%f +/- %f -  %f +/- %f\n",b[i][0],b[i][1],b2[i][0],b2[i][1]);
+
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    xpt[i] = (ptmin+ptmax)/2;
+    expt[i] = (-ptmin+ptmax)/2;
+    eff[i] = b2[i][0]/(b[i][0]-b2[i][0]*weightS);
+
+    b[i][0] = b[i][0]-b2[i][0]*weightS;
+
+    efferr[i] = TMath::Sqrt(b[i][1]*b[i][1]/b[i][0]/b[i][0] + b2[i][1]*b2[i][1]/b2[i][0]/b2[i][0])*(b2[i][0]+b2[i][1])*(1+weightS*(b2[i][0]-b2[i][1])/b[i][0])/b[i][0];//*(1-eff[i]);//der2*der2*(b[i][1]*b[i][1] - b2[i][1]*b2[i][1]));
+
+    if(TMath::Abs(efferr[i]) > 1)efferr[i]=1;
+  }
+  new TCanvas();
+  TGraphErrors *geff = new TGraphErrors(nptbin,xpt,eff,expt,efferr);
+  geff->Draw("AP");
+
+  char flag[100];
+  sprintf(flag,"");
+
+  if(isMC){
+    if(selectTrue) sprintf(flag,"true");
+    else if(!keepTrue) sprintf(flag,"back");
+  }
+
+  char flag2[100];
+  sprintf(flag2,"");
+
+  char etarange[100];
+  sprintf(etarange,"_%.1f-%.1f_",etaminkp,etamaxkp);
+
+  if(kGoodMatch)
+    sprintf(flag2,"GM");
+
+  if(bayesVsigma)
+    sprintf(flag2,"BayesVsSigma");
+
+  if(kSigma2vs3)
+    sprintf(flag2,"Sigma2vs3");
+
+  if(kOverAll)
+    sprintf(flag2,"OverAll");
+
+  if(pos){
+    if(prob >=0.2) sprintf(name,"pionPos%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
+    else sprintf(name,"pionPos%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
+  }
+  else{
+    if(prob >=0.2) sprintf(name,"pionNeg%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
+    else sprintf(name,"pionNeg%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
+  }
+
+  TFile *fout = new TFile(name,"RECREATE");
+  geff->Write();
+  fout->Close();
+
+
+  TF1 *ff = new TF1("ff","[0] - [1]*TMath::Exp([2]*x)",0,3);
+  ff->SetParameter(0,0.67);
+  ff->SetParameter(1,1.14383e+00);
+  ff->SetParameter(2,-2.29910);
+  ff->SetLineColor(4);
+  ff->SetLineColor(2);
+  ff->Draw("SAME");
+
+  TF1 *ff2 = new TF1("ff2","[0] - [1]*TMath::Exp([2]*x)",0,3);
+  ff2->SetParameter(0,0.67);
+  ff2->SetParameter(1,9.23126e-01);
+  ff2->SetParameter(2,-1.851);
+  ff2->SetLineColor(4);
+  ff2->Draw("SAME");
+}
+
+TH2F *GetHistoPip(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkp=0,Float_t pMinkn=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+
+  Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,pt+0.001,xmin[3]+0.001,pMinkp+0.001,pMinkn+0.001,(pMinkp>0.09)+0.001,kTOFmatch+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
+  Float_t x2[] = {xmax[0],etamaxkp-0.001,ptM-0.001,xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
+
+  if(kGoodMatch){
+    x[6] = 1.0001;
+    if(pMinkp > 0)
+      x2[9] = 4.9;
+      
+  }
+
+  if(kTOFmatch){
+    x[6] = 1.0001;
+  }
+
+  if(kSigma2vs3){
+    x[6] = 1.0001;
+    x2[9] = 3;
+    if(pMinkp > 0)
+      x2[9] = 2;
+  }
+
+  if(bayesVsigma){
+    if(pMinkp > 0){
+      x[4] = 0.2001;
+      x2[9] = 5;
+    }
+    else{
+      x2[9] = 3;
+    }
+
+    
+  }
+
+  if(require5sigma) x2[9] = 4.9;
+
+  AliPIDperfContainer *tmp = fContPid1;
+
+  TH2F *h = tmp->GetQA(0, x, x2);
+
+  h->GetXaxis()->SetTitle("M_{K^{0}_{s}} (GeV/#it{c}^{2})");
+  h->GetYaxis()->SetTitle("centrality [%]");
+
+  return h;
+}
+
+TH2F *GetHistoPin(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkn=0,Float_t pMinkp=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+
+  Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,xmin[2]+0.001,pt+0.001,pMinkp+0.001,pMinkn+0.001,kTOFmatch+0.001,(pMinkn>0.09)+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
+  Float_t x2[] = {xmax[0],etamaxkp-0.001,xmax[2],ptM-0.001,xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
+
+  if(kGoodMatch){
+    x[7] = 1.0001;
+    if(pMinkn > 0)
+      x2[10] = 4.9;
+      
+  }
+
+  if(kTOFmatch){
+    x[7] = 1.0001;
+  }
+
+  if(kSigma2vs3){
+    x[7] = 1.0001;
+    x2[10] = 3;
+    if(pMinkn > 0)
+      x2[10] = 2;
+  }
+ if(bayesVsigma){
+    if(pMinkn > 0){
+      x[5] = 0.2001;
+      x2[10] = 5;
+    }
+    else{
+      x2[10] = 3;
+    }    
+  }
+
+  if(require5sigma) x2[10] = 4.9;
+
+  AliPIDperfContainer *tmp = fContPid2;
+
+  TH2F *h = tmp->GetQA(0, x, x2);
+
+  h->GetXaxis()->SetTitle("M_{K^{0}_{s}} (GeV/#it{c}^{2})");
+  h->GetYaxis()->SetTitle("centrality [%]");
+
+  return h;
+}
+
+
+fit(TH1D *h,Float_t *a=NULL,char *opt="",char *opt2="",Float_t pt=1.5){
+  if(h->Integral(1,h->GetNbinsX()) < 1){
+    if(a){
+      a[0]=0.001;
+      a[1]=1;
+    }
+    return;
+  }
+  
+
+ fall->SetParameter(0,100);
+ fall->SetParameter(1,0.4971);
+ fall->SetParameter(2,0.00006);
+ fall->SetParameter(3,0.0035);
+
+ fall->SetParLimits(0,0.00001,10000);
+ fall->SetParLimits(1,0.4965,0.4985);
+ fall->SetParLimits(2,0.00005,0.001);
+ fall->SetParLimits(3,0.002,0.01);
+
+ fall->ReleaseParameter(4);
+ fall->ReleaseParameter(5);
+
+ if(selectTrue){
+   fall->FixParameter(4,0);
+   fall->FixParameter(5,0);
+ }
+
+ char name[100];
+ TF1 *ftmp=fall;
+
+ TF1 *ftmp2=new TF1(*fsign);
+ sprintf(name,"fsign%i",ifunc);
+ ftmp2->SetName(name);
+
+ TF1 *ftmp3=new TF1(*fback);
+ sprintf(name,"ftmp3%i",ifunc);
+ ftmp3->SetName(name);
+
+ ifunc++;
+
+ h->Fit(ftmp,opt,opt2,fitmin,fitmax);
+ h->Draw("ERR");
+
+ ftmp2->SetParameter(0,ftmp->GetParameter(0));
+ ftmp2->SetParameter(1,ftmp->GetParameter(1));
+ ftmp2->SetParameter(2,ftmp->GetParameter(2));
+ ftmp2->SetParameter(3,ftmp->GetParameter(3));
+ ftmp2->Draw("SAME");
+ ftmp3->SetParameter(0,ftmp->GetParameter(4));
+ ftmp3->SetParameter(1,ftmp->GetParameter(5));
+ ftmp3->Draw("SAME");
+
+ Float_t mean = ftmp->GetParameter(1);
+ Float_t sigma = 0.0044;
+
+ Float_t signI = ftmp2->Integral(mean-10*sigma,mean+10*sigma)/h->GetBinWidth(1);
+ Float_t backI = ftmp3->Integral(mean-3*sigma,mean+3*sigma)/h->GetBinWidth(1);
+
+ Float_t errI = TMath::Sqrt(ftmp->GetParError(0)*ftmp->GetParError(0)/(0.001+ftmp->GetParameter(0))/(0.001+ftmp->GetParameter(0)));
+
+ printf("signal(5 sigma) = %f +/- %f(fit) +/- %f(stat)\n",signI,errI*signI,TMath::Sqrt(signI));
+ printf("backgr(3sigma) = %f\n",backI);
+ printf("significance(3 sigma) = %f\n",signI/sqrt(signI+backI));
+
+ if(a){
+   a[0]=signI;
+   a[1]=signI*errI*signI*errI + signI;
+   a[1] = TMath::Sqrt(a[1]);
+   if(ftmp->GetNDF()) a[2] = ftmp->GetChisquare()/ftmp->GetNDF();
+
+
+   if(selectTrue){
+     a[0] = h->Integral(1,h->GetNbinsX());
+     a[1] = TMath::Sqrt(a[0]);
+   }
+ }
+}
+
+AddHisto(TH2F *h1,TH2F *h2,Float_t w){
+  Int_t nbinx = h1->GetNbinsX();
+  Int_t nbiny = h1->GetNbinsY();
+
+  for(Int_t i=1;i<=nbinx;i++){
+    for(Int_t j=1;j<=nbiny;j++){
+      Float_t val = h1->GetBinContent(i,j) + h2->GetBinContent(i,j)*w;
+      Float_t err = TMath::Min(TMath::Sqrt(val),val);
+      h1->SetBinContent(i,j,val);
+      h1->SetBinError(i,j,err);
+    }
+  }
+}
diff --git a/PWGPP/pid/doeffPr.C b/PWGPP/pid/doeffPr.C
new file mode 100644 (file)
index 0000000..1bca131
--- /dev/null
@@ -0,0 +1,500 @@
+TObject* fContPid1;
+TObject* fContPid2;
+const Int_t nBinPid = 14; // pt,eta, ptPip, ptPin, PPip, PPin, TOF3sigmaPip, TOF3sigmaPin, isPhiTrue, nsigmaPip, nsigmaPin
+// 0.985 < mass < 1.045 (60) and 0 < centrality < 100 (10)
+Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaPr*/,20/*pt+*/,20/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
+Float_t xmin[nBinPid] = {1,-0.8,0.3,0.3,0,0,-0.5,-0.5,-0.5,0,0,-TMath::Pi(),-TMath::Pi(),-TMath::Pi()/2};
+Float_t xmax[nBinPid] = {5,0.8,4.3,4.3,1,1,1.5,1.5,1.5,7.5,7.5,TMath::Pi(),TMath::Pi(),TMath::Pi()/2};
+
+TF1 *fsign;
+TF1 *fall;
+TF1 *fback;
+
+Int_t ifunc=0;
+
+Float_t fitmin = 1.08;
+Float_t fitmax = 1.155;
+
+Int_t cmin = 1;
+Int_t cmax = 8;
+
+Float_t weightS = -1.;
+
+Int_t rebinsize = 2;
+
+Int_t parplotted = 2;
+
+Bool_t isMC = kFALSE; // don't change this (is set automatically)
+Bool_t selectTrue = kTRUE; // put it to true to remove background (only for MC)
+Bool_t keepTrue = kFALSE; // put it to false to fit only background (only for MC)
+
+Bool_t kGoodMatch = kFALSE; // to check good matching
+
+Bool_t kSigma2vs3 = kFALSE; // to check good matching
+
+Bool_t require5sigma = kFALSE; // don't touch this flag
+
+Bool_t bayesVsigma = kFALSE; // only to do checks
+
+Bool_t kTOFmatch = kFALSE; // for combined PID requires TOF matching
+
+Bool_t kOverAll = kFALSE;
+
+Bool_t kLoaded=kFALSE;
+LoadLib(){
+  weightS = -1.;
+
+  if(! kLoaded){
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libPhysics.so");
+    gSystem->Load("libTree.so");
+    gSystem->Load("libMinuit.so");
+    gSystem->Load("libSTEERBase.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libAOD.so");
+    gSystem->Load("libESD.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libANALYSISalice.so");
+    gSystem->Load("libCORRFW.so");
+    gSystem->Load("libNetx.so");
+    gSystem->Load("libPWGPPpid.so");
+
+    TFile *f = new TFile("AnalysisResults.root");
+    f->ls();
+    TList *l = (TList *) f->Get("contLambdaBayes1");
+    l->ls();
+    fContPid1 = (AliPIDperfContainer *) l->At(0);
+    fContPid2 = (AliPIDperfContainer *) l->At(1);
+  }
+  kLoaded = kTRUE;
+
+  // check if MC
+  Float_t x[] = {xmin[0]+0.001,xmin[1]+0.001,xmin[2]+0.001,xmin[3]+0.001,xmin[4]+0.001,xmin[5]+0.001,xmin[6]+0.001,xmin[7]+0.001,1/*trueMC*/,xmin[9],xmin[10]};
+  Float_t x2[] = {xmax[0],xmax[1],xmax[2],xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],xmax[8],xmax[9],xmax[10]};
+
+  AliPIDperfContainer *tmp = fContPid1;
+  TH1D *h = tmp->GetQA(0, x, x2)->ProjectionX("checkMC");
+
+  if(h->GetEntries()) isMC = kTRUE;
+  else isMC=kFALSE;
+
+  if(!isMC){
+    selectTrue = kFALSE;
+    keepTrue = kTRUE;
+  }
+  else{
+    printf("MC truth found!!!!!!\nIt is MC!!!!!!");
+  }
+
+  fsign = new TF1("fsign","[0]*TMath::Voigt(x-[1],[3],[2])",fitmin,fitmax);
+  fback = new TF1("fback","pol2",fitmin,fitmax);
+  fall = new TF1("fall","[0]*TMath::Voigt(x-[1],[3],[2]) + pol2(4)",fitmin,fitmax);
+
+  fsign->SetLineColor(2);
+  fback->SetLineColor(4);
+
+  if(kSigma2vs3){
+    kGoodMatch=kFALSE;
+    kOverAll = 0;
+  }
+
+  if(bayesVsigma){
+    kOverAll = 0;
+    kGoodMatch=kFALSE;
+    kSigma2vs3=kFALSE;
+    kTOFmatch=kTRUE;
+    weightS = -0.7;
+  }
+}
+
+doeffPr(Int_t pos=1,Float_t prob=0.1,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+  LoadLib();
+
+  Int_t nptbin = binPid[2];
+  Float_t minptbin = xmin[2];
+  Float_t maxptbin = xmax[2];
+
+  if(pos == 0){
+    nptbin = binPid[3];
+    minptbin = xmin[3];
+    maxptbin = xmax[3];
+  }
+
+  if(prob > 0.1999){
+    kGoodMatch = kFALSE;
+    kSigma2vs3 = kFALSE;
+    if(! kOverAll) require5sigma = kTRUE;
+    if(!isMC) weightS = -0.95;
+  }
+
+  TCanvas *c = new TCanvas();
+  c->Divide((nptbin+1)/2,2);
+  TH2F *hh.*hh2;
+  TH1D *h,*h2;
+  char name[100];
+  Float_t b[50][3];
+
+  Double_t xx[50],yy[50];
+  Double_t exx[50],eyy[50];
+
+  for(Int_t i=0;i < nptbin;i++){
+    c->cd(i+1)->SetLogy();
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    xx[i] = (ptmin+ptmax)/2;
+    exx[i] = (-ptmin+ptmax)/2;
+
+    Float_t pp=0.1;
+    if(prob < 0.2) pp = 0.;
+    if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
+    sprintf(name,"TOF matched: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
+    hh->SetTitle(name);
+    sprintf(name,"hNoPid%i",i);
+    
+    pp=prob;
+    if(prob < 0.2) pp = 0.1;
+    if(pos) hh2=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh2=GetHistoPin(ptmin,ptmax,pp,0.0);
+    AddHisto(hh,hh2,weightS);
+
+    h = hh->ProjectionX(name,cmin,cmax);
+    h->RebinX(rebinsize);
+    h->Draw("ERR");
+    h->SetMarkerStyle(24);
+    b[i][0]=-1;
+    Int_t ntrial = 0;
+    Float_t chi2 = 10000;
+    while(ntrial < 10 && (chi2 > 20 + 1000*selectTrue)){
+      fit(h,b[i],"WW","",xx[i]);
+      c1->Update();
+//       getchar();
+      fit(h,b[i],"","",xx[i]);
+      ntrial++;
+      chi2 = b[i][2];
+      printf("chi2 = %f\n",chi2);
+      c1->Update();
+//       getchar();
+      
+    }
+
+    yy[i] = fall->GetParameter(parplotted);
+    eyy[i] = fall->GetParError(parplotted);
+  }
+
+  TGraphErrors *gpar = new TGraphErrors(nptbin,xx,yy,exx,eyy);
+  c->cd(8);
+  gpar->Draw("AP");
+  gpar->SetMarkerStyle(20);
+
+  TCanvas *c2 = new TCanvas();
+  c2->Divide((nptbin+1)/2,2);
+  Float_t b2[50][3];
+
+  for(Int_t i=0;i < nptbin;i++){
+    c2->cd(i+1);
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    Float_t pp=prob;
+    if(prob < 0.2) pp = 0.1;
+    if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
+    else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
+    sprintf(name,"P_{TOF} > 0.8: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
+    hh->SetTitle(name);
+    sprintf(name,"hPid60_%i",i);
+    h = hh->ProjectionX(name,cmin,cmax);
+    h->RebinX(rebinsize);
+    h->Draw("ERR");
+    h->SetMarkerStyle(24);
+    b2[i][0]=-1;
+    Int_t ntrial = 0;
+    Float_t chi2 = 10000;
+    while(ntrial < 40 && (chi2 > 20 + 1000*selectTrue)){
+      fit(h,b2[i],"WW","");
+      fit(h,b2[i],"","");
+      ntrial++;
+      chi2 = b2[i][2];
+      printf("chi2 = %f\n",chi2);
+    }
+    yy[i] = fall->GetParameter(parplotted);
+    eyy[i] = fall->GetParError(parplotted);
+
+  }
+
+  TGraphErrors *gpar2 = new TGraphErrors(nptbin,xx,yy,exx,eyy);
+  c2->cd(8);
+  gpar2->Draw("AP");
+  gpar2->SetMarkerStyle(20);
+  
+  Double_t xpt[50],expt[50],eff[50],efferr[50];
+  for(Int_t i=0;i<nptbin;i++){
+    printf("%f +/- %f -  %f +/- %f\n",b[i][0],b[i][1],b2[i][0],b2[i][1]);
+
+    Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
+    Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
+
+    xpt[i] = (ptmin+ptmax)/2;
+    expt[i] = (-ptmin+ptmax)/2;
+    eff[i] = b2[i][0]/(b[i][0]-b2[i][0]*weightS);
+
+    b[i][0] = b[i][0]-b2[i][0]*weightS;
+
+    efferr[i] = TMath::Sqrt(b[i][1]*b[i][1]/b[i][0]/b[i][0] + b2[i][1]*b2[i][1]/b2[i][0]/b2[i][0])*(b2[i][0]+b2[i][1])*(1+weightS*(b2[i][0]-b2[i][1])/b[i][0])/b[i][0];//*(1-eff[i]);//der2*der2*(b[i][1]*b[i][1] - b2[i][1]*b2[i][1]));
+
+    if(TMath::Abs(efferr[i]) > 1)efferr[i]=1;
+  }
+  new TCanvas();
+  TGraphErrors *geff = new TGraphErrors(nptbin,xpt,eff,expt,efferr);
+  geff->Draw("AP");
+
+  char flag[100];
+  sprintf(flag,"");
+
+  if(isMC){
+    if(selectTrue) sprintf(flag,"true");
+    else if(!keepTrue) sprintf(flag,"back");
+  }
+
+  char flag2[100];
+  sprintf(flag2,"");
+
+  char etarange[100];
+  sprintf(etarange,"_%.1f-%.1f_",etaminkp,etamaxkp);
+
+  if(kGoodMatch)
+    sprintf(flag2,"GM");
+
+  if(bayesVsigma)
+    sprintf(flag2,"BayesVsSigma");
+
+  if(kSigma2vs3)
+    sprintf(flag2,"Sigma2vs3");
+
+  if(kOverAll)
+    sprintf(flag2,"OverAll");
+
+  if(pos){
+    if(prob >=0.2) sprintf(name,"protonPos%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
+    else sprintf(name,"protonPos%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
+  }
+  else{
+    if(prob >=0.2) sprintf(name,"protonNeg%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
+    else sprintf(name,"protonNeg%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
+  }
+
+  TFile *fout = new TFile(name,"RECREATE");
+  geff->Write();
+  fout->Close();
+
+
+  TF1 *ff = new TF1("ff","[0] - [1]*TMath::Exp([2]*x)",0,3);
+  ff->SetParameter(0,0.67);
+  ff->SetParameter(1,1.14383e+00);
+  ff->SetParameter(2,-2.29910);
+  ff->SetLineColor(4);
+  ff->SetLineColor(2);
+  ff->Draw("SAME");
+
+  TF1 *ff2 = new TF1("ff2","[0] - [1]*TMath::Exp([2]*x)",0,3);
+  ff2->SetParameter(0,0.67);
+  ff2->SetParameter(1,9.23126e-01);
+  ff2->SetParameter(2,-1.851);
+  ff2->SetLineColor(4);
+  ff2->Draw("SAME");
+}
+
+TH2F *GetHistoPip(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkp=0,Float_t pMinkn=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+
+  Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,pt+0.001,xmin[3]+0.001,pMinkp+0.001,pMinkn+0.001,(pMinkp>0.09)+0.001,kTOFmatch+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
+  Float_t x2[] = {xmax[0],etamaxkp-0.001,ptM-0.001,xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
+
+  if(kGoodMatch){
+    x[6] = 1.0001;
+    if(pMinkp > 0)
+      x2[9] = 4.9;
+      
+  }
+
+  if(kTOFmatch){
+    x[6] = 1.0001;
+  }
+
+  if(kSigma2vs3){
+    x[6] = 1.0001;
+    x2[9] = 3;
+    if(pMinkp > 0)
+      x2[9] = 2;
+  }
+
+  if(bayesVsigma){
+    if(pMinkp > 0){
+      x[4] = 0.2001;
+      x2[9] = 5;
+    }
+    else{
+      x2[9] = 3;
+    }
+
+    
+  }
+
+  if(require5sigma) x2[9] = 4.9;
+
+  AliPIDperfContainer *tmp = fContPid1;
+
+  TH2F *h = tmp->GetQA(0, x, x2);
+
+  h->GetXaxis()->SetTitle("M_{#Lambda} (GeV/#it{c}^{2})");
+  h->GetYaxis()->SetTitle("centrality [%]");
+
+  return h;
+}
+
+TH2F *GetHistoPin(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkn=0,Float_t pMinkp=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
+
+  Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,xmin[2]+0.001,pt+0.001,pMinkp+0.001,pMinkn+0.001,kTOFmatch+0.001,(pMinkn>0.09)+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
+  Float_t x2[] = {xmax[0],etamaxkp-0.001,xmax[2],ptM-0.001,xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
+
+  if(kGoodMatch){
+    x[7] = 1.0001;
+    if(pMinkn > 0)
+      x2[10] = 4.9;
+      
+  }
+
+  if(kTOFmatch){
+    x[7] = 1.0001;
+  }
+
+  if(kSigma2vs3){
+    x[7] = 1.0001;
+    x2[10] = 3;
+    if(pMinkn > 0)
+      x2[10] = 2;
+  }
+ if(bayesVsigma){
+    if(pMinkn > 0){
+      x[5] = 0.2001;
+      x2[10] = 5;
+    }
+    else{
+      x2[10] = 3;
+    }    
+  }
+
+  if(require5sigma) x2[10] = 4.9;
+
+  AliPIDperfContainer *tmp = fContPid2;
+
+  TH2F *h = tmp->GetQA(0, x, x2);
+
+  h->GetXaxis()->SetTitle("M_{#Lambda} (GeV/#it{c}^{2})");
+  h->GetYaxis()->SetTitle("centrality [%]");
+
+  return h;
+}
+
+fit(TH1D *h,Float_t *a=NULL,char *opt="",char *opt2="",Float_t pt=1.5){
+  if(h->Integral(1,h->GetNbinsX()) < 1){
+    if(a){
+      a[0]=0.01;
+      a[1]=1;
+    }
+    return;
+  }
+
+
+ fall->SetParameter(0,100);
+ fall->SetParameter(1,1.115);
+ fall->SetParameter(2,0.0003);
+ fall->SetParameter(3,0.001);
+
+ fall->SetParLimits(0,0.00001,10000);
+ fall->SetParLimits(1,1.105,1.125);//1.01898 + 2.4e-04*pt-1e-03,1.01898 + 2.4e-04*pt+1e-03);
+ fall->SetParLimits(2,0.00001,0.001);
+ fall->SetParLimits(3,0.0005,0.01);
+
+ // fall->FixParameter(2,5E-4);
+
+//  fall->FixParameter(1,1.01884 + 2.9891e-04*pt);
+//  fall->FixParameter(2,0.0044);
+//  fall->FixParameter(3,7.57574e-04 + 3.85408e-04*pt);
+
+ fall->ReleaseParameter(4);
+ fall->ReleaseParameter(5);
+ fall->ReleaseParameter(6);
+
+ if(selectTrue){
+   fall->FixParameter(4,0);
+   fall->FixParameter(5,0);
+   fall->FixParameter(6,0);
+ }
+
+ char name[100];
+ TF1 *ftmp=fall;
+
+ TF1 *ftmp2=new TF1(*fsign);
+ sprintf(name,"fsign%i",ifunc);
+ ftmp2->SetName(name);
+
+ TF1 *ftmp3=new TF1(*fback);
+ sprintf(name,"ftmp3%i",ifunc);
+ ftmp3->SetName(name);
+
+ ifunc++;
+
+ h->Fit(ftmp,opt,opt2,fitmin,fitmax);
+ h->Draw("ERR");
+
+ ftmp2->SetParameter(0,ftmp->GetParameter(0));
+ ftmp2->SetParameter(1,ftmp->GetParameter(1));
+ ftmp2->SetParameter(2,ftmp->GetParameter(2));
+ ftmp2->SetParameter(3,ftmp->GetParameter(3));
+ ftmp2->Draw("SAME");
+ ftmp3->SetParameter(0,ftmp->GetParameter(4));
+ ftmp3->SetParameter(1,ftmp->GetParameter(5));
+ ftmp3->SetParameter(2,ftmp->GetParameter(6));
+ ftmp3->Draw("SAME");
+
+ Float_t mean = ftmp->GetParameter(1);
+ Float_t sigma = TMath::Abs(ftmp->GetParameter(3)) + TMath::Abs(ftmp->GetParameter(2));
+
+ Float_t signI = ftmp2->Integral(mean-10*sigma,mean+10*sigma)/h->GetBinWidth(1);
+ Float_t backI = ftmp3->Integral(mean-3*sigma,mean+3*sigma)/h->GetBinWidth(1);
+
+ Float_t errI = TMath::Sqrt(ftmp->GetParError(0)*ftmp->GetParError(0)/(0.001+ftmp->GetParameter(0))/(0.001+ftmp->GetParameter(0)));
+
+ printf("signal(5 sigma) = %f +/- %f(fit) +/- %f(stat)\n",signI,errI*signI,TMath::Sqrt(signI));
+ printf("backgr(3sigma) = %f\n",backI);
+ printf("significance(3 sigma) = %f\n",signI/sqrt(signI+backI));
+
+ if(a){
+   a[0]=signI;
+   a[1]=signI*errI*signI*errI + signI;
+   a[1] = TMath::Sqrt(a[1]);
+   if(ftmp->GetNDF()) a[2] = ftmp->GetChisquare()/ftmp->GetNDF();
+
+
+    if(selectTrue){
+      a[0] = h->Integral(1,h->GetNbinsX());
+      a[1] = TMath::Sqrt(a[0]);
+    }
+ }
+}
+
+AddHisto(TH2F *h1,TH2F *h2,Float_t w){
+  Int_t nbinx = h1->GetNbinsX();
+  Int_t nbiny = h1->GetNbinsY();
+
+  for(Int_t i=1;i<=nbinx;i++){
+    for(Int_t j=1;j<=nbiny;j++){
+      Float_t val = h1->GetBinContent(i,j) + h2->GetBinContent(i,j)*w;
+      Float_t err = TMath::Min(TMath::Sqrt(val),val);
+      h1->SetBinContent(i,j,val);
+      h1->SetBinError(i,j,err);
+    }
+  }
+}