]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx
Carlos Perez: fix typo
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskVnV0.cxx
index 129199be9738b67dd23c8f0229d8df656e9d42b2..cbe3c6fa21fce81c0055ce7e90f5722dde4db2ad 100644 (file)
 #include "AliAODMCHeader.h"
 #include "AliAODMCParticle.h"
 #include "TChain.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDVertex.h"
+#include "AliEventplane.h"
+#include "TProfile2D.h"
 
 // STL includes
 //#include <iostream>
@@ -33,6 +37,12 @@ Float_t AliAnalysisTaskVnV0::fgPsi2tpc=999.;
 Float_t AliAnalysisTaskVnV0::fgPsi3v0a=999.;
 Float_t AliAnalysisTaskVnV0::fgPsi3v0c=999.;
 Float_t AliAnalysisTaskVnV0::fgPsi3tpc=999.;
+Float_t AliAnalysisTaskVnV0::fgPsi2v0aMC=999.;
+Float_t AliAnalysisTaskVnV0::fgPsi2v0cMC=999.;
+Float_t AliAnalysisTaskVnV0::fgPsi2tpcMC=999.;
+Float_t AliAnalysisTaskVnV0::fgPsi3v0aMC=999.;
+Float_t AliAnalysisTaskVnV0::fgPsi3v0cMC=999.;
+Float_t AliAnalysisTaskVnV0::fgPsi3tpcMC=999.;
 
 //_____________________________________________________________________________
 AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
@@ -40,11 +50,15 @@ AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
   fEtaCut(0.8),   // cut on |eta| < fEtaCut
   fMinPt(0.15),   // cut on pt > fMinPt
+  fMinDistV0(0),
+  fMaxDistV0(100),
   fV2(kTRUE),
   fV3(kTRUE),
   fIsMC(kFALSE),
   fQAsw(kFALSE),
+  fIsAfter2011(kFALSE),
   fRun(-1),
+  fNcluster(70),
   fList(new TList()),
   fList2(new TList()),
   fList3(new TList()),
@@ -92,7 +106,19 @@ AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
   fContAllChargesMCCv3(NULL),
   fFillDCA(kFALSE),
   fContQApid(NULL),
-  fModulationDEDx(kFALSE)
+  fModulationDEDx(kFALSE),
+  fZvtx(0.),
+  fNK0s(0),
+  fNpiPos(0),
+  fNpiNeg(0),
+  fHKsPhi(NULL),
+  fHKsPhiEP(NULL),
+  fHK0sMass(NULL),
+  fHK0sMass2(NULL),
+  fHK0vsLambda(NULL),
+  fHctauPtEP(NULL),
+  fHctauAt1EP(NULL),
+  fCutsDaughter(NULL)
 {
   // Default constructor (should not be used)
   fList->SetName("resultsV2");
@@ -106,6 +132,13 @@ AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
   fList4->SetOwner(kTRUE); 
 
   fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
+
+  for(Int_t i=0;i < 1000;i++){
+    fPhiK0s[i] = 0.0;
+    fPtK0s[i] = 0.0;
+    fIPiPos[i] = 0;
+    fIPiNeg[i] = 0;
+  }
 }
 
 //______________________________________________________________________________
@@ -114,11 +147,15 @@ AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name):
   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
   fEtaCut(0.8),   // cut on |eta| < fEtaCut
   fMinPt(0.15),   // cut on pt > fMinPt
+  fMinDistV0(0),
+  fMaxDistV0(100),
   fV2(kTRUE),
   fV3(kTRUE),
   fIsMC(kFALSE),
   fQAsw(kFALSE),
+  fIsAfter2011(kFALSE),
   fRun(-1),
+  fNcluster(70),
   fList(new TList()),
   fList2(new TList()),
   fList3(new TList()),
@@ -166,7 +203,19 @@ AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name):
   fContAllChargesMCCv3(NULL),
   fFillDCA(kFALSE),
   fContQApid(NULL),
-  fModulationDEDx(kFALSE)
+  fModulationDEDx(kFALSE),
+  fZvtx(0.),
+  fNK0s(0),
+  fNpiPos(0),
+  fNpiNeg(0),
+  fHKsPhi(NULL),
+  fHKsPhiEP(NULL),
+  fHK0sMass(NULL),
+  fHK0sMass2(NULL),
+  fHK0vsLambda(NULL),
+  fHctauPtEP(NULL),
+  fHctauAt1EP(NULL),
+  fCutsDaughter(NULL)
 {
 
   DefineOutput(1, TList::Class());
@@ -186,6 +235,13 @@ AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name):
   fList4->SetOwner(kTRUE); 
 
   fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
+
+  for(Int_t i=0;i < 1000;i++){
+    fPhiK0s[i] = 0.0;
+    fPtK0s[i] = 0.0;
+    fIPiPos[i] = 0;
+    fIPiNeg[i] = 0;
+  }
 }
 //_____________________________________________________________________________
 AliAnalysisTaskVnV0::~AliAnalysisTaskVnV0()
@@ -214,10 +270,10 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
   // Container analyses (different steps mean different species)
   const Int_t nPtBinsTOF = 45;
   Double_t binsPtTOF[nPtBinsTOF+1] = {0., 0.05,  0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.25, 2.5, 2.75,3.0,3.25,3.5,3.75,4.0,4.5,5,5.5,6,6.5,7,8,9,10,12,15,20};
-  const Int_t nCentrTOF = 9;
+  const Int_t nCentrTOF = nCentrBin;
   const Int_t nPsiTOF = 10;  
   const Int_t nChargeBinsTOFres = 2; 
-  const Int_t nCentrTOFres = 9;
+  const Int_t nCentrTOFres = nCentrBin;
   const Int_t nProbTOFres = 4;
   const Int_t nPsiTOFres = 10;
   const Int_t nMaskPID = 3;
@@ -238,7 +294,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 
   // v2 container
   fContAllChargesV0A = new AliFlowVZEROResults("v2A",6,binsTOF);
-  fContAllChargesV0A->SetVarRange(0,-0.5,8.5); // centrality
+  fContAllChargesV0A->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fContAllChargesV0A->SetVarRange(1,-1.5,1.5);  // charge
   fContAllChargesV0A->SetVarRange(2,0.6,1.0001);// prob
   fContAllChargesV0A->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
