Luke's 1st commit of Analysis task
authorlbarnby <lbarnby@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Jul 2011 09:12:19 +0000 (09:12 +0000)
committerlbarnby <lbarnby@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Jul 2011 09:12:19 +0000 (09:12 +0000)
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.cxx [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.h [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/runLukeV0.C [new file with mode: 0644]

diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.cxx b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.cxx
new file mode 100644 (file)
index 0000000..c1ed3f6
--- /dev/null
@@ -0,0 +1,1681 @@
+ /**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliAnalysisTaskLukeV0.cxx 46301 2011-01-06 14:25:27Z agheata $ */
+
+/* AliAnalysisTaskLukeV0.cxx
+ *
+ * Task analysing lambda, antilambda & K0 spectra
+ *
+ * Based on tutorial example from offline pages
+ * Edited by Arvinder Palaha
+ * Adapted by Luke Hanratty
+ *
+ */
+
+#include "AliAnalysisTaskLukeV0.h"
+
+#include "Riostream.h"
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TList.h"
+#include "TPDGCode.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliStack.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDEvent.h"
+#include "AliESDv0.h"
+#include "AliESDInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliMCEvent.h"
+#include "AliMCVertex.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+
+ClassImp(AliAnalysisTaskLukeV0)
+
+//________________________________________________________________________
+AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0() // All data members should be initialised here
+   :AliAnalysisTaskSE(),
+    fOutputList(0),
+    fTrackCuts(0),
+       fPIDResponse(0),
+    fHistPt(0), 
+    fHistEta(0),
+       fHistLuke(0),
+       fHistBetaV0(0), 
+       fHistCosPA(0), 
+       fHistCosTheta(0), 
+       fHistCosThetaPi2(0), 
+       fHistCosThetaProton2(0), 
+       fHistDCAV0Daughters(0), 
+       fHistDecayL(0), 
+       fHistDecayLxy(0), 
+       fHistDeltaTheta(0), 
+       fHistdNV0sdT2(0), 
+       fHistEtaNTracks(0),
+       fHistEtaPTracks(0),
+       fHistEtaV0s(0),
+       fHistImpactxyN(0), 
+       fHistImpactzN(0), 
+       fHistImpactxyP(0), 
+       fHistImpactzP(0), 
+       fHistKinkIndexFalse(0), 
+       fHistKinkIndexTrue(0),
+       fHistLambdaBgRapidity(0),
+       fHistLambdaBgEta(0),
+       fHistLambdaBgPt(0),     
+       fHistLambdaSigRapidity(0),
+       fHistLambdaSigEta(0),
+       fHistLambdaSigPt(0),
+       fHistMagneticField(0),
+       fHistMcLog(0),
+       fHistMcNLambdaPrimary(0),
+       fHistMcNLambda(0),
+       fHistMcNAntilambda(0),
+       fHistMcNKshort(0),
+       fHistMcPtV0La(0),
+       fHistMcPtV0Lb(0),
+       fHistMcPtV0K0(0),
+       fHistMcSigmaPProton(0),
+       fHistMcSigmaPNProton(0),
+       fHistMcSLambdaEta(0),
+       fHistMcSLambdaDaughter(0),
+       fHistMcSLambdaPt(0),
+       fHistMcTPCTrackLength(0),
+       fHistMK0(0), 
+       fHistMK00(0), 
+       fHistMK01(0), 
+       fHistMK02(0), 
+       fHistMK03(0), 
+       fHistMK04(0), 
+       fHistMK05(0), 
+       fHistMK06(0), 
+       fHistMK07(0), 
+       fHistMK08(0), 
+       fHistMK09(0),
+       fHistMLa(0), 
+       fHistMLa0(0), 
+       fHistMLa1(0), 
+       fHistMLa2(0), 
+       fHistMLa3(0), 
+       fHistMLa4(0), 
+       fHistMLa5(0), 
+       fHistMLa6(0), 
+       fHistMLa7(0), 
+       fHistMLa8(0), 
+       fHistMLa9(0), 
+       fHistMLb(0), 
+       fHistMLb0(0), 
+       fHistMLb1(0), 
+       fHistMLb2(0), 
+       fHistMLb3(0), 
+       fHistMLb4(0), 
+       fHistMLb5(0), 
+       fHistMLb6(0), 
+       fHistMLb7(0), 
+       fHistMLb8(0), 
+       fHistMLb9(0), 
+       fHistNLambda(0), 
+       fHistNV0(0), 
+       fHistNTracks(0), 
+       fHistPtV0(0), 
+       fHistPVZ(0), 
+       fHistTauLa(0), 
+       fHistTheta(0), 
+       fHistThetaPi2(0), 
+       fHistThetaProton2(0), 
+       fHistV0Z(0), 
+       fHistZ(0),
+       fHistBetaERatio(0), 
+       fHistBetaPRatio(0), 
+       fHistBetaPXRatio(0), 
+       fHistBetheBlochITSNeg(0), 
+       fHistBetheBlochITSPos(0), 
+       fHistBetheBlochTPCNeg(0), 
+       fHistBetheBlochTPCPos(0), 
+       fHistCosPAMLa(0), 
+       fHistDCAPtPSig(0),
+       fHistDCAPtPBg(0),
+       fHistDCAPtPbarSig(0),
+       fHistDCAPtPbarBg(0),
+       fHistDecayLDCA(0), 
+       fHistDecayLMLa(0), 
+       fHistDeltaThetaMLa(0), 
+       fHistImpactxyImpactz(0), 
+       fHistImpactxyMLa(0), 
+       fHistImpactzMLa(0), 
+       fHistMcLamdaPProductionVertex(0),
+       fHistMcLamdaSProductionVertex(0),
+       fHistMcLamdaSDecayVertex(0),
+       fHistMcPMK0Pt(0),
+       fHistMcPMLaPt(0),
+       fHistMcPMLbPt(0),
+       fHistMcSLambdaDaughterPairs(0),
+       fHistMcV0MK0Pt(0),
+       fHistMcV0MLaPt(0),
+       fHistMcV0MLbPt(0),
+       fHistMcV0LamdaSProductionVertex(0),
+       fHistMcV0IDLamdaSProductionVertex(0),
+       fHistMK0PtArm(0), 
+       fHistMK0Pt(0), 
+       fHistMK0R(0), 
+       fHistMLaPtArm(0), 
+       fHistMLaPt(0), 
+       fHistMLaR(0), 
+       fHistMLbPtArm(0), 
+       fHistMLbPt(0), 
+       fHistMLbR(0),
+       fHistNegChi2PerClusterMLa(0), 
+       fHistNegTPCRefitMLa(0), 
+       fHistNTPCNegClustersMLa(0), 
+       fHistNTPCPosClustersMLa(0), 
+       fHistNV0sNTracks(0),  
+       fHistPosChi2PerClusterMLa(0), 
+       fHistPosTPCRefitMLa(0), 
+       fHistPtArm(0), 
+       fHistPtArmR(0), 
+       fHistPtV0Z(0), 
+       fHistRZ(0), 
+       fHistTauLaMLa(0), 
+       fHistXZ(0), 
+       fHistYZ(0) // The last in the above list should not have a comma after it
+{
+    // Dummy constructor ALWAYS needed for I/O.
+}
+
+//________________________________________________________________________
+AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0(const char *name) // All data members should be initialised here
+   :AliAnalysisTaskSE(name),
+    fOutputList(0),
+       fTrackCuts(0),
+       fPIDResponse(0),
+       fHistPt(0), 
+       fHistEta(0),
+       fHistLuke(0),
+       fHistBetaV0(0), 
+       fHistCosPA(0), 
+       fHistCosTheta(0), 
+       fHistCosThetaPi2(0), 
+       fHistCosThetaProton2(0), 
+       fHistDCAV0Daughters(0), 
+       fHistDecayL(0), 
+       fHistDecayLxy(0), 
+       fHistDeltaTheta(0), 
+       fHistdNV0sdT2(0), 
+       fHistEtaNTracks(0),
+       fHistEtaPTracks(0),
+       fHistEtaV0s(0),
+       fHistImpactxyN(0), 
+       fHistImpactzN(0), 
+       fHistImpactxyP(0), 
+       fHistImpactzP(0), 
+       fHistKinkIndexFalse(0), 
+       fHistKinkIndexTrue(0),
+       fHistLambdaBgRapidity(0),
+       fHistLambdaBgEta(0),
+       fHistLambdaBgPt(0),     
+       fHistLambdaSigRapidity(0),
+       fHistLambdaSigEta(0),
+       fHistLambdaSigPt(0),
+       fHistMagneticField(0),
+       fHistMcLog(0),
+       fHistMcNLambdaPrimary(0),
+       fHistMcNLambda(0),
+       fHistMcNAntilambda(0),
+       fHistMcNKshort(0),
+       fHistMcPtV0La(0),
+       fHistMcPtV0Lb(0),
+       fHistMcPtV0K0(0),
+       fHistMcSigmaPProton(0),
+       fHistMcSigmaPNProton(0),
+       fHistMcSLambdaEta(0),
+       fHistMcSLambdaDaughter(0),
+       fHistMcSLambdaPt(0),
+       fHistMcTPCTrackLength(0),
+       fHistMK0(0), 
+       fHistMK00(0), 
+       fHistMK01(0), 
+       fHistMK02(0), 
+       fHistMK03(0), 
+       fHistMK04(0), 
+       fHistMK05(0), 
+       fHistMK06(0), 
+       fHistMK07(0), 
+       fHistMK08(0), 
+       fHistMK09(0),
+       fHistMLa(0), 
+       fHistMLa0(0), 
+       fHistMLa1(0), 
+       fHistMLa2(0), 
+       fHistMLa3(0), 
+       fHistMLa4(0), 
+       fHistMLa5(0), 
+       fHistMLa6(0), 
+       fHistMLa7(0), 
+       fHistMLa8(0), 
+       fHistMLa9(0), 
+       fHistMLb(0), 
+       fHistMLb0(0), 
+       fHistMLb1(0), 
+       fHistMLb2(0), 
+       fHistMLb3(0), 
+       fHistMLb4(0), 
+       fHistMLb5(0), 
+       fHistMLb6(0), 
+       fHistMLb7(0), 
+       fHistMLb8(0), 
+       fHistMLb9(0), 
+       fHistNLambda(0), 
+       fHistNV0(0), 
+       fHistNTracks(0), 
+       fHistPtV0(0), 
+       fHistPVZ(0), 
+       fHistTauLa(0), 
+       fHistTheta(0), 
+       fHistThetaPi2(0), 
+       fHistThetaProton2(0), 
+       fHistV0Z(0), 
+       fHistZ(0),
+       fHistBetaERatio(0), 
+       fHistBetaPRatio(0), 
+       fHistBetaPXRatio(0), 
+       fHistBetheBlochITSNeg(0), 
+       fHistBetheBlochITSPos(0), 
+       fHistBetheBlochTPCNeg(0), 
+       fHistBetheBlochTPCPos(0), 
+       fHistCosPAMLa(0), 
+       fHistDCAPtPSig(0),
+       fHistDCAPtPBg(0),
+       fHistDCAPtPbarSig(0),
+       fHistDCAPtPbarBg(0),
+       fHistDecayLDCA(0), 
+       fHistDecayLMLa(0), 
+       fHistDeltaThetaMLa(0), 
+       fHistImpactxyImpactz(0), 
+       fHistImpactxyMLa(0), 
+       fHistImpactzMLa(0), 
+       fHistMcLamdaPProductionVertex(0),
+       fHistMcLamdaSProductionVertex(0),
+       fHistMcLamdaSDecayVertex(0),
+       fHistMcPMK0Pt(0),
+       fHistMcPMLaPt(0),
+       fHistMcPMLbPt(0),
+       fHistMcSLambdaDaughterPairs(0),
+       fHistMcV0MK0Pt(0),
+       fHistMcV0MLaPt(0),
+       fHistMcV0MLbPt(0),
+       fHistMcV0LamdaSProductionVertex(0),
+       fHistMcV0IDLamdaSProductionVertex(0),
+       fHistMK0PtArm(0), 
+       fHistMK0Pt(0), 
+       fHistMK0R(0), 
+       fHistMLaPtArm(0), 
+       fHistMLaPt(0), 
+       fHistMLaR(0), 
+       fHistMLbPtArm(0), 
+       fHistMLbPt(0), 
+       fHistMLbR(0),
+       fHistNegChi2PerClusterMLa(0), 
+       fHistNegTPCRefitMLa(0), 
+       fHistNTPCNegClustersMLa(0), 
+       fHistNTPCPosClustersMLa(0), 
+       fHistNV0sNTracks(0), 
+       fHistPosChi2PerClusterMLa(0), 
+       fHistPosTPCRefitMLa(0), 
+       fHistPtArm(0), 
+       fHistPtArmR(0), 
+       fHistPtV0Z(0), 
+       fHistRZ(0), 
+       fHistTauLaMLa(0), 
+       fHistXZ(0), 
+       fHistYZ(0)// The last in the above list should not have a comma after it
+{
+    // Constructor
+    // Define input and output slots here (never in the dummy constructor)
+    // Input slot #0 works with a TChain - it is connected to the default input container
+    // Output slot #1 writes into a TH1 container
+    DefineOutput(1, TList::Class());                                            // for output list
+}
+
+//________________________________________________________________________
+AliAnalysisTaskLukeV0::~AliAnalysisTaskLukeV0()
+{
+    // Destructor. Clean-up the output list, but not the histograms that are put inside
+    // (the list is owner and will clean-up these histograms). Protect in PROOF case.
+    if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+        delete fOutputList;
+    }
+    if (fTrackCuts) delete fTrackCuts;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskLukeV0::UserCreateOutputObjects()
+{
+    // Create histograms
+    // Called once (on the worker node)
+        
+    fOutputList = new TList();
+    fOutputList->SetOwner();  // IMPORTANT!
+    
+       fTrackCuts = new AliESDtrackCuts();     
+       
+       AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+       AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+       fPIDResponse = inputHandler->GetPIDResponse();
+
+    // Create histograms - original test histograms
+    Int_t ptBins = 15;
+    Float_t ptLow = 0.1, ptUp = 3.1;
+    fHistPt = new TH1F("fHistPt", "P_{T} distribution for reconstructed", ptBins, ptLow, ptUp);        
+    Int_t etaBins = 40;
+    Float_t etaLow = -2.0, etaUp = 2.0;
+    fHistEta = new TH1F("fHistEta","#eta distribution for reconstructed",etaBins, etaLow, etaUp);
+       fHistLuke = new TH1F("fHistLuke","Lukes Histogram",100, 0, 10);
+
+       // lambda plot parameters
+       int div = 96;
+       float max = 1.2;
+       float min = 1.08;
+       
+       // Create remaining histograms
+       // TH1F first
+       fHistBetaV0 = new       TH1F("fHistBetaV0","Beta of the v0 candidate; Beta; NV0s",50,0,1);
+       fHistCosPA = new        TH1F("fHistCosPA", "Cosine of Pointing Angle of V0s; Cos PA; N(v0s)",202,0.8,1.01);
+       fHistCosTheta = new     TH1F("fHistCosTheta","Cos Theta of decay in Lambda rest frame ;Cos Theta ;N V0s",180,-1,1);
+       fHistCosThetaPi2 = new  TH1F("fHistCosThetaPi2","Cos Theta of pion decay in lab frame ;Cos Theta ;N V0s",180,-1,1);
+       fHistCosThetaProton2 = new      TH1F("fHistCosThetaProton2","Cos Theta proton in lab frame ;Cos Theta ;N V0s",180,-1,1);
+       fHistDCAV0Daughters = new       TH1F("fHistDCAV0Daughters", "DCA between V0 daughters; DCA (cm); N V0s", 100, 0, 2);
+       fHistDecayL = new       TH1F("fHistDecayL", "Distance between V0 and PV; Distance(cm); N(v0s)",200,-0.1,30);
+       fHistDecayLxy = new     TH1F("fHistDecayLxy", "Distance between V0 and PV in xy plane; Distance(cm); N(v0s)",200,-0.1,30);
+       fHistDeltaTheta = new   TH1F("fHistDeltaTheta","Difference in theta of Proton and Pi in lambda rest frame; Delta Theta; NV0s",180,-7,7);
+       fHistdNV0sdT2 = new     TH1F("fHistdNV0sdT2", "No of V0s over No of Tracks squared; NV0s/T^2; N", 30, 0, 2e-3);
+       fHistEtaNTracks = new TH1F("fHistEtaNTracks","Eta for negative tracks; Eta; N",100,-5,5);
+       fHistEtaPTracks = new TH1F("fHistEtaPTracks","Eta for positive tracks; Eta; N",100,-5,5);
+       fHistEtaV0s = new TH1F("fHistEtaV0s","Eta for v0s; Eta; N",100,-5,5);
+       fHistImpactxyN = new    TH1F("fHistImpactxyN", "RSM DCA between negative particle and primary vertex in xy plane; RSM DCA (cm); N(v0s)",100,0,1);
+       fHistImpactzN = new     TH1F("fHistImpactzN", "RSM DCA between negative particle and primary vertex in z direction; RSM DCA (cm); N(v0s)",100,0,1);
+       fHistImpactxyP = new    TH1F("fHistImpactxyP", "RSM DCA between positive particle and primary vertex in xy plane; RSM DCA (cm); N(v0s)",100,0,1);
+       fHistImpactzP = new     TH1F("fHistImpactzP", "RSM DCA between positive particle and primary vertex in z direction; RSM DCA (cm); N(v0s)",100,0,1);
+       fHistKinkIndexFalse = new       TH1F("fHistKinkIndexFalse","Lambda mass of non-kink candidates; M(p#pi^{-}) (GeV/c^{2})",96,1.08,1.2);
+       fHistKinkIndexTrue = new        TH1F("fHistKinkIndexTrue","Lambda mass of kink candidates; M(p#pi^{-}) (GeV/c^{2})",96,1.08,1.2);
+       fHistLambdaBgRapidity = new     TH1F("fHistLambdaBgRapidity","Rapidity of background V0s near Lambda mass; rapidity ",100,-1.5,1.5);
+       fHistLambdaBgEta = new  TH1F("fHistLambdaBgEta","Psuedorapidity of background V0s near Lambda mass; Eta ",100,-1.5,1.5);
+       fHistLambdaBgPt = new   TH1F("fHistLambdaBgPt","Transverse momentum of background V0s near Lambda mass; Pt(GeV/c) ",150,0,15);
+       fHistLambdaSigRapidity = new    TH1F("fHistLambdaSigRapidity","Rapidity of signal V0s on Lambda mass; rapidity ",100,-1.5,1.5);
+       fHistLambdaSigEta = new TH1F("fHistLambdaSigEta","Psuedorapidity of signal V0s on Lambda mass; Eta ",100,-1.5,1.5);
+       fHistLambdaSigPt = new  TH1F("fHistLambdaSigPt","Transverse momentum of signal V0s on Lambda mass; Pt(GeV/c) ",150,0,15);       
+       fHistMagneticField = new        TH1F("fHistMagneticField", "Magnetic Field; Magnetic Field; N",100,-100,100);
+       fHistMcLog = new        TH1F("fHistMcLog", "Log: 1=P ID, 2=P nID, 3=nP ID, 4=nP nID; Code; N",21,-0.25,10.25);
+       fHistMcNLambdaPrimary = new     TH1F("fHistMcNLambdaPrimary","Number of primary lambdas in MC; NLambdas; i",6,-0.25,2.25);
+       fHistMcNLambda = new    TH1F("fHistMcNLambda","Number of lambdas in MC; NLambdas; i",31,-0.5,30);
+       fHistMcNAntilambda = new        TH1F("fHistMcNAntilambda","Number of antilambdas in MC; NAntiLambdas; i",31,-0.5,30);
+       fHistMcNKshort = new    TH1F("fHistMcNKshort","Number of K0s in MC; NKshort; i",31,-0.5,30);
+       fHistMcPtV0La = new     TH1F("fHistMcPtV0La","Pt distribution of V0s confirmed as lambdas; Pt (GeV/c); dN/dydPt",200,0,10);
+       fHistMcPtV0Lb = new     TH1F("fHistMcPtV0Lb","Pt distribution of V0s confirmed as antilambdas; Pt (GeV/c); dN/dydPt",200,0,10);
+       fHistMcPtV0K0 = new     TH1F("fHistMcPtV0K0","Pt distribution of V0s confirmed as K0s; Pt (GeV/c); dN/dydPt",200,0,10);
+       fHistMcSigmaPProton = new       TH1F("fHistMcSigmaPProton","Sigma of being a proton from TPC response (MC protons); N Sigma; N",200,0,10);
+       fHistMcSigmaPNProton = new      TH1F("fHistMcSigmaPNProton","Sigma of being a proton from TPC response (MC not protons); N Sigma; N",200,0,10); 
+       fHistMcSLambdaEta = new TH1F("fHistMcSLambdaEta","#eta distribution for MC secondary lambdas",60, -3, 3);
+       fHistMcSLambdaDaughter = new TH1F("fHistMcSLambda_Daughter","PDG code for first daughter of MC secondary lambdas",8001, -4000.5,4000.5);
+       fHistMcSLambdaPt = new TH1F("fHistMcSLambdaPt","#Pt distribution for MC secondary lambdas",500, 0, 100);
+       fHistMcTPCTrackLength = new TH1F("fHistMcTPCTrackLength","TPC track length for charged Lambda daughters",500, 0, 100);
+       fHistMK0 = new  TH1F("fHistMK0","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK00 = new TH1F("fHistMK00","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK01 = new TH1F("fHistMK01","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK02 = new TH1F("fHistMK02","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK03 = new TH1F("fHistMK03","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK04 = new TH1F("fHistMK04","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK05 = new TH1F("fHistMK05","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK06 = new TH1F("fHistMK06","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK07 = new TH1F("fHistMK07","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK08 = new TH1F("fHistMK08","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMK09 = new TH1F("fHistMK09","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
+       fHistMLa = new  TH1F("fHistMLa","Lambda Mass; M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa0 = new TH1F("fHistMLa0", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa1 = new TH1F("fHistMLa1", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa2 = new TH1F("fHistMLa2", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa3 = new TH1F("fHistMLa3", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa4 = new TH1F("fHistMLa4", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa5 = new TH1F("fHistMLa5", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa6 = new TH1F("fHistMLa6", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa7 = new TH1F("fHistMLa7", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa8 = new TH1F("fHistMLa8", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLa9 = new TH1F("fHistMLa9", "V0 Mass; M(M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb = new  TH1F("fHistMLb","AntiLambda Mass; M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb0 = new TH1F("fHistMLb0", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb1 = new TH1F("fHistMLb1", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb2 = new TH1F("fHistMLb2", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb3 = new TH1F("fHistMLb3", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb4 = new TH1F("fHistMLb4", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb5 = new TH1F("fHistMLb5", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb6 = new TH1F("fHistMLb6", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb7 = new TH1F("fHistMLb7", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb8 = new TH1F("fHistMLb8", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistMLb9 = new TH1F("fHistMLb9", "V0 Mass; M(M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
+       fHistNLambda = new      TH1F("fHistNLambda", "Number of lambda per event; N(lambda); N(events)",50,-0.5,49.5);
+       fHistNV0 = new  TH1F("fHistNV0","V0 frequency distribution; Number of V0 Candidates",1000,0,100000);
+       fHistNTracks = new      TH1F("fHistNTracks", "Track frequency distribution; Number of Tracks; N Tracks", 1000, 0, 100000);
+       fHistPtV0 = new TH1F("fHistPtV0","V0 P_{T}; P_{T} (GeV/c);dN/dP_{T} (GeV/c)^{-1}",40,0.,4.);
+       fHistPVZ = new  TH1F("fHistPVZ","Z primary; Z (cm); Counts",100,-10,10);
+       fHistTauLa = new        TH1F("fHistTauLa", "Lifetime under Lambda mass hypothesis; Lifetime(s); N(v0s)",200,0,100);
+       fHistTheta = new        TH1F("fHistTheta","Angle of decay in Lambda rest frame ;Theta ;N V0s",180,-4,4);
+       fHistThetaPi2 = new     TH1F("fHistThetaPi2","Angle of pion decay in lab frame ;Theta ;N V0s",180,-4,4);
+       fHistThetaProton2 = new TH1F("fHistThetaProton2","Angle of proton decay in lab frame;Theta ;N V0s",180,-4,4);
+       fHistV0Z = new  TH1F("fHistV0Z","Z decay; Z (cm); Counts",100,-10,10);
+       fHistZ = new    TH1F("fHistZ","Z decay - Z primary; Z (cm); Counts",100,-10,10);
+       
+       //TH2F follow
+       fHistBetaERatio = new   TH2F("fHistBetaERatio","Ratio of Energies of Pion & Proton ;Beta ;Ratio",100,0,1,100,0,1);
+       fHistBetaPRatio = new   TH2F("fHistBetaPRatio","Ratio of Momentum of Pion & Proton ;Beta ;Ratio",100,0,1,100,0,2);
+       fHistBetaPXRatio = new  TH2F("fHistBetaPXRatio","Ratio of Momentum in boost direction of Pion & Proton ;Beta ;Ratio",100,0,1,100,-10,10);
+       fHistBetheBlochITSNeg = new     TH2F("fHistBetheBlochITSNeg","-dE/dX against Momentum for negative daughter from ITS; Log10 P (GeV); -dE/dx (keV/cm ?)",1000,-1,1,1000,0,200);
+       fHistBetheBlochITSPos = new     TH2F("fHistBetheBlochITSPos","-dE/dX against Momentum for positive daughter from ITS; Log10 P (GeV); -dE/dx (keV/cm ?)",1000,-1,1,1000,0,200);
+       fHistBetheBlochTPCNeg = new     TH2F("fHistBetheBlochTPCNeg","-dE/dX against Momentum for negative daughter from TPC; Log10 P (GeV); -dE/dx (keV/cm ?)",1000,-1,1,1000,0,200);
+       fHistBetheBlochTPCPos = new     TH2F("fHistBetheBlochTPCPos","-dE/dX against Momentum for positive daughter from TPC; Log10 P (GeV); -dE/dx (keV/cm ?)",1000,-1,1,1000,0,200);
+       fHistCosPAMLa = new     TH2F("fHistCosPAMLa", "Cosine of Pointing Angle of V0s; Cos PA; M(p#pi^{-}) (GeV/c^{2})",202,0.8,1.01,96,1.08,1.2);
+       fHistDCAPtPSig = new    TH2F("fHistDCAPtPSig","DCA of proton daughter from Lambda in peak region to PV versus Pt; P_{perp} (GeV/c); DCA(cm)",200,0,10,200,0,20);
+       fHistDCAPtPBg = new     TH2F("fHistDCAPtPBg","DCA of proton daughter from Lambda in sideband region to PV versus Pt; P_{perp} (GeV/c); DCA(cm)",200,0,10,200,0,20);
+       fHistDCAPtPbarSig = new TH2F("fHistDCAPtPbarSig","DCA of antiproton daughter from antiLambda in peak region to PV versus Pt; P_{perp} (GeV/c); DCA(cm)",200,0,10,200,0,20);
+       fHistDCAPtPbarBg = new  TH2F("fHistDCAPtPbarBg","DCA of antiproton daughter from antiLambda in sideband region to PV versus Pt; P_{perp} (GeV/c); DCA(cm)",200,0,10,200,0,20);
+       fHistDecayLDCA = new    TH2F("fHistDecayLDCA", "Distance between V0 and PV against DCA between daughter tracks; Distance(cm); DCA (cm) ",301,-0.1,30,100,0,2);
+       fHistDecayLMLa = new    TH2F("fHistDecayLMLa", "Distance between V0 and PV; Distance(cm); M(p#pi^{-}) (GeV/c^{2})",301,-0.1,30,96,1.08,1.2);
+       fHistDeltaThetaMLa = new        TH2F("fHistDeltaThetaMLa"," Delta theta against effective lambda mass; Delta Theta; MLa", 180, -3, 3, 96, 1.08, 1.2);
+       fHistImpactxyImpactz = new      TH2F("fHistImpactxyImpactz", "RSM DCA between negative particle and primary vertex in xy plane; RSM DCA xy (cm); RSM DCA z (cm)",100,0,1,100,0,10);
+       fHistImpactxyMLa = new  TH2F("fHistImpactxyMLa", "RSM DCA between negative particle and primary vertex in xy plane; RSM DCA (cm); M(p#pi^{-}) (GeV/c^{2})",100,0,1,96,1.08,1.2);
+       fHistImpactzMLa = new   TH2F("fHistImpactzMLa", "RSM DCA between negative particle and primary vertex in xy plane; RSM DCA (cm); M(p#pi^{-}) (GeV/c^{2})",100,0,1,96,1.08,1.2);
+       fHistMcLamdaPProductionVertex = new     TH2F("fHistMcLamdaPProductionVertex", "Production vertex of primary Lambdas in RZ plane; Zv; Rv",100,-100,100,100,0,100);
+       fHistMcLamdaSProductionVertex = new     TH2F("fHistMcLamdaSProductionVertex", "Production vertex of secondary Lambdas in RZ plane; Zv; Rv",100,-100,100,100,0,100);
+       fHistMcLamdaSDecayVertex = new  TH2F("fHistMcLamdaSDecayVertex", "Decay vertex of secondary Lambdas in RZ plane; Zv; Rv",100,-100,100,100,0,100);
+       fHistMcPMK0Pt = new     TH2F("fHistMcPMK0Pt","Monte Carlo primary K0 Mass versus Pt; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",200,0,10,140,0.414,0.582);
+       fHistMcPMLaPt = new     TH2F("fHistMcPMLaPt","Monte Carlo primary (& sigma0) Lambda Mass versus Pt; P_{perp} (GeV/c); M(p#pi^{-}) (GeV/c^2)",200,0,10,96,1.08,1.2);
+       fHistMcPMLbPt = new     TH2F("fHistMcPMLbPt","Monte Carlo primary (& sigma0) AntiLambda Mass versus Pt; P_{perp} (GeV/c); M(#bar{p}#pi^{+}) (GeV/c^2)",200,0,10,96,1.08,1.2);
+       fHistMcSLambdaDaughterPairs = new       TH2F("fHistMcSLambdaDaughterPairs", "PDG codes of daughter particles of secondary Lambdas; Daughter1; Daughter2",21,-10.5,10.5,21,-10.5,10.5);
+       fHistMcV0MK0Pt = new    TH2F("fHistMcV0MK0Pt","Monte Carlo V0s passing cuts. K0 Mass versus Pt; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",200,0,10,140,0.414,0.582);
+       fHistMcV0MLaPt = new    TH2F("fHistMcV0MLaPt","Monte Carlo V0s passing cuts. Lambda Mass versus Pt; P_{perp} (GeV/c); Lambda Mass (GeV/c^2)",200,0,10,96,1.08,1.2);
+       fHistMcV0MLbPt = new    TH2F("fHistMcV0MLbPt","Monte Carlo V0s passing cuts. Antilambda Mass versus Pt; P_{perp} (GeV/c); Antilambda Mass (GeV/c^2)",200,0,10,96,1.08,1.2);
+       fHistMcV0LamdaSProductionVertex = new   TH2F("fHistMcV0LamdaSProductionVertex", "Production vertex of V0 secondary Lambdas in RZ plane; Zv; Rv",100,-100,100,100,0,100);
+       fHistMcV0IDLamdaSProductionVertex = new TH2F("fHistMcV0IDLamdaSProductionVertex", "Production vertex of identified V0 secondary Lambdas in RZ plane; Zv; Rv",100,-100,100,100,0,100);   
+       fHistMK0PtArm = new     TH2F("fHistMK0PtArm","K0 Mass versus PtArm; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",40,0,0.25,140,0.414,0.582);
+       fHistMK0Pt = new        TH2F("fHistMK0Pt","K0 Mass versus Pt; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",200,0,10,140,0.414,0.582);
+       fHistMK0R = new TH2F("fHistMK0R","K0 Mass versus R; R (cm); K0 Mass (GeV/c^2)",120,0,120,140,0.414,0.582);
+       fHistMLaPtArm = new     TH2F("fHistMLaPtArm","Lambda Mass versus PtArm; P_{perp} (GeV/c); M(p#pi^{-}) (GeV/c^2)",40,0,0.25,96,1.08,1.2);
+       fHistMLaPt = new        TH2F("fHistMLaPt","Lambda Mass versus Pt; P_{perp} (GeV/c); M(p#pi^{-}) (GeV/c^2)",200,0,10,96,1.08,1.2);
+       fHistMLaR = new TH2F("fHistMLaR","Lambda Mass versus R; R (cm); M(p#pi^{-}) (GeV/c^2)",120,0,120,96,1.08,1.2);
+       fHistMLbPtArm = new     TH2F("fHistMLbPtArm","AntiLambda Mass versus PtArm; P_{perp} (GeV/c); M(#bar{p}#pi^{+}) (GeV/c^2)",40,0,0.25,96,1.08,1.2);
+       fHistMLbPt = new        TH2F("fHistMLbPt","AntiLambda Mass versus Pt; P_{perp} (GeV/c); M(#bar{p}#pi^{+}) (GeV/c^2)",200,0,10,96,1.08,1.2);
+       fHistMLbR = new TH2F("fHistMLbR","AntiLambda Mass versus R; R (cm); M(#bar{p}#pi^{+}) (GeV/c^2)",120,0,120,96,1.08,1.2);
+       fHistNegChi2PerClusterMLa = new TH2F("fHistNegChi2PerClusterMLa","Chi Squared per cluster against Lambda mass (negative Track);Chi^2 per cluster; M(p#pi^{-}) (GeV/c^{2})",200,0,7,96,1.08,1.2);; 
+       fHistNegTPCRefitMLa = new       TH2F("fHistNegTPCRefitMLa","TPC refit?; TPC refit?; M(p#pi^{-}) (GeV/c^{2})",200,0,100,96,1.08,1.2);; 
+       fHistNTPCNegClustersMLa = new   TH2F("fHistNTPCNegClustersMLa","Number of TPC clusters per negative track against Lambda mass;Number of TPC clusters; M(p#pi^{-}) (GeV/c^{2})",200,0,200,96,1.08,1.2);; 
+       fHistNTPCPosClustersMLa = new   TH2F("fHistNTPCPosClustersMLa","Number of TPC clusters per positive track against Lambda mass;Number of TPC clusters; M(p#pi^{-}) (GeV/c^{2})",200,0,200,96,1.08,1.2);; 
+       fHistNV0sNTracks = new  TH2F("fHistNV0sNTracks", "Number of tracks squared against number of v0s; NTracks^2; NV0s", 1000, 0, 1e8, 1000, 0, 1.5e5);
+       fHistPosChi2PerClusterMLa = new TH2F("fHistPosChi2PerClusterMLa","Chi Squared per cluster against Lambda mass (positive Track);Chi^2 per cluster; M(p#pi^{-}) (GeV/c^{2})",200,0,7,96,1.08,1.2);; 
+       fHistPosTPCRefitMLa = new       TH2F("fHistPosTPCRefitMLa","TPC refit?; TPC refit?; M(p#pi^{-}) (GeV/c^{2})",200,0,100,96,1.08,1.2);; 
+       fHistPtArm = new        TH2F("fHistPtArm","Podolanski-Armenteros Plot; #alpha; P_{perp} (GeV/c)",40,-1,1,80,0,0.5);
+       fHistPtArmR = new       TH2F("fHistPtArmR","PtArm versus R; R decay (cm); P_{perp}",120,0,120,50,0,0.25);
+       fHistPtV0Z = new        TH2F("fHistPtV0Z","Pt of V0 vs Z position; Pt (GeV/c); Z (cm)",200,-0.1,1.9,200,-50,50);
+       fHistRZ = new   TH2F("fHistRZ","R decay versus Z decay; Z (cm); R (cm)",100,-50,50,120,0,220);
+       fHistTauLaMLa = new     TH2F("fHistTauLaMLa", "Lifetime under Lambda mass hypothesis; Lifetime(s); M(p#pi^{-}) (GeV/c^{2})",200,0,100,96,1.08,1.2);
+       fHistXZ = new   TH2F("fHistXZ","X decay versus Z decay; Z (cm); X (cm)",100,-50,50,200,-200,200);
+       fHistYZ = new   TH2F("fHistYZ","Y decay versus Z decay; Z (cm); Y (cm)",100,-50,50,200,-200,200);  
+       
+       
+        
+       // All histograms must be added to output list
+        
+       fOutputList->Add(fHistPt); 
+    fOutputList->Add(fHistEta);
+       fOutputList->Add(fHistLuke);
+       fOutputList->Add(fHistBetaV0); 
+       fOutputList->Add(fHistCosPA); 
+       fOutputList->Add(fHistCosTheta); 
+       fOutputList->Add(fHistCosThetaPi2); 
+       fOutputList->Add(fHistCosThetaProton2); 
+       fOutputList->Add(fHistDCAV0Daughters); 
+       fOutputList->Add(fHistDecayL); 
+       fOutputList->Add(fHistDecayLxy); 
+       fOutputList->Add(fHistDeltaTheta); 
+       fOutputList->Add(fHistdNV0sdT2); 
+       fOutputList->Add(fHistEtaNTracks);
+       fOutputList->Add(fHistEtaPTracks);
+       fOutputList->Add(fHistEtaV0s);
+       fOutputList->Add(fHistImpactxyN); 
+       fOutputList->Add(fHistImpactzN); 
+       fOutputList->Add(fHistImpactxyP); 
+       fOutputList->Add(fHistImpactzP); 
+       fOutputList->Add(fHistKinkIndexFalse); 
+       fOutputList->Add(fHistKinkIndexTrue);
+       fOutputList->Add(fHistLambdaBgRapidity);
+       fOutputList->Add(fHistLambdaBgEta);
+       fOutputList->Add(fHistLambdaBgPt);      
+       fOutputList->Add(fHistLambdaSigRapidity);
+       fOutputList->Add(fHistLambdaSigEta);
+       fOutputList->Add(fHistLambdaSigPt);
+       fOutputList->Add(fHistMagneticField);
+       fOutputList->Add(fHistMcLog);
+       fOutputList->Add(fHistMcNLambdaPrimary);
+       fOutputList->Add(fHistMcNLambda);
+       fOutputList->Add(fHistMcNAntilambda);
+       fOutputList->Add(fHistMcNKshort);
+       fOutputList->Add(fHistMcPtV0La);
+       fOutputList->Add(fHistMcPtV0Lb);
+       fOutputList->Add(fHistMcPtV0K0);
+       fOutputList->Add(fHistMcSigmaPProton);
+       fOutputList->Add(fHistMcSigmaPNProton);
+       fOutputList->Add(fHistMcSLambdaEta);
+       fOutputList->Add(fHistMcSLambdaDaughter);
+       fOutputList->Add(fHistMcSLambdaPt);
+       fOutputList->Add(fHistMcTPCTrackLength);
+       fOutputList->Add(fHistMK0); 
+       fOutputList->Add(fHistMK00); 
+       fOutputList->Add(fHistMK01); 
+       fOutputList->Add(fHistMK02); 
+       fOutputList->Add(fHistMK03); 
+       fOutputList->Add(fHistMK04); 
+       fOutputList->Add(fHistMK05); 
+       fOutputList->Add(fHistMK06); 
+       fOutputList->Add(fHistMK07); 
+       fOutputList->Add(fHistMK08); 
+       fOutputList->Add(fHistMK09);
+       fOutputList->Add(fHistMLa); 
+       fOutputList->Add(fHistMLa0); 
+       fOutputList->Add(fHistMLa1); 
+       fOutputList->Add(fHistMLa2); 
+       fOutputList->Add(fHistMLa3); 
+       fOutputList->Add(fHistMLa4); 
+       fOutputList->Add(fHistMLa5); 
+       fOutputList->Add(fHistMLa6); 
+       fOutputList->Add(fHistMLa7); 
+       fOutputList->Add(fHistMLa8); 
+       fOutputList->Add(fHistMLa9); 
+       fOutputList->Add(fHistMLb); 
+       fOutputList->Add(fHistMLb0); 
+       fOutputList->Add(fHistMLb1); 
+       fOutputList->Add(fHistMLb2); 
+       fOutputList->Add(fHistMLb3); 
+       fOutputList->Add(fHistMLb4); 
+       fOutputList->Add(fHistMLb5); 
+       fOutputList->Add(fHistMLb6); 
+       fOutputList->Add(fHistMLb7); 
+       fOutputList->Add(fHistMLb8); 
+       fOutputList->Add(fHistMLb9); 
+       fOutputList->Add(fHistNLambda); 
+       fOutputList->Add(fHistNV0); 
+       fOutputList->Add(fHistNTracks); 
+       fOutputList->Add(fHistPtV0); 
+       fOutputList->Add(fHistPVZ); 
+       fOutputList->Add(fHistTauLa); 
+       fOutputList->Add(fHistTheta); 
+       fOutputList->Add(fHistThetaPi2); 
+       fOutputList->Add(fHistThetaProton2); 
+       fOutputList->Add(fHistV0Z); 
+       fOutputList->Add(fHistZ);
+       fOutputList->Add(fHistBetaERatio); 
+       fOutputList->Add(fHistBetaPRatio); 
+       fOutputList->Add(fHistBetaPXRatio); 
+       fOutputList->Add(fHistBetheBlochITSNeg); 
+       fOutputList->Add(fHistBetheBlochITSPos); 
+       fOutputList->Add(fHistBetheBlochTPCNeg); 
+       fOutputList->Add(fHistBetheBlochTPCPos); 
+       fOutputList->Add(fHistCosPAMLa); 
+       fOutputList->Add(fHistDCAPtPSig);
+       fOutputList->Add(fHistDCAPtPBg);
+       fOutputList->Add(fHistDCAPtPbarSig);
+       fOutputList->Add(fHistDCAPtPbarBg);
+       fOutputList->Add(fHistDecayLDCA); 
+       fOutputList->Add(fHistDecayLMLa); 
+       fOutputList->Add(fHistDeltaThetaMLa); 
+       fOutputList->Add(fHistImpactxyImpactz); 
+       fOutputList->Add(fHistImpactxyMLa); 
+       fOutputList->Add(fHistImpactzMLa); 
+       fOutputList->Add(fHistMcLamdaPProductionVertex);
+       fOutputList->Add(fHistMcLamdaSProductionVertex);
+       fOutputList->Add(fHistMcLamdaSDecayVertex);
+       fOutputList->Add(fHistMcPMK0Pt);
+       fOutputList->Add(fHistMcPMLaPt);
+       fOutputList->Add(fHistMcPMLbPt);
+       fOutputList->Add(fHistMcSLambdaDaughterPairs);
+       fOutputList->Add(fHistMcV0MK0Pt);
+       fOutputList->Add(fHistMcV0MLaPt);
+       fOutputList->Add(fHistMcV0MLbPt);
+       fOutputList->Add(fHistMcV0LamdaSProductionVertex);
+       fOutputList->Add(fHistMcV0IDLamdaSProductionVertex);
+       fOutputList->Add(fHistMK0PtArm); 
+       fOutputList->Add(fHistMK0Pt); 
+       fOutputList->Add(fHistMK0R); 
+       fOutputList->Add(fHistMLaPtArm); 
+       fOutputList->Add(fHistMLaPt); 
+       fOutputList->Add(fHistMLaR); 
+       fOutputList->Add(fHistMLbPtArm); 
+       fOutputList->Add(fHistMLbPt); 
+       fOutputList->Add(fHistMLbR);
+       fOutputList->Add(fHistNegChi2PerClusterMLa); 
+       fOutputList->Add(fHistNegTPCRefitMLa); 
+       fOutputList->Add(fHistNTPCNegClustersMLa); 
+       fOutputList->Add(fHistNTPCPosClustersMLa); 
+       fOutputList->Add(fHistNV0sNTracks); 
+       fOutputList->Add(fHistPosChi2PerClusterMLa); 
+       fOutputList->Add(fHistPosTPCRefitMLa); 
+       fOutputList->Add(fHistPtArm); 
+       fOutputList->Add(fHistPtArmR); 
+       fOutputList->Add(fHistPtV0Z); 
+       fOutputList->Add(fHistRZ); 
+       fOutputList->Add(fHistTauLaMLa); 
+       fOutputList->Add(fHistXZ); 
+       fOutputList->Add(fHistYZ);
+       
+       
+    PostData(1, fOutputList); // Post data for ALL output slots >0 here, to get at least an empty histogram
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskLukeV0::UserExec(Option_t *) 
+{
+    // Main loop
+    // Called for each event
+       
+       // paramaters used for most cuts, to minimise editing
+       double cutCosPa(0.998), cutcTau(2);
+       double cutNImpact(-999), cutDCA(0.4);
+       double cutBetheBloch(3);
+       double cutMinNClustersTPC(70), cutMaxChi2PerClusterTPC(-999);
+       double isMonteCarlo(false);
+       double cutEta(0.8);
+       
+       //Track Cuts set here
+       if(cutMinNClustersTPC != -999)
+       {(fTrackCuts->SetMinNClustersTPC(int(cutMinNClustersTPC)));}
+       if(cutMaxChi2PerClusterTPC != -999)
+       {fTrackCuts->SetMaxChi2PerClusterTPC(cutMaxChi2PerClusterTPC);}
+       fTrackCuts->SetAcceptKinkDaughters(kFALSE); 
+       fTrackCuts->SetRequireTPCRefit(kTRUE);
+       
+        
+    // Create pointer to reconstructed event
+
+       AliVEvent *event = InputEvent();
+       if (!event) { Printf("ERROR: Could not retrieve event"); return; }
+       
+       // create pointer to event
+       AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(event);
+       if (!fESD) {
+               AliError("Cannot get the ESD event");
+               return;
+       }
+
+       /*********************************************************************/
+       // MONTE CARLO SECTION
+       // This section loops over all MC tracks
+
+       int nLambdaMC = 0;
+       int nAntilambdaMC = 0;
+       int nKshortMC = 0;
+       
+       if(isMonteCarlo) 
+       {
+
+               // If the task accesses MC info, this can be done as in the commented block below:
+               
+               // Create pointer to reconstructed event
+               AliMCEvent *mcEvent = MCEvent();
+               if (!mcEvent) 
+                       { 
+                               Printf("ERROR: Could not retrieve MC event"); 
+                               //return; 
+                       }
+               else
+                       {
+                               Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
+                       }
+                        
+               // set up a stack for use in check for primary/stable particles
+               AliStack* mcStack = mcEvent->Stack();
+               if( !mcStack ) { Printf( "Stack not available"); return; }
+               
+               AliMCVertex *mcpv = (AliMCVertex *) mcEvent->GetPrimaryVertex();
+               Double_t mcpvPos[3];
+               if (mcpv != 0)
+               {
+                       mcpv->GetXYZ(mcpvPos);
+               }
+               else 
+               {
+                       Printf("ERROR: Could not resolve MC primary vertex");
+                       return;
+               }
+               
+               //loop over all MC tracks
+               for(Int_t iMCtrack = 0; iMCtrack < mcEvent->GetNumberOfTracks(); iMCtrack++)
+               {
+                       
+                       //booleans to check if track is La, Lb, K0 and primary
+                       bool lambdaMC = false;
+                       bool antilambdaMC = false;
+                       bool kshortMC = false;
+                       bool isprimaryMC = false;
+                       
+                       AliMCParticle *mcPart = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iMCtrack));
+                       
+                       if(mcPart->PdgCode() == kLambda0)
+                       {
+                               lambdaMC = true;
+                               nLambdaMC++;
+                       }
+                       else if(mcPart->PdgCode() == kK0Short)
+                       {
+                               kshortMC = true;
+                               nKshortMC++;
+                       }
+                       else if(mcPart->PdgCode() == kLambda0Bar)
+                       {
+                               antilambdaMC = true;
+                               nAntilambdaMC++;
+                       }
+                       
+                       //if only interested in La,Lb,K0, can use this to terminate loop for other particles
+                       //if(lambdaMC == false && antilambdaMC == false && kshortMC == false)
+                       //{continue;}
+                       
+                       Int_t firstDaughterLabel = mcPart->GetFirstDaughter();
+                       Int_t lastDaughterLabel = mcPart->GetLastDaughter();
+                       Int_t motherLabel = mcPart->GetMother();
+                       
+                       if(firstDaughterLabel < 0 || lastDaughterLabel < 0)
+                       {continue;}
+                       
+                       AliMCParticle *mcFirstDaughter = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(firstDaughterLabel));
+                       AliMCParticle *mcLastDaughter = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(lastDaughterLabel));
+                       AliMCParticle *mcMother = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(motherLabel));
+               
+                       double mcRadius = TMath::Sqrt((mcPart->Xv())*(mcPart->Xv())+(mcPart->Yv())*(mcPart->Yv()));
+                       double mcFirstDaughterRadius = TMath::Sqrt((mcFirstDaughter->Xv())*(mcFirstDaughter->Xv())+(mcFirstDaughter->Yv())*(mcFirstDaughter->Yv()));
+                       double mcLastDaughterRadius = TMath::Sqrt((mcLastDaughter->Xv())*(mcLastDaughter->Xv())+(mcLastDaughter->Yv())*(mcLastDaughter->Yv()));                 
+                       
+                       double motherType = -1;
+                       if(motherLabel >= 0)
+                       {motherType = mcMother->PdgCode();}
+                       
+                       // this block of code is used to include Sigma0 decays as primary lambda/antilambda
+                       bool sigma0MC = false;
+                       if(motherType == 3212 || motherType == -3212)
+                               {
+                               if(mcEvent->IsPhysicalPrimary(motherLabel))
+                                  {sigma0MC = true;}
+                               }
+                                  
+                       
+                       if(mcEvent->IsPhysicalPrimary(iMCtrack) || sigma0MC)
+                          {
+                                  isprimaryMC = true;
+                                  if(lambdaMC)
+                                  {
+                                          fHistMcNLambdaPrimary->Fill(1);
+                                          fHistMcLamdaPProductionVertex->Fill(mcPart->Zv(),mcRadius);
+                                          
+                                          if(TMath::Abs(mcPart->Y())<=0.5)
+                                          {fHistMcPMLaPt->Fill(mcPart->Pt(),mcPart->M());}
+                                  }
+                                  if(antilambdaMC)
+                                  {
+                                          if(TMath::Abs(mcPart->Y())<=0.5)
+                                          {fHistMcPMLbPt->Fill(mcPart->Pt(),mcPart->M());}
+                                  }
+                                  if(kshortMC)
+                                  {
+                                          if(TMath::Abs(mcPart->Y())<=0.5)
+                                          {fHistMcPMK0Pt->Fill(mcPart->Pt(),mcPart->M());}
+                                  }
+                               }
+                       else 
+                               {
+                                       isprimaryMC = false;
+                                       if(lambdaMC)
+                                       {
+                                               fHistMcNLambdaPrimary->Fill(2);
+                                               
+                                               fHistMcSLambdaDaughter->Fill(mcFirstDaughter->PdgCode());
+                                               
+                                               if(motherLabel >=0)
+                                               {
+                                                       if(mcMother->PdgCode() != kLambda0)
+                                                       {fHistMcLamdaSProductionVertex->Fill(mcPart->Zv(),mcRadius);}
+                                               }
+                                               else
+                                                       {fHistMcLamdaSProductionVertex->Fill(mcPart->Zv(),mcRadius);}
+                                               
+                                               if(mcFirstDaughter->PdgCode() == kProton || mcFirstDaughter->PdgCode() == kPiMinus)
+                                                       {fHistMcLamdaSDecayVertex->Fill(mcFirstDaughter->Zv(),mcFirstDaughterRadius);}
+                                                               
+                                               fHistMcSLambdaEta->Fill(mcPart->Eta());
+                                               fHistMcSLambdaPt->Fill(mcPart->Pt());
+                                       }
+                               }
+                       
+               }
+               
+
+       } 
+
+
+       fHistMcNLambda->Fill(nLambdaMC);
+       fHistMcNAntilambda->Fill(nAntilambdaMC);
+       fHistMcNKshort->Fill(nKshortMC);
+       
+       //END OF MONTE CARLO SECTION
+       /*********************************************************************/
+       
+                       
+       
+
+       
+    // Do some fast cuts first
+    // check for good reconstructed vertex
+    if(!(fESD->GetPrimaryVertex()->GetStatus())) return;
+    // if vertex is from spd vertexZ, require more stringent cut
+    if (fESD->GetPrimaryVertex()->IsFromVertexerZ()) {
+        if (fESD->GetPrimaryVertex()->GetDispersion()>0.02 ||  fESD->GetPrimaryVertex()->GetZRes()>0.25 ) return; // bad vertex from VertexerZ
+    }
+     
+       Double_t tV0Position[3];
+       Double_t tPVPosition[3];
+       Double_t radius;
+       
+       // physics selection
+       UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+       if(!maskIsSelected)
+    {
+               //printf("Event failed physics selection\n");
+               return;
+    }
+       
+       //if additional initial cuts wanted, can set conditions here
+       bool isCut = (fESD->GetNumberOfTracks()==0);
+       if (isCut)
+       {return;}
+       
+       //gets primary vertex for the event
+       const AliESDVertex *kPv = ((AliESDEvent *)fESD)->GetPrimaryVertex();
+       if ( kPv != 0 ) 
+    {
+               tPVPosition[0] = kPv->GetXv();
+               tPVPosition[1] = kPv->GetYv();
+               tPVPosition[2] = kPv->GetZv();
+               if( tPVPosition[2] == 0. ) 
+               {
+                       //printf("WARNING: Primary vertex a Z = 0, aborting\n");
+                       return;
+               }
+    }
+       else 
+    {
+               //printf("ERROR: Primary vertex not found\n");
+               return;
+    }
+       if( !kPv->GetStatus())
+    {return;}
+       
+       
+
+       
+       int nLambda(0);
+       int nV0s(0);
+       double lambdaLabel[fESD->GetNumberOfV0s()];
+       
+       // V0 loop - runs over every v0 in the event
+       for (Int_t iV0 = 0; iV0 < fESD->GetNumberOfV0s(); iV0++) 
+    {
+               
+               AliESDv0 *v0 = fESD->GetV0(iV0);
+               if (!v0) 
+               {
+                       printf("ERROR: Could not receive v0 %d\n", iV0);
+                       continue;
+               }
+               
+               bool lambdaCandidate = true;
+               bool antilambdaCandidate = true;
+               bool kshortCandidate = true;
+               
+               // keep only events of interest for fHistMLa plots
+               
+        if (v0->GetEffMass(4,2) < 1.08 || v0->GetEffMass(4,2) > 1.2 || TMath::Abs(v0->Y(3122))>0.5 )
+               {lambdaCandidate = false;}
+        if (v0->GetEffMass(2,4) < 1.08 || v0->GetEffMass(2,4) > 1.2 || TMath::Abs(v0->Y(-3122))>0.5)
+               {antilambdaCandidate = false;}
+        if (v0->GetEffMass(2,2) < 0.414 || v0->GetEffMass(2,2) > 0.582 || TMath::Abs(v0->Y(310))>0.5)
+               {kshortCandidate = false;}
+               if (v0->GetOnFlyStatus())
+               {continue;}
+               
+               if(!isMonteCarlo) 
+               {if(lambdaCandidate == false && antilambdaCandidate == false && kshortCandidate == false)
+               {continue;}}
+               
+               
+               //gets details of the v0
+               v0->GetXYZ(tV0Position[0],tV0Position[1],tV0Position[2]);
+               radius = TMath::Sqrt(tV0Position[0]*tV0Position[0]+tV0Position[1]*tV0Position[1]);
+               
+               double decayLength = (sqrt((tV0Position[0]-tPVPosition[0])*(tV0Position[0]-tPVPosition[0])+(tV0Position[1]-tPVPosition[1])*(tV0Position[1]-tPVPosition[1])+(tV0Position[2]-tPVPosition[2])*(tV0Position[2]-tPVPosition[2])));
+               double cTauLa = decayLength*(v0->GetEffMass(4,2))/(v0->P());
+               double cTauLb = decayLength*(v0->GetEffMass(2,4))/(v0->P());
+               double cTauK0 = decayLength*(v0->GetEffMass(2,2))/(v0->P());
+               
+               Int_t indexP, indexN;
+               indexP = TMath::Abs(v0->GetPindex());
+               AliESDtrack *posTrack = ((AliESDEvent*)fESD)->GetTrack(indexP);
+               indexN = TMath::Abs(v0->GetNindex());
+               AliESDtrack *negTrack = ((AliESDEvent*)fESD)->GetTrack(indexN);
+               
+               if(!posTrack || !negTrack)
+               {continue;}
+               
+               double pTrackMomentum[3];
+               double nTrackMomentum[3];
+               double pV0x, pV0y, pV0z;
+               posTrack->GetConstrainedPxPyPz(pTrackMomentum);
+               negTrack->GetConstrainedPxPyPz(nTrackMomentum);
+               v0->GetPxPyPz(pV0x, pV0y, pV0z);
+
+               const double kMLambda = 1.115;
+               const double kMProton = 0.938;
+               const double kMPi     = 0.140;
+               
+               double pPos2 = sqrt(pTrackMomentum[0]*pTrackMomentum[0]+pTrackMomentum[1]*pTrackMomentum[1]+pTrackMomentum[2]*pTrackMomentum[2]);
+               double pNeg2 = sqrt(nTrackMomentum[0]*nTrackMomentum[0]+nTrackMomentum[1]*nTrackMomentum[1]+nTrackMomentum[2]*nTrackMomentum[2]);
+               double pV02 = sqrt(pV0x*pV0x+pV0y*pV0y+pV0z*pV0z);
+               
+               //to prevent segfaults when ratios etc taken
+               if(pV02 < 0.01 || pPos2 <0.01 || pNeg2 <0.01)
+               {continue;}
+               
+               Float_t pImpactxy(0), pImpactz(0);
+               Float_t nImpactxy(0), nImpactz(0);
+               posTrack->GetImpactParameters(pImpactxy,pImpactz);
+               negTrack->GetImpactParameters(nImpactxy,nImpactz);
+               nImpactxy = sqrt((nImpactxy*nImpactxy));
+               nImpactz  = sqrt((nImpactz *nImpactz ));
+               pImpactxy = sqrt((pImpactxy*pImpactxy));
+               pImpactz  = sqrt((pImpactz *pImpactz ));
+               
+               /*********************************************************************/
+               // Cuts are implemented here.
+               
+               if(!(fTrackCuts->IsSelected(posTrack)) || !(fTrackCuts->IsSelected(negTrack)))
+               {
+                       lambdaCandidate = false;
+                       antilambdaCandidate = false;
+                       kshortCandidate = false;
+               }
+               
+               //extra cut to account for difference between p2 & p1 data
+               if(nImpactxy < 0.1 || pImpactxy < 0.1)
+               {
+                       lambdaCandidate = false;
+                       antilambdaCandidate = false;
+                       kshortCandidate = false;
+               }
+               
+               //psuedorapidity cut
+               if(cutEta != -999)
+               {
+                       if(TMath::Abs(posTrack->Eta()) > cutEta || TMath::Abs(negTrack->Eta())  >cutEta)
+                       {
+                               lambdaCandidate = false;
+                               antilambdaCandidate = false;
+                               kshortCandidate = false;
+                       }
+               }
+               
+               //pointing angle cut
+               if(cutCosPa != -999)
+               {
+                       if (v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]) < cutCosPa)
+                       {
+                               lambdaCandidate = false;
+                               antilambdaCandidate = false;
+                               kshortCandidate = false;
+                       }
+               }
+               
+               //lifetime cut
+               if(cutcTau != -999)
+               {
+                       if(cTauLa < cutcTau)
+                       {
+                               lambdaCandidate = false;
+                       }
+                       if(cTauLb < cutcTau)
+                       {
+                               antilambdaCandidate = false;
+                       }
+                       if(cTauK0 < cutcTau)
+                       {
+                               kshortCandidate = false;
+                       }
+                       
+               }
+               
+               // Impact paramater cut (on neg particle)
+               if(cutNImpact != -999)
+               {
+                       if(nImpactxy < cutNImpact || nImpactz < cutNImpact)
+                       {
+                               lambdaCandidate = false;
+                       }
+                       if(pImpactxy < cutNImpact || pImpactz < cutNImpact)
+                       {
+                               antilambdaCandidate = false;
+                       }
+               }
+               
+
+               // DCA between daughterscut
+               if(cutDCA != -999)
+               {
+                       if(v0->GetDcaV0Daughters() > cutDCA)
+                       {
+                               lambdaCandidate = false;
+                               antilambdaCandidate = false;
+                               kshortCandidate = false;
+                       }
+               }
+               
+               // Bethe Bloch cut. Made sightly complicated as options for crude cuts still included. Should probably reduce to just 'official' cuts
+               if(cutBetheBloch != -999)
+               { 
+                       if(posTrack->GetTPCsignal() <0 || negTrack->GetTPCsignal()<0)
+                       {continue;}
+                       
+                       if(lambdaCandidate)
+                       {
+                               if(cutBetheBloch > 0)
+                               {
+                               if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kProton)) > cutBetheBloch )
+                               {lambdaCandidate = false;}
+                               if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(negTrack, AliPID::kPion)) > cutBetheBloch )
+                               {lambdaCandidate = false;}
+                               }
+                               
+                               if(cutBetheBloch == -4)  
+                               {                               
+                               if(isMonteCarlo) 
+                                       {
+                                       double beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+0.9*0.9))),2);
+                                       double gamma2 = 1.0/(1.0-beta2);
+                                       if(posTrack->GetTPCsignal() < (2.0/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
+                                       {lambdaCandidate = false;}
+                                       }
+                               else 
+                                       { 
+                                       double beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+kMProton*kMProton))),2);
+                                       double gamma2 = 1.0/(1.0-beta2);
+                                       if(posTrack->GetTPCsignal() < (2.3/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
+                                       {lambdaCandidate = false;}
+                                       }
+                               }
+                               
+                       }
+                       
+                       if(antilambdaCandidate)
+                       {
+                               if(cutBetheBloch > 0)
+                               {
+                                       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(negTrack, AliPID::kProton)) > cutBetheBloch )
+                                       {antilambdaCandidate = false;}
+                                       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kPion)) > cutBetheBloch )
+                                       {antilambdaCandidate = false;}
+                               }
+                               
+                               if(cutBetheBloch == -4)  
+                               { 
+                                       if(isMonteCarlo) 
+                                       {
+                                               double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+0.9*0.9))),2);
+                                               double gamma2 = 1.0/(1.0-beta2);
+                                               if(negTrack->GetTPCsignal() < (2.0/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
+                                               {antilambdaCandidate = false;}
+                                       }
+                                       else 
+                                       { 
+                                               double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+0.9*0.9))),2);
+                                               double gamma2 = 1.0/(1.0-beta2);
+                                               if(negTrack->GetTPCsignal() < (2.3/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
+                                               {antilambdaCandidate = false;}
+                                       }
+                               }
+                               
+                       }
+                       
+                       if(kshortCandidate)
+                       {
+                               if(cutBetheBloch > 0)
+                               {
+                                       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(negTrack, AliPID::kPion)) > cutBetheBloch )
+                                       {kshortCandidate = false;}
+                                       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kPion)) > cutBetheBloch )
+                                       {kshortCandidate = false;}
+                               }
+                               
+                               
+                               if(cutBetheBloch == -4)  
+                               { 
+                                       double par0 = 0.20;
+                                       double par1 = 4.2;
+                                       double par2 = 1000000;
+                                       
+                                       if(isMonteCarlo) 
+                                       { 
+                                               double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+par0*par0))),2);
+                                               double gamma2 = 1.0/(1.0-beta2);
+                                               if(negTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pNeg2) > -0.6)
+                                               {kshortCandidate = false;}
+                                               
+                                               beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+par0*par0))),2);
+                                               gamma2 = 1.0/(1.0-beta2);
+                                               if(posTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pNeg2) > -0.6)
+                                               {kshortCandidate = false;}
+                                       }
+                                       else 
+                                       { 
+                                               double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+par0*par0))),2);
+                                               double gamma2 = 1.0/(1.0-beta2);
+                                               if(negTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pNeg2) > -0.6)
+                                               {kshortCandidate = false;}
+                                               
+                                               beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+par0*par0))),2);
+                                               gamma2 = 1.0/(1.0-beta2);
+                                               if(posTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pPos2) > -0.6)
+                                               {kshortCandidate = false;}
+                                       }
+                               }
+                               
+                       }
+               }
+               
+               // Selection of random cuts which I've been playing with
+               /*if(nImpactxy > 3 || pImpactxy > 2)
+                {                              
+                lambdaCandidate = false;
+                }*/
+               
+               /*if(nImpactxy < 0.4 || pImpactxy < 0.4 || nImpactxy > 2.5 || pImpactxy >2.5)
+                {                              
+                antilambdaCandidate = false;
+                }      */
+               
+               //if(decayLength > 25 )
+               //{lambdaCandidate = false;}
+               
+               // Cuts to look at just lowpt region of lambdas
+               if(v0->Pt() < 0.3 || v0->Pt() > 0.7 || !lambdaCandidate)
+               {continue;}
+               
+               // cuts to just look at signal/background region of lambda mass
+               //if(!((v0->GetEffMass(4,2) >= 1.096 && v0->GetEffMass(4,2) < 1.106) || (v0->GetEffMass(4,2) >= 1.126 && v0->GetEffMass(4,2) < 1.136) ))
+               //if(!(v0->GetEffMass(4,2) >= 1.106 && v0->GetEffMass(4,2) < 1.126 ))
+               //{continue;}
+               /*********************************************************************/
+
+               /*********************************************************************/
+               // MONTE CARLO SECTION 2
+               // this section looks at the individual V0s
+               
+               bool mcLambdaCandidate = true;
+               bool mcAntilambdaCandidate = true;
+               bool mcK0Candidate = true;
+               bool isPrimaryMC2 = false;
+               bool realParticle = true;
+               
+               if(isMonteCarlo) 
+               {
+                       
+                       AliMCEvent *mcEvent = MCEvent();
+                       AliStack* mcStack = mcEvent->Stack();
+                       if( !mcStack ) { Printf( "Stack not available"); return; }
+                       
+                       TParticle *negParticle = mcStack->Particle( TMath::Abs(negTrack->GetLabel())); 
+                       TParticle *posParticle = mcStack->Particle( TMath::Abs(posTrack->GetLabel())); 
+                       
+                       Int_t negParticleMotherLabel = negParticle->GetFirstMother();
+                       Int_t posParticleMotherLabel = posParticle->GetFirstMother();
+                       
+                       if( negParticleMotherLabel == -1 || posParticleMotherLabel == -1)
+                       {
+                       realParticle = false;
+                       mcLambdaCandidate = false;
+                       mcAntilambdaCandidate = false;
+                       mcK0Candidate  =false;
+                       }
+
+                               
+                       
+
+                       
+                       if( negParticleMotherLabel != posParticleMotherLabel)
+                       {
+                               mcLambdaCandidate = false;
+                               mcAntilambdaCandidate = false;
+                               mcK0Candidate  =false;
+                       }
+                       
+                       if(realParticle == true)
+                       {
+                               AliMCParticle *mcPart2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(negParticleMotherLabel));
+                               fHistMcTPCTrackLength->Fill(mcPart2->GetNumberOfTrackReferences());
+                               Int_t firstDaughterLabel2 = mcPart2->GetFirstDaughter();
+                               Int_t lastDaughterLabel2 = mcPart2->GetLastDaughter();
+                               AliMCParticle *mcFirstDaughter2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(firstDaughterLabel2));
+                               AliMCParticle *mcLastDaughter2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(lastDaughterLabel2));
+                               fHistMcTPCTrackLength->Fill(mcFirstDaughter2->GetNumberOfTrackReferences());
+                               fHistMcTPCTrackLength->Fill(mcLastDaughter2->GetNumberOfTrackReferences());
+                               
+                               if(mcPart2->PdgCode() != kLambda0)
+                                       {mcLambdaCandidate = false;}
+                               if(mcPart2->PdgCode() != kLambda0Bar)
+                                       {mcAntilambdaCandidate = false;}
+                               if(mcPart2->PdgCode() != kK0Short)
+                                       {mcK0Candidate  =false;}
+                               
+                               if( mcEvent->IsPhysicalPrimary(negParticleMotherLabel) )
+                                       {isPrimaryMC2 = true;}
+                               
+                               Int_t motherMotherLabel2 = mcPart2->GetMother();
+                               AliMCParticle *mcMotherMother2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(motherMotherLabel2));
+                               if(motherMotherLabel2 >= 0)
+                               {
+                                       if((mcMotherMother2->PdgCode() == 3212 && mcLambdaCandidate) || (mcMotherMother2->PdgCode() == -3212 && mcAntilambdaCandidate))
+                                       {
+                                               if(mcEvent->IsPhysicalPrimary(motherMotherLabel2) )
+                                               {isPrimaryMC2 = true;}
+                                       }
+                               }
+                               
+                               
+                               
+                               double mcRadius = TMath::Sqrt((mcPart2->Xv())*(mcPart2->Xv())+(mcPart2->Yv())*(mcPart2->Yv()));
+                               
+                               if(mcLambdaCandidate) 
+                               {
+                                       lambdaLabel[iV0] = negParticleMotherLabel;
+                                       bool once = false;
+                                       for(int j = 0; j < iV0; j++)
+                                       {
+                                               if(lambdaLabel[iV0] == lambdaLabel[j])
+                                               {       
+                                                       if(!once)
+                                                               {
+                                                                       fHistMcLog->Fill(6);
+                                                                       once = true;
+                                                               }
+                                                       fHistMcLog->Fill(7);
+                                               }       
+                                       }
+                               }
+                               
+                               if(mcLambdaCandidate && lambdaCandidate)
+                               {
+                                       fHistMcPtV0La->Fill(v0->Pt());
+                                       fHistMcV0MLaPt->Fill(v0->Pt(),v0->GetEffMass(4,2)); 
+                               }
+                               if(mcAntilambdaCandidate && antilambdaCandidate)
+                               {
+                                       fHistMcPtV0Lb->Fill(v0->Pt());
+                                       fHistMcV0MLbPt->Fill(v0->Pt(),v0->GetEffMass(2,4));
+                               }
+                               if(mcK0Candidate && kshortCandidate)
+                               {
+                                       fHistMcPtV0K0->Fill(v0->Pt());
+                                       fHistMcV0MK0Pt->Fill(v0->Pt(),v0->GetEffMass(2,2));
+                               }
+                               
+                               if(mcLambdaCandidate && isPrimaryMC2 && lambdaCandidate)
+                                       {fHistMcLog->Fill(1);}
+                               if(mcLambdaCandidate && isPrimaryMC2 && !lambdaCandidate)
+                                       {fHistMcLog->Fill(2);}
+                               if(mcLambdaCandidate && !isPrimaryMC2 && lambdaCandidate)
+                               {
+                                       fHistMcLog->Fill(3);
+                                       fHistMcV0IDLamdaSProductionVertex->Fill(mcPart2->Zv(),mcRadius);
+                                       fHistMcV0LamdaSProductionVertex->Fill(mcPart2->Zv(),mcRadius);
+                               }
+                               if(mcLambdaCandidate && !isPrimaryMC2 && !lambdaCandidate)
+                               {
+                                       fHistMcLog->Fill(4);
+                                       fHistMcV0LamdaSProductionVertex->Fill(mcPart2->Zv(),mcRadius);
+                               }
+
+                                       
+                       }
+                       
+
+                       
+                       AliMCParticle *posParticle2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(posTrack->GetLabel())));
+                       AliMCParticle *negParticle2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(negTrack->GetLabel())));
+                       
+                       if(posParticle2->PdgCode() == kProton)
+                       {
+                               fHistMcSigmaPProton->Fill(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kProton)));
+                       }
+                       if(posParticle2->PdgCode() != kProton)
+                       {
+                               fHistMcSigmaPNProton->Fill(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kProton)));
+                       }
+               }
+       
+               // END OF MONTE CARLO SECTION 2
+               /*********************************************************************/
+
+               //remove all non-candidates
+               if(lambdaCandidate == false && antilambdaCandidate == false && kshortCandidate == false)
+               {continue;}
+               
+               
+               //count number of valid v0s
+               nV0s+=1;
+                       
+               /*********************************************************************/
+               //This section fills histograms
+               
+               if(lambdaCandidate)
+               {
+                       if((v0->GetEffMass(4,2) >= 1.096 && v0->GetEffMass(4,2) < 1.106) || (v0->GetEffMass(4,2) >= 1.126 && v0->GetEffMass(4,2) < 1.136) )
+                       {fHistDCAPtPBg->Fill(v0->Pt(),nImpactxy);}
+                       if(v0->GetEffMass(4,2) >= 1.106 && v0->GetEffMass(4,2) < 1.126 )
+                       {fHistDCAPtPSig->Fill(v0->Pt(),nImpactxy);}
+               }
+               
+               if(antilambdaCandidate)
+               {
+                       if((v0->GetEffMass(2,4) >= 1.096 && v0->GetEffMass(2,4) < 1.106) || (v0->GetEffMass(2,4) >= 1.126 && v0->GetEffMass(2,4) < 1.136) )
+                       {fHistDCAPtPbarBg->Fill(v0->Pt(),pImpactxy);}
+                       if(v0->GetEffMass(2,4) >= 1.106 && v0->GetEffMass(2,4) < 1.126 )
+                       {fHistDCAPtPbarSig->Fill(v0->Pt(),pImpactxy);}
+               }
+                       
+               
+               
+               fHistEtaPTracks->Fill(posTrack->Eta());
+               fHistEtaNTracks->Fill(negTrack->Eta());
+               fHistEtaV0s->Fill(v0->Eta());
+               
+               
+               if(lambdaCandidate)
+               {
+                       double laMass = v0->GetEffMass(4,2);
+                       
+                       if((laMass > 1.096 && laMass < 1.106) || (laMass > 1.126 && laMass < 1.136) )
+                       {
+                       fHistLambdaBgRapidity->Fill(v0->Y(3122));
+                       fHistLambdaBgEta->Fill(v0->Eta());
+                       fHistLambdaBgPt->Fill(v0->Pt());
+                       }
+                       if(laMass > 1.106 && laMass < 1.126 )
+                       {
+                       fHistLambdaSigRapidity->Fill(v0->Y(3122));
+                       fHistLambdaSigEta->Fill(v0->Eta());
+                       fHistLambdaSigPt->Fill(v0->Pt());
+                       }
+                       
+                       fHistNTPCPosClustersMLa->Fill((posTrack->GetTPCNcls()),(v0->GetEffMass(4,2)));
+                       fHistNTPCNegClustersMLa->Fill((negTrack->GetTPCNcls()),(v0->GetEffMass(4,2)));
+                       fHistPosChi2PerClusterMLa->Fill((posTrack->GetTPCchi2())/(posTrack->GetTPCNcls()),(v0->GetEffMass(4,2)));
+                       fHistNegChi2PerClusterMLa->Fill((negTrack->GetTPCchi2())/(negTrack->GetTPCNcls()),(v0->GetEffMass(4,2)));
+                       if((posTrack->GetKinkIndex(0) <= 0) && (negTrack->GetKinkIndex(0) <= 0))
+                       {fHistKinkIndexFalse->Fill((v0->GetEffMass(4,2)));}
+                       else
+                       {fHistKinkIndexTrue->Fill((v0->GetEffMass(4,2)));}
+                       UInt_t status = posTrack->GetStatus();
+                       fHistPosTPCRefitMLa->Fill((status&AliESDtrack::kTPCrefit),(v0->GetEffMass(4,2)));
+                       status = negTrack->GetStatus();
+                       fHistNegTPCRefitMLa->Fill((status&AliESDtrack::kTPCrefit),(v0->GetEffMass(4,2)));                                
+               }
+               
+               fHistDCAV0Daughters->Fill(v0->GetDcaV0Daughters());
+               fHistDecayLDCA->Fill(decayLength,v0->GetDcaV0Daughters());
+               
+               fHistCosPA->Fill(v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]));
+               fHistDecayL->Fill(decayLength);
+               fHistDecayLxy->Fill(sqrt((tV0Position[0]-tPVPosition[0])*(tV0Position[0]-tPVPosition[0])+(tV0Position[1]-tPVPosition[1])*(tV0Position[1]-tPVPosition[1])));
+               fHistTauLa->Fill(cTauLa);
+               if(lambdaCandidate)
+               {
+                       fHistCosPAMLa->Fill(v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]),(v0->GetEffMass(4,2)));
+                       fHistDecayLMLa->Fill(decayLength,(v0->GetEffMass(4,2)));
+                       fHistTauLaMLa->Fill(cTauLa,(v0->GetEffMass(4,2)));
+               }
+               
+
+               {fHistBetheBlochTPCPos->Fill(TMath::Log10(pPos2),posTrack->GetTPCsignal());     }
+               {fHistBetheBlochTPCNeg->Fill(TMath::Log10(pNeg2),negTrack->GetTPCsignal());}
+               
+               fHistBetheBlochITSPos->Fill(TMath::Log10(pPos2),posTrack->GetITSsignal());
+               fHistBetheBlochITSNeg->Fill(TMath::Log10(pNeg2),negTrack->GetITSsignal());
+               
+               fHistImpactxyN->Fill(nImpactxy);
+               fHistImpactzN->Fill(nImpactz);
+               fHistImpactxyP->Fill(pImpactxy);
+               fHistImpactzP->Fill(pImpactz);
+               
+               if(lambdaCandidate)
+               {
+                       fHistImpactxyMLa->Fill(nImpactxy,v0->GetEffMass(4,2));
+                       fHistImpactzMLa->Fill(nImpactz,v0->GetEffMass(4,2));
+               }
+               fHistImpactxyImpactz->Fill(nImpactxy,nImpactz);
+               
+               fHistV0Z->Fill(tV0Position[2]);
+               fHistZ->Fill(tV0Position[2]-tPVPosition[2]);
+               fHistPtArmR->Fill(radius,v0->PtArmV0());
+               
+               fHistRZ->Fill(tV0Position[2],radius);
+               fHistPtV0Z->Fill(v0->Pt(),tV0Position[2]);
+               
+               fHistPtArm->Fill(v0->AlphaV0(),v0->PtArmV0());
+               fHistXZ->Fill(tV0Position[2],tV0Position[0]);
+               fHistYZ->Fill(tV0Position[2],tV0Position[1]);
+               fHistPtArmR->Fill(radius,v0->PtArmV0());
+               fHistPtV0->Fill(v0->Pt());
+               
+               //effective mass histograms
+               
+               //sets assumed particle type of pos/neg daughters. 
+               // 0 = electron, 1 = Muon, 2 = pion, 3 = kaon, 4 = proton.
+               int dPos = 0;
+               int dNeg = 0;
+               
+               //    v0->ChangeMassHypothesis(kLambda0);
+               dPos = 4;
+               dNeg = 2;
+               if(v0->GetEffMass(dPos,dNeg) > 1.11 && v0->GetEffMass(dPos,dNeg) < 1.13)
+               {
+                       if(!(v0->GetOnFlyStatus()))
+                       {
+                               nLambda++;
+                       }
+               }
+               if(lambdaCandidate)
+               {
+                       fHistMLa->Fill(v0->GetEffMass(dPos,dNeg));
+                       fHistMLaR->Fill(radius,v0->GetEffMass(dPos,dNeg));
+                       fHistMLaPtArm->Fill(v0->PtArmV0(),v0->GetEffMass(dPos,dNeg));   
+                       fHistMLaPt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));   
+               }
+               
+               //    v0->ChangeMassHypothesis(kK0Short);
+               dPos = 2;
+               dNeg = 2;
+               if(kshortCandidate)
+               {
+                       fHistMK0->Fill(v0->GetEffMass(dPos,dNeg));
+                       fHistMK0R->Fill(radius,v0->GetEffMass(dPos,dNeg));
+                       fHistMK0PtArm->Fill(v0->PtArmV0(),v0->GetEffMass(dPos,dNeg));
+                       fHistMK0Pt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
+               }
+               //    v0->ChangeMassHypothesis(kLambda0Bar);
+               dPos = 2;
+               dNeg = 4;
+               if(antilambdaCandidate)
+               {
+                       fHistMLb->Fill(v0->GetEffMass(dPos,dNeg));
+                       fHistMLbR->Fill(radius,v0->GetEffMass(dPos,dNeg));
+                       fHistMLbPtArm->Fill(v0->PtArmV0(),v0->GetEffMass(dPos,dNeg));
+                       fHistMLbPt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
+               }
+               
+               float binMin = 0.20;
+               float binWidth = 0.35;
+               
+               //   v0->ChangeMassHypothesis(kLambda0);
+               
+               if(lambdaCandidate)
+               {
+                       dPos = 4;
+                       dNeg = 2;
+                       if((v0->Pt()) < binMin || (v0->Pt())>= (binMin + 10*binWidth))
+                               continue;
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 1*binWidth))    
+                               fHistMLa0->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin  && (v0->Pt())< (binMin + 2*binWidth))    
+                               fHistMLa1->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 3*binWidth))    
+                               fHistMLa2->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 4*binWidth))    
+                               fHistMLa3->Fill(v0->GetEffMass(dPos,dNeg));    
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 5*binWidth))    
+                               fHistMLa4->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 6*binWidth))    
+                               fHistMLa5->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin  && (v0->Pt())< (binMin + 7*binWidth))    
+                               fHistMLa6->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 8*binWidth))    
+                               fHistMLa7->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 9*binWidth))    
+                               fHistMLa8->Fill(v0->GetEffMass(dPos,dNeg));    
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 10*binWidth))    
+                               fHistMLa9->Fill(v0->GetEffMass(dPos,dNeg));
+               }
+               //    v0->ChangeMassHypothesis(kLambda0Bar);
+               
+               if(antilambdaCandidate)
+               {
+                       dPos = 2;
+                       dNeg = 4;
+                       if((v0->Pt()) < binMin || (v0->Pt())>= (binMin + 10*binWidth))
+                               continue;
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 1*binWidth))    
+                               fHistMLb0->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin  && (v0->Pt())< (binMin + 2*binWidth))    
+                               fHistMLb1->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 3*binWidth))    
+                               fHistMLb2->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 4*binWidth))    
+                               fHistMLb3->Fill(v0->GetEffMass(dPos,dNeg));    
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 5*binWidth))    
+                               fHistMLb4->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 6*binWidth))    
+                               fHistMLb5->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin  && (v0->Pt())< (binMin + 7*binWidth))    
+                               fHistMLb6->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 8*binWidth))    
+                               fHistMLb7->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 9*binWidth))    
+                               fHistMLb8->Fill(v0->GetEffMass(dPos,dNeg));    
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 10*binWidth))    
+                               fHistMLb9->Fill(v0->GetEffMass(dPos,dNeg));
+               }
+               
+               if(kshortCandidate)
+               {
+                       dPos = 2;
+                       dNeg = 2;
+                       if((v0->Pt()) < binMin || (v0->Pt())>= (binMin + 10*binWidth))
+                               continue;
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 1*binWidth))    
+                               fHistMK00->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin  && (v0->Pt())< (binMin + 2*binWidth))    
+                               fHistMK01->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 3*binWidth))    
+                               fHistMK02->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 4*binWidth))    
+                               fHistMK03->Fill(v0->GetEffMass(dPos,dNeg));    
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 5*binWidth))    
+                               fHistMK04->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 6*binWidth))    
+                               fHistMK05->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin  && (v0->Pt())< (binMin + 7*binWidth))    
+                               fHistMK06->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 8*binWidth))    
+                               fHistMK07->Fill(v0->GetEffMass(dPos,dNeg));
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 9*binWidth))    
+                               fHistMK08->Fill(v0->GetEffMass(dPos,dNeg));    
+                       else if((v0->Pt()) >= binMin && (v0->Pt())< (binMin + 10*binWidth))    
+                               fHistMK09->Fill(v0->GetEffMass(dPos,dNeg));
+               }
+               
+    }
+       
+       // track loop, runs over all tracks in event
+       int nTracks = 0;
+       for (Int_t iTrack=0;  iTrack < fESD->GetNumberOfTracks(); iTrack++)
+       {
+               
+               AliESDtrack* track=fESD->GetTrack(iTrack);
+               //Float_t xy=0;
+               //Float_t z=0;
+               //track->GetImpactParameters(xy,z);
+               //&&(xy<1.0)&&(z<1.0)) 
+               if ( !fTrackCuts->IsSelected(track)) 
+               {continue;}
+               if (cutEta != -999 && TMath::Abs(track->Eta()) > cutEta)
+               {continue;}
+               
+               nTracks += 1;
+               fHistEta->Fill(track->Eta());
+               
+       }
+       
+       fHistPVZ->Fill(tPVPosition[2]);
+       fHistNV0->Fill(nV0s);
+       fHistNTracks->Fill(nTracks);    
+       fHistNV0sNTracks->Fill(nTracks*nTracks,nV0s); 
+       fHistdNV0sdT2->Fill( (1.0*nV0s)/(nTracks*nTracks));
+       fHistMagneticField->Fill(fESD->GetMagneticField()); 
+       fHistNLambda->Fill(nLambda);    
+       
+       
+       fHistLuke->Fill(3); // test histogram
+    // NEW HISTO should be filled before this point, as PostData puts the
+    // information for this iteration of the UserExec in the container
+       PostData(1, fOutputList);
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskLukeV0::Terminate(Option_t *) 
+{
+    // Draw result to screen, or perform fitting, normalizations
+    // Called once at the end of the query
+        
+    fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+    if(!fOutputList) { Printf("ERROR: could not retrieve TList fOutputList"); return; }
+      
+ }
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.h b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.h
new file mode 100644 (file)
index 0000000..7082a87
--- /dev/null
@@ -0,0 +1,192 @@
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliAnalysisTaskLukeV0.h 45956 2010-12-10 12:55:37Z agheata $ */
+/* AliAnalysisTaskLukeV0.h
+ *
+ * Header file for task AliAnalysisTaskLukeV0.cxx
+ * The task is designed to extract Lambda, Antilambda & K0s spectra from PbPb collision data/MC 
+ *
+ * Based on tutorial example from offline pages
+ * Edited by Arvinder Palaha
+ * Further adapted by Luke Hanratty
+ */
+
+#ifndef ALIANALYSISTASKLukeV0_H
+#define ALIANALYSISTASKLukeV0_H
+
+class TH1F;
+class TH2F;
+class TList;
+class AliESDtrackCuts;
+class AliPIDResponse;
+
+#ifndef ALIANALYSISTASKSE_H
+#include "AliAnalysisTaskSE.h"
+#endif
+
+class AliAnalysisTaskLukeV0 : public AliAnalysisTaskSE {
+ public:
+    AliAnalysisTaskLukeV0();
+    AliAnalysisTaskLukeV0(const char *name);
+    virtual ~AliAnalysisTaskLukeV0();
+    
+    virtual void     UserCreateOutputObjects();
+    virtual void     UserExec(Option_t *option);
+    virtual void     Terminate(Option_t *);
+    
+ private:
+    TList           *fOutputList;                      // Output list
+    AliESDtrackCuts *fTrackCuts;                       // Track cuts
+       AliPIDResponse  *fPIDResponse;                  // PID
+    TH1F            *fHistPt;                          // Pt spectrum
+    TH1F            *fHistEta;                         // pseudorapidity spectrum
+       TH1F                    *fHistLuke;                             // test added histogram
+       TH1F                    *fHistBetaV0;                   // beta of V0s
+       TH1F                    *fHistCosPA;                    // Cosine of pointing angle
+       TH1F                    *fHistCosTheta;                 // Cosine of decay angle in 'lambda' rest frame
+       TH1F                    *fHistCosThetaPi2;              //      cosine of 'pion' decay angle in lab frame
+       TH1F                    *fHistCosThetaProton2;  //      cosine of 'proton' decay angle in lab frame
+       TH1F                    *fHistDCAV0Daughters;   //      DCA between the v0 daughters
+       TH1F                    *fHistDecayL;                   //      3d decay length of V0
+       TH1F                    *fHistDecayLxy;                 //      2d decay length of V0 in xy plane
+       TH1F                    *fHistDeltaTheta;               //      Difference between theta of 'pion' & 'proton' in 'lambda' rest frame
+       TH1F                    *fHistdNV0sdT2;                 //      Ratio of number of v0s to number of tracks squared
+       TH1F                    *fHistEtaNTracks;               //      Eta of negative tracks
+       TH1F                    *fHistEtaPTracks;               //      Eta of positive tracks
+       TH1F                    *fHistEtaV0s;                   //      Eta of v0s
+       TH1F                    *fHistImpactxyN;                //      Impact paramater of negative daughter in xy plane
+       TH1F                    *fHistImpactzN;                 //  Impact paramater of negative daughter in z direction
+       TH1F                    *fHistImpactxyP;                //  Impact paramater of positive daughter in xy plane
+       TH1F                    *fHistImpactzP;                 //  Impact paramater of positive daughter in z direction
+       TH1F                    *fHistKinkIndexFalse;   //      Lambda mass for non-kink candidates
+       TH1F                    *fHistKinkIndexTrue;    //      Lambda mass for kink candidates
+       TH1F                    *fHistLambdaBgRapidity; //      Rapididity of V0s in sidebands of lambda mass
+       TH1F                    *fHistLambdaBgEta;              //      Eta of V0s in sidebands of lambda mass
+       TH1F                    *fHistLambdaBgPt;               //      Transverse momentum of V0s in sidebands of lambda mass
+       TH1F                    *fHistLambdaSigRapidity;//      Rapididity of V0s in peak of lambda mass
+       TH1F                    *fHistLambdaSigEta;             //      Eta of V0s in peak of lambda mass
+       TH1F                    *fHistLambdaSigPt;              //      Transverse momentum of V0s in peak of lambda mass
+       TH1F                    *fHistMagneticField;    //      Magnetic field of events
+       TH1F                    *fHistMcLog;                    //      Log of monte-carlo related information. Code in histogram.
+       TH1F                    *fHistMcNLambdaPrimary; //      Number of primary lambdas in MC event
+       TH1F                    *fHistMcNLambda;                //      Number of lambdas in MC event
+       TH1F                    *fHistMcNAntilambda;    //      Number of antilambdas in MC event
+       TH1F                    *fHistMcNKshort;                //      Number of K0s in MC event
+       TH1F                    *fHistMcPtV0La;                 //      Transverse Momentum distribution of lambdas in MC data
+       TH1F                    *fHistMcPtV0Lb;                 //      Transverse Momentum distribution of antilambdas in MC data
+       TH1F                    *fHistMcPtV0K0;                 //      Transverse Momentum distribution of K0s in MC data
+       TH1F                    *fHistMcSigmaPProton;   //      TPC response sigma of being a proton for MC protons
+       TH1F                    *fHistMcSigmaPNProton;  //      TPC response sigma of being a proton for MC non-protons
+       TH1F                    *fHistMcSLambdaEta;             //      Eta distribution of secondary lambdas in MC
+       TH1F                    *fHistMcSLambdaDaughter;//      Daughter PDG code of Secondary lambdas in MC
+       TH1F                    *fHistMcSLambdaPt;              //      Transverse momentum distribution of secondary lambdas in MC
+       TH1F                    *fHistMcTPCTrackLength; //      TPC track length in MC - DOESNT WORK
+       TH1F                    *fHistMK0;                              //      Reconstructed K0 mass for all V0s
+       TH1F                    *fHistMK00;                             //      Reconstructed K0 mass satisfying chosen conditions, typically momentum range 
+       TH1F                    *fHistMK01;                     //      "
+       TH1F                    *fHistMK02;                     //      "
+       TH1F                    *fHistMK03;                     //      "
+       TH1F                    *fHistMK04;                     //      "
+       TH1F                    *fHistMK05;                     //      "
+       TH1F                    *fHistMK06;                     //      "
+       TH1F                    *fHistMK07;                     //      "
+       TH1F                    *fHistMK08;                     //      "
+       TH1F                    *fHistMK09;                     //      "
+       TH1F                    *fHistMLa;                              //      Reconstructed Lambda mass for all V0s
+       TH1F                    *fHistMLa0;                     //      Reconstructed Lambda mass satisfying chosen conditions, typically momentum range
+       TH1F                    *fHistMLa1;                     //      "
+       TH1F                    *fHistMLa2;                     //      "
+       TH1F                    *fHistMLa3;                     //      "
+       TH1F                    *fHistMLa4;                     //      "
+       TH1F                    *fHistMLa5;                     //      "
+       TH1F                    *fHistMLa6;                     //      "
+       TH1F                    *fHistMLa7;                     //      "
+       TH1F                    *fHistMLa8;                     //      "
+       TH1F                    *fHistMLa9;                     //      "
+       TH1F                    *fHistMLb;                              //      Reconstructed Antilambda mass for all V0s
+       TH1F                    *fHistMLb0;                     //      Reconstructed Antilambda mass satisfying chosen conditions, typically momentum range
+       TH1F                    *fHistMLb1;                     //      "
+       TH1F                    *fHistMLb2;                     //      "
+       TH1F                    *fHistMLb3;                     //      "
+       TH1F                    *fHistMLb4;                     //      "
+       TH1F                    *fHistMLb5;                     //      "
+       TH1F                    *fHistMLb6;                     //      "
+       TH1F                    *fHistMLb7;                     //      "
+       TH1F                    *fHistMLb8;                     //      "
+       TH1F                    *fHistMLb9;                     //      "
+       TH1F                    *fHistNLambda;                  //      Number of lambda candidates in peak region per event
+       TH1F                    *fHistNV0;                              //      Number of V0s per event
+       TH1F                    *fHistNTracks;                  //      Number of tracks per event
+       TH1F                    *fHistPtV0;                             //      Transverse momentum distribution of V0s
+       TH1F                    *fHistPVZ;                              //      Distribution of Z coordinate of primary vertex
+       TH1F                    *fHistTauLa;                    //      Distribution of 'lambda' lifetime*c
+       TH1F                    *fHistTheta;                    //      Angle of decay in 'lambda' rest frame
+       TH1F                    *fHistThetaPi2;                 //      Angle of decay of 'pion' in rest frame
+       TH1F                    *fHistThetaProton2;             //      Angle of decay of 'proton' in rest frame
+       TH1F                    *fHistV0Z;                              //      Z coordinate of V0 decay
+       TH1F                    *fHistZ;                                //      Distance in Z between primary vertex and v0 decay
+       TH2F                    *fHistBetaERatio;               //      Beta of V0 vs ratio of energies between 'pion' & 'proton'
+       TH2F                    *fHistBetaPRatio;               //      Beta of V0 vs ratio of momenta between 'pion' & 'proton'
+       TH2F                    *fHistBetaPXRatio;              //      Beta of V0 vs ratio of x-momenta between 'pion' & 'proton'
+       TH2F                    *fHistBetheBlochITSNeg; //      ITS response for negative daughters vs momentum
+       TH2F                    *fHistBetheBlochITSPos; //      ITS response for positive daughters vs momentum
+       TH2F                    *fHistBetheBlochTPCNeg; //      TPC response for negative daughters vs momentum
+       TH2F                    *fHistBetheBlochTPCPos; //      TPC response for positive daughters vs momentum
+       TH2F                    *fHistCosPAMLa;                 //      Cosine of pointing angle vs reconstructed lambda mass
+       TH2F                    *fHistDCAPtPSig;                //      DCA vs momenum for 'proton' from 'lambda' for V0s in lambda mass peak
+       TH2F                    *fHistDCAPtPBg;                 //      DCA vs momenum for 'proton' from 'lambda' for V0s in lambda mass sidebands
+       TH2F                    *fHistDCAPtPbarSig;             //      DCA vs momenum for 'antiproton' from 'antilambda' for V0s in antilambda mass peak
+       TH2F                    *fHistDCAPtPbarBg;              //      DCA vs momenum for 'antiproton' from 'antilambda' for V0s in antilambda mass sidebands
+       TH2F                    *fHistDecayLDCA;                //      Decay length of V0 vs DCA
+       TH2F                    *fHistDecayLMLa;                //      Decay length of v0 vs reconstructed 'lambda' mass
+       TH2F                    *fHistDeltaThetaMLa;    //      Difference between theta of 'pion' & 'proton' in 'lambda' rest frame vs reconstructed 'lambda' mass
+       TH2F                    *fHistImpactxyImpactz;  //      Impact paramater in xy vs z
+       TH2F                    *fHistImpactxyMLa;              //      Impact paramater in xy vs reconstructed lambda mass
+       TH2F                    *fHistImpactzMLa;               //  Impact paramater in z vs reconstructed lambda mass
+       TH2F                    *fHistMcLamdaPProductionVertex;         //      Production vertex of primary lambdas in MC 
+       TH2F                    *fHistMcLamdaSProductionVertex;         //      Production vertex of secondary lambdas in MC 
+       TH2F                    *fHistMcLamdaSDecayVertex;                      //      Decay vertex of secondary lambdas in MC 
+       TH2F                    *fHistMcPMK0Pt;                 //      Transverse momentum distribution vs reconstructed K0 mass of primary K0s in MC
+       TH2F                    *fHistMcPMLaPt;                 //      Transverse momentum distribution vs reconstructed Lambda mass of primary Lambda in MC
+       TH2F                    *fHistMcPMLbPt;                 //      Transverse momentum distribution vs reconstructed Antilambd mass of primary Antilambda in MC
+       TH2F                    *fHistMcSLambdaDaughterPairs;           //      PDG codes of daughters of Secondary lambdas in MC
+       TH2F                    *fHistMcV0MK0Pt;                //      Reconstructed K0 mass vs transverse momentum of V0s passing cuts in MC
+       TH2F                    *fHistMcV0MLaPt;                //      Reconstructed Lambda mass vs transverse momentum of V0s passing cuts in MC
+       TH2F                    *fHistMcV0MLbPt;                //      Reconstructed Antilambda mass vs transverse momentum of V0s passing cuts in MC
+       TH2F                    *fHistMcV0LamdaSProductionVertex;       //      Production vertex of V0 secondary Lambdas in RZ plane in MC
+       TH2F                    *fHistMcV0IDLamdaSProductionVertex;     //  Production vertex of identified V0 secondary Lambdas in RZ plane in MC
+       TH2F                    *fHistMK0PtArm;                 //      Mass of 'K0' vs Armenteros transverse momentum
+       TH2F                    *fHistMK0Pt;                    //      Mass of 'K0' vs transverse momentum
+       TH2F                    *fHistMK0R;                             //      Mass of 'K0' vs radius
+       TH2F                    *fHistMLaPtArm;                 //      Mass of 'Lambda' vs Armenteros transverse momentum
+       TH2F                    *fHistMLaPt;                    //      Mass of 'Lambda' vs transverse momentum
+       TH2F                    *fHistMLaR;                             //      Mass of 'Lambda' vs radius
+       TH2F                    *fHistMLbPtArm;                 //      Mass of 'Antilambda' vs Armenteros transverse momentum
+       TH2F                    *fHistMLbPt;                    //      Mass of 'Antilambda' vs transverse momentum
+       TH2F                    *fHistMLbR;                             //      Mass of 'Antilambda' vs radius  
+       TH2F                    *fHistNegChi2PerClusterMLa;                     //      Chi2 per cluster of negative daughter vs reconstructed 'lambda' mass
+       TH2F                    *fHistNegTPCRefitMLa;   //      TPC refit t/f of negative daughter vs reconstructed 'lambda' mass
+       TH2F                    *fHistNTPCNegClustersMLa;                       //      Number of TPC clusters  of negative daughter vs reconstructed 'lambda' mass
+       TH2F                    *fHistNTPCPosClustersMLa;                       // Number of TPC clusters of positive daughter vs reconstructed 'lambda' mass
+       TH2F                    *fHistNV0sNTracks;              //      Number of V0s vs number of tracks per event
+       TH2F                    *fHistPosChi2PerClusterMLa;                     //      Chi2 per cluster of positive daughter vs reconstructed 'lambda' mass
+       TH2F                    *fHistPosTPCRefitMLa;   //      TPC refit t/f of positive daughter vs reconstructed 'lambda' mass
+       TH2F                    *fHistPtArm;                    //      Podolanski-Armenteros plot of V0s
+       TH2F                    *fHistPtArmR;                   //      PtArm vs radius of v0
+       TH2F                    *fHistPtV0Z;                    //      Transverse momentum vs Z coordinate of v0 decay
+       TH2F                    *fHistRZ;                               //      Radius vs z of V0 decay
+       TH2F                    *fHistTauLaMLa;                 //      'lambda' lifetime vs reconstructed 'lambda' mass
+       TH2F                    *fHistXZ;                               //      X vs z of V0 decay
+       TH2F                    *fHistYZ;                               //      Y vs Z of V0 decay
+       
+    // NEW HISTO to be declared here
+    
+    AliAnalysisTaskLukeV0(const AliAnalysisTaskLukeV0&); // not implemented
+    AliAnalysisTaskLukeV0& operator=(const AliAnalysisTaskLukeV0&); // not implemented
+    
+    ClassDef(AliAnalysisTaskLukeV0, 1); // example of analysis
+};
+
+#endif
+
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/runLukeV0.C b/PWG2/SPECTRA/LambdaK0PbPb/runLukeV0.C
new file mode 100644 (file)
index 0000000..4222612
--- /dev/null
@@ -0,0 +1,272 @@
+// runLukeV0.C
+//
+// Run macro for AliAnalysisTaskLukeV0.cxx/.h with physics selection & PID
+//
+// sample proof datasets: 
+// "/alice/sim/LHC10f6a_000126432" (pp MC) 
+// "/alice/data/LHC10c_000120821_p1" (pp) 
+//
+// "/alice/data/LHC10h_000139172_p2" (PbPb pass 2)
+// "/alice/data/LHC10h_000138150_p1" (PbPb pass 1) 
+//
+// Original author: Arvinder Palaha
+// Adapted by: Luke Hanratty
+//
+
+#define myRunType "proof" // local, proof or grid
+#define myGridMode "full" // full, test, offline, submit or terminate
+#define mybMCtruth 0 // 0 or 1; MCEvent handler is on or off
+#define mybMCphyssel 0 // 0 = real data, 1 = MC
+#define myNEntries 5000 //num //2000 // local and proof mode only; 1234567890 = all
+#define myFirstEntry 0//8100000 // local and proof mode only - first event
+#define myProofDataset "/alice/data/LHC10h_000139172_p2" //set path
+#define myProofCluster "hanratty@alice-caf.cern.ch" // set proof cluster
+#define myTaskName "luke_task" // set name of this task
+#define myAPIVersion "V1.1x"   // set version of API
+#define myROOTVersion "v5-28-00e"   // set version of Root
+#define myAliROOTVersion "v4-21-27-AN"   // set version of AliRoot
+#define myLocalFiles "LocalFiles" // should be a txt file with a list of local files
+#define myNproofWorkers 0 // can limit the maximum number of workers. 0 is default 40.
+class AliAnalysisGrid;
+
+//______________________________________________________________________________
+void runLukeV0(
+                const char* runtype = myRunType,
+                const char* gridmode = myGridMode,
+                const bool bMCtruth = mybMCtruth,
+                const bool bMCphyssel = mybMCphyssel,
+                const Long64_t nentries = myNEntries,
+                const Long64_t firstentry = myFirstEntry,
+                const char* proofdataset = myProofDataset,
+                const char* proofcluster = myProofCluster,
+                const char* taskname = myTaskName
+                )
+{
+    // check run type
+    if(runtype != "local" && runtype != "proof" && runtype != "grid"){
+        Printf("\n\tIncorrect run option, check first argument of run macro");
+        Printf("\tint runtype = local, proof or grid\n");
+        return;
+    }
+    Printf("%s analysis chosen",runtype);
+  
+    // load libraries
+    gSystem->Load("libCore.so");        
+    gSystem->Load("libGeom.so");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libPhysics.so");
+    gSystem->Load("libTree.so");
+    gSystem->Load("libSTEERBase.so");
+    gSystem->Load("libESD.so");
+    gSystem->Load("libAOD.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libANALYSISalice.so");
+  
+    // add aliroot include path
+    gROOT->ProcessLine(Form(".include %s/include",gSystem->ExpandPathName("$ALICE_ROOT")));
+    gROOT->SetStyle("Plain");
+        
+    // analysis manager
+    AliAnalysisManager* mgr = new AliAnalysisManager(taskname);
+    
+    // create the alien handler and attach it to the manager
+    AliAnalysisGrid *plugin = CreateAlienHandler(taskname, gridmode, proofcluster, proofdataset); 
+    mgr->SetGridHandler(plugin);
+    
+    AliVEventHandler* esdH = new AliESDInputHandler();
+    mgr->SetInputEventHandler(esdH);
+        
+    // mc event handler
+    if(bMCtruth) {
+        AliMCEventHandler* mchandler = new AliMCEventHandler();
+        // Not reading track references
+        mchandler->SetReadTR(kFALSE);
+        mgr->SetMCtruthEventHandler(mchandler);
+    }   
+
+    // === Physics Selection Task ===
+    //
+    // In SelectCollisionCandidate(), default is kMB, so the task UserExec() 
+    // function is only called for these events.
+    // Options are:
+    //    kMB             Minimum Bias trigger
+    //    kMBNoTRD        Minimum bias trigger where the TRD is not read out
+    //    kMUON           Muon trigger
+    //    kHighMult       High-Multiplicity Trigger
+    //    kUserDefined    For manually defined trigger selection
+    //
+    // Multiple options possible with the standard AND/OR operators && and ||
+    // These all have the usual offline SPD or V0 selections performed.
+    //
+    // With a pointer to the physics selection object using physSelTask->GetPhysicsSelection(),
+    // one can manually set the selected and background classes using:
+    //    AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL")
+    //    AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");
+    //
+    // One can also specify multiple classes at once, or require a class to NOT
+    // trigger, for e.g.
+    //    AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");
+    //
+    // NOTE that manually setting the physics selection overrides the standard
+    // selection, so it must be done in completeness.
+    //
+    // ALTERNATIVELY, one can make the physics selection inside the task
+    // UserExec().
+    // For this case, comment out the task->SelectCol.... line, 
+    // and see AliBasicTask.cxx UserExec() function for details on this.
+
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+    AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(bMCphyssel);
+    if(!physSelTask) { Printf("no physSelTask"); return; }
+    //AliPhysicsSelection *physSel = physSelTask->GetPhysicsSelection();
+    //physSel->AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL");// #3119 #769");
+    
+       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+    AliAnalysisTask *PIDTask = AddTaskPIDResponse(bMCtruth,kTRUE);
+    if(!PIDTask) { Printf("no PIDtask"); return; }
+       
+    // create task
+    gROOT->LoadMacro("AliAnalysisTaskLukeV0.cxx++g");
+    AliAnalysisTaskSE* task = new AliAnalysisTaskLukeV0(taskname);
+    task->SelectCollisionCandidates(AliVEvent::kMB); // if physics selection performed in UserExec(), this line should be commented
+    mgr->AddTask(task);
+    
+    // set output root file name for different analysis
+    TString outfilename = Form("routput.%s.root",runtype);
+  
+    // create containers for input/output
+    AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+    AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("coutput1", TList::Class(), AliAnalysisManager::kOutputContainer, outfilename);
+        
+    // connect input/output
+    mgr->ConnectInput(task, 0, cinput);
+    mgr->ConnectOutput(task, 1, coutput1);
+        
+    // enable debug printouts
+    mgr->SetDebugLevel(2);
+    if (!mgr->InitAnalysis()) return;
+    mgr->PrintStatus();
+  
+    // start analysis
+    Printf("Starting Analysis....");
+    mgr->StartAnalysis(runtype,nentries,firstentry);
+}
+
+//______________________________________________________________________________
+AliAnalysisGrid* CreateAlienHandler(const char *taskname, const char *gridmode, const char *proofcluster, const char *proofdataset)
+{
+    AliAnalysisAlien *plugin = new AliAnalysisAlien();
+    // Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+    plugin->SetRunMode(gridmode);
+
+    // Set versions of used packages
+    plugin->SetAPIVersion(myAPIVersion);
+    plugin->SetROOTVersion(myROOTVersion);
+    plugin->SetAliROOTVersion(myAliROOTVersion);
+       
+    // Declare input data to be processed.
+
+    // Method 1: Create automatically XML collections using alien 'find' command.
+    // Define production directory LFN
+    plugin->SetGridDataDir("/alice/data/2010/LHC10b");
+    // On real reconstructed data:
+    // plugin->SetGridDataDir("/alice/data/2009/LHC09d");
+    // Set data search pattern
+    //plugin->SetDataPattern("*ESDs.root"); // THIS CHOOSES ALL PASSES
+    // Data pattern for reconstructed data
+    plugin->SetDataPattern("*ESDs/pass2/*ESDs.root"); // CHECK LATEST PASS OF DATA SET IN ALIENSH
+    plugin->SetRunPrefix("000");   // real data
+    // ...then add run numbers to be considered
+    plugin->AddRunNumber(115514);
+    //plugin->SetRunRange(114917,115322);
+    plugin->SetNrunsPerMaster(1);
+    plugin->SetOutputToRunNo();
+    // comment out the next line when using the "terminate" option, unless
+    // you want separate merged files for each run
+    plugin->SetMergeViaJDL();
+
+    // Method 2: Declare existing data files (raw collections, xml collections, root file)
+    // If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir())
+    // XML collections added via this method can be combined with the first method if
+    // the content is compatible (using or not tags)
+    //   plugin->AddDataFile("tag.xml");
+    //   plugin->AddDataFile("/alice/data/2008/LHC08c/000057657/raw/Run57657.Merged.RAW.tag.root");
+
+    // Define alien work directory where all files will be copied. Relative to alien $HOME.
+    plugin->SetGridWorkingDir(taskname);
+
+    // Declare alien output directory. Relative to working directory.
+    plugin->SetGridOutputDir("out"); // In this case will be $HOME/taskname/out
+
+    // Declare the analysis source files names separated by blancs. To be compiled runtime
+    // using ACLiC on the worker nodes.
+    plugin->SetAnalysisSource("AliAnalysisTaskLukeV0.cxx");
+
+    // Declare all libraries (other than the default ones for the framework. These will be
+    // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here.
+    plugin->SetAdditionalLibs("AliAnalysisTaskLukeV0.h AliAnalysisTaskLukeV0.cxx");
+
+    // Declare the output file names separated by blancs.
+    // (can be like: file.root or file.root@ALICE::Niham::File)
+    // To only save certain files, use SetDefaultOutputs(kFALSE), and then
+    // SetOutputFiles("list.root other.filename") to choose which files to save
+    plugin->SetDefaultOutputs();
+    //plugin->SetOutputFiles("list.root");
+
+    // Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+    plugin->SetAnalysisMacro(Form("%s.C",taskname));
+
+    // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+    plugin->SetSplitMaxInputFileNumber(100);
+
+    // Optionally modify the executable name (default analysis.sh)
+    plugin->SetExecutable(Form("%s.sh",taskname));
+
+    // set number of test files to use in "test" mode
+    plugin->SetNtestFiles(10);
+
+    // Optionally resubmit threshold.
+    plugin->SetMasterResubmitThreshold(90);
+
+    // Optionally set time to live (default 30000 sec)
+    plugin->SetTTL(30000);
+
+    // Optionally set input format (default xml-single)
+    plugin->SetInputFormat("xml-single");
+
+    // Optionally modify the name of the generated JDL (default analysis.jdl)
+    plugin->SetJDLName(Form("%s.jdl",taskname));
+
+    // Optionally modify job price (default 1)
+    plugin->SetPrice(1);      
+
+    // Optionally modify split mode (default 'se')    
+    plugin->SetSplitMode("se");
+    
+    //----------------------------------------------------------
+    //---      PROOF MODE SPECIFIC SETTINGS         ------------
+    //---------------------------------------------------------- 
+    // Proof cluster
+    plugin->SetProofCluster(proofcluster);
+    // Dataset to be used   
+    plugin->SetProofDataSet(proofdataset);
+    // May need to reset proof. Supported modes: 0-no reset, 1-soft, 2-hard
+    plugin->SetProofReset(0);
+    // May limit number of workers
+    plugin->SetNproofWorkers(myNproofWorkers);
+    // May limit the number of workers per slave
+    plugin->SetNproofWorkersPerSlave(1);   
+    // May use a specific version of root installed in proof
+    plugin->SetRootVersionForProof("current");
+    // May set the aliroot mode. Check http://aaf.cern.ch/node/83 
+    plugin->SetAliRootMode("default"); // Loads AF libs by default
+    // May request ClearPackages (individual ClearPackage not supported)
+    plugin->SetClearPackages(kFALSE);
+    // Plugin test mode works only providing a file containing test file locations, used in "local" mode also
+    plugin->SetFileForTestMode(myLocalFiles); // file should contain path name to a local directory containg *ESDs.root etc
+    // Request connection to alien upon connection to grid
+    plugin->SetProofConnectGrid(kFALSE);
+
+    return plugin;
+}
+