@@ -259,9 +315,13 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
   if(fV2) fContAllChargesV0A->AddSpecies("t",nPtBinsTOF,binsPtTOF);
   if(fV2) fContAllChargesV0A->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
   if(fV2) fContAllChargesV0A->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0A->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0A->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0A->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0A->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 
   fContAllChargesV0C = new AliFlowVZEROResults("v2C",6,binsTOF);
-  fContAllChargesV0C->SetVarRange(0,-0.5,8.5); // centrality
+  fContAllChargesV0C->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fContAllChargesV0C->SetVarRange(1,-1.5,1.5);  // charge
   fContAllChargesV0C->SetVarRange(2,0.6,1.0001);// prob
   fContAllChargesV0C->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
@@ -282,13 +342,21 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
   if(fV2) fContAllChargesV0C->AddSpecies("t",nPtBinsTOF,binsPtTOF);
   if(fV2) fContAllChargesV0C->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
   if(fV2) fContAllChargesV0C->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0C->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0C->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0C->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
+  if(fV2) fContAllChargesV0C->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 
   fList->Add(fContAllChargesV0A);
   fList->Add(fContAllChargesV0C);
 
+  fHctauPtEP = new TProfile2D("hctauPtEP","K^{0}_{s} decay length;p_{T} (GeV/#it{c});#Delta#phi (rad)",40,0,5,10,-TMath::Pi(),TMath::Pi());
+  fHctauAt1EP = new TH2F("hctauAt1EP","K^{0}_{s} decay length at 1 GeV/#it{c};c#tau (cm);#Delta#phi (rad)",50,0,50,10,-TMath::Pi(),TMath::Pi());
+  // added at the end
+
   if(fIsMC && fV2){
     fContAllChargesMC = new AliFlowVZEROResults("v2mc",5,binsTOFmc);
-    fContAllChargesMC->SetVarRange(0,-0.5,8.5); // centrality
+    fContAllChargesMC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
     fContAllChargesMC->SetVarRange(1,-1.5,1.5);  // charge
     fContAllChargesMC->SetVarRange(2,0.6,1.0001);// prob
     fContAllChargesMC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
@@ -307,7 +375,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
     fList3->Add(fContAllChargesMC); 
 
     fContAllChargesMCA = new AliFlowVZEROResults("v2mcA",5,binsTOFmcPureMC);
-    fContAllChargesMCA->SetVarRange(0,-0.5,8.5); // centrality
+    fContAllChargesMCA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
     fContAllChargesMCA->SetVarRange(1,-1.5,1.5);  // charge
     fContAllChargesMCA->SetVarRange(2,0.6,1.0001);// prob
     fContAllChargesMCA->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
@@ -326,7 +394,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
     fList3->Add(fContAllChargesMCA); 
 
     fContAllChargesMCC = new AliFlowVZEROResults("v2mcC",5,binsTOFmcPureMC);
-    fContAllChargesMCC->SetVarRange(0,-0.5,8.5); // centrality
+    fContAllChargesMCC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
     fContAllChargesMCC->SetVarRange(1,-1.5,1.5);  // charge
     fContAllChargesMCC->SetVarRange(2,0.6,1.0001);// prob
     fContAllChargesMCC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
@@ -347,7 +415,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 
   // v3 container
   fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",6,binsTOF);
-  fContAllChargesV0Av3->SetVarRange(0,-0.5,8.5); // centrality
+  fContAllChargesV0Av3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fContAllChargesV0Av3->SetVarRange(1,-1.5,1.5);  // charge
   fContAllChargesV0Av3->SetVarRange(2,0.6,1.0001);// prob
   fContAllChargesV0Av3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
@@ -368,9 +436,13 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
   if(fV3) fContAllChargesV0Av3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
   if(fV3) fContAllChargesV0Av3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
   if(fV3) fContAllChargesV0Av3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Av3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Av3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Av3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Av3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 
   fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",6,binsTOF);
-  fContAllChargesV0Cv3->SetVarRange(0,-0.5,8.5); // centrality
+  fContAllChargesV0Cv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fContAllChargesV0Cv3->SetVarRange(1,-1.5,1.5);  // charge
   fContAllChargesV0Cv3->SetVarRange(2,0.6,1.0001);// prob
   fContAllChargesV0Cv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
@@ -391,13 +463,17 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
   if(fV3) fContAllChargesV0Cv3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
   if(fV3) fContAllChargesV0Cv3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
   if(fV3) fContAllChargesV0Cv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Cv3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Cv3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Cv3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
+  if(fV3) fContAllChargesV0Cv3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
 
   fList2->Add(fContAllChargesV0Av3);
   fList2->Add(fContAllChargesV0Cv3);
 
   if(fIsMC && fV3){
     fContAllChargesMCAv3 = new AliFlowVZEROResults("v3mcA",5,binsTOFmcPureMC);
-    fContAllChargesMCAv3->SetVarRange(0,-0.5,8.5); // centrality
+    fContAllChargesMCAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
     fContAllChargesMCAv3->SetVarRange(1,-1.5,1.5);  // charge
     fContAllChargesMCAv3->SetVarRange(2,0.6,1.0001);// prob
     fContAllChargesMCAv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
@@ -416,7 +492,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
     fList3->Add(fContAllChargesMCAv3); 
 
     fContAllChargesMCCv3 = new AliFlowVZEROResults("v3mcC",5,binsTOFmcPureMC);
-    fContAllChargesMCCv3->SetVarRange(0,-0.5,8.5); // centrality
+    fContAllChargesMCCv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
     fContAllChargesMCCv3->SetVarRange(1,-1.5,1.5);  // charge
     fContAllChargesMCCv3->SetVarRange(2,0.6,1.0001);// prob
     fContAllChargesMCCv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
@@ -437,18 +513,18 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 
   // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
   // v2
-  fHResTPCv0A2 = new TProfile("hResTPCv0A2","",9,0,9);
-  fHResTPCv0C2 = new TProfile("hResTPCv0C2","",9,0,9);
-  fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",9,0,9);
+  fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
+  fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
+  fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
 
   fList->Add(fHResTPCv0A2);
   fList->Add(fHResTPCv0C2);
   fList->Add(fHResv0Cv0A2);
 
   // v3
-  fHResTPCv0A3 = new TProfile("hResTPCv0A3","",9,0,9);
-  fHResTPCv0C3 = new TProfile("hResTPCv0C3","",9,0,9);
-  fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",9,0,9);
+  fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
+  fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
+  fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
 
   fList2->Add(fHResTPCv0A3);
   fList2->Add(fHResTPCv0C3);
@@ -456,17 +532,17 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 
   // MC as in the dataEP resolution (but using MC tracks)
   if(fIsMC && fV3){
-    fHResMA2 = new TProfile("hResMA2","",9,0,9);
-    fHResMC2 = new TProfile("hResMC2","",9,0,9);
-    fHResAC2 = new TProfile("hResAC2","",9,0,9);
+    fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
+    fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
+    fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
     fList3->Add(fHResMA2); 
     fList3->Add(fHResMC2); 
     fList3->Add(fHResAC2); 
   }
   if(fIsMC && fV3){
-    fHResMA3 = new TProfile("hResMA3","",9,0,9);
-    fHResMC3 = new TProfile("hResMC3","",9,0,9);
-    fHResAC3 = new TProfile("hResAC3","",9,0,9);
+    fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
+    fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
+    fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
     fList3->Add(fHResMA3); 
     fList3->Add(fHResMC3); 
     fList3->Add(fHResAC3); 
@@ -475,12 +551,12 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 
   // V0A and V0C event plane distributions
   //v2 
-  fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
-  fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
+  fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
+  fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
 
   //v3
-  fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
-  fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",9,0,9,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
+  fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
+  fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
 
   // QA container
   // v2
@@ -503,7 +579,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 
 
   fQA = new AliFlowVZEROQA("v2AQA",5,binsQA);
-  fQA->SetVarRange(0,-0.5,8.5); // centrality
+  fQA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fQA->SetVarRange(1,0,7);  // pt
   fQA->SetVarRange(2,0.,1.0001);// prob
   fQA->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
@@ -522,7 +598,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 //   fQA->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 
   fQA2 = new AliFlowVZEROQA("v2CQA",5,binsQA);
-  fQA2->SetVarRange(0,-0.5,8.5); // centrality
+  fQA2->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fQA2->SetVarRange(1,0,7);  // pt
   fQA2->SetVarRange(2,0.,1.0001);// prob
   fQA2->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
@@ -541,7 +617,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 //   fQA2->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 
   fQAv3 = new AliFlowVZEROQA("v3AQA",5,binsQAv3);
-  fQAv3->SetVarRange(0,-0.5,8.5); // centrality
+  fQAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fQAv3->SetVarRange(1,0,7);  // pt
   fQAv3->SetVarRange(2,0.,1.0001);// prob
   fQAv3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
@@ -560,7 +636,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
 //   fQAv3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
 
   fQA2v3 = new AliFlowVZEROQA("v3CQA",5,binsQAv3);
-  fQA2v3->SetVarRange(0,-0.5,8.5); // centrality
+  fQA2v3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fQA2v3->SetVarRange(1,0,7);  // pt
   fQA2v3->SetVarRange(2,0.,1.0001);// prob
   fQA2v3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
@@ -659,7 +735,7 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
     nsigmaQA[i] = -10 + 20.0*i/nbinsigma;
   }
   fContQApid = new AliFlowVZEROResults("qaPID",nBinQApid,binQApid);
-  fContQApid->SetVarRange(0,-0.5,8.5); // centrality
+  fContQApid->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
   fContQApid->SetVarRange(1,0,20);  // charge
   fContQApid->SetVarName(0,"centrality");
   fContQApid->SetVarName(1,"p_{t}");
@@ -673,6 +749,22 @@ void AliAnalysisTaskVnV0::UserCreateOutputObjects()
   
   printf("Output creation ok!!\n\n\n\n");
 
+  fList->Add(fHctauPtEP);
+  fList->Add(fHctauAt1EP);
+
+  fHKsPhi = new TH2D("hKsPhi","K^{0}_{s} #phi distributuion;v_{z} (cm);#phi (rad)",20,-10,10,20,0,2*TMath::Pi());
+  fList->Add(fHKsPhi);
+  fHKsPhiEP = new TH2D("hKsPhiEP","EP V0C #phi distributuion;v_{z} (cm);#phi (rad)",20,-10,10,20,0,2*TMath::Pi());
+  fList->Add(fHKsPhiEP);
+
+  fHK0sMass = new TH2D("hK0sMass","K^{0}_{s} mass;p_{T} (GeV/#it{c});mass (GeV/#it{c}^{2})",20,0,5,400,0,1);
+  fList->Add(fHK0sMass);
+  fHK0sMass2 = new TH2D("hK0sMass2","K^{0}_{s} mass using secondary vertex;p_{T} (GeV/#it{c});mass (GeV/#it{c}^{2})",20,0,5,400,0,1);
+  fList->Add(fHK0sMass2);
+
+  fHK0vsLambda= new TH2D("hK0vsLambda",";K^{0} mass;#Lambda mass",100,0,1,100,0.5,1.5);
+  fList->Add(fHK0vsLambda);
+
   // Post output data.
   if(fV2) PostData(1, fList);
   if(fV3) PostData(2, fList2);
@@ -694,6 +786,12 @@ void AliAnalysisTaskVnV0::UserExec(Option_t *)
     fgPsi3v0a=999.;
     fgPsi3v0c=999.;
     fgPsi3tpc=999.;
+    fgPsi2v0aMC=999.;
+    fgPsi2v0cMC=999.;
+    fgPsi2tpcMC=999.;
+    fgPsi3v0aMC=999.;
+    fgPsi3v0cMC=999.;
+    fgPsi3tpcMC=999.;
 
     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
     if(!fOutputAOD){
@@ -706,11 +804,11 @@ void AliAnalysisTaskVnV0::UserExec(Option_t *)
 
     if(run != fRun){
        // Load the calibrations run dependent
-       OpenInfoCalbration(run);
-       fRun=run;
+      if(! fIsAfter2011) OpenInfoCalbration(run);
+      fRun=run;
     }
 
-    Float_t zvtx = GetVertex(fOutputAOD);
+    fZvtx = GetVertex(fOutputAOD);
 
 
 
@@ -743,17 +841,19 @@ void AliAnalysisTaskVnV0::UserExec(Option_t *)
 
 */
 
-    if (TMath::Abs(zvtx) < fVtxCut) {
+    if (TMath::Abs(fZvtx) < fVtxCut) {
       //Centrality
       Float_t v0Centr  = -10.;
       Float_t trkCentr  = -10.;
       AliCentrality *centrality = fOutputAOD->GetCentrality();
       if (centrality){
-       v0Centr  = centrality->GetCentralityPercentile("V0M");
-       trkCentr = centrality->GetCentralityPercentile("TRK"); 
+//     printf("v0centr = %f -- tpccnetr%f\n",centrality->GetCentralityPercentile("V0M"),centrality->GetCentralityPercentile("TRK"));
+       v0Centr  = centrality->GetCentralityPercentile("V0M");
+       trkCentr = centrality->GetCentralityPercentile("TRK"); 
+       //changed
       }
 
-      if(TMath::Abs(v0Centr - trkCentr) < 5.0){ // consistency cut on centrality selection
+      if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr > 0){ // consistency cut on centrality selection
        fPID->SetDetResponse(fOutputAOD, v0Centr); // Set the PID object for each event!!!!
        Analyze(fOutputAOD,v0Centr); // Do analysis!!!!
 
@@ -767,6 +867,11 @@ void AliAnalysisTaskVnV0::UserExec(Option_t *)
 //________________________________________________________________________
 void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
 {      
+
+  Int_t nusedForK0s=0;
+  AliAODTrack *usedForK0s1[1000];
+  AliAODTrack *usedForK0s2[1000];
+
   Float_t mass[8] = {5.10998909999999971e-04, 1.05658000000000002e-01, 1.39570000000000000e-01, 4.93676999999999977e-01, 9.38271999999999995e-01,1.87783699999999998,2.81740199999999996,1.40805449999999999};
   
   // Event plane resolution for v2
@@ -774,8 +879,9 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
                         0.446480,0.612705,0.712222,0.736200,0.697907,0.610114,0.481009,0.327402,0.182277};// V0C vs. centrality
   
   Int_t iC = -1;    
-  if (v0Centr >0 && v0Centr < 80){ // analysis only for 0-80% centrality classes
-
+  if (v0Centr < 80){ // analysis only for 0-80% centrality classes
+    //  if (v0Centr >0 && v0Centr < 80){ // analysis only for 0-80% centrality classes
+    // changed
     fgIsPsiComputed = kTRUE;
 
     // centrality bins
@@ -788,7 +894,37 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
     else if(v0Centr < 60) iC = 6;
     else if(v0Centr < 70) iC = 7;
     else iC = 8;
-    
+
+    Int_t iCcal = iC;
+
+/*
+    if(nCentrBin==16){
+       iC = v0Centr/5;
+       if(iC >= nCentrBin) iC = nCentrBin-1;
+    }
+
+    // centrality bins
+    // changed
+    if(v0Centr < 10 + 10./9) iC = 0;
+    else if(v0Centr < 10 + 20./9) iC = 1;
+    else if(v0Centr < 10 + 30./9) iC = 2;
+    else if(v0Centr < 10 + 40./9) iC = 3;
+    else if(v0Centr < 10 + 50./9) iC = 4;
+    else if(v0Centr < 10 + 60./9) iC = 5;
+    else if(v0Centr < 10 + 70./9) iC = 6;
+    else if(v0Centr < 10 + 90./9) iC = 7;
+    else if(v0Centr < 10 + 100./9) iC = 8;
+    else if(v0Centr < 10 + 110./9) iC = 9;
+    else if(v0Centr < 10 + 120./9) iC = 10;
+    else if(v0Centr < 10 + 130./9) iC = 11;
+    else if(v0Centr < 10 + 140./9) iC = 12;
+    else if(v0Centr < 10 + 150./9) iC = 13;
+    else if(v0Centr < 10 + 160./9) iC = 14;
+    else if(v0Centr < 10 + 170./9) iC = 15;
+    else iC = 16;
+    if(iC >= nCentrBin) iC= nCentrBin - 1;
+*/
+
     //reset Q vector info      
     Double_t Qxa2 = 0, Qya2 = 0;
     Double_t Qxc2 = 0, Qyc2 = 0;
@@ -847,6 +983,9 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
            fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
            fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
            fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
+            fgPsi2v0aMC = EvPlaneMCV2[0];
+            fgPsi2v0cMC = EvPlaneMCV2[1];
+            fgPsi2tpcMC = EvPlaneMCV2[2];
          }
          if(fV3){
            EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
@@ -855,6 +994,9 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
            fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
            fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
            fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
+            fgPsi3v0aMC = EvPlaneMCV3[0];
+            fgPsi3v0cMC = EvPlaneMCV3[1];
+            fgPsi3tpcMC = EvPlaneMCV3[2];
          }
 
          // flow A and C side
@@ -924,43 +1066,101 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
       }
     }
 
+    // TPC EP needed for resolution studies (TPC subevent)
+    Double_t Qx2 = 0, Qy2 = 0;
+    Double_t Qx3 = 0, Qy3 = 0;
+
+    for(Int_t iT = 0; iT < nAODTracks; iT++) {
+      
+      AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
+      
+      if (!aodTrack){
+       continue;
+      }
+      
+      Bool_t trkFlag = aodTrack->TestFilterBit(1);
+
+      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(fOutputAOD->GetPrimaryVertex(), fOutputAOD->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());
+      
+    }
+    
+    evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
+    evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
+
+    fgPsi2tpc = evPlAng2;
+    fgPsi3tpc = evPlAng3;
+
+    SelectK0s();
+
     //V0 info    
     AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
 
     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
       Float_t multv0 = aodV0->GetMultiplicity(iv0);
-      if (iv0 < 32){ // V0C
-       Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
-       Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
-       Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
-       Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
-      } else {       // V0A
-       Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
-       Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
-       Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
-       Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+
+      if(! fIsAfter2011){
+       if(! fIsMC){
+         if (iv0 < 32){ // V0C
+           Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+           Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+           Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+           Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+         } else {       // V0A
+           Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+           Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+           Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+           Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+         }
+       }
+       else{
+         if (iv0 < 32){ // V0C
+           Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+           Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+           Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+           Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+         } else {       // V0A
+           Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+           Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+           Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+           Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+         }
+       }
       }
     }
 
     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
-    Double_t Qxamean2 = fMeanQ[iC][1][0];
-    Double_t Qxarms2  = fWidthQ[iC][1][0];
-    Double_t Qyamean2 = fMeanQ[iC][1][1];
-    Double_t Qyarms2  = fWidthQ[iC][1][1];
-    Double_t Qxamean3 = fMeanQv3[iC][1][0];
-    Double_t Qxarms3  = fWidthQv3[iC][1][0];
-    Double_t Qyamean3 = fMeanQv3[iC][1][1];
-    Double_t Qyarms3  = fWidthQv3[iC][1][1];
+    Double_t Qxamean2 = fMeanQ[iCcal][1][0];
+    Double_t Qxarms2  = fWidthQ[iCcal][1][0];
+    Double_t Qyamean2 = fMeanQ[iCcal][1][1];
+    Double_t Qyarms2  = fWidthQ[iCcal][1][1];
+    Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
+    Double_t Qxarms3  = fWidthQv3[iCcal][1][0];
+    Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
+    Double_t Qyarms3  = fWidthQv3[iCcal][1][1];
     
-    Double_t Qxcmean2 = fMeanQ[iC][0][0];
-    Double_t Qxcrms2  = fWidthQ[iC][0][0];
-    Double_t Qycmean2 = fMeanQ[iC][0][1];
-    Double_t Qycrms2  = fWidthQ[iC][0][1];     
-    Double_t Qxcmean3 = fMeanQv3[iC][0][0];
-    Double_t Qxcrms3  = fWidthQv3[iC][0][0];
-    Double_t Qycmean3 = fMeanQv3[iC][0][1];
-    Double_t Qycrms3  = fWidthQv3[iC][0][1];   
+    Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
+    Double_t Qxcrms2  = fWidthQ[iCcal][0][0];
+    Double_t Qycmean2 = fMeanQ[iCcal][0][1];
+    Double_t Qycrms2  = fWidthQ[iCcal][0][1];  
+    Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
+    Double_t Qxcrms3  = fWidthQv3[iCcal][0][0];
+    Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
+    Double_t Qycrms3  = fWidthQv3[iCcal][0][1];        
     
     Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
     Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
@@ -971,29 +1171,48 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
     Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
     Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
        
-    evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
-    evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
-    evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
-    evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
+    if(! fIsAfter2011){
+      if(! fIsMC){
+       evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
+       evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
+       evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
+       evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
+      }
+      else{
+       evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
+       evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
+       evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
+       evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
+      }
+    }
+    else{
+      AliEventplane *ep =  aodEvent->GetEventplane();
+      evPlAngV0ACor2 = ep->GetEventplane("V0A", aodEvent, 2);
+      evPlAngV0CCor2 = ep->GetEventplane("V0C", aodEvent, 2);
+      evPlAngV0ACor3 = ep->GetEventplane("V0A", aodEvent, 3);
+      evPlAngV0CCor3 = ep->GetEventplane("V0C", aodEvent, 3);
+    }
+
 
     fgPsi2v0a = evPlAngV0ACor2;
     fgPsi2v0c = evPlAngV0CCor2;
     fgPsi3v0a = evPlAngV0ACor3;
     fgPsi3v0c = evPlAngV0CCor3;
-                                
+
+    fHKsPhiEP->Fill(fZvtx,fgPsi2v0c);                           
     //loop track and get pid
     for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
       AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
        
       if (!aodTrack){
-       aodTrack->Delete();
        continue;
       }
       
       Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
-      if(fFillDCA) trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
+      if(fFillDCA) 
+       trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
 
-      if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < 70) || !trkFlag){
+      if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
        continue;
       }
 
@@ -1004,9 +1223,9 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
            
       if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
        continue;
-           
-      if(fFillDCA && TMath::Abs(b[0]) > 3.0 && TMath::Abs(b[1]) > 3)
-       continue;
+      
+      if(fFillDCA && (TMath::Abs(b[0]) > 3.0 || TMath::Abs(b[1]) > 3))
+       continue;
       
       // re-map the container in an array to do the analysis for V0A and V0C within a loop
       Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
@@ -1052,7 +1271,7 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
 
        Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
        Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
-           
+
        fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
        Float_t dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
        Float_t *probRead = fPID->GetProb();
@@ -1148,44 +1367,44 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
          if(TMath::Abs(nsigmaTPC[2])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[2])<5))){ //pi
            xQA[2] = prob[2];
            xQA3[2] = xQA[2];
-           if(fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
-           if(fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
+           if(fQAsw && fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
+           if(fQAsw && fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
          }
          if(TMath::Abs(nsigmaTPC[3])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[3])<5))){ //K
            xQA[2] = prob[3];
            xQA3[2] = xQA[2];
-           if(fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
-//         if(fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);    
+           if(fQAsw && fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
+//         if(fQAsw && fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);   
          }
          if(TMath::Abs(nsigmaTPC[4])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[4])<5))){//p
            xQA[2] = prob[4];
            xQA3[2] = xQA[2];
-           if(fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
-//         if(fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);    
+           if(fQAsw && fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
+//         if(fQAsw && fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);   
          }
          if(TMath::Abs(nsigmaTPC[0])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[0])<5))){//e
            xQA[2] = prob[0];
            xQA3[2] = xQA[2];
-//         if(fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
-//         if(fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);    
+//         if(fQAsw && fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
+//         if(fQAsw && fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);   
          }
          if(TMath::Abs(nsigmaTPC[5])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[5])<5))){//d
            xQA[2] = prob[5];
            xQA3[2] = xQA[2];
-           //    if(fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
-           //    if(fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);      
+           //    if(fQAsw && fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
+           //    if(fQAsw && fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);     
          }
          if(TMath::Abs(nsigmaTPC[6])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[6])<5))){//t
            xQA[2] = prob[6];
            xQA3[2] = xQA[2];
-           //    if(fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
-           //    if(fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);      
+           //    if(fQAsw && fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
+           //    if(fQAsw && fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);     
          }
          if(TMath::Abs(nsigmaTPC[7])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[7])<5))){//He3
            xQA[2] = prob[7];
            xQA3[2] = xQA[2];
-           //    if(fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
-           //    if(fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);      
+           //    if(fQAsw && fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
+           //    if(fQAsw && fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);     
          }
        }
 
@@ -1352,51 +1571,221 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
       } // end side loop
     } // end track loop
 
-    // Fill EP distribution histograms
-    if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
-    if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
-    
-    if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
-    if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
+    // my V0 loop
+    for(Int_t imy=0;imy<fNK0s;imy++){
+      Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
+      Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
+      
+      AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};        
+      AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
+     
+      for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
+       Float_t x[6] = {iC,-1/*my K0s are negative for convention*/,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
+       Float_t x3[6] = {iC,-1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
 
-    // TPC EP needed for resolution studies (TPC subevent)
-    Double_t Qx2 = 0, Qy2 = 0;
-    Double_t Qx3 = 0, Qy3 = 0;
+       Float_t v2V0 = TMath::Cos(2*(fPhiK0s[imy] - evPlAngV0[iV0]));
+       Float_t v3V0 = TMath::Cos(3*(fPhiK0s[imy] - evPlAngV0v3[iV0]));
+       if(fV2) contV0[iV0]->Fill(9,fPtK0s[imy],v2V0,x);
+       if(fV3) contV0v3[iV0]->Fill(9,fPtK0s[imy],v3V0,x3);
+      }
+    }
 
-    for(Int_t iT = 0; iT < nAODTracks; iT++) {
-      
-      AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
-      
-      if (!aodTrack){
-       aodTrack->Delete();
-       continue;
+    // V0 loop
+    Int_t nV0s = fOutputAOD->GetNumberOfV0s();
+    AliAODv0 *myV0;
+    Double_t dQT, dALPHA, dPT, 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()) > fEtaCut) continue; // skipping low momentum
+      Int_t pass = PassesAODCuts(myV0,fOutputAOD,0);
+      if(pass) {
+       dMASS = myV0->MassK0Short();
+       pass = 3;
+       fHK0sMass2->Fill(myV0->Pt(),dMASS);
+      }
+      if(TMath::Abs(dMASS-0.497)/0.005 > 3){
+       pass = PassesAODCuts(myV0,fOutputAOD,1);
+       if(pass) dMASS = myV0->MassLambda();
+       if(pass==2) dMASS = myV0->MassAntiLambda();
       }
-      
-      Bool_t trkFlag = aodTrack->TestFilterBit(1);
 
-      if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70)  || !trkFlag) 
-       continue;
+      if(pass){// 1 lambda, 2 antilambda, 3=K0s
+       dPT=myV0->Pt();
+       dQT=myV0->PtArmV0();
+       dALPHA=myV0->AlphaV0();
+
+       Int_t iPos, iNeg;
+       AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0);
+       if(iT->Charge()>0) {
+         iPos = 0; iNeg = 1;
+       } else {
+         iPos = 1; iNeg = 0;
+       }
+
+       // check if one of the daugthers was already used
+       if(pass == 3 && TMath::Abs(dMASS-0.497)/0.005 < 1){
+         fHKsPhi->Fill(fZvtx, myV0->Phi());
+       }
+
+       if(pass == 1000){ // disable
+         Bool_t used = kFALSE;
+         for(Int_t ii=0;ii<nusedForK0s;ii++){
+           if(myV0->GetDaughter(iNeg) == usedForK0s1[ii] || myV0->GetDaughter(iPos) == usedForK0s2[ii]){
+             used = kTRUE;
+           }
+         }
+         if((!used) && nusedForK0s < 1000){
+           nusedForK0s++;
+           usedForK0s1[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iNeg);
+           usedForK0s2[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iPos);
+           printf("accepted\n");
+         }
+         else{
+           dMASS = 0;
+           printf("rejected\n");
+         }
+       }
+
+       iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
+       AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
+
+       // re-map the container in an array to do the analysis for V0A and V0C within a loop
+       Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
+       AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
        
-      Double_t b[2] = {-99., -99.};
-      Double_t bCov[3] = {-99., -99., -99.};
-      if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
-       continue;
+       Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
+       AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
+
+       for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
+        
+         if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
+
+         Float_t v2V0 = TMath::Cos(2*(myV0->Phi() - evPlAngV0[iV0]));
+         Float_t v3V0 = TMath::Cos(3*(myV0->Phi() - evPlAngV0v3[iV0]));
+         
+         Float_t deltaphi = myV0->Phi()- evPlAngV0[iV0];
+         if(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
+         if(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
+
+         Float_t x[6] = {iC,1,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
+         Float_t x3[6] = {iC,1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
+         
+         Float_t decaylength = myV0->DecayLengthXY(fOutputAOD->GetPrimaryVertex());
+         //      printf("decay length = %f\n",decaylength);
+
+         if(pass==2){ // anti-lambda charge = -1
+           x[1] = -1;
+           x3[1] = -1;
+         }
+
+         if(decaylength < fMinDistV0) pass = 0;          
+         if(decaylength > fMaxDistV0) pass = 0;          
+
+         Float_t nsigma = 0;
+         if(pass < 3)
+           nsigma = TMath::Abs(dMASS-1.116)/0.0016;
+         else if(pass == 3)
+           nsigma = TMath::Abs(dMASS-0.497)/0.005;
+
+         if(nsigma < 1)
+           x[2] = 0.95;
+         else if(nsigma < 2)
+           x[2] = 0.85;
+         else if(nsigma < 3)
+           x[2] = 0.75;
+         else if(nsigma < 4)
+           x[2] = 0.65;
+         else
+           x[2] = 0.5;
+                   
+         x3[2] = x[2];
+
+         // Fill Container for lambda and Ks
+         if(fV2 && pass == 3 && x[2] > 0.6){
+           contV0[iV0]->Fill(9,myV0->Pt(),v2V0,x);
+           fHctauPtEP->Fill(myV0->Pt(),deltaphi,decaylength);//ciao
+           if(myV0->Pt() < 1.1 && myV0->Pt() > 0.9) fHctauAt1EP->Fill(decaylength,deltaphi);
+         }
+         if(fV3 && pass == 3 && x[2] > 0.6) contV0v3[iV0]->Fill(9,myV0->Pt(),v3V0,x3);
+         if(fV2 && pass < 3 && x[2] > 0.6) contV0[iV0]->Fill(10,myV0->Pt(),v2V0,x);
+         if(fV3 && pass < 3 && x[2] > 0.6) contV0v3[iV0]->Fill(10,myV0->Pt(),v3V0,x3);
+
+         if(pass < 3){ // lambda
+           AliAODTrack* aodTrack = iT;
+           if(pass==2) aodTrack=jT;
+
+           v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
+           v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
+
+           fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
+           Float_t *probRead = fPID->GetProb();
+           Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
+           Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis 
+
+            if(prob[4] < 0.61) prob[4] = 0.61;
            
-      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());
-      
-    }
-    
-    evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
-    evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
+           Float_t xdec[6] = {iC,aodTrack->Charge(),prob[4],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
+           Float_t xdec3[6] = {iC,aodTrack->Charge(),prob[4],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
+
+           // Fill Container for (anti)proton from lambda
+           if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
+             if(fV2) contV0[iV0]->Fill(11,aodTrack->Pt(),v2V0,xdec);
+             if(fV3) contV0v3[iV0]->Fill(11,aodTrack->Pt(),v3V0,xdec3);
+           }
+         }
+          else if(pass == 3){
+            AliAODTrack* aodTrack = iT;
+
+            v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
+            v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
+
+            fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
+            Float_t *probRead = fPID->GetProb();
+            Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
+            Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
+
+            if(prob[2] < 0.61) prob[2] = 0.61;
+
+            Float_t xdec[6] = {iC,aodTrack->Charge(),prob[2],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to 
+            Float_t xdec3[6] = {iC,aodTrack->Charge(),prob[2],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to 
+
+            if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
+              if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdec);
+              if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdec3);
+            }
+           
+            aodTrack = jT;
+            v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
+            v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
+
+            fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
+            Float_t *probRead2 = fPID->GetProb();
+            Float_t prob2[8] = {probRead2[0],probRead2[1],probRead2[2],probRead2[3],probRead2[4],probRead2[5],probRead2[6],probRead2[7]};
+            Float_t tofMismProb2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
+
+            if(prob2[2] < 0.61) prob2[2] = 0.61;
+
+            Float_t xdecB[6] = {iC,aodTrack->Charge(),prob2[2],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5,0}; // to
+            Float_t xdecB3[6] = {iC,aodTrack->Charge(),prob2[2],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5,0}; // to
+
+            if(nsigma < 2 && xdecB[2] > 0.6  && TMath::Abs(aodTrack->Eta()) < 0.8){
+              if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdecB);
+              if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdecB3);
+            }
+          }
+       }
+       
+      }
+    } // end loop on V0
 
-    fgPsi2tpc = evPlAng2;
-    fgPsi3tpc = evPlAng3;
+
+    // Fill EP distribution histograms
+    if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
+    if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
+    
+    if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
+    if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
 
     // Fill histograms needed for resolution evaluation
     if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
@@ -1408,6 +1797,13 @@ void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
     if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
   }
   
+
+
+  // clean track array
+  for(Int_t i=0;i < nusedForK0s;i++){
+    usedForK0s1[i] = NULL;
+    usedForK0s2[i] = NULL;
+  }
 }
 
 //_____________________________________________________________________________
@@ -1460,7 +1856,7 @@ void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){
 
     for(Int_t iside=0;iside<2;iside++){
        for(Int_t icoord=0;icoord<2;icoord++){
-           for(Int_t i=0;i  < nCentrBin;i++){
+           for(Int_t i=0;i  < 9;i++){
                char namecont[100];
                if(iside==0 && icoord==0)
                  snprintf(namecont,100,"hQxc2_%i",i);
@@ -1511,3 +1907,316 @@ void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){
        }
     }
 }
+//=======================================================================
+Int_t AliAnalysisTaskVnV0::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
+{
+  Int_t set = 0;
+  Float_t fV0Cuts[9];
+  // defines cuts to be used
+  // fV0Cuts[9] dl dca ctp d0 d0d0 qt minEta maxEta PID
+  switch(set) {
+  case(0): // No cuts
+    fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
+    fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
+    fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 0;
+    break;
+  case(1): // Tight cuts
+    fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
+    fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
+    fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 0;
+    break;
+  case(2): // Tight cuts + PID
+    fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
+    fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
+    fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 1;
+    break;
+  case(3): // No cuts + PID
+    fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
+    fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
+    fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 1;
+    break;
+  }
+
+  // daughter cuts
+  if(! fCutsDaughter){
+    fCutsDaughter = new AliESDtrackCuts(Form("daughter_cuts_%s","ESD") );
+    fCutsDaughter->SetPtRange(0.2,10.0);
+    fCutsDaughter->SetEtaRange(-0.8, 0.8 );
+    fCutsDaughter->SetMinNClustersTPC(80);
+    fCutsDaughter->SetMaxChi2PerClusterTPC(4.0);
+    fCutsDaughter->SetRequireTPCRefit(kTRUE);
+    fCutsDaughter->SetAcceptKinkDaughters(kFALSE);
+  }
+
+  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
+  AliESDtrack ieT( iT );
+  ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
+  ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
+  ieT.SetTPCPointsF( iT->GetTPCNclsF() );
+  ieT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
+  if (!fCutsDaughter->IsSelected( &ieT ) ) return 0;
+
+  jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
+  AliESDtrack jeT( jT );
+  jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
+  jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
+  jeT.SetTPCPointsF( jT->GetTPCNclsF() );
+  jeT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
+  if (!fCutsDaughter->IsSelected( &jeT ) ) 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 );
+  Double_t dDCA=myV0->DcaV0Daughters();
+  Double_t dCTP=myV0->CosPointingAngle( pvertex );
+  Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
+  Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
+  Double_t dD0D0=dD0P*dD0M;
+  Double_t dQT=myV0->PtArmV0();
+  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
+//   Double_t dPT=myV0->Pt();
+  Double_t dETA=myV0->Eta();
+  Int_t passes = 1;
+  if(dDL  <fV0Cuts[0]) passes = 0;
+  if(dDCA >fV0Cuts[1]) passes = 0;
+  if(dCTP <fV0Cuts[2]) passes = 0;
+  if(TMath::Abs(dD0P) <fV0Cuts[3]) passes = 0;
+  if(TMath::Abs(dD0M) <fV0Cuts[3]) passes = 0;
+  if(dD0D0>fV0Cuts[4]) passes = 0;
+  if(dETA <fV0Cuts[6]) passes = 0;
+  if(dETA >fV0Cuts[7]) passes = 0;
+  if(specie==0) if(dQT<fV0Cuts[5]) passes = 0;
+  if(specie==1&&passes==1&&dALPHA<0) passes = 2; // antilambda
+
+
+//   if(jT->Pt() < 0.5*myV0->Pt() || iT->Pt() < 0.5*myV0->Pt()) passes = 0;
+
+
+  // additional cut
+//   if(!(iT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
+//   if(!(jT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
+
+//   if(!(iT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
+//   if(!(jT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
+
+//  if(!(iT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
+//  if(!(jT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
+
+  Bool_t trkFlag = iT->TestFilterBit(1); // TPC only tracks (4,global track)
+  Bool_t trkFlag2 = jT->TestFilterBit(1); // TPC only tracks (4,global track)
+
+  if(!trkFlag) passes = 0;
+  if(!trkFlag2) passes = 0;
+
+  if(passes&&fV0Cuts[8]) {
+
+    Double_t dedxExp[8];
+    fPID->ComputeProb(iT,tAOD); // compute Bayesian probabilities
+    Float_t nsigmaTPC[8];
+
+    Int_t tofMatch=0;
+    Int_t tofMatch2=0;
+
+    if(iT->GetDetPid()){ // check the PID object is available
+      for(Int_t iS=0;iS < 8;iS++){
+       dedxExp[iS] = fPID->GetExpDeDx(iT,iS);
+       nsigmaTPC[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
+      }
+    }
+    else{
+      for(Int_t iS=0;iS < 8;iS++)
+       nsigmaTPC[iS] = 10;
+    }
+
+
+    if(fPID->GetCurrentMask(1)) // if TOF is present
+      tofMatch = 1;
+
+//     Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis 
+
+    Float_t *probRead = fPID->GetProb();
+    Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
+
+    fPID->ComputeProb(jT,tAOD); // compute Bayesian probabilities
+    Float_t nsigmaTPC2[8];
+    if(jT->GetDetPid()){ // check the PID object is available
+      for(Int_t iS=0;iS < 8;iS++){
+       dedxExp[iS] = fPID->GetExpDeDx(jT,iS);
+       nsigmaTPC2[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
+      }
+    }
+    else{
+      for(Int_t iS=0;iS < 8;iS++)
+       nsigmaTPC2[iS] = 10;
+    }
+
+    if(fPID->GetCurrentMask(1)) // if TOF is present
+      tofMatch2 = 1;
+
+//     Float_t tofMismProbMC2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis 
+
+    probRead = fPID->GetProb();
+    Float_t prob2[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
+    
+    if(jT->GetTPCNcls() < fNcluster) passes = 0;
+    else if(iT->GetTPCNcls() < fNcluster) passes = 0;
+
+//     if(! (tofMatch && tofMatch2)) passes = 0; 
+
+    /*
+    Float_t dMASS = myV0->MassK0Short();
+    Float_t nsigmaMass = TMath::Abs(dMASS-0.497)/0.005;
+    if(specie == 0 && TMath::Abs(nsigmaMass) < 1 && myV0->Pt() > 1) printf("candidate i=(pt=%f-phi=%f-tof=%i) j=(pt=%f-phi=%f-tof=%i) \n",iT->Pt(),iT->Phi(),tofMatch,jT->Pt(),jT->Phi(),tofMatch2);
+    */
+
+    switch(specie) {
+    case 0: // K0 PID
+      if(0){
+       if( ((jT->GetTPCmomentum()<15) &&
+            (TMath::Abs(nsigmaTPC2[2])>3.)) || prob2[2] < 0.9)
+         passes = 0;
+       if( ((iT->GetTPCmomentum()<15) &&
+            (TMath::Abs(nsigmaTPC[2])>3.))|| prob[2] < 0.9 )
+         passes = 0;
+      }
+      break;
+    case 1: // Lambda PID  i==pos j ==neg
+      if(passes==1) {
+       if( (iT->GetTPCmomentum()<15) &&
+           (TMath::Abs(nsigmaTPC[4])>3.) )
+         passes = 0;
+       if( (jT->GetTPCmomentum()<15) &&
+           (TMath::Abs(nsigmaTPC2[2])>3.) )
+         passes = 0;
+      }
+      if(passes==2) {
+       if( (iT->GetTPCmomentum()<15) &&
+           (TMath::Abs(nsigmaTPC[2])>3.) )
+         passes = 0;
+       if( (jT->GetTPCmomentum()<15) &&
+           (TMath::Abs(nsigmaTPC2[4])>3.) )
+         passes = 0;
+      }
+      break;
+    }
+  }
+  return passes;
+}
+//=======================================================================
+void AliAnalysisTaskVnV0::SelectK0s(){
+  fNK0s=0;
+  fNpiPos=0;
+  fNpiNeg=0;
+
+  if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAng2,1.0); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
+
+  // 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(1); // TPC only tracks
+//    trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
+
+    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.5/aodTrack->Pt()) continue;
+
+    fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
+    Float_t *probRead = fPID->GetProb();
+    Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
+ //    Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis 
+
+    Int_t charge = aodTrack->Charge();
+    if(prob[2] > 0.9){
+      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.2&& length > 1 && length < 25){
+       fHK0sMass->Fill(pt,mass);
+       
+       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);
+
+       fHK0vsLambda->Fill(mass,TMath::Min(massaL,massaAL));
+
+       if(TMath::Abs(mass-0.497)/0.005 < 1 && massaL > 1.15 && massaAL > 1.15){
+         fPhiK0s[fNK0s] = phi;
+         fPtK0s[fNK0s] = pt;
+         fNK0s++;
+       }
+      }
+    }
+  }
+}