Initial commit
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Jan 2011 11:52:40 +0000 (11:52 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Jan 2011 11:52:40 +0000 (11:52 +0000)
PWG2/SPECTRA/LambdaK0PbPb/AddTaskLambdaK0PbPb.C [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisCentralitySelector.cxx [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisCentralitySelector.h [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.cxx [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.h [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AnalysisStrange.C [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/CreateAlienHandler.C [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/run.C [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/run.sh [new file with mode: 0755]

diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AddTaskLambdaK0PbPb.C b/PWG2/SPECTRA/LambdaK0PbPb/AddTaskLambdaK0PbPb.C
new file mode 100644 (file)
index 0000000..e66e729
--- /dev/null
@@ -0,0 +1,72 @@
+AliAnalysisTaskPerformanceStrange ** AddTaskLambdaK0PbPb(const char * outfilename, AliAnalysisCentralitySelector * centr, Int_t &nbin) {
+
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskPhysicsSelection", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AddTaskPhysicsSelection", "This task requires an input event handler");
+    return NULL;
+  }
+  TString inputDataType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  
+  if (inputDataType != "ESD") {
+    Printf("ERROR! This task can only run on ESDs!");
+  }
+
+  // Configure analysis
+  //===========================================================================
+    
+  Int_t binMin = 0; // FIXME: settable? different percentiles?
+  Int_t binMax = 10;
+  nbin = binMax - binMin + 1;
+
+  AliAnalysisTaskPerformanceStrange ** task = new AliAnalysisTaskPerformanceStrange*[nbin];
+  Int_t itask = -1;
+  for(Int_t ibin = binMin; ibin <= binMax; ibin++){  
+    itask++;
+
+    task[itask] = new AliAnalysisTaskPerformanceStrange("TaskLambdaK0");
+    cout << "Booking " << ibin << "  "<< itask << " " << task[itask] <<endl;
+    
+    mgr->AddTask(task[itask]);
+  
+    // // Set Cuts
+    // if (!esdTrackCuts)
+    //   {
+    //     printf("ERROR: esdTrackCuts could not be created\n");
+    //     return;
+    //   }  
+    // task->SetTrackCuts(esdTrackCuts);
+    
+    // set centrality
+    AliAnalysisCentralitySelector * centrBin = (AliAnalysisCentralitySelector*) centr->Clone();
+    centrBin->SetCentralityBin(ibin);
+    task[itask]->SetCentralitySelector(centrBin);
+
+    TString outfilenameCentr = outfilename;
+    outfilenameCentr.ReplaceAll(".root",Form("_%2.2d.root",ibin));
+
+    AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer();
+    AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(Form("clambdak0Histo_%2.2d",ibin), TList::Class(),AliAnalysisManager::kOutputContainer, outfilenameCentr.Data());
+    AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(Form("clambdak0Centr_%2.2d",ibin), AliAnalysisCentralitySelector::Class(),AliAnalysisManager::kOutputContainer, outfilenameCentr.Data());
+
+    mgr->ConnectInput (task[itask], 0, mgr->GetCommonInputContainer());
+    mgr->ConnectOutput(task[itask],1,coutput1);
+    mgr->ConnectOutput(task[itask],2,coutput2);
+
+  }
+  // TODO:
+  // IO into folders in a file?
+
+  // Set I/O
+
+  return task;
+}   
+
+
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisCentralitySelector.cxx b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisCentralitySelector.cxx
new file mode 100644 (file)
index 0000000..64e62a2
--- /dev/null
@@ -0,0 +1,104 @@
+// AliAnalysisMultPbCentralitySelector 
+// Interface class to centrality estimators for the PbPb
+// track-multiplicity analysis
+// Michele Floris, CERN
+
+#include "AliAnalysisCentralitySelector.h"
+#include "AliESDtrackCuts.h"
+#include "AliCentrality.h"
+#include "AliESDEvent.h"
+#include "AliLog.h"
+#include "AliESDVZERO.h"
+#include <iostream>
+#include "AliMultiplicity.h"
+
+using namespace std;
+
+
+
+ClassImp(AliAnalysisCentralitySelector)
+
+Bool_t AliAnalysisCentralitySelector::IsCentralityBinSelected(AliESDEvent* aEsd, AliESDtrackCuts * trackCuts) {
+
+  // Centrality selection
+  // On MC cuts on the number of good tracks,
+  // On data cuts using AliESDCentrality and the cut requested in ntracks
+
+  //  cout << "Tracks " << trackCuts->CountAcceptedTracks(aEsd) << endl;
+  ///  cout << "CENTRALITY " << fUseV0CutRange << " " << fUseMultRange << " " << fMultMin << " " << fMultMax << endl;
+  
+  if (fUseMultRange && fUseV0CutRange && fUseSPDOuterRange) {
+    AliFatal(Form("Cannot use multiple estimators at once: fUseMultRange [%d], fUseV0CutRange[%d], fUseSPDOuterRange[%d]!!",
+                 fUseMultRange , fUseV0CutRange , fUseSPDOuterRange)); 
+  }
+
+  if (fUseV0CutRange) {
+
+    Float_t multV0=0;
+    AliESDVZERO* esdV0 = aEsd->GetVZEROData();
+    Float_t multV0A=esdV0->GetMTotV0A();
+    Float_t multV0C=esdV0->GetMTotV0C();
+    multV0 = multV0A+multV0C;
+    
+    if (fIsMC) multV0 = 0.85871 * multV0; // FIXME: still valid?
+    if (multV0 < fMultMin) return kFALSE;
+    if (multV0 > fMultMax) return kFALSE;
+    //    cout << "ok" << endl;
+
+  } 
+  else if (fUseSPDOuterRange) {
+
+    const AliMultiplicity * mult = aEsd->GetMultiplicity();
+    Float_t outerLayerSPD = mult->GetNumberOfITSClusters(1);  
+    
+    if (outerLayerSPD < fMultMin) return kFALSE;
+    if (outerLayerSPD > fMultMax) return kFALSE;
+    //    cout << "ok" << endl;
+
+  }
+ else if(fIsMC || fUseMultRange) {
+    if(!trackCuts){
+      AliFatal("Track cuts object is invalid");
+    }
+    Float_t ntracks = trackCuts->CountAcceptedTracks(aEsd);
+    //    cout << "Hey! " << fCentrBin << " " << ntracks << " " << fMultMin <<" - " << fMultMax << endl;
+    
+    if (fCentrBin == -1 && !fUseMultRange) return kTRUE;
+    if (ntracks < fMultMin) return kFALSE;
+    if (ntracks > fMultMax) return kFALSE;                                                    
+  } else {
+
+   AliCentrality *centrality = (AliCentrality*) aEsd->GetCentrality(); // FIXME: change to alicentrality?
+    if(!centrality && !fUseMultRange) {
+      AliFatal("Centrality object not available"); 
+    }
+    else {
+      Int_t centrBin = centrality->GetCentralityClass5(fCentrEstimator.Data()) ;    
+      if (centrBin != fCentrBin && fCentrBin != -1 && !fUseMultRange) return kFALSE;
+    }
+  }
+
+  //  cout << "Selected" << endl;
+  
+
+  return kTRUE;
+
+}
+
+void AliAnalysisCentralitySelector::Print(Option_t* option ) const {
+  // Print some information
+
+  Printf("AliAnalysisCentralitySelector [%s]", option);
+  Printf(" - Centrality estimator [%s]",fCentrEstimator.Data());
+  Printf(" - Centrality bin       [%d]",fCentrBin);
+  if ( fUseMultRange ) {
+    Printf ("Using multiplicity range [%1.1f - %1.1f]",fMultMin,fMultMax);
+  }
+  if ( fUseV0CutRange ) {
+    Printf ("Using V0 range [%1.1f - %1.1f]",fMultMin,fMultMax);
+  }
+  if ( fIsMC ) {    
+    Printf("Running on Monte Carlo, actual cut was on tracks multiplicity [%1.1f - %1.1f]",fMultMin,fMultMax);    
+  } 
+  
+}
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisCentralitySelector.h b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisCentralitySelector.h
new file mode 100644 (file)
index 0000000..6604746
--- /dev/null
@@ -0,0 +1,71 @@
+#ifndef ALIANALYSISCENTRALITYSELECTOR_H
+#define ALIANALYSISCENTRALITYSELECTOR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                    AliAnalysisCentralitySelector
+// 
+// This class selects collision candidates from data runs, applying selection cuts on triggers 
+// and background rejection based on the content of the ESD
+//
+// Author: Michele Floris, CERN
+//-------------------------------------------------------------------------
+
+#include <AliAnalysisCuts.h>
+#include <AliLog.h>
+
+#define VERBOSE_STAT
+
+class AliESDEvent;
+class TH2F;
+class TH1F;
+class TCollection;
+class AliTriggerAnalysis;
+class AliAnalysisTaskSE;
+class AliESDtrackCuts;
+
+class AliAnalysisCentralitySelector : public AliAnalysisCuts
+{
+public:
+
+  AliAnalysisCentralitySelector() : fIsMC (0), fCentrEstimator(""), fCentrBin(-1), fMultMin(0), fMultMax(1000000), fUseMultRange(kFALSE), fUseV0CutRange(kFALSE), fUseSPDOuterRange(kFALSE) {;}
+  virtual ~AliAnalysisCentralitySelector(){}
+    
+  // AliAnalysisCuts interface
+  virtual UInt_t GetSelectionMask(const TObject* obj) { return (UInt_t) IsCentralityBinSelected((AliESDEvent*) obj, NULL); }
+  virtual Bool_t IsSelected(TList*) { AliFatal("Not implemented"); return kFALSE; }
+  virtual Bool_t IsSelected(TObject* obj)  {return (UInt_t) IsCentralityBinSelected ( (AliESDEvent*) obj, NULL);}
+    
+  Bool_t IsCentralityBinSelected(AliESDEvent* aEsd, AliESDtrackCuts * trackCuts);
+    
+  void SetIsMC(Bool_t flag = kTRUE, Int_t multMin = 0, Int_t multMax=10000) { fIsMC = flag; fMultMin = multMin; fMultMax = multMax; }
+  void SetMultRange(Int_t multMin = 0, Int_t multMax=10000) { fMultMin = multMin; fMultMax = multMax; }
+  void SetUseMultRange(Bool_t flag = kTRUE) {fUseMultRange = flag;}
+  void SetUseV0Range(Bool_t flag = kTRUE) {fUseV0CutRange = flag;}
+  void SetUseSPDOuterRange(Bool_t flag = kTRUE) {fUseSPDOuterRange = flag;}
+  void SetCentralityEstimator(const char * estimator) { fCentrEstimator = estimator; }
+  void SetCentralityBin(Int_t bin) { fCentrBin = bin; }
+  virtual void Print(Option_t* option = "") const ;
+  virtual Long64_t Merge(TCollection* list){list->GetEntries();return 0;}
+  
+protected:
+  Bool_t fIsMC;             // flag if MC is analyzed
+  TString fCentrEstimator;  // Centrality estimator for AliCentrality
+  Int_t   fCentrBin; // centrality bin to be selected
+  Float_t fMultMin ; // Minimum multiplicity, because on MC we cut on tracks rather than on the estimator . Also used for other estimators
+  Float_t fMultMax ; // Maximum multiplicity, because on MC we cut on tracks rather than on the estimator . Also used for other estimators
+  Bool_t fUseMultRange; // if true, use track bins rather than multiplicity estimator
+  Bool_t fUseV0CutRange; // if true, use v0 range rather than multiplicity estimator
+  Bool_t fUseCorrV0; // linearized V0
+  Bool_t fUseSPDOuterRange; // if true, use SPD outer cluster range rather than multiplicity estimator
+
+  ClassDef(AliAnalysisCentralitySelector, 1)
+    
+  private:
+  AliAnalysisCentralitySelector(const AliAnalysisCentralitySelector&); // not implemented
+  AliAnalysisCentralitySelector& operator=(const AliAnalysisCentralitySelector&); // not implemented
+};
+
+#endif
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.cxx b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.cxx
new file mode 100644 (file)
index 0000000..9f1a9d7
--- /dev/null
@@ -0,0 +1,3108 @@
+/**************************************************************************
+ * Copyright(c) 1998-2009, 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//             AliAnalysisTaskPerformanceStrange class
+//    This task is for a performance study of V0 identification.
+//                It works with MC info and ESD tree.
+//                 Author: H.Ricaud, H.Ricaud@gsi.de
+//-----------------------------------------------------------------
+
+#include <Riostream.h>
+
+#include <stdio.h>
+#include <iostream>
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+//#include "TH3F.h"
+#include "TF1.h"
+#include "TList.h"
+#include "TMath.h"
+#include "TCanvas.h"
+
+#include "AliAnalysisManager.h"
+
+#include "AliPhysicsSelection.h"
+#include "AliBackgroundSelection.h"
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrack.h"
+#include "AliESDv0.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDpid.h"
+#include "AliMultiplicity.h"
+
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODv0.h"
+#include "AliAODMCHeader.h"
+#include "AliAODInputHandler.h"
+
+//#include "AliV0vertexer.h"
+
+#include "AliAODMCParticle.h"
+
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliGenEventHeader.h"
+
+#include "AliLog.h"
+
+#include "AliKFVertex.h"
+#include "AliVertexerTracks.h"
+
+#include "AliAnalysisTaskPerformanceStrange.h"
+#include "AliAnalysisCentralitySelector.h"
+
+
+ClassImp(AliAnalysisTaskPerformanceStrange)
+
+
+//________________________________________________________________________
+AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange()
+: AliAnalysisTaskSE("TaskStrange"), fAnalysisMC(0), fAnalysisType("infoType"),  fCollidingSystems(0), fUsePID("infoPID"), fUseCut("infoCut"), fESD(0), fListHist(0), fCentrSelector(0),
+
+    fHistMCPrimaryVertexX(0),
+    fHistMCPrimaryVertexY(0),
+    fHistMCPrimaryVertexZ(0),
+    fHistMCMultiplicityPrimary(0),
+    fHistMCMultiplicityTracks(0),
+    fHistMCtracksProdRadiusK0s(0),
+    fHistMCtracksProdRadiusLambda(0),
+    fHistMCtracksProdRadiusAntiLambda(0),
+    fHistMCtracksDecayRadiusK0s(0),
+    fHistMCtracksDecayRadiusLambda(0),
+    fHistMCtracksDecayRadiusAntiLambda(0),
+    fHistMCPtAllK0s(0),
+    fHistMCPtAllLambda(0),
+    fHistMCPtAllAntiLambda(0),
+    fHistMCProdRadiusK0s(0),
+    fHistMCProdRadiusLambda(0),
+    fHistMCProdRadiusAntiLambda(0),
+    fHistMCRapK0s(0),
+    fHistMCRapInPtRangeK0s(0),
+    fHistMCRapLambda(0),
+    fHistMCRapInPtRangeLambda(0),
+    fHistMCRapAntiLambda(0),
+    fHistMCRapInPtRangeAntiLambda(0),
+    fHistMCRapXi(0),
+    fHistMCRapInPtRangeXi(0),
+    fHistMCRapPhi(0),
+    fHistMCRapInPtRangePhi(0),
+////////////////////////////////////////
+    fHistMCPtK0s(0),
+
+
+    fHistMCPtLambda(0),
+///////////////////////////////////////////
+
+    fHistMCPtLambdaFromSigma(0),
+    fHistMCPtAntiLambdaFromSigma(0),
+    fHistNTimesRecK0s(0),
+//    fHistNTimesRecK0sMI(0),
+    fHistNTimesRecLambda(0),
+//    fHistNTimesRecLambdaMI(0),
+    fHistNTimesRecAntiLambda(0),
+//    fHistNTimesRecAntiLambdaMI(0),
+    fHistNTimesRecK0sVsPt(0),
+//    fHistNTimesRecK0sVsPtMI(0),
+    fHistNTimesRecLambdaVsPt(0),
+//    fHistNTimesRecLambdaVsPtMI(0),
+    fHistNTimesRecAntiLambdaVsPt(0),
+//    fHistNTimesRecAntiLambdaVsPtMI(0),
+    fHistNumberEvents(0),
+    fHistTrackPerEvent(0),
+    fHistTrackletPerEvent(0),
+    fHistMCDaughterTrack(0),
+    fHistSPDPrimaryVertexZ(0),
+    fHistPrimaryVertexX(0),
+    fHistPrimaryVertexY(0),
+    fHistPrimaryVertexZ(0),
+    fHistPrimaryVertexResX(0),
+    fHistPrimaryVertexResY(0),
+    fHistPrimaryVertexResZ(0),
+    fHistPrimaryVertexPosXV0events(0), 
+    fHistPrimaryVertexPosYV0events(0), 
+    fHistPrimaryVertexPosZV0events(0),
+    fHistDaughterPt(0),
+    fHistDcaPosToPrimVertex(0),
+    fHistDcaNegToPrimVertex(0),
+    fHistDcaPosToPrimVertexZoom(0),
+    fHistDcaNegToPrimVertexZoom(0),
+    fHistRadiusV0(0),
+    fHistDecayLengthV0(0),
+    fHistDcaV0Daughters(0),
+    fHistChi2(0),
+    fHistCosPointAngle(0),
+    fHistCosPointAngleZoom(0),
+    fHistProdRadius(0),
+//    fHistProdRadiusMI(0),
+    fHistV0Multiplicity(0),
+//    fHistV0MultiplicityMI(0),
+    fHistChi2KFBeforeCutK0s(0), 
+    fHistChi2KFBeforeCutLambda(0), 
+    fHistChi2KFBeforeCutAntiLambda(0),
+    fHistChi2KFAfterCutK0s(0), 
+    fHistChi2KFAfterCutLambda(0), 
+    fHistChi2KFAfterCutAntiLambda(0),
+    fHistMassK0(0),
+//    fHistMassK0MI(0),
+    fHistMassLambda(0),
+//    fHistMassLambdaMI(0),
+    fHistMassAntiLambda(0),
+//    fHistMassAntiLambdaMI(0),
+    fHistMassVsRadiusK0(0),
+//    fHistMassVsRadiusK0MI(0),
+    fHistMassVsRadiusLambda(0),
+//    fHistMassVsRadiusLambdaMI(0),
+    fHistMassVsRadiusAntiLambda(0),
+//    fHistMassVsRadiusAntiLambdaMI(0),
+
+///////////////////////////////////////
+    fHistPtVsMassK0(0),
+//    fHistPtVsMassK0MI(0),
+    fHistPtVsMassLambda(0),
+////////////////////////////////////////
+       fHistPzPtBeforeK0s(0),
+       fHistPzPtAfterK0s(0),
+       fHistPzPtBeforeLambda(0),
+       fHistPzPtAfterLambda(0),
+
+//////////////////////////////////////////
+//    fHistPtVsMassLambdaMI(0),
+//    fHistPtVsMassAntiLambda(0),
+//    fHistPtVsMassAntiLambdaMI(0),
+
+
+
+//    fHistMultVsPtVsMassK0(0),
+//    fHistMultVsPtVsMassLambda(0),
+//    fHistMultVsPtVsMassAntiLambda(0),
+
+
+    fHistArmenterosPodolanski(0),
+//    fHistArmenterosPodolanskiMI(0),
+    fHistNsigmaPosPionAntiLambda(0),
+    fHistNsigmaNegProtonAntiLambda(0),
+    fHistNsigmaPosProtonLambda(0),
+    fHistNsigmaNegPionLambda(0),
+    fHistNsigmaPosPionK0(0),
+    fHistNsigmaNegPionK0(0),
+    fHistAsMcRapK0(0),
+//    fHistAsMcRapK0MI(0),
+    fHistAsMcRapLambda(0),
+//    fHistAsMcRapLambdaMI(0),
+    fHistAsMcRapAntiLambda(0),
+//    fHistAsMcRapAntiLambdaMI(0),
+
+/////////////////////////////////////////////
+    fHistAsMcPtK0(0),
+
+//    fHistAsMcPtK0MI(0),
+    fHistAsMcPtLambda(0),
+//////////////////////////////////////////////
+
+
+//    fHistAsMcPtAntiLambdaMI(0),
+    fHistAsMcPtZoomK0(0),
+//    fHistAsMcPtZoomK0MI(0),
+    fHistAsMcPtZoomLambda(0),
+//    fHistAsMcPtZoomLambdaMI(0),
+    fHistAsMcProdRadiusK0(0),
+//    fHistAsMcProdRadiusK0MI(0),
+    fHistAsMcProdRadiusLambda(0),
+//    fHistAsMcProdRadiusLambdaMI(0),
+    fHistAsMcProdRadiusAntiLambda(0),
+//    fHistAsMcProdRadiusAntiLambdaMI(0),
+    fHistAsMcProdRadiusXvsYK0s(0),
+//    fHistAsMcProdRadiusXvsYK0sMI(0),
+    fHistAsMcProdRadiusXvsYLambda(0),
+//    fHistAsMcProdRadiusXvsYLambdaMI(0),
+    fHistAsMcProdRadiusXvsYAntiLambda(0),
+//    fHistAsMcProdRadiusXvsYAntiLambdaMI(0),
+    fHistPidMcMassK0(0),
+//    fHistPidMcMassK0MI(0),
+    fHistPidMcMassLambda(0),
+//    fHistPidMcMassLambdaMI(0),
+    fHistPidMcMassAntiLambda(0),
+//    fHistPidMcMassAntiLambdaMI(0),
+    fHistAsMcMassK0(0),
+//    fHistAsMcMassK0MI(0),
+    fHistAsMcMassLambda(0),
+//    fHistAsMcMassLambdaMI(0),
+    fHistAsMcMassAntiLambda(0),
+//    fHistAsMcMassAntiLambdaMI(0),
+    fHistAsMcPtVsMassK0(0),
+//    fHistAsMcPtVsMassK0MI(0),
+    fHistAsMcPtVsMassLambda(0),
+//    fHistAsMcPtVsMassLambdaMI(0),
+    fHistAsMcPtVsMassAntiLambda(0),
+//    fHistAsMcPtVsMassAntiLambdaMI(0),
+    fHistAsMcMassVsRadiusK0(0),
+//    fHistAsMcMassVsRadiusK0MI(0),
+    fHistAsMcMassVsRadiusLambda(0),
+//    fHistAsMcMassVsRadiusLambdaMI(0),
+    fHistAsMcMassVsRadiusAntiLambda(0),
+//    fHistAsMcMassVsRadiusAntiLambdaMI(0),
+    fHistAsMcResxK0(0),
+    fHistAsMcResyK0(0),
+    fHistAsMcReszK0(0),
+    fHistAsMcResrVsRadiusK0(0),
+    fHistAsMcReszVsRadiusK0(0),
+//    fHistAsMcResxK0MI(0),
+//    fHistAsMcResyK0MI(0),
+ //   fHistAsMcReszK0MI(0),
+//    fHistAsMcResrVsRadiusK0MI(0),
+//    fHistAsMcReszVsRadiusK0MI(0),
+    fHistAsMcResxLambda(0),
+    fHistAsMcResyLambda(0),
+    fHistAsMcReszLambda(0),
+    fHistAsMcResrVsRadiusLambda(0),
+    fHistAsMcReszVsRadiusLambda(0),
+//    fHistAsMcResxLambdaMI(0),
+//    fHistAsMcResyLambdaMI(0),
+//    fHistAsMcReszLambdaMI(0),
+//    fHistAsMcResrVsRadiusLambdaMI(0),
+//    fHistAsMcReszVsRadiusLambdaMI(0),
+    fHistAsMcResxAntiLambda(0),
+    fHistAsMcResyAntiLambda(0),
+    fHistAsMcReszAntiLambda(0),
+    fHistAsMcResrVsRadiusAntiLambda(0),
+    fHistAsMcReszVsRadiusAntiLambda(0),
+//    fHistAsMcResxAntiLambdaMI(0),
+//    fHistAsMcResyAntiLambdaMI(0),
+//    fHistAsMcReszAntiLambdaMI(0),
+//    fHistAsMcResrVsRadiusAntiLambdaMI(0),
+//    fHistAsMcReszVsRadiusAntiLambdaMI(0),
+    fHistAsMcResPtK0(0),
+//    fHistAsMcResPtK0MI(0),
+    fHistAsMcResPtLambda(0),
+//    fHistAsMcResPtLambdaMI(0),
+    fHistAsMcResPtAntiLambda(0),
+//    fHistAsMcResPtAntiLambdaMI(0),
+    fHistAsMcResPtVsRapK0(0),
+//    fHistAsMcResPtVsRapK0MI(0),
+    fHistAsMcResPtVsRapLambda(0),
+//    fHistAsMcResPtVsRapLambdaMI(0),
+    fHistAsMcResPtVsRapAntiLambda(0),
+//    fHistAsMcResPtVsRapAntiLambdaMI(0),
+    fHistAsMcResPtVsPtK0(0),
+//    fHistAsMcResPtVsPtK0MI(0),
+    fHistAsMcResPtVsPtLambda(0),
+//    fHistAsMcResPtVsPtLambdaMI(0),
+    fHistAsMcResPtVsPtAntiLambda(0),
+//    fHistAsMcResPtVsPtAntiLambdaMI(0),
+    fHistAsMcMotherPdgCodeK0s(0),
+//    fHistAsMcMotherPdgCodeK0sMI(0),
+    fHistAsMcMotherPdgCodeLambda(0),
+//    fHistAsMcMotherPdgCodeLambdaMI(0),
+    fHistAsMcMotherPdgCodeAntiLambda(0),
+//    fHistAsMcMotherPdgCodeAntiLambdaMI(0),
+    fHistAsMcPtLambdaFromSigma(0),
+//    fHistAsMcPtLambdaFromSigmaMI(0),
+    fHistAsMcPtAntiLambdaFromSigma(0),
+//    fHistAsMcPtAntiLambdaFromSigmaMI(0),
+    fHistAsMcSecondaryPtVsRapK0s(0),
+//    fHistAsMcSecondaryPtVsRapK0sMI(0),
+    fHistAsMcSecondaryPtVsRapLambda(0),
+//    fHistAsMcSecondaryPtVsRapLambdaMI(0),
+    fHistAsMcSecondaryPtVsRapAntiLambda(0),
+//    fHistAsMcSecondaryPtVsRapAntiLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusK0s(0),
+//    fHistAsMcSecondaryProdRadiusK0sMI(0),
+    fHistAsMcSecondaryProdRadiusLambda(0),
+//    fHistAsMcSecondaryProdRadiusLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusAntiLambda(0),
+//    fHistAsMcSecondaryProdRadiusAntiLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusXvsYK0s(0),
+//    fHistAsMcSecondaryProdRadiusXvsYK0sMI(0),
+    fHistAsMcSecondaryProdRadiusXvsYLambda(0),
+//    fHistAsMcSecondaryProdRadiusXvsYLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusXvsYAntiLambda(0),
+//    fHistAsMcSecondaryProdRadiusXvsYAntiLambdaMI(0),
+    fHistAsMcSecondaryMotherPdgCodeK0s(0),
+//    fHistAsMcSecondaryMotherPdgCodeK0sMI(0),
+    fHistAsMcSecondaryMotherPdgCodeLambda(0),
+//    fHistAsMcSecondaryMotherPdgCodeLambdaMI(0),
+    fHistAsMcSecondaryMotherPdgCodeAntiLambda(0),
+//    fHistAsMcSecondaryMotherPdgCodeAntiLambdaMI(0),
+    fHistAsMcSecondaryPtLambdaFromSigma(0),
+//    fHistAsMcSecondaryPtLambdaFromSigmaMI(0),
+    fHistAsMcSecondaryPtAntiLambdaFromSigma(0)
+//    fHistAsMcSecondaryPtAntiLambdaFromSigmaMI(0)
+    
+{
+  // Constructor
+
+  // New V0 cuts
+/*  fCuts[0]=33;    // max allowed chi2
+  fCuts[1]=0.05;  // min allowed impact parameter for the 1st daughter
+  fCuts[2]=0.05;  // min allowed impact parameter for the 2nd daughter
+  fCuts[3]=0.5;   // max allowed DCA between the daughter tracks
+  fCuts[4]=0.00;  // max allowed cosine of V0's pointing angle
+  fCuts[5]=0.2;   // min radius of the fiducial volume
+  fCuts[6]=100;   // max radius of the fiducial volume
+*/
+}
+
+
+
+
+
+//________________________________________________________________________
+AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange(const char *name)
+  : AliAnalysisTaskSE(name), fAnalysisMC(0), fAnalysisType("infoType"), fCollidingSystems(0), fUsePID("infoPID"), fUseCut("infocut"), fESD(0), fListHist(),fCentrSelector(0),
+    fHistMCPrimaryVertexX(0),
+    fHistMCPrimaryVertexY(0),
+    fHistMCPrimaryVertexZ(0),
+    fHistMCMultiplicityPrimary(0),
+    fHistMCMultiplicityTracks(0),
+    fHistMCtracksProdRadiusK0s(0),
+    fHistMCtracksProdRadiusLambda(0),
+    fHistMCtracksProdRadiusAntiLambda(0),
+    fHistMCtracksDecayRadiusK0s(0),
+    fHistMCtracksDecayRadiusLambda(0),
+    fHistMCtracksDecayRadiusAntiLambda(0),
+    fHistMCPtAllK0s(0),
+    fHistMCPtAllLambda(0),
+    fHistMCPtAllAntiLambda(0),
+    fHistMCProdRadiusK0s(0),
+    fHistMCProdRadiusLambda(0),
+    fHistMCProdRadiusAntiLambda(0),
+    fHistMCRapK0s(0),
+    fHistMCRapInPtRangeK0s(0),
+    fHistMCRapLambda(0),
+    fHistMCRapInPtRangeLambda(0),
+    fHistMCRapAntiLambda(0),
+    fHistMCRapInPtRangeAntiLambda(0),
+    fHistMCRapXi(0),
+    fHistMCRapInPtRangeXi(0),
+    fHistMCRapPhi(0),
+    fHistMCRapInPtRangePhi(0),
+
+////////////////////////////////////////////////
+    fHistMCPtK0s(0),
+
+    fHistMCPtLambda(0),
+/////////////////////////////////////////////////
+
+    fHistMCPtLambdaFromSigma(0),
+    fHistMCPtAntiLambdaFromSigma(0),
+    fHistNTimesRecK0s(0),
+//    fHistNTimesRecK0sMI(0),
+    fHistNTimesRecLambda(0),
+//    fHistNTimesRecLambdaMI(0),
+    fHistNTimesRecAntiLambda(0),
+//    fHistNTimesRecAntiLambdaMI(0),
+    fHistNTimesRecK0sVsPt(0),
+//    fHistNTimesRecK0sVsPtMI(0),
+    fHistNTimesRecLambdaVsPt(0),
+//    fHistNTimesRecLambdaVsPtMI(0),
+    fHistNTimesRecAntiLambdaVsPt(0),
+//    fHistNTimesRecAntiLambdaVsPtMI(0),
+    fHistNumberEvents(0),
+    fHistTrackPerEvent(0),
+    fHistTrackletPerEvent(0),
+    fHistMCDaughterTrack(0),
+    fHistSPDPrimaryVertexZ(0),
+    fHistPrimaryVertexX(0),
+    fHistPrimaryVertexY(0),
+    fHistPrimaryVertexZ(0),
+    fHistPrimaryVertexResX(0),
+    fHistPrimaryVertexResY(0),
+    fHistPrimaryVertexResZ(0),
+    fHistPrimaryVertexPosXV0events(0), 
+    fHistPrimaryVertexPosYV0events(0), 
+    fHistPrimaryVertexPosZV0events(0),
+    fHistDaughterPt(0),
+    fHistDcaPosToPrimVertex(0),
+    fHistDcaNegToPrimVertex(0),
+    fHistDcaPosToPrimVertexZoom(0),
+    fHistDcaNegToPrimVertexZoom(0),
+    fHistRadiusV0(0),
+    fHistDecayLengthV0(0),
+    fHistDcaV0Daughters(0),
+    fHistChi2(0),
+    fHistCosPointAngle(0),
+    fHistCosPointAngleZoom(0),
+    fHistProdRadius(0),
+//    fHistProdRadiusMI(0),
+    fHistV0Multiplicity(0),
+//    fHistV0MultiplicityMI(0),
+    fHistChi2KFBeforeCutK0s(0), 
+    fHistChi2KFBeforeCutLambda(0), 
+    fHistChi2KFBeforeCutAntiLambda(0),
+    fHistChi2KFAfterCutK0s(0), 
+    fHistChi2KFAfterCutLambda(0), 
+    fHistChi2KFAfterCutAntiLambda(0),
+    fHistMassK0(0),
+//    fHistMassK0MI(0),
+    fHistMassLambda(0),
+//    fHistMassLambdaMI(0),
+    fHistMassAntiLambda(0),
+//    fHistMassAntiLambdaMI(0),
+    fHistMassVsRadiusK0(0),
+//    fHistMassVsRadiusK0MI(0),
+    fHistMassVsRadiusLambda(0),
+//    fHistMassVsRadiusLambdaMI(0),
+    fHistMassVsRadiusAntiLambda(0),
+//    fHistMassVsRadiusAntiLambdaMI(0),
+
+/////////////////////////////////////////////
+    fHistPtVsMassK0(0),
+//    fHistPtVsMassK0MI(0),
+    fHistPtVsMassLambda(0),
+//    fHistPtVsMassLambdaMI(0),
+///////////////////////////////////////////////////
+       fHistPzPtBeforeK0s(0),
+       fHistPzPtAfterK0s(0),
+       fHistPzPtBeforeLambda(0),
+       fHistPzPtAfterLambda(0),
+
+//////////////////////////////////////////
+
+//    fHistMultVsPtVsMassK0(0),
+//    fHistMultVsPtVsMassLambda(0),
+//    fHistMultVsPtVsMassAntiLambda(0),
+
+
+    fHistArmenterosPodolanski(0),
+//    fHistArmenterosPodolanskiMI(0),
+    fHistNsigmaPosPionAntiLambda(0),
+    fHistNsigmaNegProtonAntiLambda(0),
+    fHistNsigmaPosProtonLambda(0),
+    fHistNsigmaNegPionLambda(0),
+    fHistNsigmaPosPionK0(0),
+    fHistNsigmaNegPionK0(0),
+    fHistAsMcRapK0(0),
+//    fHistAsMcRapK0MI(0),
+    fHistAsMcRapLambda(0),
+//    fHistAsMcRapLambdaMI(0),
+    fHistAsMcRapAntiLambda(0),
+//    fHistAsMcRapAntiLambdaMI(0),
+
+///////////////////////////////////
+    fHistAsMcPtK0(0),
+
+
+//    fHistAsMcPtK0MI(0),
+    fHistAsMcPtLambda(0),
+
+/////////////////////////////////////
+
+
+//    fHistAsMcPtAntiLambdaMI(0),
+    fHistAsMcPtZoomK0(0),
+//    fHistAsMcPtZoomK0MI(0),
+    fHistAsMcPtZoomLambda(0),
+//    fHistAsMcPtZoomLambdaMI(0),
+    fHistAsMcProdRadiusK0(0),
+//    fHistAsMcProdRadiusK0MI(0),
+    fHistAsMcProdRadiusLambda(0),
+//    fHistAsMcProdRadiusLambdaMI(0),
+    fHistAsMcProdRadiusAntiLambda(0),
+//    fHistAsMcProdRadiusAntiLambdaMI(0),
+    fHistAsMcProdRadiusXvsYK0s(0),
+//    fHistAsMcProdRadiusXvsYK0sMI(0),
+    fHistAsMcProdRadiusXvsYLambda(0),
+//    fHistAsMcProdRadiusXvsYLambdaMI(0),
+    fHistAsMcProdRadiusXvsYAntiLambda(0),
+//    fHistAsMcProdRadiusXvsYAntiLambdaMI(0),
+    fHistPidMcMassK0(0),
+//    fHistPidMcMassK0MI(0),
+    fHistPidMcMassLambda(0),
+//    fHistPidMcMassLambdaMI(0),
+    fHistPidMcMassAntiLambda(0),
+//    fHistPidMcMassAntiLambdaMI(0),
+    fHistAsMcMassK0(0),
+//    fHistAsMcMassK0MI(0),
+    fHistAsMcMassLambda(0),
+//    fHistAsMcMassLambdaMI(0),
+    fHistAsMcMassAntiLambda(0),
+//    fHistAsMcMassAntiLambdaMI(0),
+    fHistAsMcPtVsMassK0(0),
+//    fHistAsMcPtVsMassK0MI(0),
+    fHistAsMcPtVsMassLambda(0),
+//    fHistAsMcPtVsMassLambdaMI(0),
+    fHistAsMcPtVsMassAntiLambda(0),
+//    fHistAsMcPtVsMassAntiLambdaMI(0),
+    fHistAsMcMassVsRadiusK0(0),
+//    fHistAsMcMassVsRadiusK0MI(0),
+    fHistAsMcMassVsRadiusLambda(0),
+//    fHistAsMcMassVsRadiusLambdaMI(0),
+    fHistAsMcMassVsRadiusAntiLambda(0),
+//    fHistAsMcMassVsRadiusAntiLambdaMI(0),
+    fHistAsMcResxK0(0),
+    fHistAsMcResyK0(0),
+    fHistAsMcReszK0(0),
+    fHistAsMcResrVsRadiusK0(0),
+    fHistAsMcReszVsRadiusK0(0),
+//    fHistAsMcResxK0MI(0),
+//    fHistAsMcResyK0MI(0),
+ //   fHistAsMcReszK0MI(0),
+//    fHistAsMcResrVsRadiusK0MI(0),
+//    fHistAsMcReszVsRadiusK0MI(0),
+    fHistAsMcResxLambda(0),
+    fHistAsMcResyLambda(0),
+    fHistAsMcReszLambda(0),
+    fHistAsMcResrVsRadiusLambda(0),
+    fHistAsMcReszVsRadiusLambda(0),
+//    fHistAsMcResxLambdaMI(0),
+//    fHistAsMcResyLambdaMI(0),
+//    fHistAsMcReszLambdaMI(0),
+//    fHistAsMcResrVsRadiusLambdaMI(0),
+//    fHistAsMcReszVsRadiusLambdaMI(0),
+    fHistAsMcResxAntiLambda(0),
+    fHistAsMcResyAntiLambda(0),
+    fHistAsMcReszAntiLambda(0),
+    fHistAsMcResrVsRadiusAntiLambda(0),
+    fHistAsMcReszVsRadiusAntiLambda(0),
+//    fHistAsMcResxAntiLambdaMI(0),
+//    fHistAsMcResyAntiLambdaMI(0),
+//    fHistAsMcReszAntiLambdaMI(0),
+//    fHistAsMcResrVsRadiusAntiLambdaMI(0),
+//    fHistAsMcReszVsRadiusAntiLambdaMI(0),
+    fHistAsMcResPtK0(0),
+//    fHistAsMcResPtK0MI(0),
+    fHistAsMcResPtLambda(0),
+//    fHistAsMcResPtLambdaMI(0),
+    fHistAsMcResPtAntiLambda(0),
+//    fHistAsMcResPtAntiLambdaMI(0),
+    fHistAsMcResPtVsRapK0(0),
+//    fHistAsMcResPtVsRapK0MI(0),
+    fHistAsMcResPtVsRapLambda(0),
+//    fHistAsMcResPtVsRapLambdaMI(0),
+    fHistAsMcResPtVsRapAntiLambda(0),
+//    fHistAsMcResPtVsRapAntiLambdaMI(0),
+    fHistAsMcResPtVsPtK0(0),
+//    fHistAsMcResPtVsPtK0MI(0),
+    fHistAsMcResPtVsPtLambda(0),
+//    fHistAsMcResPtVsPtLambdaMI(0),
+    fHistAsMcResPtVsPtAntiLambda(0),
+//    fHistAsMcResPtVsPtAntiLambdaMI(0),
+    fHistAsMcMotherPdgCodeK0s(0),
+//    fHistAsMcMotherPdgCodeK0sMI(0),
+    fHistAsMcMotherPdgCodeLambda(0),
+//    fHistAsMcMotherPdgCodeLambdaMI(0),
+    fHistAsMcMotherPdgCodeAntiLambda(0),
+//    fHistAsMcMotherPdgCodeAntiLambdaMI(0),
+    fHistAsMcPtLambdaFromSigma(0),
+//    fHistAsMcPtLambdaFromSigmaMI(0),
+    fHistAsMcPtAntiLambdaFromSigma(0),
+//    fHistAsMcPtAntiLambdaFromSigmaMI(0),
+    fHistAsMcSecondaryPtVsRapK0s(0),
+//    fHistAsMcSecondaryPtVsRapK0sMI(0),
+    fHistAsMcSecondaryPtVsRapLambda(0),
+//    fHistAsMcSecondaryPtVsRapLambdaMI(0),
+    fHistAsMcSecondaryPtVsRapAntiLambda(0),
+//    fHistAsMcSecondaryPtVsRapAntiLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusK0s(0),
+//    fHistAsMcSecondaryProdRadiusK0sMI(0),
+    fHistAsMcSecondaryProdRadiusLambda(0),
+//    fHistAsMcSecondaryProdRadiusLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusAntiLambda(0),
+//    fHistAsMcSecondaryProdRadiusAntiLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusXvsYK0s(0),
+//    fHistAsMcSecondaryProdRadiusXvsYK0sMI(0),
+    fHistAsMcSecondaryProdRadiusXvsYLambda(0),
+//    fHistAsMcSecondaryProdRadiusXvsYLambdaMI(0),
+    fHistAsMcSecondaryProdRadiusXvsYAntiLambda(0),
+//    fHistAsMcSecondaryProdRadiusXvsYAntiLambdaMI(0),
+    fHistAsMcSecondaryMotherPdgCodeK0s(0),
+//    fHistAsMcSecondaryMotherPdgCodeK0sMI(0),
+    fHistAsMcSecondaryMotherPdgCodeLambda(0),
+//    fHistAsMcSecondaryMotherPdgCodeLambdaMI(0),
+    fHistAsMcSecondaryMotherPdgCodeAntiLambda(0),
+//    fHistAsMcSecondaryMotherPdgCodeAntiLambdaMI(0),
+    fHistAsMcSecondaryPtLambdaFromSigma(0),
+//    fHistAsMcSecondaryPtLambdaFromSigmaMI(0),
+    fHistAsMcSecondaryPtAntiLambdaFromSigma(0)
+//    fHistAsMcSecondaryPtAntiLambdaFromSigmaMI(0)
+    
+{
+  // Constructor
+
+  //New V0 cuts
+/*  fCuts[0]=33;    // max allowed chi2
+  fCuts[1]=0.05;  // min allowed impact parameter for the 1st daughter
+  fCuts[2]=0.05;  // min allowed impact parameter for the 2nd daughter
+  fCuts[3]=0.5;   // max allowed DCA between the daughter tracks
+  fCuts[4]=0.00;  // max allowed cosine of V0's pointing angle
+  fCuts[5]=0.2;   // min radius of the fiducial volume
+  fCuts[6]=100;   // max radius of the fiducial volume
+*/
+  // Define output slots only here
+  // Output slot #1 writes into a TList container
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, AliAnalysisCentralitySelector::Class());
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPerformanceStrange::UserCreateOutputObjects() 
+{
+
+  //******************
+  // Create histograms
+  //*******************
+  fListHist = new TList();
+  fListHist->SetOwner();
+
+  // Bo: tbd: condition before allocation (i.e. if (!fHistMCMultiplicityPrimary){...} for each histo...
+
+
+
+  //***************
+  // MC histograms
+  //***************
+  // Primary Vertex:
+  fHistMCPrimaryVertexX          = new TH1F("h1MCPrimaryVertexX", "MC Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
+  fListHist->Add(fHistMCPrimaryVertexX);
+
+  fHistMCPrimaryVertexY          = new TH1F("h1MCPrimaryVertexY", "MC Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
+  fListHist->Add(fHistMCPrimaryVertexY);
+
+  fHistMCPrimaryVertexZ          = new TH1F("h1MCPrimaryVertexZ", "MC Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
+  fListHist->Add(fHistMCPrimaryVertexZ);
+  
+  // Multiplicity
+  fHistMCMultiplicityPrimary           = new TH1F("h1MCMultiplicityPrimary", "MC Primary Particles;NPrimary;Count", 201, -0.5, 200.5);
+  fListHist->Add(fHistMCMultiplicityPrimary);
+
+  fHistMCMultiplicityTracks            = new TH1F("h1MCMultiplicityTracks", "MC Tracks;Ntracks;Count", 201, -0.5, 200.5);
+  fListHist->Add(fHistMCMultiplicityTracks);
+
+  // Production Radius of non-primary particles:
+  fHistMCtracksProdRadiusK0s           = new TH2F("h2MCtracksProdRadiusK0s","Non-primary MC K^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+  fListHist->Add(fHistMCtracksProdRadiusK0s);
+
+  fHistMCtracksProdRadiusLambda        = new TH2F("h2MCtracksProdRadiusLambda","Non-primary MC #Lambda^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+  fListHist->Add(fHistMCtracksProdRadiusLambda);
+
+  fHistMCtracksProdRadiusAntiLambda    = new TH2F("h2MCtracksProdRadiusAntiLambda","Non-primary MC #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+  fListHist->Add(fHistMCtracksProdRadiusAntiLambda);
+
+  // Decay Radius of non-primary particles:
+  fHistMCtracksDecayRadiusK0s          = new TH1F("h1MCtracksDecayRadiusK0s","Non-primary MC K^{0} Decay Radius;r (cm)",101,-1,100);
+  fListHist->Add(fHistMCtracksDecayRadiusK0s);
+
+  fHistMCtracksDecayRadiusLambda       = new TH1F("h1MCtracksDecayRadiusLambda","Non-primary MC #Lambda^{0} Decay Radius;r (cm)",101,-1,100);
+  fListHist->Add(fHistMCtracksDecayRadiusLambda);
+
+  fHistMCtracksDecayRadiusAntiLambda   = new TH1F("h1MCtracksDecayRadiusAntiLambda","Non-primary #bar{#Lambda}^{0} Decay Radius;r (cm)",100,1,101);
+  fListHist->Add(fHistMCtracksDecayRadiusAntiLambda);
+
+  // Pt Distribution of non-primary particles:
+  fHistMCPtAllK0s                      = new TH1F("h1MCPtAllK0s", "Non-primary MC K^{0};p_{t} (GeV/c);Counts",240,0,12);
+  fListHist->Add(fHistMCPtAllK0s);
+
+  fHistMCPtAllLambda                   = new TH1F("h1MCPtAllLambda", "Non-primary MC #Lambda^{0};p_{t} (GeV/c);Counts",240,0,12);
+  fListHist->Add(fHistMCPtAllLambda);
+
+  fHistMCPtAllAntiLambda               = new TH1F("h1MCPtAllAntiLambda", "Non-primary MC #bar{#Lambda}^{0};p_{t} (GeV/c);Counts",240,0,12);
+  fListHist->Add(fHistMCPtAllAntiLambda);
+
+  // Production Radius
+  fHistMCProdRadiusK0s                 = new TH1F("h1MCProdRadiusK0s", "MC K^{0} Production Radius;r (cm);Count", 400, -2, 2);
+  fListHist->Add(fHistMCProdRadiusK0s);
+
+  fHistMCProdRadiusLambda              = new TH1F("h1MCProdRadiusLambda", "MC #Lambda^{0} Production Radius;r (cm);Count", 400, -2, 2);
+  fListHist->Add(fHistMCProdRadiusLambda);
+
+   fHistMCProdRadiusAntiLambda         = new TH1F("h1MCProdRadiusAntiLambda", "MC #bar{#Lambda}^{0} Production Radius;r (cm);Count", 400, -2, 2);
+  fListHist->Add(fHistMCProdRadiusAntiLambda);
+
+
+  // Rapidity distribution:
+  fHistMCRapK0s                 = new TH1F("h1MCRapK0s", "K^{0};y",160,-4,4);
+  fListHist->Add(fHistMCRapK0s);
+
+  fHistMCRapInPtRangeK0s        = new TH1F("h1MCRapInPtRangeK0s", "K^{0};y",160,-4,4);
+  fListHist->Add(fHistMCRapInPtRangeK0s);
+
+  fHistMCRapLambda              = new TH1F("h1MCRapLambda", "#Lambda;y",160,-4,4);
+  fListHist->Add(fHistMCRapLambda);
+
+  fHistMCRapInPtRangeLambda     = new TH1F("h1MCRapInPtRangeLambda", "#Lambda;y",160,-4,4);
+  fListHist->Add(fHistMCRapInPtRangeLambda);
+
+  fHistMCRapAntiLambda          = new TH1F("h1MCRapAntiLambda", "#bar{#Lambda};y",160,-4,4);
+  fListHist->Add(fHistMCRapAntiLambda);
+
+  fHistMCRapInPtRangeAntiLambda = new TH1F("h1MCRapInPtRangeAntiLambda", "#bar{#Lambda};y",160,-4,4);
+  fListHist->Add(fHistMCRapInPtRangeAntiLambda);
+
+  fHistMCRapXi                  = new TH1F("h1MCRapXi", "Xi;y",160,-4,4);
+  fListHist->Add(fHistMCRapXi);
+
+  fHistMCRapInPtRangeXi         = new TH1F("h1MCRapInPtRangeXi", "Xi;y",160,-4,4);
+  fListHist->Add(fHistMCRapInPtRangeXi);
+
+  fHistMCRapPhi                  = new TH1F("h1MCRapPhi", "Phi;y",160,-4,4);
+  fListHist->Add(fHistMCRapPhi);
+
+  fHistMCRapInPtRangePhi         = new TH1F("h1MCRapInPtRangePhi", "Phi;y",160,-4,4);
+  fListHist->Add(fHistMCRapInPtRangePhi);
+
+
+  // Pt distribution:
+  fHistMCPtK0s               = new TH1F("h1MCPtK0s", "K^{0};p_{t} (GeV/c)",240,0,12);
+  fListHist->Add(fHistMCPtK0s);
+
+
+
+
+  fHistMCPtLambda            = new TH1F("h1MCPtLambda", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
+  fListHist->Add(fHistMCPtLambda);
+
+
+
+
+  // Pt distribution of Lambda coming from Sigma decay
+  fHistMCPtLambdaFromSigma      = new TH1F("h1MCPtLambdaFromSigma", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
+  fListHist->Add(fHistMCPtLambdaFromSigma);
+
+  fHistMCPtAntiLambdaFromSigma  = new TH1F("h1MCPtAntiLambdaFromSigma", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
+  fListHist->Add(fHistMCPtAntiLambdaFromSigma);
+  // Multiple reconstruction studies
+  fHistNTimesRecK0s             = new TH1F("h1NTimesRecK0s","number of times a K0s is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
+  fListHist->Add(fHistNTimesRecK0s);
+//  fHistNTimesRecK0sMI           = new TH1F("h1NTimesRecK0sMI","number of times a K0s MI is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
+//  fListHist->Add(fHistNTimesRecK0sMI);
+
+  fHistNTimesRecLambda          = new TH1F("h1NTimesRecLambda","number of times a Lambda is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
+  fListHist->Add(fHistNTimesRecLambda);
+//  fHistNTimesRecLambdaMI        = new TH1F("h1NTimesRecLambdaMI","number of times a Lambda MI is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
+//  fListHist->Add(fHistNTimesRecLambdaMI);
+
+  fHistNTimesRecAntiLambda      = new TH1F("h1NTimesRecAntiLambda","number of times an AntiLambda is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
+  fListHist->Add(fHistNTimesRecAntiLambda);
+//  fHistNTimesRecAntiLambdaMI    = new TH1F("h1NTimesRecAntiLambdaMI","number of times an AntiLambda  MI is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
+//  fListHist->Add(fHistNTimesRecAntiLambdaMI);
+
+  fHistNTimesRecK0sVsPt         = new TH2F("h2NTimesRecK0sVsPt","NTimes versus Pt, K^{0} in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
+  fListHist->Add(fHistNTimesRecK0sVsPt);
+//  fHistNTimesRecK0sVsPtMI       = new TH2F("h2NTimesRecK0sVsPtMI","NTimes versus Pt, K^{0}, on-the-fly finder, in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
+//  fListHist->Add(fHistNTimesRecK0sVsPtMI);
+
+  fHistNTimesRecLambdaVsPt      = new TH2F("h2NTimesRecLambdaVsPt","NTimes versus Pt, #Lambda^{0} in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
+  fListHist->Add(fHistNTimesRecLambdaVsPt);
+//  fHistNTimesRecLambdaVsPtMI    = new TH2F("h2NTimesRecLambdaVsPtMI","NTimes versus Pt, #Lambda^{0} on-the-fly finder in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
+//  fListHist->Add(fHistNTimesRecLambdaVsPtMI);
+
+  fHistNTimesRecAntiLambdaVsPt  = new TH2F("h2NTimesRecAntiLambdaVsPt","NTimes versus Pt, #bar{#Lambda}^{0} in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
+  fListHist->Add(fHistNTimesRecAntiLambdaVsPt);
+//  fHistNTimesRecAntiLambdaVsPtMI= new TH2F("h2NTimesRecAntiLambdaVsPtMI","NTimes versus Pt, #bar{#Lambda}^{0}, on-the-fly finder in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
+//  fListHist->Add(fHistNTimesRecAntiLambdaVsPtMI);
+
+  
+
+  //***********************************
+  // Reconstructed particles histograms
+  //***********************************
+
+  // Number of events;
+  fHistNumberEvents           = new TH1F("h1NumberEvents", "Number of events; index;Number of Events",10,0,10);
+  fListHist->Add(fHistNumberEvents);
+
+  // multiplicity
+  fHistTrackPerEvent           = new TH1F("h1TrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",1000,0,1000);
+  fListHist->Add(fHistTrackPerEvent);
+
+  fHistTrackletPerEvent       = new TH1F("h1TrackletPerEvent", "Number of tracklets;Number of tracklets per events;Number of events",1000,0,1000);
+  fListHist->Add(fHistTrackletPerEvent);
+
+  fHistMCDaughterTrack         = new TH1F("h1MCDaughterTrack","Distribution of mc id for daughters;id tags;Counts",15,0,15);
+  fListHist->Add(fHistMCDaughterTrack);
+
+  // Primary Vertex:
+  fHistSPDPrimaryVertexZ          = new TH1F("h1SPDPrimaryVertexZ", "SPD Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
+  fListHist->Add(fHistSPDPrimaryVertexZ);
+
+  fHistPrimaryVertexX          = new TH1F("h1PrimaryVertexX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
+  fListHist->Add(fHistPrimaryVertexX);
+
+  fHistPrimaryVertexY          = new TH1F("h1PrimaryVertexY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
+  fListHist->Add(fHistPrimaryVertexY);
+
+  fHistPrimaryVertexZ          = new TH1F("h1PrimaryVertexZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
+  fListHist->Add(fHistPrimaryVertexZ);
+
+
+  // Primary vertex resolution
+  fHistPrimaryVertexResX          = new TH1F("h1PrimaryVertexResX", "Primary Vertex Resolution X;Primary Vertex Resolution X (cm);Events",100,-0.25,0.25);
+  fListHist->Add(fHistPrimaryVertexResX);
+
+  fHistPrimaryVertexResY          = new TH1F("h1PrimaryVertexResY", "Primary Vertex Resolution Y;Primary Vertex Resolution Y (cm);Events",100,-0.25,0.25);
+  fListHist->Add(fHistPrimaryVertexResY);
+
+  fHistPrimaryVertexResZ          = new TH1F("h1PrimaryVertexResZ", "Primary Vertex Resolution Z;Primary Vertex Resolution Z (cm);Events",200,-0.25,0.25);
+  fListHist->Add(fHistPrimaryVertexResZ);
+  
+
+  // Primary Vertex in events with V0 candidates:
+  fHistPrimaryVertexPosXV0events       = new TH1F("h1PrimaryVertexPosXV0events", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
+  fListHist->Add(fHistPrimaryVertexPosXV0events);
+  fHistPrimaryVertexPosYV0events       = new TH1F("h1PrimaryVertexPosYV0events", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
+  fListHist->Add(fHistPrimaryVertexPosYV0events);
+  fHistPrimaryVertexPosZV0events       = new TH1F("h1PrimaryVertexPosZV0events", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20.0,20.0);
+  fListHist->Add(fHistPrimaryVertexPosZV0events);
+
+  // Daughters Pt:
+  fHistDaughterPt              = new TH2F("h2DaughterPt", "Daughter Pt;Positive Daughter Pt; Negative Daughter Pt",200,0,2,200,0,2);
+  fListHist->Add(fHistDaughterPt);
+
+  // Cut checks:
+  fHistDcaPosToPrimVertex      = new TH2F("h2DcaPosToPrimVertex", "Positive V0 daughter;dca(cm);Status",500,0,5,2,-0.5,1.5);
+  fListHist->Add(fHistDcaPosToPrimVertex);
+
+  fHistDcaNegToPrimVertex      = new TH2F("h2DcaNegToPrimVertex", "Negative V0 daughter;dca(cm);Status",500,0,5,2,-0.5,1.5);
+  fListHist->Add(fHistDcaNegToPrimVertex);
+
+  fHistDcaPosToPrimVertexZoom  = new TH2F("h2DcaPosToPrimVertexZoom", "Positive V0 daughter;dca(cm);Status",100,0,0.1,2,-0.5,1.5);
+  fListHist->Add(fHistDcaPosToPrimVertexZoom);
+
+  fHistDcaNegToPrimVertexZoom  = new TH2F("h2DcaNegToPrimVertexZoom", "Negative V0 daughter;dca(cm);Status",100,0,0.1,2,-0.5,1.5);
+  fListHist->Add(fHistDcaNegToPrimVertexZoom);
+
+  fHistRadiusV0                = new TH2F("h2RadiusV0", "Radius;Radius(cm);Status",1200,0,120,2,-0.5,1.5);
+  fListHist->Add(fHistRadiusV0);
+
+  fHistDecayLengthV0           = new TH2F("h2DecayLengthV0", "V0s decay Length;decay length(cm);Status", 240, 0, 120,2,-0.5,1.5);
+  fListHist->Add(fHistDecayLengthV0);
+
+  fHistDcaV0Daughters          = new TH2F("h2DcaV0Daughters", "DCA between daughters;dca(cm);Status", 160, 0, 4,2,-0.5,1.5);
+  fListHist->Add(fHistDcaV0Daughters);
+
+  fHistChi2                    = new TH2F("h2Chi2", "V0s chi2;chi2;Status", 33, 0, 33,2,-0.5,1.5);
+  fListHist->Add(fHistChi2);
+
+  fHistCosPointAngle           = new TH2F("h2CosPointAngle", "Cosine of V0's pointing angle", 100,0,1,2,-0.5,1.5);
+  fListHist->Add(fHistCosPointAngle);
+
+  fHistCosPointAngleZoom       = new TH2F("h2CosPointAngleZoom", "Cosine of V0's pointing angle", 100,0.9,1,2,-0.5,1.5);
+  fListHist->Add(fHistCosPointAngleZoom);
+
+  fHistProdRadius              = new TH2F("h2ProdRadius", "Production position;x (cm);y (cm)", 100,-50,50,100,-50,50);
+  fListHist->Add(fHistProdRadius);
+
+//  fHistProdRadiusMI            = new TH2F("h2ProdRadiusMI", "Production position, V0s MI;x (cm);y (cm)", 100,-50,50,100,-50,50);
+//  fListHist->Add(fHistProdRadiusMI);
+
+  // V0 Multiplicity
+  if (!fHistV0Multiplicity) {
+    if (fCollidingSystems)
+      fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of Offline V0s;Events", 200, 0, 40000);
+    else
+      fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of Offline V0s;Events", 10, 0, 10); 
+    fListHist->Add(fHistV0Multiplicity);
+  }
+
+/*  if (!fHistV0MultiplicityMI) {
+    if (fCollidingSystems)
+      fHistV0MultiplicityMI = new TH1F("fHistV0MultiplicityMI", "Multiplicity distribution;Number of On-the-fly V0s;Events", 200, 0, 40000);
+    else
+      fHistV0MultiplicityMI = new TH1F("fHistV0MultiplicityMI", "Multiplicity distribution;Number of On-the-fly V0s;Events", 10, 0, 10); 
+    fListHist->Add(fHistV0MultiplicityMI);
+  }
+*/
+  // AliKF Chi2
+  fHistChi2KFBeforeCutK0s               = new TH2F("h1Chi2KFBeforeCutK0s", "K^{0}  candidates;#Chi^{2});Counts", 250, 0, 50, 2,-0.5,1.5);
+  fListHist->Add(fHistChi2KFBeforeCutK0s);
+  fHistChi2KFBeforeCutLambda            = new TH2F("h1Chi2KFBeforeCutLambda", "#Lambda^{0}  candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
+  fListHist->Add(fHistChi2KFBeforeCutLambda);
+  fHistChi2KFBeforeCutAntiLambda        = new TH2F("h1Chi2KFBeforeCutAntiLambda", "#bar{#Lambda}^{0}  candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
+  fListHist->Add(fHistChi2KFBeforeCutAntiLambda);
+
+  fHistChi2KFAfterCutK0s               = new TH2F("h1Chi2KFAfterCutK0s", "K^{0}  candidates;#Chi^{2});Counts", 250, 0, 50, 2,-0.5,1.5);
+  fListHist->Add(fHistChi2KFAfterCutK0s);
+  fHistChi2KFAfterCutLambda            = new TH2F("h1Chi2KFAfterCutLambda", "#Lambda^{0}  candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
+  fListHist->Add(fHistChi2KFAfterCutLambda);
+  fHistChi2KFAfterCutAntiLambda        = new TH2F("h1Chi2KFAfterCutAntiLambda", "#bar{#Lambda}^{0} candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
+  fListHist->Add(fHistChi2KFAfterCutAntiLambda);
+
+  // Mass:
+  fHistMassK0                   = new TH1F("h1MassK0", "K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
+  fListHist->Add(fHistMassK0);
+//  fHistMassK0MI                 = new TH1F("h1MassK0MI", "K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
+//  fListHist->Add(fHistMassK0MI);
+
+  fHistMassLambda               = new TH1F("h1MassLambda", "#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+  fListHist->Add(fHistMassLambda);
+//  fHistMassLambdaMI             = new TH1F("h1MassLambdaMI", "#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+//  fListHist->Add(fHistMassLambdaMI);
+
+  fHistMassAntiLambda           = new TH1F("h1MassAntiLambda", "#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+  fListHist->Add(fHistMassAntiLambda);
+//  fHistMassAntiLambdaMI         = new TH1F("h1MassAntiLambdaMI", "#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+//  fListHist->Add(fHistMassAntiLambdaMI);
+
+  // invariant mass vs radius
+  const Double_t radius[10] = {0.0,2.5,2.9,3.9,7.6,15.0,23.9,37.8,42.8,100.0};
+  Int_t lNbinRadius        = 9;
+  Int_t lNbinInvMassLambda = 300;
+
+  fHistMassVsRadiusK0           = new TH2F("h2MassVsRadiusK0", "K^{0} candidates;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 200, 0.4, 0.6);
+  fListHist->Add(fHistMassVsRadiusK0);
+
+//  fHistMassVsRadiusK0MI         = new TH2F("h2MassVsRadiusK0MI", "K^{0} MI candidates;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 200, 0.4, 0.6);
+//  fListHist->Add(fHistMassVsRadiusK0MI);
+  
+  fHistMassVsRadiusLambda       = new TH2F("h2MassVsRadiusLambda", "#Lambda candidates;radius (cm);M(p#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 140, 1.06, 1.2);
+  fListHist->Add(fHistMassVsRadiusLambda);
+
+//  fHistMassVsRadiusLambdaMI     = new TH2F("h2MassVsRadiusLambdaMI", "#Lambda MI candidates;radius (cm);M(p#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 140, 1.06, 1.2);
+//  fListHist->Add(fHistMassVsRadiusLambdaMI);
+
+  fHistMassVsRadiusAntiLambda   = new TH2F("h2MassVsRadiusAntiLambda", "#bar{#Lambda} candidates;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",lNbinRadius,radius, 140, 1.06, 1.2);
+  fListHist->Add(fHistMassVsRadiusAntiLambda);
+
+//  fHistMassVsRadiusAntiLambdaMI = new TH2F("h2MassVsRadiusAntiLambdaMI", "#bar{#Lambda} candidates;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",lNbinRadius,radius, 140, 1.06, 1.2);
+//  fListHist->Add(fHistMassVsRadiusAntiLambdaMI);
+
+  // Pt Vs Mass
+  fHistPtVsMassK0               = new TH2F("h2PtVsMassK0","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
+  fListHist->Add(fHistPtVsMassK0);
+//  fHistPtVsMassK0MI             = new TH2F("h2PtVsMassK0MI","K^{0} MIcandidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
+//  fListHist->Add(fHistPtVsMassK0MI);
+
+  fHistPtVsMassLambda           = new TH2F("h2PtVsMassLambda","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
+  fListHist->Add(fHistPtVsMassLambda);
+//  fHistPtVsMassLambdaMI         = new TH2F("h2PtVsMassLambdaMI","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
+//  fListHist->Add(fHistPtVsMassLambdaMI);
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+fHistPzPtBeforeK0s     = new TH1F("h1PzPtBeforeK0s","K0s; Abs(pz/pt); count",1000,0,10);
+  fListHist->Add(fHistPzPtBeforeK0s);
+
+fHistPzPtAfterK0s      = new TH1F("h1PzPtAfterK0s","K0s; Abs(pz/pt); count",1000,0,10);
+  fListHist->Add(fHistPzPtAfterK0s);
+
+fHistPzPtBeforeLambda  = new TH1F("h1PzPtBeforeLambda","#Lambda^{0}; Abs(pz/pt); count",1000,0,10);
+  fListHist->Add(fHistPzPtBeforeLambda);
+
+fHistPzPtAfterLambda   = new TH1F("h1PzPtAfterLambda","#Lambda^{0}; Abs(pz/pt); count",1000,0,10);
+  fListHist->Add(fHistPzPtAfterLambda);
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+  fHistArmenterosPodolanski     = new TH2F("h2ArmenterosPodolanski","Armenteros-Podolanski phase space;#alpha;p_{t} arm",100,-1.0,1.0,50,0,0.5);
+//  fHistArmenterosPodolanskiMI   = new TH2F("h2ArmenterosPodolanskiMI","Armenteros-Podolanski phase space;#alpha;p_{t} arm",100,-1.0,1.0,50,0,0.5);
+
+
+  //PID
+  fHistNsigmaPosPionAntiLambda   = new TH1F("h1NsigmaPosPionAntiLambda", "Positive daughter of Antilambda;NsigmaPion;Counts",25,0,5);
+  fListHist->Add(fHistNsigmaPosPionAntiLambda);
+
+  fHistNsigmaNegProtonAntiLambda = new TH1F("h1NsigmaNegProtonAntiLambda", "Negative daughter of Antilambda;NsigmaProton;Counts",25,0,5);
+  fListHist->Add(fHistNsigmaNegProtonAntiLambda);
+  
+  fHistNsigmaPosProtonLambda     = new TH1F("h1NsigmaPosProtonLambda", "Positive daughter of Lambda;NsigmaProton;Counts",25,0,5); 
+  fListHist->Add(fHistNsigmaPosProtonLambda);
+  
+  fHistNsigmaNegPionLambda       = new TH1F("h1NsigmaNegPionLambda", "Negative daughter of Lambda;NsigmaPion;Counts",25,0,5);
+  fListHist->Add(fHistNsigmaNegPionLambda);
+  
+  fHistNsigmaPosPionK0           = new TH1F("h1NsigmaPosPionK0", "Positive daughter of K0s;NsigmaPion;Counts",25,0,5);
+  fListHist->Add(fHistNsigmaPosPionK0);
+  
+  fHistNsigmaNegPionK0           = new TH1F("h1NsigmaNegPionK0", "Negative daughter of K0s;NsigmaPion;Counts",25,0,5);
+  fListHist->Add(fHistNsigmaNegPionK0);
+
+
+  //********************************
+  // Associated particles histograms
+  //********************************
+
+  // Rap distribution
+  fHistAsMcRapK0                = new TH1F("h1AsMcRapK0", "K^{0} associated;eta;Counts", 60, -1.5, 1.5);
+  fListHist->Add(fHistAsMcRapK0);
+//  fHistAsMcRapK0MI              = new TH1F("h1AsMcRapK0MI", "K^{0} associated;eta;Counts", 60, -1.5, 1.5);
+//  fListHist->Add(fHistAsMcRapK0MI);
+
+  fHistAsMcRapLambda            = new TH1F("h1AsMcRapLambda", "#Lambda^{0} associated;eta;Counts", 60, -1.5, 1.5);
+  fListHist->Add(fHistAsMcRapLambda);
+//  fHistAsMcRapLambdaMI          = new TH1F("h1AsMcRapLambdaMI", "#Lambda^{0} associated;eta;Counts", 60, -1.5, 1.5);
+//  fListHist->Add(fHistAsMcRapLambdaMI);
+
+  fHistAsMcRapAntiLambda        = new TH1F("h1AsMcRapAntiLambda", "#bar{#Lambda}^{0} associated;eta;Counts", 60, -1.5, 1.5);
+  fListHist->Add(fHistAsMcRapAntiLambda);
+ // fHistAsMcRapAntiLambdaMI      = new TH1F("h1AsMcRapAntiLambdaMI", "#bar{#Lambda}^{0} associated;eta;Counts", 60, -1.5, 1.5);
+//  fListHist->Add(fHistAsMcRapAntiLambdaMI);
+
+
+  //Pt distribution
+  fHistAsMcPtK0                = new TH1F("h1AsMcPtK0", "K^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
+  fListHist->Add(fHistAsMcPtK0);
+
+
+//  fHistAsMcPtK0MI              = new TH1F("h1AsMcPtK0MI", "K^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
+//  fListHist->Add(fHistAsMcPtK0MI);
+
+  fHistAsMcPtLambda            = new TH1F("h1AsMcPtLambda", "#Lambda^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
+  fListHist->Add(fHistAsMcPtLambda);
+
+
+
+//  fHistAsMcPtAntiLambdaMI      = new TH1F("h1AsMcPtAntiLambdaMI", "#bar{#Lambda}^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
+//  fListHist->Add(fHistAsMcPtAntiLambdaMI);
+
+  fHistAsMcPtZoomK0            = new TH1F("h1AsMcPtZoomK0", "K^{0} candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
+  fListHist->Add(fHistAsMcPtZoomK0);
+//  fHistAsMcPtZoomK0MI          = new TH1F("h1AsMcPtZoomK0MI", "K^{0} MI candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
+//  fListHist->Add(fHistAsMcPtZoomK0MI);
+
+  fHistAsMcPtZoomLambda        = new TH1F("h1AsMcPtZoomLambda", "#Lambda^{0} candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
+  fListHist->Add(fHistAsMcPtZoomLambda);
+//  fHistAsMcPtZoomLambdaMI      = new TH1F("h1AsMcPtZoomLambdaMI", "#Lambda^{0} MI candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
+//  fListHist->Add(fHistAsMcPtZoomLambdaMI);
+
+
+  // Radius distribution
+  fHistAsMcProdRadiusK0               = new TH1F("h1AsMcProdRadiusK0", "K^{0} associated;r (cm);Counts", 500, 0, 100);
+  fListHist->Add(fHistAsMcProdRadiusK0);
+//  fHistAsMcProdRadiusK0MI             = new TH1F("h1AsMcProdRadiusK0MI", "K^{0} associated;r (cm);Counts", 500, 0, 100);
+//  fListHist->Add(fHistAsMcProdRadiusK0MI);
+
+  fHistAsMcProdRadiusLambda           = new TH1F("h1AsMcProdRadiusLambda", "#Lambda^{0} associated;r (cm);Counts", 500, 0, 100);
+  fListHist->Add(fHistAsMcProdRadiusLambda);
+//  fHistAsMcProdRadiusLambdaMI         = new TH1F("h1AsMcProdRadiusLambdaMI", "#Lambda^{0} associated;r (cm);Counts", 500, 0, 100);
+//  fListHist->Add(fHistAsMcProdRadiusLambdaMI);
+
+  fHistAsMcProdRadiusAntiLambda       = new TH1F("h1AsMcProdRadiusAntiLambda", "#bar{#Lambda}^{0} associated;r (cm);Counts", 500, 0, 100);
+  fListHist->Add(fHistAsMcProdRadiusAntiLambda);
+//  fHistAsMcProdRadiusAntiLambdaMI     = new TH1F("h1AsMcProdRadiusAntiLambdaMI", "#bar{#Lambda}^{0} associated;r (cm);Counts", 500, 0, 100);
+//  fListHist->Add(fHistAsMcProdRadiusAntiLambdaMI);
+
+  fHistAsMcProdRadiusXvsYK0s          = new TH2F("h2AsMcProdRadiusXvsYK0s","Associated Secondary K^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+  fListHist->Add(fHistAsMcProdRadiusXvsYK0s);
+//  fHistAsMcProdRadiusXvsYK0sMI        = new TH2F("h2AsMcProdRadiusXvsYK0sMI","Associated Secondary K^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+//  fListHist->Add(fHistAsMcProdRadiusXvsYK0sMI);
+
+  fHistAsMcProdRadiusXvsYLambda       = new TH2F("h2AsMcProdRadiusXvsYLambda","Associated Secondary #Lambda^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+  fListHist->Add(fHistAsMcProdRadiusXvsYLambda);
+//  fHistAsMcProdRadiusXvsYLambdaMI     = new TH2F("h2AsMcProdRadiusXvsYLambdaMI","Associated Secondary #Lambda^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+//  fListHist->Add(fHistAsMcProdRadiusXvsYLambdaMI);
+
+  fHistAsMcProdRadiusXvsYAntiLambda   = new TH2F("h2AsMcProdRadiusXvsYAntiLambda","Associated Secondary #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+  fListHist->Add(fHistAsMcProdRadiusXvsYAntiLambda);
+//  fHistAsMcProdRadiusXvsYAntiLambdaMI = new TH2F("h2AsMcProdRadiusXvsYAntiLambdaMI","Associated Secondary #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
+//  fListHist->Add(fHistAsMcProdRadiusXvsYAntiLambdaMI);
+
+
+
+  // Mass
+  fHistPidMcMassK0             = new TH1F("h1PidMcMassK0", "K^{0} MC PId checked;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
+  fListHist->Add(fHistPidMcMassK0);
+//  fHistPidMcMassK0MI           = new TH1F("h1PidMcMassK0MI", "K^{0} MC PId checked;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
+//  fListHist->Add(fHistPidMcMassK0MI);
+
+  fHistPidMcMassLambda         = new TH1F("h1PidMcMassLambda", "#Lambda^{0} MC PId checked;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+  fListHist->Add(fHistPidMcMassLambda);
+//  fHistPidMcMassLambdaMI       = new TH1F("h1PidMcMassLambdaMI", "#Lambda^{0} MC PId checked;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+//  fListHist->Add(fHistPidMcMassLambdaMI);
+  
+  fHistPidMcMassAntiLambda     = new TH1F("h1PidMcMassAntiLambda", "#bar{#Lambda}^{0} MC PId checked;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+  fListHist->Add(fHistPidMcMassAntiLambda);
+//  fHistPidMcMassAntiLambdaMI   = new TH1F("h1PidMcMassAntiLambdaMI", "#bar{#Lambda}^{0} MC PId checked;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+//  fListHist->Add(fHistPidMcMassAntiLambdaMI);
+
+  fHistAsMcMassK0              = new TH1F("h1AsMcMassK0", "K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
+  fListHist->Add(fHistAsMcMassK0);
+//  fHistAsMcMassK0MI            = new TH1F("h1AsMcMassK0MI", "K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
+//  fListHist->Add(fHistAsMcMassK0MI);
+  
+  fHistAsMcMassLambda          = new TH1F("h1AsMcMassLambda", "#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+  fListHist->Add(fHistAsMcMassLambda);
+//  fHistAsMcMassLambdaMI        = new TH1F("h1AsMcMassLambdaMI", "#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+//  fListHist->Add(fHistAsMcMassLambdaMI);
+
+  fHistAsMcMassAntiLambda      = new TH1F("h1AsMcMassAntiLambda", "#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+  fListHist->Add(fHistAsMcMassAntiLambda);
+//  fHistAsMcMassAntiLambdaMI    = new TH1F("h1AsMcMassAntiLambdaMI", "#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
+//  fListHist->Add(fHistAsMcMassAntiLambdaMI);
+
+  //Pt versus Mass
+  fHistAsMcPtVsMassK0               = new TH2F("h2AsMcPtVsMassK0","K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
+  fListHist->Add(fHistAsMcPtVsMassK0);
+//  fHistAsMcPtVsMassK0MI             = new TH2F("h2AsMcPtVsMassK0MI","K^{0} MIassociated;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
+//  fListHist->Add(fHistAsMcPtVsMassK0MI);
+
+  fHistAsMcPtVsMassLambda           = new TH2F("h2AsMcPtVsMassLambda","#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
+  fListHist->Add(fHistAsMcPtVsMassLambda);
+//  fHistAsMcPtVsMassLambdaMI         = new TH2F("h2AsMcPtVsMassLambdaMI","#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
+//  fListHist->Add(fHistAsMcPtVsMassLambdaMI);
+
+  fHistAsMcPtVsMassAntiLambda       = new TH2F("h2AsMcPtVsMassAntiLambda","#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
+  fListHist->Add(fHistAsMcPtVsMassAntiLambda);
+//  fHistAsMcPtVsMassAntiLambdaMI     = new TH2F("h2AsMcPtVsMassAntiLambdaMI","#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
+//  fListHist->Add(fHistAsMcPtVsMassAntiLambdaMI);
+
+
+  // invariant mass vs radius
+  fHistAsMcMassVsRadiusK0             = new TH2F("h2AsMcMassVsRadiusK0", "K^{0} associated;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 500, 0.47, 0.52);
+  fListHist->Add(fHistAsMcMassVsRadiusK0);
+
+//  fHistAsMcMassVsRadiusK0MI           = new TH2F("h2AsMcMassVsRadiusK0MI", "K^{0} MI associated;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 500, 0.47, 0.52);
+//  fListHist->Add(fHistAsMcMassVsRadiusK0MI);
+  
+  fHistAsMcMassVsRadiusLambda         = new TH2F("h2AsMcMassVsRadiusLambda", "#Lambda associated;radius (cm);M(p#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, lNbinInvMassLambda, 1.10, 1.13);
+  fListHist->Add(fHistAsMcMassVsRadiusLambda);
+
+//  fHistAsMcMassVsRadiusLambdaMI       = new TH2F("h2AsMcMassVsRadiusLambdaMI", "#Lambda MI associated;radius (cm);M(p#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, lNbinInvMassLambda, 1.10, 1.13);
+//  fListHist->Add(fHistAsMcMassVsRadiusLambdaMI);
+
+  fHistAsMcMassVsRadiusAntiLambda     = new TH2F("h2AsMcMassVsRadiusAntiLambda", "#bar{#Lambda} associated;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",lNbinRadius,radius,lNbinInvMassLambda , 1.10, 1.13);
+  fListHist->Add(fHistAsMcMassVsRadiusAntiLambda);
+  
+//  fHistAsMcMassVsRadiusAntiLambdaMI   = new TH2F("h2AsMcMassVsRadiusAntiLambdaMI", "#bar{#Lambda} MI associated;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",lNbinRadius,radius,lNbinInvMassLambda , 1.10, 1.13);
+//  fListHist->Add(fHistAsMcMassVsRadiusAntiLambdaMI);
+    
+
+  // Position Resolution
+  fHistAsMcResxK0                     = new TH1F("h1AsMcResxK0", "K^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResxK0);
+  fHistAsMcResyK0                     = new TH1F("h1AsMcResyK0", "K^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResyK0);
+  fHistAsMcReszK0                     = new TH1F("h1AsMcReszK0", "K^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcReszK0);
+  fHistAsMcResrVsRadiusK0             = new TH2F("h2AsMcResrVsRadiusK0", "K^{0} associated;Radius (cm);#Delta r (cm)",200,0.0,50., 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResrVsRadiusK0);
+  fHistAsMcReszVsRadiusK0             = new TH2F("h2AsMcReszVsRadiusK0", "K^{0} associated;Radius (cm);#Delta z (cm)",200,0.0,50.0, 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcReszVsRadiusK0);
+
+//  fHistAsMcResxK0MI                   = new TH1F("h1AsMcResxK0MI", "K^{0} MI associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResxK0MI);
+//  fHistAsMcResyK0MI                   = new TH1F("h1AsMcResyK0MI", "K^{0} MI associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResyK0MI);
+//  fHistAsMcReszK0MI                   = new TH1F("h1AsMcReszK0MI", "K^{0} MI associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcReszK0MI);
+//  fHistAsMcResrVsRadiusK0MI           = new TH2F("h2AsMcResrVsRadiusK0MI", "K^{0} MI associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResrVsRadiusK0MI);
+//  fHistAsMcReszVsRadiusK0MI           = new TH2F("h2AsMcReszVsRadiusK0MI", "K^{0} MI associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcReszVsRadiusK0MI);
+
+  fHistAsMcResxLambda                 = new TH1F("h1AsMcResxLambda", "#Lambda^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResxLambda);
+  fHistAsMcResyLambda                 = new TH1F("h1AsMcResyLambda", "#Lambda^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResyLambda);
+  fHistAsMcReszLambda                 = new TH1F("h1AsMcReszLambda", "#Lambda^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcReszLambda);
+  fHistAsMcResrVsRadiusLambda         = new TH2F("h2AsMcResrVsRadiusLambda", "#Lambda^{0} associated;Radius (cm);#Delta r (cm)",200,0.0,50.0, 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResrVsRadiusLambda);
+  fHistAsMcReszVsRadiusLambda         = new TH2F("h2AsMcReszVsRadiusLambda", "#Lambda^{0} associated;Radius (cm);#Delta z (cm)",200,0.0,50.0, 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcReszVsRadiusLambda);
+
+//  fHistAsMcResxLambdaMI               = new TH1F("h1AsMcResxLambdaMI", "#Lambda^{0} MI associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResxLambdaMI);
+//  fHistAsMcResyLambdaMI               = new TH1F("h1AsMcResyLambdaMI", "#Lambda^{0} MI associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResyLambdaMI);
+//  fHistAsMcReszLambdaMI               = new TH1F("h1AsMcReszLambdaMI", "#Lambda^{0} MI associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcReszLambdaMI);
+//  fHistAsMcResrVsRadiusLambdaMI       = new TH2F("h2AsMcResrVsRadiusLambdaMI", "#Lambda^{0} MI associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResrVsRadiusLambdaMI);
+//  fHistAsMcReszVsRadiusLambdaMI       = new TH2F("h2AsMcReszVsRadiusLambdaMI", "#Lambda^{0} MI associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcReszVsRadiusLambdaMI);
+
+  fHistAsMcResxAntiLambda             = new TH1F("h1AsMcResxAntiLambda", "#bar{#Lambda}^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResxAntiLambda);
+  fHistAsMcResyAntiLambda             = new TH1F("h1AsMcResyAntiLambda", "#bar{#Lambda}^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResyAntiLambda);
+  fHistAsMcReszAntiLambda             = new TH1F("h1AsMcReszAntiLambda", "#bar{#Lambda}^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcReszAntiLambda);
+  fHistAsMcResrVsRadiusAntiLambda     = new TH2F("h2AsMcResrVsRadiusAntiLambda", "#bar{#Lambda}^{0} associated;Radius (cm);#Delta r (cm)",200,0.0,50.0, 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcResrVsRadiusAntiLambda);
+  fHistAsMcReszVsRadiusAntiLambda     = new TH2F("h2AsMcReszVsRadiusAntiLambda", "#bar{#Lambda}^{0} associated;Radius (cm);#Delta z (cm)",200,0.0,50.0, 50, -0.25, 0.25);
+  fListHist->Add(fHistAsMcReszVsRadiusAntiLambda);
+
+//  fHistAsMcResxAntiLambdaMI           = new TH1F("h1AsMcResxAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResxAntiLambdaMI);
+//  fHistAsMcResyAntiLambdaMI           = new TH1F("h1AsMcResyAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResyAntiLambdaMI);
+//  fHistAsMcReszAntiLambdaMI           = new TH1F("h1AsMcReszAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcReszAntiLambdaMI);
+//  fHistAsMcResrVsRadiusAntiLambdaMI   = new TH2F("h2AsMcResrVsRadiusAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcResrVsRadiusAntiLambdaMI);
+//  fHistAsMcReszVsRadiusAntiLambdaMI   = new TH2F("h2AsMcReszVsRadiusAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
+//  fListHist->Add(fHistAsMcReszVsRadiusAntiLambdaMI);
+
+  // Pt Resolution
+  fHistAsMcResPtK0                   = new TH1F("h1AsMcResPtK0","Pt Resolution K^{0};#Delta Pt;Counts",200,-1,1);
+  fListHist->Add(fHistAsMcResPtK0);
+//  fHistAsMcResPtK0MI                 = new TH1F("h1AsMcResPtK0MI","Pt Resolution K^{0} MI;#Delta Pt;Counts",200,-1,1);
+//  fListHist->Add(fHistAsMcResPtK0MI);
+  
+  fHistAsMcResPtLambda               = new TH1F("h1AsMcResPtLambda","Pt Resolution #Lambda^{0};#Delta Pt;Counts",200,-1,1);
+  fListHist->Add(fHistAsMcResPtLambda);
+//  fHistAsMcResPtLambdaMI             = new TH1F("h1AsMcResPtLambdaMI","Pt Resolution #Lambda^{0} MI;#Delta Pt;Counts",200,-1,1);
+//  fListHist->Add(fHistAsMcResPtLambdaMI);
+
+  fHistAsMcResPtAntiLambda           = new TH1F("h1AsMcResPtAntiLambda","Pt Resolution #bar{#Lambda}^{0};#Delta Pt;Counts",200,-1,1);
+  fListHist->Add(fHistAsMcResPtAntiLambda);
+//  fHistAsMcResPtAntiLambdaMI         = new TH1F("h1AsMcResPtAntiLambdaMI","Pt Resolution #bar{#Lambda}^{0} MI;#Delta Pt;Counts",200,-1,1);
+///  fListHist->Add(fHistAsMcResPtAntiLambdaMI);
+
+
+  fHistAsMcResPtVsRapK0              = new TH2F("h2AsMcResPtVsRapK0","Pt Resolution K^{0};#Delta Pt;Rap",200,-1,1,20,-1,1);
+  fListHist->Add(fHistAsMcResPtVsRapK0);
+//  fHistAsMcResPtVsRapK0MI            = new TH2F("h2AsMcResPtVsRapK0MI","Pt Resolution K^{0} MI;#Delta Pt;Rap",200,-1,1,20,-1,1);
+//  fListHist->Add(fHistAsMcResPtVsRapK0MI);
+  
+  fHistAsMcResPtVsRapLambda          = new TH2F("h2AsMcResPtVsRapLambda","Pt Resolution #Lambda^{0};#Delta Pt;Rap",200,-1,1,20,-1,1);
+  fListHist->Add(fHistAsMcResPtVsRapLambda);
+//  fHistAsMcResPtVsRapLambdaMI        = new TH2F("h2AsMcResPtVsRapLambdaMI","Pt Resolution #Lambda^{0} MI;#Delta Pt;Rap",200,-1,1,20,-1,1);
+//  fListHist->Add(fHistAsMcResPtVsRapLambdaMI);
+
+  fHistAsMcResPtVsRapAntiLambda      = new TH2F("h2AsMcResPtVsRapAntiLambda","Pt Resolution #bar{#Lambda}^{0};#Delta Pt;Rap",200,-1,1,20,-1,1);
+  fListHist->Add(fHistAsMcResPtVsRapAntiLambda);
+//  fHistAsMcResPtVsRapAntiLambdaMI    = new TH2F("h2AsMcResPtVsRapAntiLambdaMI","Pt Resolution #bar{#Lambda}^{0} MI;#Delta Pt;Rap",200,-1,1,20,-1,1);
+//  fListHist->Add(fHistAsMcResPtVsRapAntiLambdaMI);
+
+
+  fHistAsMcResPtVsPtK0               = new TH2F("h2AsMcResPtVsPtK0","Pt Resolution K^{0};#Delta Pt;Pt",600,-0.15,0.15,240,0,12);
+  fListHist->Add(fHistAsMcResPtVsPtK0);
+//  fHistAsMcResPtVsPtK0MI             = new TH2F("h2AsMcResPtVsPtK0MI","Pt Resolution K^{0} MI;#Delta Pt;Pt",600,-0.15,0.15,240,0,12);
+//  fListHist->Add(fHistAsMcResPtVsPtK0MI);
+    
+  fHistAsMcResPtVsPtLambda           = new TH2F("h2AsMcResPtVsPtLambda","Pt Resolution #Lambda^{0};#Delta Pt;Pt",600,-0.15,0.15,240,0,12);
+  fListHist->Add(fHistAsMcResPtVsPtLambda);
+//  fHistAsMcResPtVsPtLambdaMI         = new TH2F("h2AsMcResPtVsPtLambdaMI","Pt Resolution #Lambda^{0} MI;#Delta Pt;Pt",600,-0.15,0.15,240,0,12);
+//  fListHist->Add(fHistAsMcResPtVsPtLambdaMI);
+
+  fHistAsMcResPtVsPtAntiLambda       = new TH2F("h2AsMcResPtVsPtAntiLambda","Pt Resolution #bar{#Lambda}^{0};#Delta Pt;Pt",300,-0.15,0.15,240,0,12);
+  fListHist->Add(fHistAsMcResPtVsPtAntiLambda);
+//  fHistAsMcResPtVsPtAntiLambdaMI     = new TH2F("h2AsMcResPtVsPtAntiLambdaMI","Pt Resolution #bar{#Lambda}^{0} MI;#Delta Pt;Pt",300,-0.15,0.15,240,0,12);
+//  fListHist->Add(fHistAsMcResPtVsPtAntiLambdaMI);
+
+
+  // pdgcode of mother
+  fHistAsMcMotherPdgCodeK0s           = new TH1F("h1AsMcMotherPdgCodeK0s","Mother of Associated K^{0};mother;counts",11,0,11);
+  fListHist->Add(fHistAsMcMotherPdgCodeK0s);
+//  fHistAsMcMotherPdgCodeK0sMI         = new TH1F("h1AsMcMotherPdgCodeK0sMI","Mother of Associated K^{0} MI;mother;counts",11,0,11);
+//  fListHist->Add(fHistAsMcMotherPdgCodeK0sMI);
+
+  fHistAsMcMotherPdgCodeLambda        = new TH1F("h1AsMcMotherPdgCodeLambda","Mother of Associated #Lambda^{0};mother;counts",11,0,11);
+  fListHist->Add(fHistAsMcMotherPdgCodeLambda);
+//  fHistAsMcMotherPdgCodeLambdaMI      = new TH1F("h1AsMcMotherPdgCodeLambdaMI","Mother of Associated #Lambda^{0} MI;mother;counts",11,0,11);
+//  fListHist->Add(fHistAsMcMotherPdgCodeLambdaMI);
+
+  fHistAsMcMotherPdgCodeAntiLambda    = new TH1F("h1AsMcMotherPdgCodeAntiLambda","Mother of Associated #bar{#Lambda}^{0};mother;counts",11,0,11);
+  fListHist->Add(fHistAsMcMotherPdgCodeAntiLambda);
+//  fHistAsMcMotherPdgCodeAntiLambdaMI  = new TH1F("h1AsMcMotherPdgCodeAntiLambdaMI","Mother of Associated #bar{Lambda}^{0} MI;mother;counts",11,0,11);
+//  fListHist->Add(fHistAsMcMotherPdgCodeAntiLambdaMI);
+
+
+  // Pt distribution Lambda from Sigma
+  fHistAsMcPtLambdaFromSigma          = new TH1F("h1AsMcPtLambdaFromSigma","#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+  fListHist->Add(fHistAsMcPtLambdaFromSigma);
+//  fHistAsMcPtLambdaFromSigmaMI        = new TH1F("h1AsMcPtLambdaFromSigmaMI","#Lambda^{0} MI associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+//  fListHist->Add(fHistAsMcPtLambdaFromSigmaMI);
+
+  fHistAsMcPtAntiLambdaFromSigma      = new TH1F("h1AsMcPtAntiLambdaFromSigma","#bar{#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+  fListHist->Add(fHistAsMcPtAntiLambdaFromSigma);
+//  fHistAsMcPtAntiLambdaFromSigmaMI    = new TH1F("h1AsMcPtAntiLambdaFromSigmaMI","#bar{#Lambda}^{0} MI associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+//  fListHist->Add(fHistAsMcPtAntiLambdaFromSigmaMI);
+
+
+  // Associated secondary particles:
+  // Pt and rapidity distribution
+  fHistAsMcSecondaryPtVsRapK0s          = new TH2F("h2AsMcSecondaryPtVsRapK0s", "K^{0} associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
+  fListHist->Add(fHistAsMcSecondaryPtVsRapK0s);
+//  fHistAsMcSecondaryPtVsRapK0sMI        = new TH2F("h2AsMcSecondaryPtVsRapK0sMI", "K^{0} MI associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
+//  fListHist->Add(fHistAsMcSecondaryPtVsRapK0sMI);
+
+  fHistAsMcSecondaryPtVsRapLambda       = new TH2F("h2AsMcSecondaryPtVsRapLambda", "#Lambda^{0} associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
+  fListHist->Add(fHistAsMcSecondaryPtVsRapLambda);
+//  fHistAsMcSecondaryPtVsRapLambdaMI     = new TH2F("h2AsMcSecondaryPtVsRapLambdaMI", "#Lambda^{0} MI associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
+//  fListHist->Add(fHistAsMcSecondaryPtVsRapLambdaMI);
+
+  fHistAsMcSecondaryPtVsRapAntiLambda   = new TH2F("h2AsMcSecondaryPtVsRapAntiLambda", "#bar{#Lambda}^{0} associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
+  fListHist->Add(fHistAsMcSecondaryPtVsRapAntiLambda);
+//  fHistAsMcSecondaryPtVsRapAntiLambdaMI = new TH2F("h2AsMcSecondaryPtVsRapAntiLambdaMI", "#bar{#Lambda}^{0} MI associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
+//  fListHist->Add(fHistAsMcSecondaryPtVsRapAntiLambdaMI);
+
+  // Production radius
+  fHistAsMcSecondaryProdRadiusK0s              = new TH1F("h1AsMcSecondaryProdRadiusK0s", "K^{0} Production Radius;r (cm);Count", 170, -2, 15);
+  fListHist->Add(fHistAsMcSecondaryProdRadiusK0s);
+//  fHistAsMcSecondaryProdRadiusK0sMI            = new TH1F("h1AsMcSecondaryProdRadiusK0sMI", "K^{0} MI Production Radius;r (cm);Count", 170, -2, 15);
+//  fListHist->Add(fHistAsMcSecondaryProdRadiusK0sMI);
+
+  fHistAsMcSecondaryProdRadiusLambda           = new TH1F("h1AsMcSecondaryProdRadiusLambda", "#Lambda^{0} Production Radius;r (cm);Count", 170, -2, 15);
+  fListHist->Add(fHistAsMcSecondaryProdRadiusLambda);
+//  fHistAsMcSecondaryProdRadiusLambdaMI         = new TH1F("h1AsMcSecondaryProdRadiusLambdaMI", "#Lambda^{0} MI Production Radius;r (cm);Count", 170, -2, 15);
+//  fListHist->Add(fHistAsMcSecondaryProdRadiusLambdaMI);
+
+  fHistAsMcSecondaryProdRadiusAntiLambda       = new TH1F("h1AsMcSecondaryProdRadiusAntiLambda", "#bar{#Lambda}^{0} Production Radius;r (cm);Count", 170, -2, 15);
+  fListHist->Add(fHistAsMcSecondaryProdRadiusAntiLambda);  
+//  fHistAsMcSecondaryProdRadiusAntiLambdaMI     = new TH1F("h1AsMcSecondaryProdRadiusAntiLambdaMI", "#bar{#Lambda}^{0} MI Production Radius;r (cm);Count", 170, -2, 15);
+//  fListHist->Add(fHistAsMcSecondaryProdRadiusAntiLambdaMI);
+
+  fHistAsMcSecondaryProdRadiusXvsYK0s          = new TH2F("h2AsMcSecondaryProdRadiusXvsYK0s","Associated Secondary K^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
+  fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYK0s);
+//  fHistAsMcSecondaryProdRadiusXvsYK0sMI        = new TH2F("h2AsMcSecondaryProdRadiusXvsYK0sMI","Associated Secondary K^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
+//  fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYK0sMI);
+
+  fHistAsMcSecondaryProdRadiusXvsYLambda       = new TH2F("h2AsMcSecondaryProdRadiusXvsYLambda","Associated Secondary #Lambda^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
+  fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYLambda);
+//  fHistAsMcSecondaryProdRadiusXvsYLambdaMI     = new TH2F("h2AsMcSecondaryProdRadiusXvsYLambdaMI","Associated Secondary #Lambda^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
+//  fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYLambdaMI);
+
+  fHistAsMcSecondaryProdRadiusXvsYAntiLambda   = new TH2F("h2AsMcSecondaryProdRadiusXvsYAntiLambda","Associated Secondary #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
+  fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYAntiLambda);
+//  fHistAsMcSecondaryProdRadiusXvsYAntiLambdaMI = new TH2F("h2AsMcSecondaryProdRadiusXvsYAntiLambdaMI","Associated Secondary #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
+//  fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYAntiLambdaMI);
+
+  fHistAsMcSecondaryMotherPdgCodeK0s           = new TH1F("h1AsMcSecondaryMotherPdgCodeK0s","Mother of Associated Secondary K^{0};mother;counts",11,0,11);
+  fListHist->Add(fHistAsMcSecondaryMotherPdgCodeK0s);
+//  fHistAsMcSecondaryMotherPdgCodeK0sMI         = new TH1F("h1AsMcSecondaryMotherPdgCodeK0sMI","Mother of Associated Secondary K^{0} MI;mother;counts",11,0,11);
+//  fListHist->Add(fHistAsMcSecondaryMotherPdgCodeK0sMI);
+
+  fHistAsMcSecondaryMotherPdgCodeLambda        = new TH1F("h1AsMcSecondaryMotherPdgCodeLambda","Mother of Associated Secondary #Lambda^{0};mother;counts",11,0,11);
+  fListHist->Add(fHistAsMcSecondaryMotherPdgCodeLambda);
+//  fHistAsMcSecondaryMotherPdgCodeLambdaMI      = new TH1F("h1AsMcSecondaryMotherPdgCodeLambdaMI","Mother of Associated Secondary #Lambda^{0} MI;mother;counts",11,0,11);
+//  fListHist->Add(fHistAsMcSecondaryMotherPdgCodeLambdaMI);
+
+  fHistAsMcSecondaryMotherPdgCodeAntiLambda    = new TH1F("h1AsMcSecondaryMotherPdgCodeAntiLambda","Mother of Associated Secondary #bar{#Lambda}^{0};mother;counts",11,0,11);
+  fListHist->Add(fHistAsMcSecondaryMotherPdgCodeAntiLambda);
+//  fHistAsMcSecondaryMotherPdgCodeAntiLambdaMI  = new TH1F("h1AsMcSecondaryMotherPdgCodeAntiLambdaMI","Mother of Associated Secondary #bar{Lambda}^{0} MI;mother;counts",11,0,11);
+//  fListHist->Add(fHistAsMcSecondaryMotherPdgCodeAntiLambdaMI);
+
+  // Pt distribution Lambda from Sigma
+  fHistAsMcSecondaryPtLambdaFromSigma          = new TH1F("h1AsMcSecondaryPtLambdaFromSigma","#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+  fListHist->Add(fHistAsMcSecondaryPtLambdaFromSigma);
+//  fHistAsMcSecondaryPtLambdaFromSigmaMI        = new TH1F("h1AsMcSecondaryPtLambdaFromSigmaMI","#Lambda^{0} MI associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+//  fListHist->Add(fHistAsMcSecondaryPtLambdaFromSigmaMI);
+
+  fHistAsMcSecondaryPtAntiLambdaFromSigma      = new TH1F("h1AsMcSecondaryPtAntiLambdaFromSigma","#bar{#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+  fListHist->Add(fHistAsMcSecondaryPtAntiLambdaFromSigma);
+//  fHistAsMcSecondaryPtAntiLambdaFromSigmaMI    = new TH1F("h1AsMcSecondaryPtAntiLambdaFromSigmaMI","#bar{#Lambda}^{0} MI associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
+//  fListHist->Add(fHistAsMcSecondaryPtAntiLambdaFromSigmaMI);
+
+  PostData(1, fListHist);
+  PostData(2, fCentrSelector);
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPerformanceStrange::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+  AliStack* stack = NULL;
+  TClonesArray *mcArray = NULL;
+  TArrayF mcPrimaryVtx;
+
+  fESD=(AliESDEvent *)InputEvent();
+
+  if (!fESD) {
+    Printf("ERROR: fESD not available");
+    return;
+  }
+
+
+  AliVEvent* lEvent = InputEvent();
+  
+  if (!lEvent) {
+    Printf("ERROR: Event not available");
+    return;
+  }
+
+  fHistNumberEvents->Fill(0.5);
+
+  //******************
+  // Trigger Selection ! Warning Works only for ESD, add protection in case of AOD loop
+  //******************
+
+  Bool_t isSelected = 
+    (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() 
+     & AliVEvent::kMB);
+  if (!isSelected) return;
+  
+  // Centrality selection
+  static AliESDtrackCuts * trackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); // FIXME: make it a data member
+  Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,trackCuts); // FIXME esd track cuts if needed  
+  if(!isCentralitySelected) return;
+  // FIXME: add to hist number events another entry for centrality.
+
+ // Done by the AliPhysicsSelection Task ! Only the selected events are passed to this task
+
+  fHistNumberEvents->Fill(1.5);
+
+
+  //*************************
+  //End track multiplicity
+  //*************************
+
+
+  // Remove Events with no tracks
+  //if (!(fESD->GetNumberOfTracks()))  return;
+
+  fHistNumberEvents->Fill(2.5);
+  fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks());
+
+  //*************************************
+  // Cut used:
+  //*************************************
+      
+  // Cut Rapidity:
+  Double_t lCutRap  = 0.75;
+
+  // Cut AliKF Chi2 for Reconstructed particles
+  Double_t cutChi2KF  = 1E3;
+
+  // If PID is used:
+  Double_t lLimitPPID    = 0.7;
+  Float_t cutNSigmaLowP  = 1E3;
+  Float_t cutNSigmaHighP = 1E3;
+  if (fUsePID.Contains("withPID")) {
+    cutNSigmaLowP  = 5.0;
+    cutNSigmaHighP = 3.0;
+  }
+
+
+  // Cut Daughters pt (GeV/c):
+  Double_t cutMinPtDaughter = 0.160;
+
+  // Cut primary vertex:
+  Double_t cutPrimVertex = 10.0;
+
+  // Min number of TPC clusters:
+  Int_t nbMinTPCclusters = 80;
+
+  //*******************
+  // PID parameters:
+  //*******************
+      
+  Double_t fAlephParameters[5] = {0,0,0,0,0,};
+
+  fAlephParameters[0] = 0.0283086;
+  fAlephParameters[1] = 2.63394e+01;
+  fAlephParameters[2] = 5.04114e-11;
+  fAlephParameters[3] = 2.12543e+00;
+  fAlephParameters[4] = 4.88663e+00; 
+
+
+  //*******************
+  // Access MC:
+  //*******************
+  if (fAnalysisMC) {
+    if(fAnalysisType == "ESD") {
+      AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+      if (!eventHandler) {
+       Printf("ERROR: Could not retrieve MC event handler");
+       return;
+      }    
+      AliMCEvent* mcEvent = eventHandler->MCEvent();
+      if (!mcEvent) {
+       Printf("ERROR: Could not retrieve MC event");
+       return;
+      }    
+      stack = mcEvent->Stack();
+      if (!stack) {
+       Printf("ERROR: Could not retrieve stack");
+       return;
+      }
+      
+      AliGenEventHeader* mcHeader=mcEvent->GenEventHeader();
+      if(!mcHeader) return;
+      mcHeader->PrimaryVertex(mcPrimaryVtx);
+      
+    }
+    
+    else if(fAnalysisType == "AOD") {
+      
+      // load MC particles
+      mcArray = (TClonesArray*)fESD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+      if(!mcArray) {
+       Printf("strange analysis::UserExec: MC particles branch not found!\n");
+       return;
+      }
+      
+      // load MC header
+      AliAODMCHeader *mcHeader = 
+       (AliAODMCHeader*)fESD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+      if(!mcHeader) {
+       Printf("strange analysis::UserExec: MC header branch not found!\n");
+       return;
+      }
+    }
+
+    // PID parameters for MC simulations:
+    fAlephParameters[0] = 2.15898e+00/50.;
+    fAlephParameters[1] = 1.75295e+01;
+    fAlephParameters[2] = 3.40030e-09;
+    fAlephParameters[3] = 1.96178e+00;
+    fAlephParameters[4] = 3.91720e+00; 
+  }
+
+
+  //**********************************************
+  // MC loop
+  //**********************************************
+
+  Double_t lmcPrimVtxR      = 0;
+
+  Int_t lNbMCPrimary        = 0;
+  Int_t lNbMCPart           = 0;
+
+  Int_t lPdgcodeCurrentPart = 0;
+  Double_t lRapCurrentPart  = 0;
+  Double_t lPtCurrentPart   = 0;
+  
+  Int_t lComeFromSigma      = 0;
+
+  
+  // Production Radius
+  Double_t mcPosX     = 0.0,  mcPosY      = 0.0,  mcPosZ      = 0.0;
+  Double_t mcPosR     = 0.0;
+
+  // Decay Radius
+  Double_t mcDecayPosX = 0, mcDecayPosY = 0, mcDecayPosR = 0;
+
+  // current mc particle 's mother
+  Int_t iCurrentMother  = 0, lPdgCurrentMother    = 0;
+  Bool_t lCurrentMotherIsPrimary;
+
+  // current mc particles 's daughter:
+  Int_t lPdgCurrentDaughter0 = 0, lPdgCurrentDaughter1 = 0; 
+
+  // variables for multiple reconstruction studies:
+  Int_t id0           = 0, id1          = 0;
+  //Int_t lLabelTrackN  = 0, lLabelTrackP = 0;
+  //Int_t lPartNMother  = 0, lPartPMother = 0;
+  //Int_t lPartPMotherPDGcode      = 0;
+  Int_t lNtimesReconstructedK0s   = 0, lNtimesReconstructedLambda   = 0, lNtimesReconstructedAntiLambda   = 0;
+ // Int_t lNtimesReconstructedK0sMI = 0, lNtimesReconstructedLambdaMI = 0, lNtimesReconstructedAntiLambdaMI = 0;
+
+  //****************************
+  // Start loop over MC particles
+  if (fAnalysisMC) {
+
+    // Primary vertex
+    fHistMCPrimaryVertexX->Fill(mcPrimaryVtx.At(0));
+    fHistMCPrimaryVertexY->Fill(mcPrimaryVtx.At(1));
+    fHistMCPrimaryVertexZ->Fill(mcPrimaryVtx.At(2));
+    
+    lmcPrimVtxR = TMath::Sqrt(mcPrimaryVtx.At(0)*mcPrimaryVtx.At(0)+mcPrimaryVtx.At(1)*mcPrimaryVtx.At(1));
+  
+
+    if(fAnalysisType == "ESD") {
+      
+      lNbMCPrimary = stack->GetNprimary();
+      lNbMCPart    = stack->GetNtrack();
+      
+      fHistMCMultiplicityPrimary->Fill(lNbMCPrimary);
+      fHistMCMultiplicityTracks->Fill(lNbMCPart);
+      
+      
+      for (Int_t iMc = 0; iMc < (stack->GetNtrack()); iMc++) {  
+       TParticle *p0 = stack->Particle(iMc);
+       if (!p0) {
+         //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
+         continue;
+       }
+       lPdgcodeCurrentPart = p0->GetPdgCode();
+       
+       // Keep only K0s, Lambda and AntiLambda, Xi and Phi:
+       if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) && (lPdgcodeCurrentPart != 3312 ) && (lPdgcodeCurrentPart != -3312) && (lPdgcodeCurrentPart != -333) ) continue;
+       
+       lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
+       //lEtaCurrentPart   = p0->Eta();
+       lPtCurrentPart    = p0->Pt();
+
+       iCurrentMother    = p0->GetFirstMother();
+
+//     lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();
+       if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}      
+
+       mcPosX = p0->Vx();
+       mcPosY = p0->Vy();
+       mcPosZ = p0->Vz();
+       mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
+       
+       id0  = p0->GetDaughter(0);
+       id1  = p0->GetDaughter(1);
+
+       // Decay Radius and Production Radius
+       if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
+         TParticle *pDaughter0 = stack->Particle(id0);
+         TParticle *pDaughter1 = stack->Particle(id1);
+         lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
+         lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
+         
+         mcDecayPosX = pDaughter0->Vx();
+         mcDecayPosY = pDaughter0->Vy();
+         mcDecayPosR = TMath::Sqrt(mcDecayPosX*mcDecayPosX+mcDecayPosY*mcDecayPosY);
+       }
+       else  {
+         //FIXME: shouldn't this be a fatal?
+         //Printf("ERROR: particle with label %d and/or %d not found in stack (mc loop)", id0,id1);
+         mcDecayPosR = -1.0;
+       }
+       
+       if (lPdgcodeCurrentPart==310)   {
+         fHistMCtracksProdRadiusK0s->Fill(mcPosX,mcPosY);
+         fHistMCtracksDecayRadiusK0s->Fill(mcDecayPosR);
+         if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllK0s->Fill(lPtCurrentPart);
+       }
+       else if (lPdgcodeCurrentPart==3122)  {
+         fHistMCtracksProdRadiusLambda->Fill(mcPosX,mcPosY);
+         fHistMCtracksDecayRadiusLambda->Fill(mcDecayPosR);
+         if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllLambda->Fill(lPtCurrentPart);
+       }
+       else if (lPdgcodeCurrentPart==-3122) {
+         fHistMCtracksProdRadiusAntiLambda->Fill(mcPosX,mcPosY);
+         fHistMCtracksDecayRadiusAntiLambda->Fill(mcDecayPosR);
+         if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiLambda->Fill(lPtCurrentPart);
+       }
+       
+         // FIXME: not sure if I understand this: is it correct?
+       if ( ( ( TMath::Abs(lPdgCurrentMother) == 3212)  ||
+              ( TMath::Abs(lPdgCurrentMother) == 3224)  ||
+              ( TMath::Abs(lPdgCurrentMother) == 3214)  ||
+              ( TMath::Abs(lPdgCurrentMother) == 3114) )
+            && ( iCurrentMother <= lNbMCPrimary )
+            ) lComeFromSigma = 1;
+       else lComeFromSigma = 0;
+       
+       //*********************************************
+       // Now keep only primary particles   
+       if ( ( iMc > lNbMCPrimary ) && (!lComeFromSigma) ) continue;
+
+       //********************************************
+       //check if V0 is reconstructed several times  
+     
+       lNtimesReconstructedK0s   = 0; lNtimesReconstructedLambda   = 0; lNtimesReconstructedAntiLambda   = 0;
+       //lNtimesReconstructedK0sMI = 0; lNtimesReconstructedLambdaMI = 0; lNtimesReconstructedAntiLambdaMI = 0;
+
+       //for (Int_t jV0 = 0; jV0 < fESD->GetNumberOfV0s(); jV0++) {
+       
+       //lLabelTrackN  = 0; lLabelTrackP = 0;
+       //lPartNMother  = 0; lPartPMother = 0;
+       
+       //AliESDv0    *vertexESD = ((AliESDEvent*)fESD)->GetV0(jV0);
+       //if (!vertexESD) continue;
+       
+       //AliESDtrack *trackNESD = ((AliESDEvent*)fESD)->GetTrack(TMath::Abs(vertexESD->GetNindex()));
+       //lLabelTrackN = (UInt_t)TMath::Abs(trackNESD->GetLabel());
+       //if (lLabelTrackN!=id0 && lLabelTrackN!=id1) continue;
+       
+       //AliESDtrack *trackPESD = ((AliESDEvent*)fESD)->GetTrack(TMath::Abs(vertexESD->GetPindex()));
+       //lLabelTrackP = (UInt_t)TMath::Abs(trackPESD->GetLabel());
+       //if (lLabelTrackP!=id0 && lLabelTrackP!=id1) continue;
+       
+       //TParticle   *lPartNESD = stack->Particle(lLabelTrackN);
+       //TParticle   *lPartPESD = stack->Particle(lLabelTrackP);
+       //lPartNMother = lPartNESD->GetFirstMother();
+       //lPartPMother = lPartPESD->GetFirstMother();
+       
+       //lPartPMotherPDGcode = stack->Particle(lPartPMother)->GetPdgCode();
+       
+       //switch (vertexESD->GetOnFlyStatus()){
+       
+       //case 0 : 
+       //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0s++;
+       //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambda++;
+       //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambda++;
+       //break;
+       
+       //case 1 :
+       //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0sMI++;
+       //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambdaMI++;
+       //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambdaMI++;
+       //break;
+       
+       //}     
+       //} // end loop over reconstructed V0s inside MC loop
+;
+        // Rap distribution
+        if (lPdgcodeCurrentPart==310) {
+         fHistMCRapK0s->Fill(lRapCurrentPart);
+         if (lPtCurrentPart < 0.2 && lPtCurrentPart < 3.0)
+           fHistMCRapInPtRangeK0s->Fill(lRapCurrentPart);
+       }
+
+       if (lPdgcodeCurrentPart==3122) {
+         fHistMCRapLambda->Fill(lRapCurrentPart);
+         if (lPtCurrentPart < 0.6 && lPtCurrentPart < 3.5)
+           fHistMCRapInPtRangeLambda->Fill(lRapCurrentPart);
+       }
+
+       if (lPdgcodeCurrentPart==-3122) {
+         fHistMCRapAntiLambda->Fill(lRapCurrentPart);
+         if (lPtCurrentPart < 0.6 && lPtCurrentPart < 3.5)
+           fHistMCRapInPtRangeAntiLambda->Fill(lRapCurrentPart);
+       }
+
+       if (lPdgcodeCurrentPart==3312 || lPdgcodeCurrentPart==-3312) {
+         fHistMCRapXi->Fill(lRapCurrentPart);
+         if (lPtCurrentPart < 0.6 && lPtCurrentPart < 3.0)
+           fHistMCRapInPtRangeXi->Fill(lRapCurrentPart);
+       }
+
+       if (lPdgcodeCurrentPart==333) {
+         fHistMCRapPhi->Fill(lRapCurrentPart);
+         if (lPtCurrentPart < 0.7 && lPtCurrentPart < 3.0)
+           fHistMCRapInPtRangePhi->Fill(lRapCurrentPart);
+       }
+       // Rapidity Cut
+       if (TMath::Abs(lRapCurrentPart) > lCutRap) continue;
+       if (lPdgcodeCurrentPart==310) {
+         fHistMCProdRadiusK0s->Fill(mcPosR);
+
+         fHistMCPtK0s->Fill(lPtCurrentPart);
+
+
+
+         fHistNTimesRecK0s->Fill(lNtimesReconstructedK0s);
+       //  fHistNTimesRecK0sMI->Fill(lNtimesReconstructedK0s);
+         fHistNTimesRecK0sVsPt->Fill(lPtCurrentPart,lNtimesReconstructedK0s);
+       //  fHistNTimesRecK0sVsPtMI->Fill(lPtCurrentPart,lNtimesReconstructedK0sMI);
+       }
+       else 
+       if (lPdgcodeCurrentPart==3122) {
+         fHistMCProdRadiusLambda->Fill(mcPosR);
+
+         fHistMCPtLambda->Fill(lPtCurrentPart);          
+
+
+         fHistNTimesRecLambda->Fill(lNtimesReconstructedLambda);
+       //  fHistNTimesRecLambdaMI->Fill(lNtimesReconstructedLambdaMI);
+         fHistNTimesRecLambdaVsPt->Fill(lPtCurrentPart,lNtimesReconstructedLambda);
+        // fHistNTimesRecLambdaVsPtMI->Fill(lPtCurrentPart,lNtimesReconstructedLambdaMI);
+         if (lComeFromSigma) fHistMCPtLambdaFromSigma->Fill(lPtCurrentPart);
+
+         //printf("found Lambda MC pT=%e\n",lPtCurrentPart);
+         //printf("found Lambda MC Plabel=%d PPDGcode=%d Nlabel=%d NPDGcode=%d\n\n",id0,lPdgCurrentDaughter0,id1,lPdgCurrentDaughter1); 
+         
+       }
+
+       
+      } // end loop ESD MC
+      
+    } // end ESD condition
+
+    else if(fAnalysisType == "AOD") {
+      lNbMCPart = mcArray->GetEntriesFast();
+      lNbMCPrimary = 0;
+      
+      fHistMCMultiplicityTracks->Fill(lNbMCPart);
+      
+      for (Int_t iMc = 0; iMc < lNbMCPart; iMc++) {  
+       
+       // Primary vertex TO DO !!
+       //
+       
+       AliAODMCParticle *mcAODPart = (AliAODMCParticle*)mcArray->At(iMc);
+       if (!mcAODPart) {
+         //Printf("Strange analysis task (mc loop): particle with label %d not found", iMc);
+         continue;
+       }
+       lPdgcodeCurrentPart = mcAODPart->GetPdgCode();
+       if (mcAODPart->IsPhysicalPrimary()) {lNbMCPrimary = lNbMCPrimary +1;}
+       
+       // Keep only K0s, Lambda and AntiLambda:
+       if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) ) continue;
+       
+       //lEtaCurrentPart   = mcAODPart->Eta();
+       lRapCurrentPart   = mcAODPart->Y();
+       lPtCurrentPart    = mcAODPart->Pt();
+       iCurrentMother    = mcAODPart->GetMother();
+       lPdgCurrentMother = ((AliAODMCParticle*)mcArray->At(iCurrentMother))->GetPdgCode();
+       lCurrentMotherIsPrimary = ((AliAODMCParticle*)mcArray->At(iCurrentMother))->IsPhysicalPrimary();
+       
+       mcPosX = mcAODPart->Xv();
+       mcPosY = mcAODPart->Yv();
+       mcPosZ = mcAODPart->Zv();
+       mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
+       
+       id0  = mcAODPart->GetDaughter(0);
+       id1  = mcAODPart->GetDaughter(1);
+       
+       // Decay Radius and Production Radius
+       if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
+         AliAODMCParticle *mcAODDaughter1 = (AliAODMCParticle*)mcArray->At(id1);
+         if (!mcAODPart) {
+           //Printf("Strange analysis task (mc loop): daughter not found");
+           continue;
+         }
+         mcDecayPosX = mcAODDaughter1->Xv();
+         mcDecayPosY = mcAODDaughter1->Yv();
+         mcDecayPosR = TMath::Sqrt(mcDecayPosX*mcDecayPosX+mcDecayPosY*mcDecayPosY);
+       }
+       else  {
+         //Printf("ERROR: particle with label %d and/or %d not found in stack (mc loop)", id0,id1);
+         mcDecayPosR = -1.0;
+       }
+       
+       if (lPdgcodeCurrentPart==310)   {
+         fHistMCtracksProdRadiusK0s->Fill(mcPosX,mcPosY);
+         fHistMCtracksDecayRadiusK0s->Fill(mcDecayPosR);
+         if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllK0s->Fill(lPtCurrentPart);
+       }
+       else if (lPdgcodeCurrentPart==3122)  {
+         fHistMCtracksProdRadiusLambda->Fill(mcPosX,mcPosY);
+         fHistMCtracksDecayRadiusLambda->Fill(mcDecayPosR);
+         if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllLambda->Fill(lPtCurrentPart);
+       }
+       else if (lPdgcodeCurrentPart==-3122) {
+         fHistMCtracksProdRadiusAntiLambda->Fill(mcPosX,mcPosY);
+         fHistMCtracksDecayRadiusAntiLambda->Fill(mcDecayPosR);
+         if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiLambda->Fill(lPtCurrentPart);
+       }
+       
+       if ( ( ( TMath::Abs(lPdgCurrentMother) == 3212)  ||
+              ( TMath::Abs(lPdgCurrentMother) == 3224)  ||
+              ( TMath::Abs(lPdgCurrentMother) == 3214)  ||
+              ( TMath::Abs(lPdgCurrentMother) == 3114) )
+            && (lCurrentMotherIsPrimary)
+            ) lComeFromSigma = 1;
+       else lComeFromSigma = 0;
+      
+       //*********************************************
+        // Now keep only primary particles 
+        
+       // FIX IT !!!!    iMC is not defined !!!! FIX IT also in ESD/AOD loop !!
+       if ( ( iMc > lNbMCPrimary ) && (!lComeFromSigma) ) continue;
+
+       //********************************************
+       // check if V0 is reconstructed several times  
+       
+       //lNtimesReconstructedK0s   = 0; lNtimesReconstructedLambda   = 0; lNtimesReconstructedAntiLambda   = 0;
+        //lNtimesReconstructedK0sMI = 0; lNtimesReconstructedLambdaMI = 0; lNtimesReconstructedAntiLambdaMI = 0;
+           
+        //for (Int_t jV0 = 0; jV0 < fESD->GetNumberOfV0s(); jV0++) {
+       
+       //lLabelTrackN  = 0; lLabelTrackP = 0;
+       //lPartNMother  = 0; lPartPMother = 0;
+
+       //AliAODv0    *vertexAOD= ((AliAODEvent*)fESD)->GetV0(jV0);
+       //if (!vertexAOD) continue;
+       //printf("enter!!");
+       //AliVParticle  *trackP  = ((AliVEvent*)fESD)->GetTrack(vertexAOD->GetPosID());
+       //if (!trackP) continue;
+       //lLabelTrackP = TMath::Abs(trackP->GetLabel());
+       //if (lLabelTrackP!=id0 && lLabelTrackP!=id1) continue;
+       
+       //AliVParticle  *trackN  = ((AliVEvent*)fESD)->GetTrack(vertexAOD->GetNegID());
+       //if (!trackN) continue;
+       //lLabelTrackN = TMath::Abs(trackN->GetLabel());
+       //if (lLabelTrackN!=id0 && lLabelTrackN!=id1) continue;
+       
+       //AliAODMCParticle *lPartNAOD = (AliAODMCParticle*)mcArray->At(lLabelTrackN);
+       //if (!lPartNAOD) continue;
+       //AliAODMCParticle *lPartPAOD = (AliAODMCParticle*)mcArray->At(lLabelTrackP);
+       //if (!lPartPAOD) continue;
+       
+       //lPartNMother = lPartNAOD->GetMother();
+       //lPartPMother = lPartPAOD->GetMother();
+
+       //lPartPMotherPDGcode = ((AliAODMCParticle*)mcArray->At(lPartPMother))->GetPdgCode();
+       
+       //switch (vertexAOD->GetOnFlyStatus()){
+         
+       //case 0 : 
+         //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0s++;
+         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambda++;
+         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambda++;
+         //break;
+         
+       //case 1 :
+         //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0sMI++;
+         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambdaMI++;
+         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambdaMI++;
+         //break;
+         
+       ///}    
+      //} // end loop over reconstructed V0s inside MC loop
+    
+       if (TMath::Abs(lRapCurrentPart) > lCutRap) continue;
+       if (lPdgcodeCurrentPart==310) {
+         fHistMCProdRadiusK0s->Fill(mcPosR);
+         fHistMCPtK0s->Fill(lPtCurrentPart);
+         fHistNTimesRecK0s->Fill(lNtimesReconstructedK0s);
+       //  fHistNTimesRecK0sMI->Fill(lNtimesReconstructedK0s);
+         fHistNTimesRecK0sVsPt->Fill(lPtCurrentPart,lNtimesReconstructedK0s);
+       //  fHistNTimesRecK0sVsPtMI->Fill(lPtCurrentPart,lNtimesReconstructedK0sMI);
+       }
+       else if (lPdgcodeCurrentPart==3122) {
+         fHistMCProdRadiusLambda->Fill(mcPosR);
+         fHistMCPtLambda->Fill(lPtCurrentPart);
+         fHistNTimesRecLambda->Fill(lNtimesReconstructedLambda);
+       //  fHistNTimesRecLambdaMI->Fill(lNtimesReconstructedLambdaMI);
+         fHistNTimesRecLambdaVsPt->Fill(lPtCurrentPart,lNtimesReconstructedLambda);
+       //  fHistNTimesRecLambdaVsPtMI->Fill(lPtCurrentPart,lNtimesReconstructedLambdaMI);
+         if (lComeFromSigma) fHistMCPtLambdaFromSigma->Fill(lPtCurrentPart);
+       }
+       else if (lPdgcodeCurrentPart==-3122) {
+         fHistMCProdRadiusAntiLambda->Fill(mcPosR);
+         //fHistMCPtAntiLambda->Fill(lPtCurrentPart);
+         fHistNTimesRecAntiLambda->Fill(lNtimesReconstructedAntiLambda);
+       //  fHistNTimesRecAntiLambdaMI->Fill(lNtimesReconstructedAntiLambdaMI);
+         fHistNTimesRecAntiLambdaVsPt->Fill(lPtCurrentPart,lNtimesReconstructedAntiLambda);
+       //  fHistNTimesRecAntiLambdaVsPtMI->Fill(lPtCurrentPart,lNtimesReconstructedAntiLambdaMI);
+         if (lComeFromSigma) fHistMCPtAntiLambdaFromSigma->Fill(lPtCurrentPart);
+       }
+
+      } // end loop over AODMC particles 
+      fHistMCMultiplicityPrimary->Fill(lNbMCPrimary);
+      
+    } // end AOD condition
+
+  } // End Loop over MC condition
+
+  
+
+
+
+  //************************************
+  // ESD or AOD loop 
+  //************************************
+
+  Double_t lMagneticField = 999;
+
+  //Multiplcity:
+  Int_t    nv0sTot= 0, nv0s = 0;
+//  Int_t nv0sMI =0;   
+  // Variables:
+  Double_t  lV0Position[3];
+  Double_t lDcaPosToPrimVertex = 0;
+  Double_t lDcaNegToPrimVertex = 0;
+  Double_t lDcaV0Daughters     = 0;
+  Double_t lV0cosPointAngle    = 0;
+  Double_t lChi2V0             = 0;
+  Double_t lV0DecayLength      = 0;
+  Double_t lV0Radius           = 0;
+  Double_t lDcaV0ToPrimVertex  = 0;
+  
+  Int_t    lOnFlyStatus        = 0;
+  //Float_t   tdcaPosToPrimVertexXYZ[2], tdcaNegToPrimVertexXYZ[2]; // ..[0] = Impact parameter in XY plane and ..[1] = Impact parameter in Z            
+  //Double_t  tdcaDaughterToPrimVertex[2];                          // ..[0] = Pos and ..[1] = Neg
+
+  
+
+  Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
+  Double_t lPtK0s      = 0, lPtLambda      = 0, lPtAntiLambda      = 0;
+  Double_t lRapK0s     = 0, lRapLambda     = 0, lRapAntiLambda     = 0;
+  Double_t lEtaK0s     = 0, lEtaLambda     = 0, lEtaAntiLambda     = 0;
+  Double_t lAlphaV0      = 0, lPtArmV0       = 0;
+
+  Double_t lPzK0s      = 0, lPzLambda      = 0, lPzAntiLambda      = 0;
+
+
+  Double_t lV0Eta = 999;
+  
+  // to study Associated V0s:
+  Int_t    lIndexTrackPos       = 0, lIndexTrackNeg         = 0;
+  UInt_t   lLabelTrackPos       = 0, lLabelTrackNeg         = 0;
+  Int_t    lCheckPIdK0Short     = 0, lCheckMcK0Short        = 0;
+  Int_t    lCheckPIdLambda      = 0, lCheckMcLambda         = 0;
+  Int_t    lCheckPIdAntiLambda  = 0, lCheckMcAntiLambda     = 0;
+  Int_t    lCheckSecondaryK0s   = 0, lCheckSecondaryLambda  = 0, lCheckSecondaryAntiLambda  = 0;
+  Int_t    lCheckGamma          = 0;
+  Double_t mcPosMotherX         = 0, mcPosMotherY           = 0, mcPosMotherZ  = 0;
+  Double_t mcPosMotherR         = 0;
+  Double_t mcMotherPt           = 0;
+
+  Int_t lIndexPosMother        = 0;
+  Int_t lIndexNegMother        = 0;
+  Int_t lIndexMotherOfMother   = 0;
+  Int_t lPDGCodePosDaughter    = 0;
+  Int_t lPDGCodeNegDaughter    = 0;
+  Int_t lPdgcodeMother         = 0;
+  Int_t lPdgcodeMotherOfMother = 0;
+
+  // Reconstructed position
+  Double_t rcPosXK0s        = 0,  rcPosYK0s        = 0, rcPosZK0s        = 0;
+  Double_t rcPosRK0s        = 0;
+  Double_t rcPosXLambda     = 0,  rcPosYLambda     = 0, rcPosZLambda     = 0;
+  Double_t rcPosRLambda     = 0;
+  Double_t rcPosXAntiLambda = 0,  rcPosYAntiLambda = 0, rcPosZAntiLambda = 0;
+  Double_t rcPosRAntiLambda = 0;
+
+  // Pt resolution
+  Double_t deltaPtK0s  = 0, deltaPtLambda  = 0, deltaPtAntiLambda  = 0;
+
+  // Daughters
+  AliESDtrack  *myTrackPos  = NULL;
+  AliESDtrack  *myTrackNeg  = NULL;
+  AliVParticle *lVPartPos   = NULL;
+  AliVParticle *lVPartNeg   = NULL;
+
+  // Daughters' momentum:
+  Double_t  lMomPos[3] = {999,999,999};
+  Double_t  lMomNeg[3] = {999,999,999};
+  Double_t  lPtPos = 999, lPtNeg = 999;
+  Double_t  lPPos = 999, lPNeg = 999;
+
+  // Inner Wall parameters:
+  Double_t  lMomInnerWallPos =999, lMomInnerWallNeg = 999;
+
+  // AliKF Chi2 and Armenteros variables
+  Double_t lChi2KFK0s  = 0, lChi2KFLambda = 0,  lChi2KFAntiLambda = 0;
+  Double_t lAlphaV0K0s = 0, lAlphaV0Lambda = 0,  lAlphaV0AntiLambda = 0;
+  Double_t lPtArmV0K0s = 0, lPtArmV0Lambda = 0,  lPtArmV0AntiLambda = 0;
+  Double_t lQlPos   = 0, lQlNeg   = 0;
+
+
+  // PID
+  Float_t nSigmaPosPion   = 0;
+  Float_t nSigmaNegPion   = 0;
+
+  Float_t nSigmaPosProton = 0;
+  Float_t nSigmaNegProton = 0;
+
+  Int_t lCheckPIDK0sPosDaughter        = 0, lCheckPIDK0sNegDaughter        = 0;
+  Int_t lCheckPIDLambdaPosDaughter     = 0, lCheckPIDLambdaNegDaughter     = 0;
+  Int_t lCheckPIDAntiLambdaPosDaughter = 0, lCheckPIDAntiLambdaNegDaughter = 0;
+
+  
+  
+  //***********************
+  // Primary Vertex cuts &
+  // Magnetic field and Quality tracks cuts 
+
+  Double_t  lPrimaryVtxPosition[3];
+  Double_t  lPrimaryVtxCov[6];
+  Double_t  lPrimaryVtxChi2 = 999;
+  Double_t  lResPrimaryVtxX = 999;
+  Double_t  lResPrimaryVtxY = 999;
+  Double_t  lResPrimaryVtxZ = 999;
+     
+  AliAODVertex *myPrimaryVertex = NULL;
+  //const AliVVertex *mySPDPrimaryVertex = NULL;
+
+  AliESDtrackCuts *myTracksCuts = NULL;
+     
+  const AliMultiplicity *myMultiplicty = ((AliESDEvent*)fESD)->GetMultiplicity();
+
+  if(fAnalysisType == "ESD") {  
+  
+    // Best Primary Vertex:  
+    const AliESDVertex *myBestPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertex();
+    myBestPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertex();
+    if (!myBestPrimaryVertex) return;
+    if (!myBestPrimaryVertex->GetStatus()) return;
+    fHistNumberEvents->Fill(3.5);
+    myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
+    myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
+    if ( ( TMath::Abs(lPrimaryVtxPosition[2]) ) > cutPrimVertex) return ;
+    fHistNumberEvents->Fill(4.5);    
+    lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
+    lResPrimaryVtxX = myBestPrimaryVertex->GetXRes();
+    lResPrimaryVtxY = myBestPrimaryVertex->GetYRes();
+    lResPrimaryVtxZ = myBestPrimaryVertex->GetZRes();
+
+    // remove TPC-only primary vertex : retain only events with tracking + SPD vertex
+    const AliESDVertex *mySPDPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertexSPD();
+    if (!mySPDPrimaryVertex) return;
+    fHistSPDPrimaryVertexZ->Fill(mySPDPrimaryVertex->GetZ());
+    const AliESDVertex *myPrimaryVertexTracking = ((AliESDEvent*)fESD)->GetPrimaryVertexTracks();
+    if (!myPrimaryVertexTracking) return;
+    if (!mySPDPrimaryVertex->GetStatus() && !myPrimaryVertexTracking->GetStatus() ) return;
+    fHistNumberEvents->Fill(5.5);
+
+    
+    myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
+    if (!myPrimaryVertex) return;
+
+
+     // Number of Tracklets:
+    //const AliMultiplicity *myMultiplicty = ((AliESDEvent*)fESD)->GetMultiplicity();
+    //if (myMultiplicty->GetNumberOfTracklets() < 10) return;
+    fHistTrackletPerEvent->Fill(myMultiplicty->GetNumberOfTracklets());
+
+    lMagneticField = ((AliESDEvent*)fESD)->GetMagneticField();
+
+    myTracksCuts = new AliESDtrackCuts();
+    // require TPC refit
+    myTracksCuts->SetRequireTPCRefit(kTRUE);
+    // minimum number of clusters in TPC
+    myTracksCuts->SetMinNClustersTPC(nbMinTPCclusters);
+
+  }
+  
+  else if(fAnalysisType == "AOD") {
+    printf("enter AOD!!");
+    myPrimaryVertex = ((AliAODEvent*)fESD)->GetPrimaryVertex();
+    if (!myPrimaryVertex) return;
+
+    lPrimaryVtxPosition[0] = myPrimaryVertex->GetX();
+    lPrimaryVtxPosition[1] = myPrimaryVertex->GetY();
+    lPrimaryVtxPosition[2] = myPrimaryVertex->GetZ();
+
+    // Cut on SPD vertex and fill histo Nevents: FIX it !
+
+     // Tracks cuts FIX IT !
+
+    // FIX it !!!
+    lMagneticField = 999;   
+
+  }
+
+  fHistPrimaryVertexX->Fill(lPrimaryVtxPosition[0]);
+  fHistPrimaryVertexY->Fill(lPrimaryVtxPosition[1]);
+  fHistPrimaryVertexZ->Fill(lPrimaryVtxPosition[2]);
+  //Double_t lrcPrimVtxR = TMath::Sqrt(lPrimaryVtxPosition[0]*lPrimaryVtxPosition[0]+lPrimaryVtxPosition[0]*lPrimaryVtxPosition[0]);
+
+  fHistPrimaryVertexResX->Fill(lResPrimaryVtxX);
+  fHistPrimaryVertexResY->Fill(lResPrimaryVtxY);
+  fHistPrimaryVertexResZ->Fill(lResPrimaryVtxZ);
+
+
+  //***********************
+  // AliKF Primary Vertex
+
+  AliKFVertex primaryVtxKF( *myPrimaryVertex );
+  AliKFParticle::SetField(lMagneticField);
+
+
+  //************************************
+  // PID
+
+  AliESDpid *fESDpid = new AliESDpid(); // FIXME delete
+  fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]); 
+      
+
+
+  //***Rerun the V0 finder
+
+//  fESD->ResetV0s();
+//  AliV0vertexer v0Vertexer;
+//  v0Vertexer.SetCuts(fCuts);
+//  v0Vertexer.Tracks2V0vertices(fESD);
+  
+  //*************************
+  // V0 loop
+      
+  nv0sTot = fESD->GetNumberOfV0s();
+  if (!nv0sTot) fHistNumberEvents->Fill(6.5);
+
+  for (Int_t iV0 = 0; iV0 < nv0sTot; iV0++) {
+    
+    // ALiKF
+    AliKFParticle* negPiKF = NULL;
+    AliKFParticle* posPiKF = NULL;
+    AliKFParticle* posPKF  = NULL;
+    AliKFParticle* negAPKF = NULL;
+
+    
+    lIndexPosMother     = 0; lIndexNegMother     = 0; lIndexMotherOfMother       = 0;
+    lCheckPIdK0Short    = 0; lCheckMcK0Short     = 0; lCheckSecondaryK0s         = 0;
+    lCheckPIdLambda     = 0; lCheckMcLambda      = 0; lCheckSecondaryLambda      = 0;
+    lCheckPIdAntiLambda = 0; lCheckMcAntiLambda  = 0; lCheckSecondaryAntiLambda  = 0;       
+    lComeFromSigma      = -1;
+    
+    
+    if(fAnalysisType == "ESD") {
+
+
+      AliESDv0 *v0 = ((AliESDEvent*)fESD)->GetV0(iV0);
+      if (!v0) continue;
+      
+      // Primary vertex:
+      fHistPrimaryVertexPosXV0events->Fill(lPrimaryVtxPosition[0]);
+      fHistPrimaryVertexPosYV0events->Fill(lPrimaryVtxPosition[1]);
+      fHistPrimaryVertexPosZV0events->Fill(lPrimaryVtxPosition[2]);
+      
+      // V0's Daughters
+      lIndexTrackPos = TMath::Abs(v0->GetPindex());
+      lIndexTrackNeg = TMath::Abs(v0->GetNindex());
+      AliESDtrack *myTrackPosTest = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
+      AliESDtrack *myTrackNegTest = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
+      if (!myTrackPosTest || !myTrackNegTest) {
+       Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
+       continue;
+      }
+      // Remove like-sign
+      if ( myTrackPosTest->GetSign() == myTrackNegTest->GetSign()){
+       continue;
+      } 
+     
+      // VO's main characteristics to check the reconstruction cuts
+      lOnFlyStatus       = v0->GetOnFlyStatus();
+      lChi2V0            = v0->GetChi2V0();
+      lDcaV0Daughters    = v0->GetDcaV0Daughters();
+      lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
+      lV0cosPointAngle   = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1], lPrimaryVtxPosition[2]);
+
+      v0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
+
+      lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
+      lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
+                                  TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
+                                  TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
+
+
+      if( myTrackPosTest->GetSign() ==1){
+       
+       myTrackPos = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
+       myTrackNeg = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
+
+       // Daughters' momentum;
+       v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
+       v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
+
+       if (negPiKF) delete negPiKF; negPiKF=NULL;
+       if (posPiKF) delete posPiKF; posPiKF=NULL;
+       if (posPKF)  delete posPKF;  posPKF=NULL;
+       if (negAPKF) delete negAPKF; negAPKF=NULL;
+       
+       negPiKF = new AliKFParticle( *(v0->GetParamN()) ,-211);
+       posPiKF = new AliKFParticle( *(v0->GetParamP()) ,211);
+       posPKF  = new AliKFParticle( *(v0->GetParamP()) ,2212);
+       negAPKF = new AliKFParticle( *(v0->GetParamN()) ,-2212);
+       
+      }
+           
+      if( myTrackPosTest->GetSign() ==-1){
+       
+       myTrackPos = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
+       myTrackNeg = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
+
+       // Daughters' momentum;
+       v0->GetPPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
+       v0->GetNPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
+
+       if (negPiKF) delete negPiKF; negPiKF=NULL;
+       if (posPiKF) delete posPiKF; posPiKF=NULL;
+       if (posPKF)  delete posPKF;  posPKF=NULL;
+       if (negAPKF) delete negAPKF; negAPKF=NULL;
+       
+       negPiKF = new AliKFParticle( *(v0->GetParamP()) ,-211);
+       posPiKF = new AliKFParticle( *(v0->GetParamN()) ,211);
+       posPKF  = new AliKFParticle( *(v0->GetParamN()) ,2212);
+       negAPKF = new AliKFParticle( *(v0->GetParamP()) ,-2212);
+
+      }
+
+      lLabelTrackPos = (UInt_t)TMath::Abs(myTrackPos->GetLabel());
+      lLabelTrackNeg = (UInt_t)TMath::Abs(myTrackNeg->GetLabel());
+
+      // Daughters Pt and P:
+      lPtPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1]);
+      lPtNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1]);
+
+      lPPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1] + lMomPos[2]*lMomPos[2]);
+      lPNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1] + lMomNeg[2]*lMomNeg[2]);
+
+      // Inner Wall parameter:
+      const AliExternalTrackParam *myInnerWallTrackPos = myTrackPos->GetInnerParam(); 
+      if(myInnerWallTrackPos) lMomInnerWallPos = myInnerWallTrackPos->GetP(); 
+      const AliExternalTrackParam *myInnerWallTrackNeg = myTrackNeg->GetInnerParam(); 
+      if(myInnerWallTrackNeg) lMomInnerWallNeg = myInnerWallTrackNeg->GetP(); 
+             
+      // DCA between daughter and Primary Vertex:
+      if (myTrackPos) lDcaPosToPrimVertex = TMath::Abs(myTrackPos->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lMagneticField) );
+      
+      if (myTrackNeg) lDcaNegToPrimVertex = TMath::Abs(myTrackNeg->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lMagneticField) );
+      
+      // Quality tracks cuts:
+      if ( !(myTracksCuts->IsSelected(myTrackPos)) || !(myTracksCuts->IsSelected(myTrackNeg)) ) continue;
+
+      // Armenteros variables:
+      lAlphaV0      =  v0->AlphaV0();
+      lPtArmV0      =  v0->PtArmV0();
+
+      // Pseudorapidity:
+      lV0Eta = v0->Eta();
+      
+      // PID
+      if (fUsePID.Contains("withPID")) {
+       nSigmaPosPion   = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackPos,AliPID::kPion));
+       
+       nSigmaNegPion   = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackNeg,AliPID::kPion));
+
+       nSigmaPosProton = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackPos,AliPID::kProton));
+       
+       nSigmaNegProton = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackNeg,AliPID::kProton));
+      }
+      else {
+       nSigmaPosPion = 0; nSigmaNegPion =0; nSigmaPosProton = 0; nSigmaNegProton= 0;
+      }
+      
+      
+      
+      // Monte-Carlo particle associated to reconstructed particles: 
+      if (fAnalysisMC) {
+       //if (lLabelTrackPos < 0 || lLabelTrackNeg < 0) continue;
+       TParticle  *lMCESDPartPos  = stack->Particle(lLabelTrackPos);
+       if(!lMCESDPartPos) { 
+         Printf("no MC particle for positive and/or negative daughter\n");
+         continue;
+       }
+       TParticle  *lMCESDPartNeg  = stack->Particle(lLabelTrackNeg);
+       if (!lMCESDPartNeg) continue;
+       lPDGCodePosDaughter = lMCESDPartPos->GetPdgCode();
+       lPDGCodeNegDaughter = lMCESDPartNeg->GetPdgCode();
+       lIndexPosMother = lMCESDPartPos->GetFirstMother();
+       lIndexNegMother = lMCESDPartNeg->GetFirstMother();
+       
+       if (lIndexPosMother == -1) continue;
+       TParticle  *lMCESDMother    = stack->Particle(lIndexPosMother);
+       if (!lMCESDMother) continue;
+       lPdgcodeMother         = lMCESDMother->GetPdgCode();
+       lIndexMotherOfMother   = lMCESDMother->GetFirstMother();
+       if (lIndexMotherOfMother ==-1) lPdgcodeMotherOfMother = 0;
+       else {
+         TParticle  *lMCESDMotherOfMother    = stack->Particle(lIndexMotherOfMother);
+         if (!lMCESDMotherOfMother) continue;
+         lPdgcodeMotherOfMother = lMCESDMotherOfMother->GetPdgCode();
+       }
+       
+       mcPosX = lMCESDPartPos->Vx();
+       mcPosY = lMCESDPartPos->Vy();
+       mcPosZ = lMCESDPartPos->Vz();
+       mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
+       mcPosMotherX = lMCESDMother->Vx();
+       mcPosMotherY = lMCESDMother->Vy();
+       mcPosMotherZ = lMCESDMother->Vz();
+       mcPosMotherR = TMath::Sqrt(mcPosMotherX*mcPosMotherX+mcPosMotherY*mcPosMotherY);
+       
+       mcMotherPt   = lMCESDMother->Pt();
+      }
+
+    } // end ESD condition
+
+
+    
+    else if(fAnalysisType == "AOD") { 
+
+      AliAODv0     *myAODv0 = ((AliAODEvent*)fESD)->GetV0(iV0);
+      if (!myAODv0) continue;
+
+      // Primary vertex:
+      fHistPrimaryVertexPosXV0events->Fill(lPrimaryVtxPosition[0]);
+      fHistPrimaryVertexPosYV0events->Fill(lPrimaryVtxPosition[1]);
+      fHistPrimaryVertexPosZV0events->Fill(lPrimaryVtxPosition[2]);
+      
+
+      //Multiplicity:
+      if(!lOnFlyStatus) nv0s++;
+//      else  if(lOnFlyStatus) nv0sMI++;
+
+      // V0's Daughters
+      lIndexTrackPos = TMath::Abs(myAODv0->GetPosID());
+      lIndexTrackNeg = TMath::Abs(myAODv0->GetNegID());
+      
+      AliVParticle  *lVPartPosTest  = ((AliVEvent*)fESD)->GetTrack(lIndexTrackPos);
+      AliVParticle  *lVPartNegTest  = ((AliVEvent*)fESD)->GetTrack(lIndexTrackNeg);
+      //AliAODTrack  *lVPartPos  = ((AliAODEvent*)fESD)->GetTrack(lIndexTrackPos);
+      //AliAODTrack  *lVPartNeg  = ((AliAODEvent*)fESD)->GetTrack(lIndexTrackNeg);
+
+      if (!lVPartPosTest ||(!lVPartNegTest )) {
+       Printf("strange analysis::UserExec:: Could not retreive one of the daughter track\n");
+       continue;
+      }
+
+      // Quality cuts:
+      // TO DO !!!!!!!
+
+      // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
+      //if( !(lVPartPosTest->GetStatus() & AliAODTrack::kTPCrefit)) continue;      
+      //if( !(lVPartNegTest->GetStatus() & AliAODTrack::kTPCrefit)) continue;
+      
+
+      lDcaPosToPrimVertex = myAODv0->DcaPosToPrimVertex();     
+      lDcaNegToPrimVertex = myAODv0->DcaNegToPrimVertex();
+      lOnFlyStatus        = myAODv0->GetOnFlyStatus();
+      lChi2V0             = myAODv0->Chi2V0();
+      lDcaV0Daughters     = myAODv0->DcaV0Daughters();
+      lDcaV0ToPrimVertex  = myAODv0->DcaV0ToPrimVertex();
+      lV0DecayLength      = myAODv0->DecayLengthV0(lPrimaryVtxPosition);
+      lV0cosPointAngle    = myAODv0->CosPointingAngle(lPrimaryVtxPosition);
+      lV0Radius           = myAODv0->RadiusV0();
+
+      if( lVPartPosTest->Charge() ==1){
+       
+       lVPartPos = ((AliVEvent*)fESD)->GetTrack(lIndexTrackPos);
+       lVPartNeg = ((AliVEvent*)fESD)->GetTrack(lIndexTrackNeg);
+       
+       
+       if (negPiKF) delete negPiKF; negPiKF=NULL;
+       if (posPiKF) delete posPiKF; posPiKF=NULL;
+       if (posPKF)  delete posPKF; posPKF=NULL;
+       if (negAPKF) delete negAPKF; negAPKF=NULL;
+       
+       //negPiKF = new AliKFParticle( *(myAODv0->GetParamN()) ,-211);
+       //posPiKF = new AliKFParticle( *(myAODv0->GetParamP()) ,211);
+       //posPKF  = new AliKFParticle( *(myAODv0->GetParamP()) ,2212);
+       //negAPKF = new AliKFParticle( *(myAODv0->GetParamN()) ,-2212);
+       // TO DO !!!!!!
+       negPiKF = NULL;
+       posPiKF = NULL;
+       posPKF  = NULL;
+       negAPKF = NULL;
+       
+      }
+           
+      if( lVPartPosTest->Charge() ==-1){
+       
+       lVPartPos = ((AliVEvent*)fESD)->GetTrack(lIndexTrackNeg);
+       lVPartNeg = ((AliVEvent*)fESD)->GetTrack(lIndexTrackPos);
+       
+       if (negPiKF) delete negPiKF; negPiKF=NULL;
+       if (posPiKF) delete posPiKF; posPiKF=NULL;
+       if (posPKF)  delete posPKF; posPKF=NULL;
+       if (negAPKF) delete negAPKF; negAPKF=NULL;
+       
+       //negPiKF = new AliKFParticle( *(myAODv0->GetParamP()) ,-211);
+       //posPiKF = new AliKFParticle( *(myAODv0->GetParamN()) ,211);
+       //posPKF  = new AliKFParticle( *(myAODv0->GetParamN()) ,2212);
+       //negAPKF = new AliKFParticle( *(myAODv0->GetParamP()) ,-2212);
+       negPiKF = NULL;
+       posPiKF = NULL;
+       posPKF  = NULL;
+       negAPKF = NULL;
+      }
+
+      lLabelTrackPos  = TMath::Abs(lVPartPos->GetLabel());
+      lLabelTrackNeg  = TMath::Abs(lVPartNeg->GetLabel());
+      
+      // Armenteros variables:
+      lAlphaV0   = myAODv0->AlphaV0();
+      lPtArmV0   = myAODv0->PtArmV0();
+
+      // Pseudorapidity:
+      lV0Eta = myAODv0->PseudoRapV0();
+      
+      // PID not accessible with AOD !
+      nSigmaPosPion = 0; nSigmaNegPion =0; nSigmaPosProton = 0; nSigmaNegProton= 0;
+
+      
+      // Monte-Carlo particle associated to reconstructed particles:  
+      if (fAnalysisMC) {
+       AliAODMCParticle *lMCAODPartPos = (AliAODMCParticle*)mcArray->At(lLabelTrackPos);
+       if (!lMCAODPartPos) continue;
+       AliAODMCParticle *lMCAODPartNeg = (AliAODMCParticle*)mcArray->At(lLabelTrackNeg);
+       if(!lMCAODPartNeg) { 
+        // Printf("strange analysis::UserExec:no MC particle for negative daughter\n");
+         continue;
+       }
+       lPDGCodePosDaughter = lMCAODPartPos->GetPdgCode();
+       lPDGCodeNegDaughter = lMCAODPartNeg->GetPdgCode();
+       lIndexPosMother = lMCAODPartPos->GetMother();
+       lIndexNegMother = lMCAODPartNeg->GetMother();
+       
+       AliAODMCParticle *lMCAODMother = (AliAODMCParticle*)mcArray->At(lIndexPosMother);
+       lPdgcodeMother = lMCAODMother->GetPdgCode();
+       lIndexMotherOfMother  = lMCAODMother->GetMother();
+       if (lIndexMotherOfMother ==-1) lPdgcodeMotherOfMother = 0;
+       else {
+         lPdgcodeMotherOfMother   = ((AliAODMCParticle*)mcArray->At(lIndexMotherOfMother))->GetPdgCode();
+       }
+       
+       mcPosX = lMCAODPartPos->Xv();
+       mcPosY = lMCAODPartPos->Yv();
+       mcPosZ = lMCAODPartPos->Zv();
+       mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
+       mcPosMotherX = lMCAODMother->Xv();
+       mcPosMotherY = lMCAODMother->Yv();
+       mcPosMotherZ = lMCAODMother->Zv();
+       mcPosMotherR = TMath::Sqrt(mcPosMotherX*mcPosMotherX+mcPosMotherY*mcPosMotherY);
+       mcMotherPt   = lMCAODMother->Pt();
+      }
+            
+    } // end AOD condition
+    
+    
+    // Multiplicity:
+    if(!lOnFlyStatus) nv0s++;
+//    else  if(lOnFlyStatus) nv0sMI++;
+
+    // Daughter momentum cut: ! FIX it in case of AOD !
+    if ( (lPtPos  < cutMinPtDaughter ) ||
+         (lPtNeg  < cutMinPtDaughter )
+        ) continue;
+    
+    AliKFParticle v0K0sKF;
+    v0K0sKF+=(*negPiKF);
+    v0K0sKF+=(*posPiKF);
+    v0K0sKF.SetProductionVertex(primaryVtxKF);
+    
+    AliKFParticle v0LambdaKF;
+    v0LambdaKF+=(*negPiKF);
+    v0LambdaKF+=(*posPKF);     
+    v0LambdaKF.SetProductionVertex(primaryVtxKF);
+    
+    AliKFParticle v0AntiLambdaKF;
+    v0AntiLambdaKF+=(*posPiKF);
+    v0AntiLambdaKF+=(*negAPKF);
+    v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
+    
+    // Invariant mass
+    lInvMassK0s        = v0K0sKF.GetMass();
+    lInvMassLambda     = v0LambdaKF.GetMass();
+    lInvMassAntiLambda = v0AntiLambdaKF.GetMass();
+    
+    // Rapidity:
+    lRapK0s        = 0.5*TMath::Log((v0K0sKF.GetE()+v0K0sKF.GetPz())/(v0K0sKF.GetE()-v0K0sKF.GetPz()+1.e-13));
+    lRapLambda     = 0.5*TMath::Log((v0LambdaKF.GetE()+v0LambdaKF.GetPz())/(v0LambdaKF.GetE()-v0LambdaKF.GetPz()+1.e-13));
+    lRapAntiLambda = 0.5*TMath::Log((v0AntiLambdaKF.GetE()+v0AntiLambdaKF.GetPz())/(v0AntiLambdaKF.GetE()-v0AntiLambdaKF.GetPz()+1.e-13));
+
+    // Pseudo-rapidity
+    lEtaK0s     = v0K0sKF.GetEta();
+    lEtaLambda  = v0LambdaKF.GetEta();
+    lEtaAntiLambda  = v0AntiLambdaKF.GetEta();
+
+    // Pz:
+    lPzK0s        = v0K0sKF.GetPz();
+    lPzLambda     = v0LambdaKF.GetPz();
+    lPzAntiLambda = v0AntiLambdaKF.GetPz();
+    
+    // Pt:
+    lPtK0s        = v0K0sKF.GetPt();
+    lPtLambda     = v0LambdaKF.GetPt();
+    lPtAntiLambda = v0AntiLambdaKF.GetPt();
+
+if (lPtK0s==0) continue;
+if (lPtLambda==0) continue;
+
+    // Pt Resolution
+    deltaPtK0s        = (lPtK0s - mcMotherPt)/mcMotherPt;
+    deltaPtLambda     = (lPtLambda - mcMotherPt)/mcMotherPt;
+    deltaPtAntiLambda = (lPtAntiLambda - mcMotherPt)/mcMotherPt;
+
+    // KF Chi2
+    lChi2KFK0s        = v0K0sKF.GetChi2();
+    lChi2KFLambda     = v0LambdaKF.GetChi2();
+    lChi2KFAntiLambda = v0AntiLambdaKF.GetChi2();
+    
+    // Reconstructed Position
+    rcPosXK0s = v0K0sKF.GetX();
+    rcPosYK0s = v0K0sKF.GetY(); 
+    rcPosZK0s = v0K0sKF.GetZ();
+    rcPosRK0s = TMath::Sqrt(rcPosXK0s*rcPosXK0s+rcPosYK0s*rcPosYK0s);
+
+    rcPosXLambda = v0LambdaKF.GetX(); 
+    rcPosYLambda = v0LambdaKF.GetY(); 
+    rcPosZLambda = v0LambdaKF.GetZ();
+    rcPosRLambda = TMath::Sqrt(rcPosXLambda*rcPosXLambda+rcPosYLambda*rcPosYLambda); 
+
+    rcPosXAntiLambda = v0AntiLambdaKF.GetX();
+    rcPosYAntiLambda = v0AntiLambdaKF.GetY(); 
+    rcPosZAntiLambda = v0AntiLambdaKF.GetZ();
+    rcPosRAntiLambda = TMath::Sqrt(rcPosXAntiLambda*rcPosXAntiLambda+rcPosYAntiLambda*rcPosYAntiLambda); 
+
+    TVector3 momPos(lMomPos[0],lMomPos[1],lMomPos[2]);
+    TVector3 momNeg(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
+    TVector3 momTotK0s(v0K0sKF.GetPx(),v0K0sKF.GetPy(),v0K0sKF.GetPz());
+    TVector3 momTotLambda(v0LambdaKF.GetPx(),v0LambdaKF.GetPy(),v0LambdaKF.GetPz());
+    TVector3 momTotAntiLambda(v0AntiLambdaKF.GetPx(),v0AntiLambdaKF.GetPy(),v0AntiLambdaKF.GetPz());
+    
+    lQlPos = momPos.Dot(momTotK0s)/momTotK0s.Mag();
+    lQlNeg = momNeg.Dot(momTotK0s)/momTotK0s.Mag();
+    lAlphaV0K0s = 1.-2./(1.+lQlPos/lQlNeg);
+    lQlPos = momPos.Dot(momTotLambda)/momTotLambda.Mag();
+    lQlNeg = momNeg.Dot(momTotLambda)/momTotLambda.Mag();
+    lAlphaV0Lambda = 1.-2./(1.+lQlPos/lQlNeg);
+    lQlPos = momPos.Dot(momTotAntiLambda)/momTotAntiLambda.Mag();
+    lQlNeg = momNeg.Dot(momTotAntiLambda)/momTotAntiLambda.Mag();
+    lAlphaV0AntiLambda = 1.-2./(1.+lQlPos/lQlNeg);
+    
+    lPtArmV0K0s = momPos.Perp(momTotK0s);
+    lPtArmV0Lambda = momPos.Perp(momTotLambda);
+    lPtArmV0AntiLambda = momPos.Perp(momTotAntiLambda);
+    
+    // Look for associated particles:
+    if (fAnalysisMC) {
+      if( (lIndexPosMother==-1) || (lIndexNegMother==-1) ) {
+       fHistMCDaughterTrack->Fill(1);
+      }
+      
+      else if( ( (lPDGCodePosDaughter==+211) && (lPDGCodeNegDaughter==-211) )    
+              ) {
+       lCheckPIdK0Short    = 1;
+       fHistMCDaughterTrack->Fill(3);
+       if ( (lIndexPosMother==lIndexNegMother) &&
+            (lPdgcodeMother==310) ) {
+         if (lIndexPosMother <= lNbMCPrimary) lCheckMcK0Short  = 1;
+         else lCheckSecondaryK0s = 1;
+       }
+      }
+      else if( ( (lPDGCodePosDaughter==+2212) && (lPDGCodeNegDaughter==-211)  )  
+              ) {
+       lCheckPIdLambda     = 1;
+       fHistMCDaughterTrack->Fill(5);
+       if ( (lIndexPosMother==lIndexNegMother) &&
+            (lPdgcodeMother==3122)  ){
+         if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
+              ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
+              ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
+              ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
+              ) lComeFromSigma = 1;
+         else lComeFromSigma = 0; 
+         if ( (lIndexPosMother <= lNbMCPrimary) || 
+            ( ( lIndexPosMother > lNbMCPrimary) && (lComeFromSigma) )
+            ) lCheckMcLambda  = 1; 
+         else lCheckSecondaryLambda    = 1;
+       }
+      }
+      else if( ( (lPDGCodePosDaughter==211)   && (lPDGCodeNegDaughter==-2212) )             
+              ) {
+       lCheckPIdAntiLambda = 1;
+       fHistMCDaughterTrack->Fill(7);
+       if ( (lIndexPosMother==lIndexNegMother) &&
+            (lPdgcodeMother==-3122) ) {
+         if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
+              ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
+              ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
+              ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
+              ) lComeFromSigma = 1;
+         else lComeFromSigma = 0;  
+         if ( (lIndexPosMother <= lNbMCPrimary) || 
+            ( ( lIndexPosMother > lNbMCPrimary) && (lComeFromSigma) )
+            ) lCheckMcAntiLambda  = 1;
+         else lCheckSecondaryAntiLambda = 1;
+       }
+      }
+      
+      // Gamma conversion
+      else if ( (lPDGCodePosDaughter==11) &&
+               (lPDGCodeNegDaughter==-11) &&
+               (lPdgcodeMother==22 ) )
+       lCheckGamma = 1;
+    } // end "look for associated particles  
+   
+    
+    // Cuts:
+/*    if (fUseCut.Contains("yes")) {
+      if ( (lDcaPosToPrimVertex < 0.036 ) ||
+          (lDcaNegToPrimVertex < 0.036 ) ||
+          (lDcaV0Daughters     > 0.5   ) ||
+          (lV0cosPointAngle    < 0.999 ) 
+          )    
+       continue;
+    }
+*/
+
+/*
+      if ( (lDcaV0Daughters     > 0.3   ) ||
+          (lV0cosPointAngle    < 0.998 )
+
+          )    continue;
+*/
+    // PID condition:
+    lCheckPIDK0sPosDaughter        = 0, lCheckPIDK0sNegDaughter        = 0;
+    lCheckPIDLambdaPosDaughter     = 0, lCheckPIDLambdaNegDaughter     = 0;
+    lCheckPIDAntiLambdaPosDaughter = 0, lCheckPIDAntiLambdaNegDaughter = 0;
+
+    if (lMomInnerWallPos < lLimitPPID) {
+      if (nSigmaPosPion < cutNSigmaLowP)   {
+       lCheckPIDK0sPosDaughter        = 1;
+       lCheckPIDAntiLambdaPosDaughter = 1;
+      }
+      if (nSigmaPosProton < cutNSigmaLowP) lCheckPIDLambdaPosDaughter    = 1;      
+    }
+
+    else if (lMomInnerWallPos > lLimitPPID) {
+      if (nSigmaPosPion < cutNSigmaHighP)   {
+       lCheckPIDK0sPosDaughter        = 1;
+       lCheckPIDAntiLambdaPosDaughter = 1;
+      }
+      if (nSigmaPosProton < cutNSigmaHighP) lCheckPIDLambdaPosDaughter    = 1;
+    }
+
+    if (lMomInnerWallNeg < lLimitPPID) {
+      if (nSigmaNegPion < cutNSigmaLowP)    {
+       lCheckPIDK0sNegDaughter       = 1;
+       lCheckPIDLambdaNegDaughter    = 1;
+      }
+      if (nSigmaNegProton < cutNSigmaLowP)  lCheckPIDAntiLambdaNegDaughter = 1;
+      
+    }
+    else if (lMomInnerWallNeg > lLimitPPID) {
+      if (nSigmaNegPion < cutNSigmaHighP)   {
+       lCheckPIDK0sNegDaughter       = 1;
+       lCheckPIDLambdaNegDaughter    = 1;
+      }
+      if (nSigmaNegProton < cutNSigmaHighP) lCheckPIDAntiLambdaNegDaughter = 1;
+    }
+    
+    
+    //*****************************
+    // filling histograms
+
+    fHistDcaPosToPrimVertex->Fill(lDcaPosToPrimVertex,lOnFlyStatus);
+    fHistDcaNegToPrimVertex->Fill(lDcaNegToPrimVertex,lOnFlyStatus);
+    fHistDcaPosToPrimVertexZoom->Fill(lDcaPosToPrimVertex,lOnFlyStatus);
+    fHistDcaNegToPrimVertexZoom->Fill(lDcaNegToPrimVertex,lOnFlyStatus);
+    fHistRadiusV0->Fill(lV0Radius,lOnFlyStatus);
+    fHistDecayLengthV0->Fill(lV0DecayLength,lOnFlyStatus);
+    fHistDcaV0Daughters->Fill(lDcaV0Daughters,lOnFlyStatus);
+    fHistChi2->Fill(lChi2V0,lOnFlyStatus);
+    fHistCosPointAngle->Fill(lV0cosPointAngle,lOnFlyStatus);
+    if (lV0cosPointAngle >= 0.9) fHistCosPointAngleZoom->Fill(lV0cosPointAngle,lOnFlyStatus);
+    fHistChi2KFBeforeCutK0s->Fill(lChi2KFK0s,lOnFlyStatus); 
+    fHistChi2KFBeforeCutLambda->Fill(lChi2KFLambda,lOnFlyStatus);
+    fHistChi2KFBeforeCutAntiLambda->Fill(lChi2KFAntiLambda,lOnFlyStatus);
+
+
+    // Histo versus Rap and armenteros plot
+    if (!lOnFlyStatus){
+      if (lCheckMcK0Short) fHistAsMcRapK0->Fill(lRapK0s);
+      if (lCheckMcLambda) fHistAsMcRapLambda->Fill(lRapLambda);
+      if (lCheckMcAntiLambda) fHistAsMcRapLambda->Fill(lRapAntiLambda);
+//      fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0);
+//      fHistDaughterPt->Fill(lPtPos,lPtNeg);
+    }
+/*    else {
+      if (lCheckMcK0Short) fHistAsMcRapK0MI->Fill(lRapK0s);
+      if (lCheckMcLambda) fHistAsMcRapLambdaMI->Fill(lRapLambda);
+      if (lCheckMcAntiLambda) fHistAsMcRapLambdaMI->Fill(lRapAntiLambda);
+      fHistArmenterosPodolanskiMI->Fill(lAlphaV0,lPtArmV0);
+    }*/
+    
+    
+    // K0s associated histograms in |rap| < lCutRap:
+
+////////////////////////////
+    if ( lCheckPIDK0sPosDaughter && lCheckPIDK0sNegDaughter
+        && (lChi2KFK0s < cutChi2KF)) fHistPzPtBeforeK0s->Fill(TMath::Abs(lPzK0s/lPtK0s));
+/////////////////////////////
+
+    if ( lCheckPIDK0sPosDaughter && lCheckPIDK0sNegDaughter
+        && (lChi2KFK0s < cutChi2KF) && (TMath::Abs(lPzK0s/lPtK0s)<0.7) )     {
+
+
+
+       fHistPzPtAfterK0s->Fill(TMath::Abs(lPzK0s/lPtK0s));
+
+
+      
+      fHistChi2KFAfterCutK0s->Fill(lChi2KFK0s,lOnFlyStatus);
+
+      if (TMath::Abs(lRapK0s) < lCutRap) {
+
+       fHistNsigmaPosPionK0->Fill(nSigmaPosPion);
+       fHistNsigmaNegPionK0->Fill(nSigmaNegPion);
+       
+       switch (lOnFlyStatus){
+       case 0 : 
+         fHistMassK0->Fill(lInvMassK0s);
+         fHistMassVsRadiusK0->Fill(rcPosRK0s,lInvMassK0s);
+        fHistPtVsMassK0->Fill(lInvMassK0s,lPtK0s);
+
+
+//       fHistMultVsPtVsMassK0->Fill(multiplicity ,lInvMassK0s,lPtK0s);
+         if(lCheckPIdK0Short) fHistPidMcMassK0->Fill(lInvMassK0s);
+         if(lCheckMcK0Short) {
+           fHistAsMcMassK0->Fill(lInvMassK0s);
+           fHistAsMcPtK0->Fill(lPtK0s);
+
+
+           fHistAsMcPtVsMassK0->Fill(lInvMassK0s,lPtK0s);
+           if (lPtK0s <= 1) fHistAsMcPtZoomK0->Fill(lPtK0s);
+           fHistAsMcMassVsRadiusK0->Fill(rcPosRK0s,lInvMassK0s);
+           fHistAsMcResxK0->Fill(rcPosXK0s-mcPosX);
+           fHistAsMcResyK0->Fill(rcPosYK0s-mcPosY);
+           fHistAsMcReszK0->Fill(rcPosZK0s-mcPosZ);
+           fHistAsMcResrVsRadiusK0->Fill(rcPosRK0s,rcPosRK0s-mcPosR);
+           fHistAsMcReszVsRadiusK0->Fill(rcPosZK0s,rcPosZK0s-mcPosZ);
+           fHistAsMcProdRadiusK0->Fill(mcPosMotherR);
+           fHistAsMcProdRadiusXvsYK0s->Fill(mcPosMotherX,mcPosMotherY);
+           fHistAsMcResPtK0->Fill(deltaPtK0s);
+           fHistAsMcResPtVsRapK0->Fill(deltaPtK0s,lRapK0s);
+           fHistAsMcResPtVsPtK0->Fill(deltaPtK0s,lPtK0s);
+         }
+         else if (lCheckSecondaryK0s) {
+           fHistAsMcSecondaryPtVsRapK0s->Fill(lPtK0s,lRapK0s);
+           fHistAsMcSecondaryProdRadiusK0s->Fill(mcPosMotherR);
+           fHistAsMcSecondaryProdRadiusXvsYK0s->Fill(mcPosMotherX,mcPosMotherY);
+           switch (lPdgcodeMotherOfMother) {
+           case 130   : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(0.5);break; // K0L
+           case 321   : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(1.5);break; // K+
+           case -321  : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(2.5);break; // K-
+           case -3122 : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(3.5);break; //AntiLambda
+           default    : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(6.5);break;
+           }
+         }
+         break;
+         
+/*     case 1 :
+         fHistMassK0MI->Fill(lInvMassK0s);
+         fHistMassVsRadiusK0MI->Fill(rcPosRK0s,lInvMassK0s);
+         fHistPtVsMassK0MI->Fill(lInvMassK0s,lPtK0s);
+         if(lCheckPIdK0Short) fHistPidMcMassK0MI->Fill(lInvMassK0s);
+         if(lCheckMcK0Short) {
+           fHistAsMcMassK0MI->Fill(lInvMassK0s);
+           fHistAsMcPtK0MI->Fill(lPtK0s);
+           fHistAsMcPtVsMassK0MI->Fill(lInvMassK0s,lPtK0s);
+           if (lPtK0s <= 1) fHistAsMcPtZoomK0MI->Fill(lPtK0s);
+           fHistAsMcMassVsRadiusK0MI->Fill(rcPosRK0s,lInvMassK0s);
+           fHistAsMcResxK0MI->Fill(rcPosXK0s-mcPosX);
+           fHistAsMcResyK0MI->Fill(rcPosYK0s-mcPosY);
+           fHistAsMcReszK0MI->Fill(rcPosZK0s-mcPosZ);
+           fHistAsMcResrVsRadiusK0MI->Fill(rcPosRK0s,rcPosRK0s-mcPosR);
+           fHistAsMcReszVsRadiusK0MI->Fill(rcPosZK0s,rcPosZK0s-mcPosZ);
+           fHistAsMcProdRadiusK0MI->Fill(mcPosMotherR);
+           fHistAsMcProdRadiusXvsYK0sMI->Fill(mcPosMotherX,mcPosMotherY);
+           fHistAsMcResPtK0MI->Fill(deltaPtK0s);
+           fHistAsMcResPtVsRapK0MI->Fill(deltaPtK0s,lRapK0s);
+           fHistAsMcResPtVsPtK0MI->Fill(deltaPtK0s,lPtK0s);
+         }
+         else if (lCheckSecondaryK0s) {
+           fHistAsMcSecondaryPtVsRapK0sMI->Fill(lPtK0s,lRapK0s);
+           fHistAsMcSecondaryProdRadiusK0sMI->Fill(mcPosMotherR); 
+           fHistAsMcSecondaryProdRadiusXvsYK0sMI->Fill(mcPosMotherX,mcPosMotherY);
+           switch (lPdgcodeMotherOfMother) {
+           case 130   : fHistAsMcSecondaryMotherPdgCodeK0sMI->Fill(0.5);break; // K0L
+           case 321   : fHistAsMcSecondaryMotherPdgCodeK0sMI->Fill(1.5);break; // K+
+           case -321  : fHistAsMcSecondaryMotherPdgCodeK0sMI->Fill(2.5);break; // K-
+           case -3122 : fHistAsMcSecondaryMotherPdgCodeK0sMI->Fill(3.5);break; //AntiLambda
+           default    : fHistAsMcSecondaryMotherPdgCodeK0sMI->Fill(6.5);break;
+           }
+         }
+         break;        */
+       }
+      } // end rapidity condition
+    } // end nsigma condition
+    
+
+    // Associated Lambda histograms in |rap| < lCutRap
+
+////////////////////////len koly kontrole Abs(Pz/Pt)
+    if ( lCheckPIDLambdaPosDaughter && lCheckPIDLambdaNegDaughter
+        && (lChi2KFLambda < cutChi2KF)) fHistPzPtBeforeLambda->Fill(TMath::Abs(lPzLambda/lPtLambda)); 
+////////////////////////
+
+    if ( lCheckPIDLambdaPosDaughter && lCheckPIDLambdaNegDaughter
+        && (lChi2KFLambda < cutChi2KF) && (TMath::Abs(lPzLambda/lPtLambda)<0.7) )  {
+
+
+
+       fHistPzPtAfterLambda->Fill(TMath::Abs(lPzLambda/lPtLambda));      
+
+      fHistChi2KFAfterCutLambda->Fill(lChi2KFLambda,lOnFlyStatus);
+
+      if (TMath::Abs(lRapLambda) < lCutRap) {
+
+       fHistNsigmaPosProtonLambda->Fill(nSigmaPosProton);
+       fHistNsigmaNegPionLambda->Fill(nSigmaNegPion);
+       switch (lOnFlyStatus){
+       case 0 : 
+         fHistMassLambda->Fill(lInvMassLambda);
+         fHistMassVsRadiusLambda->Fill(rcPosRLambda,lInvMassLambda);
+        fHistPtVsMassLambda->Fill(lInvMassLambda,lPtLambda);
+
+
+
+
+//          fHistMultVsPtVsMassLambda->Fill(multiplicity ,lInvMassLambda,lPtLambda);
+         if(lCheckPIdLambda) fHistPidMcMassLambda->Fill(lInvMassLambda);
+         
+         if(lCheckMcLambda) {
+           fHistAsMcMassLambda->Fill(lInvMassLambda);
+           fHistAsMcPtLambda->Fill(lPtLambda);
+
+
+           fHistAsMcPtVsMassLambda->Fill(lInvMassLambda,lPtLambda);
+           if (lPtLambda <= 1) fHistAsMcPtZoomLambda->Fill(lPtLambda);
+           fHistAsMcMassVsRadiusLambda->Fill(rcPosRLambda,lInvMassLambda);
+           fHistAsMcResxLambda->Fill(rcPosXLambda-mcPosX);
+           fHistAsMcResyLambda->Fill(rcPosYLambda-mcPosY);
+           fHistAsMcReszLambda->Fill(rcPosZLambda-mcPosZ);
+           fHistAsMcResrVsRadiusLambda->Fill(rcPosRLambda,rcPosRLambda-mcPosR);
+           fHistAsMcReszVsRadiusLambda->Fill(rcPosZLambda,rcPosZLambda-mcPosZ);
+           fHistAsMcProdRadiusLambda->Fill(mcPosMotherR);
+           fHistAsMcProdRadiusXvsYLambda->Fill(mcPosMotherX,mcPosMotherY);
+           fHistAsMcResPtLambda->Fill(deltaPtLambda);
+           fHistAsMcResPtVsRapLambda->Fill(deltaPtLambda,lRapLambda);
+           fHistAsMcResPtVsPtLambda->Fill(deltaPtLambda,lPtLambda);
+           if (lComeFromSigma) fHistAsMcPtLambdaFromSigma->Fill(lPtLambda);
+           switch (lPdgcodeMotherOfMother) {
+           case 3222 : fHistAsMcMotherPdgCodeLambda->Fill(0.5); break; // Sigma +
+           case 3212 : fHistAsMcMotherPdgCodeLambda->Fill(1.5); break; // Sigma 0
+           case 3112 : fHistAsMcMotherPdgCodeLambda->Fill(2.5); break;// Sigma -
+           case 3224 : fHistAsMcMotherPdgCodeLambda->Fill(3.5); break;// Sigma(1385) +
+           case 3214 : fHistAsMcMotherPdgCodeLambda->Fill(4.5); break;// Sigma(1385) 0
+           case 3114 : fHistAsMcMotherPdgCodeLambda->Fill(5.5); break;// Sigma(1385) -
+           case 3322 : fHistAsMcMotherPdgCodeLambda->Fill(6.5); break; // Xi 0
+           case 3312 : fHistAsMcMotherPdgCodeLambda->Fill(7.5); break; // Xi -
+           case 3334 : fHistAsMcMotherPdgCodeLambda->Fill(8.5); break; // Omega
+           case -1   : fHistAsMcMotherPdgCodeLambda->Fill(9.5); break;
+           default   : fHistAsMcMotherPdgCodeLambda->Fill(10.5);break; 
+           }
+
+           
+           //printf("found Lambda RC dcaPos=%e dcaNeg=%e dcaDau=%e cosP=%e pT=%e mass=%e\n",lDcaPosToPrimVertex ,lDcaNegToPrimVertex ,lDcaV0Daughters,lV0cosPointAngle,lPtLambda,lInvMassLambda);
+           //printf("found Lambda RC Pindex=%d  Nindex=%d  Plabel=%d  Nlabel=%d\n\n",lIndexTrackPos,lIndexTrackNeg,lLabelTrackPos,lLabelTrackNeg);
+           
+         }
+         
+         else if (lCheckSecondaryLambda) {
+           fHistAsMcSecondaryPtVsRapLambda->Fill(lPtLambda,lRapLambda);
+           fHistAsMcSecondaryProdRadiusLambda->Fill(mcPosMotherR); 
+           fHistAsMcSecondaryProdRadiusXvsYLambda->Fill(mcPosMotherX,mcPosMotherY);
+           if (lComeFromSigma) fHistAsMcSecondaryPtLambdaFromSigma->Fill(lPtLambda);
+           printf(" lPdgcodeMotherOfMother= %d",lPdgcodeMotherOfMother);
+           switch (lPdgcodeMotherOfMother) {
+           case 3222 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(0.5); break;// Sigma +
+           case 3212 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(1.5); break;// Sigma 0
+           case 3112 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(2.5); break;// Sigma -
+           case 3224 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(3.5); break;// Sigma(1385) +
+           case 3214 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(4.5); break;// Sigma(1385) 0
+           case 3114 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(5.5); break;// Sigma(1385) -
+           case 3322 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(6.5); break; // Xi 0
+           case 3312 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(7.5); break; // Xi -
+           case 3334 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(8.5); break; // Omega
+           case -1   : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(9.5); break;
+           default   : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(10.5);break;
+           }
+         }
+         break;
+         
+/*     case 1 :
+         fHistMassLambdaMI->Fill(lInvMassLambda);
+         fHistMassVsRadiusLambdaMI->Fill(rcPosRLambda,lInvMassLambda);
+         fHistPtVsMassLambdaMI->Fill(lInvMassLambda,lPtLambda);
+         if(lCheckPIdLambda) fHistPidMcMassLambdaMI->Fill(lInvMassLambda);
+         
+         if(lCheckMcLambda) {
+           fHistAsMcMassLambdaMI->Fill(lInvMassLambda);
+           fHistAsMcPtLambdaMI->Fill(lPtLambda);
+           fHistAsMcPtVsMassLambdaMI->Fill(lInvMassLambda,lPtLambda);
+           fHistAsMcMassVsRadiusLambdaMI->Fill(rcPosRLambda,lInvMassLambda);
+           fHistAsMcResxLambdaMI->Fill(rcPosXLambda-mcPosX);
+           fHistAsMcResyLambdaMI->Fill(rcPosYLambda-mcPosY);
+           fHistAsMcReszLambdaMI->Fill(rcPosZLambda-mcPosZ);
+           fHistAsMcResrVsRadiusLambdaMI->Fill(rcPosRLambda,rcPosRLambda-mcPosR);
+           fHistAsMcReszVsRadiusLambdaMI->Fill(rcPosZLambda,rcPosZLambda-mcPosZ);
+           fHistAsMcProdRadiusLambdaMI->Fill(mcPosMotherR);
+           fHistAsMcProdRadiusXvsYLambdaMI->Fill(mcPosMotherX,mcPosMotherY);
+           fHistAsMcResPtLambdaMI->Fill(deltaPtLambda);
+           fHistAsMcResPtVsRapLambdaMI->Fill(deltaPtLambda,lRapLambda);
+           fHistAsMcResPtVsPtLambdaMI->Fill(deltaPtLambda,lPtLambda);
+           if (lComeFromSigma) fHistAsMcPtLambdaFromSigmaMI->Fill(lPtLambda);
+           switch (lPdgcodeMotherOfMother) {
+           case 3222 : fHistAsMcMotherPdgCodeLambdaMI->Fill(0.5); break; // Sigma +
+           case 3212 : fHistAsMcMotherPdgCodeLambdaMI->Fill(1.5); break; // Sigma 0
+           case 3112 : fHistAsMcMotherPdgCodeLambdaMI->Fill(2.5); break;// Sigma -
+           case 3224 : fHistAsMcMotherPdgCodeLambdaMI->Fill(3.5); break;// Sigma(1385) +
+           case 3214 : fHistAsMcMotherPdgCodeLambdaMI->Fill(4.5); break;// Sigma(1385) 0
+           case 3114 : fHistAsMcMotherPdgCodeLambdaMI->Fill(5.5); break;// Sigma(1385) -
+           case 3322 : fHistAsMcMotherPdgCodeLambdaMI->Fill(6.5);break; // Xi 0
+           case 3312 : fHistAsMcMotherPdgCodeLambdaMI->Fill(7.5);break; // Xi -
+           case 3334 : fHistAsMcMotherPdgCodeLambdaMI->Fill(8.5);break; // Omega
+           case -1   : fHistAsMcMotherPdgCodeLambdaMI->Fill(9.5);break;
+           default   : fHistAsMcMotherPdgCodeLambdaMI->Fill(10.5);break;
+           }
+         }
+         else if (lCheckSecondaryLambda) {
+           fHistAsMcSecondaryPtVsRapLambdaMI->Fill(lPtLambda,lRapLambda);
+           fHistAsMcSecondaryProdRadiusLambdaMI->Fill(mcPosMotherR); 
+           fHistAsMcSecondaryProdRadiusXvsYLambdaMI->Fill(mcPosMotherX,mcPosMotherY);
+           if (lComeFromSigma) fHistAsMcSecondaryPtLambdaFromSigmaMI->Fill(lPtLambda);
+           switch (lPdgcodeMotherOfMother) {
+           case 3222 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(0.5); break;// Sigma +
+           case 3212 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(1.5); break;// Sigma 0
+           case 3112 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(2.5); break;// Sigma -
+           case 3224 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(3.5); break;// Sigma(1385) +
+           case 3214 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(4.5); break;// Sigma(1385) 0
+           case 3114 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(5.5); break;// Sigma(1385) -
+           case 3322 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(6.5); break; // Xi 0
+           case 3312 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(7.5); break; // Xi -
+           case 3334 : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(8.5); break; // Omega
+           case -1   : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(9.5); break;
+           default   : fHistAsMcSecondaryMotherPdgCodeLambdaMI->Fill(10.5);break;
+           }
+         }
+         break;        
+       */
+       }
+      } // end rapidity condition
+    } //end nsigma condition - lambda
+
+
+    if (negPiKF) delete negPiKF; negPiKF= NULL;
+    if (posPiKF) delete posPiKF; posPiKF= NULL;
+    if (posPKF)  delete posPKF;  posPKF = NULL;
+    if (negAPKF) delete negAPKF; negAPKF= NULL;
+    
+  } // end V0 loop
+
+  fHistV0Multiplicity->Fill(nv0s);
+//  fHistV0MultiplicityMI->Fill(nv0sMI);
+
+  if (fAnalysisType == "ESD") { if(myPrimaryVertex) delete myPrimaryVertex; }
+
+  if(myTracksCuts) delete myTracksCuts;
+  
+  // Post output data
+  PostData(1, fListHist);
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskPerformanceStrange::Terminate(Option_t *) 
+{/*
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  TList *cRetrievedList = 0x0;
+  cRetrievedList = (TList*)GetOutputData(1);
+  
+  if(!cRetrievedList){
+    AliWarning("ERROR - AliAnalysisTaskPerformanceStrange: output data container list not available\n"); return;
+  }
+  
+  
+  fHistV0Multiplicity = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0Multiplicity"));
+  if (!fHistV0Multiplicity) {
+    Printf("ERROR: fHistV0Multiplicity not available");
+    return;
+  }
+
+  fHistV0MultiplicityMI = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0MultiplicityMI"));
+  if (!fHistV0MultiplicityMI) {
+    Printf("ERROR: fHistV0MultiplicityMI not available");
+    return;
+  }
+
+  TCanvas *canPerformanceStrange = new TCanvas("AliAnalysisTaskCheckV0","Multiplicity",10,10,510,510);
+  canPerformanceStrange->Divide(2,1);
+  if (fHistV0Multiplicity->GetMaximum() > 0.) canPerformanceStrange->cd(1)->SetLogy();
+  fHistV0Multiplicity->SetMarkerStyle(25);
+  fHistV0Multiplicity->DrawCopy("E");
+  if (fHistV0MultiplicityMI->GetMaximum() > 0.) canPerformanceStrange->cd(2)->SetLogy();
+  fHistV0MultiplicityMI->SetMarkerStyle(24);
+  fHistV0MultiplicityMI->DrawCopy("E");
+  
+
+*/ 
+}
+
+//----------------------------------------------------------------------------
+
+Double_t AliAnalysisTaskPerformanceStrange::MyRapidity(Double_t rE, Double_t rPz) const
+{
+  // Local calculation for rapidity
+  return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
+} 
+//----------------------------------------------------------------------------
+
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.h b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.h
new file mode 100644 (file)
index 0000000..ddb7789
--- /dev/null
@@ -0,0 +1,384 @@
+#ifndef ALIANALYSISTASKPERFORMANCESTRANGE_H
+#define ALIANALYSISTASKPERFORMANCESTRANGE_H
+
+/*  See cxx source for full Copyright notice */
+
+//-----------------------------------------------------------------
+//             AliAnalysisTaskPerformanceSrange class
+//    This task is for a performance study of V0 identification.
+//                It works with MC info and ESD tree.
+//                 Author: H.Ricaud, H.Ricaud@gsi.de
+//-----------------------------------------------------------------
+
+class TString;
+class TList;
+class TH1F;
+class TH2F;
+class AliAnalysisCentralitySelector;
+//class TH3F;
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskPerformanceStrange : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskPerformanceStrange();
+  AliAnalysisTaskPerformanceStrange(const char *name);
+  virtual ~AliAnalysisTaskPerformanceStrange() {}
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  void   SetCollidingSystems(Bool_t collidingSystems = 0) {fCollidingSystems = collidingSystems;}
+  void   SetAnalysisMC(Bool_t analysisMC) {fAnalysisMC = analysisMC;}
+  void   SetAnalysisType(const char* analysisType) {fAnalysisType = analysisType;}
+  void   SetUsePID(const char* usePID) {fUsePID = usePID;}
+  void   SetAnalysisCut(const char* useCut) {fUseCut = useCut;}
+  void   SetCentralitySelector(AliAnalysisCentralitySelector * centr) { fCentrSelector = centr;}
+  Double_t MyRapidity(Double_t rE, Double_t rPz) const;
+ private:
+  Double_t fCuts[7];                            //! V0 finding cuts
+  Bool_t       fAnalysisMC;                     //  0->No MC or 1->MC analysis
+  TString      fAnalysisType;                   //  "ESD" or "AOD"
+  Bool_t       fCollidingSystems;               //  Colliding systems 0/1 for pp/PbPb  
+  TString      fUsePID;                         //  "withPID" or "noPID"
+  TString      fUseCut;                         //  "yes" or "no"
+
+  AliESDEvent *fESD;                            //! ESD object
+
+  TList               *fListHist;              //! Output List
+
+  AliAnalysisCentralitySelector * fCentrSelector; // Centrality selector, used to 
+
+  // MC histograms
+  TH1F        *fHistMCPrimaryVertexX;      //! Histo
+  TH1F        *fHistMCPrimaryVertexY;      //! Histo
+  TH1F        *fHistMCPrimaryVertexZ;      //! Histo
+
+  TH1F        *fHistMCMultiplicityPrimary;       //! Histo
+  TH1F        *fHistMCMultiplicityTracks;       //! Histo
+  
+  TH2F        *fHistMCtracksProdRadiusK0s;       //! Histo
+  TH2F        *fHistMCtracksProdRadiusLambda;       //! Histo
+  TH2F        *fHistMCtracksProdRadiusAntiLambda;       //! Histo
+
+  TH1F        *fHistMCtracksDecayRadiusK0s;       //! Histo
+  TH1F        *fHistMCtracksDecayRadiusLambda;       //! Histo
+  TH1F        *fHistMCtracksDecayRadiusAntiLambda;       //! Histo
+
+  TH1F        *fHistMCPtAllK0s;       //! Histo
+  TH1F        *fHistMCPtAllLambda;       //! Histo
+  TH1F        *fHistMCPtAllAntiLambda;       //! Histo
+
+  TH1F        *fHistMCProdRadiusK0s;       //! Histo
+  TH1F        *fHistMCProdRadiusLambda;       //! Histo
+  TH1F        *fHistMCProdRadiusAntiLambda;       //! Histo
+
+  TH1F        *fHistMCRapK0s;                 //! Histo
+  TH1F        *fHistMCRapInPtRangeK0s;        //! Histo
+  TH1F        *fHistMCRapLambda;              //! Histo
+  TH1F        *fHistMCRapInPtRangeLambda;     //! Histo
+  TH1F        *fHistMCRapAntiLambda;          //! Histo
+  TH1F        *fHistMCRapInPtRangeAntiLambda; //! Histo
+  TH1F        *fHistMCRapXi;                  //! Histo
+  TH1F        *fHistMCRapInPtRangeXi;         //! Histo
+  TH1F        *fHistMCRapPhi;                 //! Histo
+  TH1F        *fHistMCRapInPtRangePhi;        //! Histo
+////////////////////////////////////////////////////////// 
+  TH1F        *fHistMCPtK0s;       //! Histo
+
+
+  TH1F        *fHistMCPtLambda;       //! Histo
+//////////////////////////////////////////////////////////
+
+  TH1F        *fHistMCPtLambdaFromSigma;       //! Histo
+  TH1F        *fHistMCPtAntiLambdaFromSigma;       //! Histo
+
+  TH1F        *fHistNTimesRecK0s;       //! Histo
+//  TH1F        *fHistNTimesRecK0sMI;       //! Histo
+  TH1F        *fHistNTimesRecLambda;       //! Histo
+//  TH1F        *fHistNTimesRecLambdaMI;       //! Histo
+  TH1F        *fHistNTimesRecAntiLambda;       //! Histo
+//  TH1F        *fHistNTimesRecAntiLambdaMI;       //! Histo
+
+  TH2F        *fHistNTimesRecK0sVsPt;       //! Histo
+//  TH2F        *fHistNTimesRecK0sVsPtMI;       //! Histo
+  TH2F        *fHistNTimesRecLambdaVsPt;       //! Histo
+//  TH2F        *fHistNTimesRecLambdaVsPtMI;       //! Histo
+  TH2F        *fHistNTimesRecAntiLambdaVsPt;       //! Histo
+//  TH2F        *fHistNTimesRecAntiLambdaVsPtMI;       //! Histo
+
+  // ESD histograms
+  TH1F        *fHistNumberEvents;        //! Histo
+  TH1F        *fHistTrackPerEvent;       //! Histo
+  TH1F        *fHistTrackletPerEvent;   //! Histo
+  TH1F        *fHistMCDaughterTrack;       //! Histo
+
+  TH1F        *fHistSPDPrimaryVertexZ;       //! Histo
+
+  TH1F        *fHistPrimaryVertexX;       //! Histo
+  TH1F        *fHistPrimaryVertexY;       //! Histo
+  TH1F        *fHistPrimaryVertexZ;       //! Histo
+
+  TH1F        *fHistPrimaryVertexResX;       //! Histo
+  TH1F        *fHistPrimaryVertexResY;       //! Histo
+  TH1F        *fHistPrimaryVertexResZ;       //! Histo
+
+  TH1F        *fHistPrimaryVertexPosXV0events;  //! Primary vertex position in X in events with V0 candidates
+  TH1F        *fHistPrimaryVertexPosYV0events;  //! Primary vertex position in Y in events with V0 candidates
+  TH1F        *fHistPrimaryVertexPosZV0events;  //! Primary vertex position in Z in events with V0 candidates
+
+  TH2F        *fHistDaughterPt;               //! Histo
+
+  TH2F        *fHistDcaPosToPrimVertex;       //! Histo
+  TH2F        *fHistDcaNegToPrimVertex;       //! Histo
+  TH2F        *fHistDcaPosToPrimVertexZoom;       //! Histo
+  TH2F        *fHistDcaNegToPrimVertexZoom;       //! Histo
+  TH2F        *fHistRadiusV0;       //! Histo
+  TH2F        *fHistDecayLengthV0;       //! Histo
+  TH2F        *fHistDcaV0Daughters;       //! Histo
+  TH2F        *fHistChi2;       //! Histo
+  TH2F        *fHistCosPointAngle;       //! Histo
+  TH2F        *fHistCosPointAngleZoom;       //! Histo
+  TH2F        *fHistProdRadius;       //! Histo
+
+//  TH2F        *fHistProdRadiusMI;       //! Histo
+
+  TH1F        *fHistV0Multiplicity;  //! Histo
+//  TH1F        *fHistV0MultiplicityMI; //! Histo
+
+  TH2F        *fHistChi2KFBeforeCutK0s;   //! Histo 
+  TH2F        *fHistChi2KFBeforeCutLambda;   //! Histo 
+  TH2F        *fHistChi2KFBeforeCutAntiLambda;   //! Histo  
+  TH2F        *fHistChi2KFAfterCutK0s;   //! Histo 
+  TH2F        *fHistChi2KFAfterCutLambda;   //! Histo 
+  TH2F        *fHistChi2KFAfterCutAntiLambda;   //! Histo
+
+  TH1F        *fHistMassK0;       //! Histo
+//  TH1F        *fHistMassK0MI;       //! Histo
+  TH1F        *fHistMassLambda;       //! Histo
+//  TH1F        *fHistMassLambdaMI;       //! Histo
+  TH1F        *fHistMassAntiLambda;       //! Histo
+//  TH1F        *fHistMassAntiLambdaMI;       //! Histo
+
+  TH2F        *fHistMassVsRadiusK0;       //! Histo
+//  TH2F        *fHistMassVsRadiusK0MI;       //! Histo
+  TH2F        *fHistMassVsRadiusLambda;       //! Histo
+//  TH2F        *fHistMassVsRadiusLambdaMI;       //! Histo
+  TH2F        *fHistMassVsRadiusAntiLambda;       //! Histo
+//  TH2F        *fHistMassVsRadiusAntiLambdaMI;       //! Histo
+
+////////////////////////////////////////////////////////////////////////////
+  TH2F        *fHistPtVsMassK0;       //! Histo
+//  TH2F        *fHistPtVsMassK0MI;       //! Histo
+  TH2F        *fHistPtVsMassLambda;       //! Histo
+/////////////////////////////////////kontrola pre |pz/pt|//////////////
+  TH1F        *fHistPzPtBeforeK0s;
+  TH1F        *fHistPzPtAfterK0s;
+  TH1F        *fHistPzPtBeforeLambda;
+  TH1F        *fHistPzPtAfterLambda;
+///////////////////////////////////////////////////////////////////////////////////
+
+//  TH2F        *fHistPtVsMassLambdaMI;       //! Histo
+//  TH2F        *fHistPtVsMassAntiLambda;       //! Histo
+//  TH2F        *fHistPtVsMassAntiLambdaMI;       //! Histo
+
+
+/////////////////////////////////////////////
+
+////////////////////////////////////////////////////////
+//  TH3F        *fHistMultVsPtVsMassK0;
+//  TH3F        *fHistMultVsPtVsMassLambda;
+//  TH3F        *fHistMultVsPtVsMassAntiLambda;
+////////////////////////////////////////////////////////
+  TH2F        *fHistArmenterosPodolanski;       //! Histo
+//  TH2F        *fHistArmenterosPodolanskiMI;       //! Histo
+
+  //PID
+  TH1F        *fHistNsigmaPosPionAntiLambda;    //! Histo
+  TH1F        *fHistNsigmaNegProtonAntiLambda;   //! Histo
+  TH1F        *fHistNsigmaPosProtonLambda;        //! Histo
+  TH1F        *fHistNsigmaNegPionLambda;           //! Histo
+  TH1F        *fHistNsigmaPosPionK0;                //! Histo
+  TH1F        *fHistNsigmaNegPionK0;                 //! Histo
+
+  // Associated particles histograms
+  TH1F        *fHistAsMcRapK0;       //! Histo
+//  TH1F        *fHistAsMcRapK0MI;       //! Histo
+  TH1F        *fHistAsMcRapLambda;       //! Histo
+//  TH1F        *fHistAsMcRapLambdaMI;       //! Histo
+  TH1F        *fHistAsMcRapAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcRapAntiLambdaMI;       //! Histo
+
+////////////////////////////////////////////////////////////////////
+  TH1F        *fHistAsMcPtK0;       //! Histo
+
+
+//  TH1F        *fHistAsMcPtK0MI;       //! Histo
+  TH1F        *fHistAsMcPtLambda;       //! Histo
+/////////////////////////////////////////////////////////////////////
+
+//  TH1F        *fHistAsMcPtAntiLambdaMI;       //! Histo
+  TH1F        *fHistAsMcPtZoomK0;       //! Histo
+//  TH1F        *fHistAsMcPtZoomK0MI;       //! Histo
+  TH1F        *fHistAsMcPtZoomLambda;       //! Histo
+//    TH1F        *fHistAsMcPtZoomLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcProdRadiusK0;       //! Histo
+//  TH1F        *fHistAsMcProdRadiusK0MI;       //! Histo
+  TH1F        *fHistAsMcProdRadiusLambda;       //! Histo
+//  TH1F        *fHistAsMcProdRadiusLambdaMI;       //! Histo
+  TH1F        *fHistAsMcProdRadiusAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcProdRadiusAntiLambdaMI;       //! Histo
+
+  TH2F        *fHistAsMcProdRadiusXvsYK0s;       //! Histo
+//  TH2F        *fHistAsMcProdRadiusXvsYK0sMI;       //! Histo
+  TH2F        *fHistAsMcProdRadiusXvsYLambda;       //! Histo
+//  TH2F        *fHistAsMcProdRadiusXvsYLambdaMI;       //! Histo
+  TH2F        *fHistAsMcProdRadiusXvsYAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcProdRadiusXvsYAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistPidMcMassK0;       //! Histo
+//  TH1F        *fHistPidMcMassK0MI;       //! Histo
+  TH1F        *fHistPidMcMassLambda;       //! Histo
+//  TH1F        *fHistPidMcMassLambdaMI;       //! Histo
+  TH1F        *fHistPidMcMassAntiLambda;       //! Histo
+//  TH1F        *fHistPidMcMassAntiLambdaMI;       //! Histo
+  TH1F        *fHistAsMcMassK0;       //! Histo
+//  TH1F        *fHistAsMcMassK0MI;       //! Histo
+  TH1F        *fHistAsMcMassLambda;       //! Histo
+//  TH1F        *fHistAsMcMassLambdaMI;       //! Histo
+  TH1F        *fHistAsMcMassAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcMassAntiLambdaMI;       //! Histo
+
+  TH2F        *fHistAsMcPtVsMassK0;       //! Histo
+//  TH2F        *fHistAsMcPtVsMassK0MI;       //! Histo
+  TH2F        *fHistAsMcPtVsMassLambda;       //! Histo
+//  TH2F        *fHistAsMcPtVsMassLambdaMI;       //! Histo
+  TH2F        *fHistAsMcPtVsMassAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcPtVsMassAntiLambdaMI;       //! Histo
+
+
+  TH2F        *fHistAsMcMassVsRadiusK0;       //! Histo
+//  TH2F        *fHistAsMcMassVsRadiusK0MI;       //! Histo
+  TH2F        *fHistAsMcMassVsRadiusLambda;       //! Histo
+//  TH2F        *fHistAsMcMassVsRadiusLambdaMI;       //! Histo
+  TH2F        *fHistAsMcMassVsRadiusAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcMassVsRadiusAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcResxK0;       //! Histo
+  TH1F        *fHistAsMcResyK0;       //! Histo
+  TH1F        *fHistAsMcReszK0;       //! Histo
+
+  TH2F        *fHistAsMcResrVsRadiusK0;       //! Histo
+  TH2F        *fHistAsMcReszVsRadiusK0;       //! Histo
+
+//  TH1F        *fHistAsMcResxK0MI;       //! Histo
+//  TH1F        *fHistAsMcResyK0MI;       //! Histo
+//  TH1F        *fHistAsMcReszK0MI;       //! Histo
+
+//  TH2F        *fHistAsMcResrVsRadiusK0MI;       //! Histo
+//  TH2F        *fHistAsMcReszVsRadiusK0MI;       //! Histo
+
+  TH1F        *fHistAsMcResxLambda;       //! Histo
+  TH1F        *fHistAsMcResyLambda;       //! Histo
+  TH1F        *fHistAsMcReszLambda;       //! Histo
+
+  TH2F        *fHistAsMcResrVsRadiusLambda;       //! Histo
+  TH2F        *fHistAsMcReszVsRadiusLambda;       //! Histo
+    
+//  TH1F        *fHistAsMcResxLambdaMI;       //! Histo
+//  TH1F        *fHistAsMcResyLambdaMI;       //! Histo
+//  TH1F        *fHistAsMcReszLambdaMI;       //! Histo
+
+//  TH2F        *fHistAsMcResrVsRadiusLambdaMI;       //! Histo
+//  TH2F        *fHistAsMcReszVsRadiusLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcResxAntiLambda;       //! Histo
+  TH1F        *fHistAsMcResyAntiLambda;       //! Histo
+  TH1F        *fHistAsMcReszAntiLambda;       //! Histo
+
+  TH2F        *fHistAsMcResrVsRadiusAntiLambda;       //! Histo
+  TH2F        *fHistAsMcReszVsRadiusAntiLambda;       //! Histo
+    
+//  TH1F        *fHistAsMcResxAntiLambdaMI;       //! Histo
+//  TH1F        *fHistAsMcResyAntiLambdaMI;       //! Histo
+//  TH1F        *fHistAsMcReszAntiLambdaMI;       //! Histo
+
+//  TH2F        *fHistAsMcResrVsRadiusAntiLambdaMI;       //! Histo
+//  TH2F        *fHistAsMcReszVsRadiusAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcResPtK0;       //! Histo
+//  TH1F        *fHistAsMcResPtK0MI;       //! Histo
+  TH1F        *fHistAsMcResPtLambda;       //! Histo
+//  TH1F        *fHistAsMcResPtLambdaMI;       //! Histo
+  TH1F        *fHistAsMcResPtAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcResPtAntiLambdaMI;       //! Histo
+
+  TH2F        *fHistAsMcResPtVsRapK0;       //! Histo
+//  TH2F        *fHistAsMcResPtVsRapK0MI;       //! Histo
+  TH2F        *fHistAsMcResPtVsRapLambda;       //! Histo
+//  TH2F        *fHistAsMcResPtVsRapLambdaMI;       //! Histo
+  TH2F        *fHistAsMcResPtVsRapAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcResPtVsRapAntiLambdaMI;       //! Histo
+  TH2F        *fHistAsMcResPtVsPtK0;       //! Histo
+ // TH2F        *fHistAsMcResPtVsPtK0MI;       //! Histo
+  TH2F        *fHistAsMcResPtVsPtLambda;       //! Histo
+//  TH2F        *fHistAsMcResPtVsPtLambdaMI;       //! Histo
+  TH2F        *fHistAsMcResPtVsPtAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcResPtVsPtAntiLambdaMI;       //! Histo
+  
+
+  TH1F        *fHistAsMcMotherPdgCodeK0s;       //! Histo
+//  TH1F        *fHistAsMcMotherPdgCodeK0sMI;       //! Histo
+  TH1F        *fHistAsMcMotherPdgCodeLambda;       //! Histo
+//  TH1F        *fHistAsMcMotherPdgCodeLambdaMI;       //! Histo
+  TH1F        *fHistAsMcMotherPdgCodeAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcMotherPdgCodeAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcPtLambdaFromSigma;       //! Histo
+//  TH1F        *fHistAsMcPtLambdaFromSigmaMI;       //! Histo
+  TH1F        *fHistAsMcPtAntiLambdaFromSigma;       //! Histo
+//  TH1F        *fHistAsMcPtAntiLambdaFromSigmaMI;       //! Histo
+
+  // Associated secondary particles:
+  TH2F        *fHistAsMcSecondaryPtVsRapK0s;       //! Histo
+//  TH2F        *fHistAsMcSecondaryPtVsRapK0sMI;       //! Histo
+  TH2F        *fHistAsMcSecondaryPtVsRapLambda;       //! Histo
+//  TH2F        *fHistAsMcSecondaryPtVsRapLambdaMI;       //! Histo
+  TH2F        *fHistAsMcSecondaryPtVsRapAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcSecondaryPtVsRapAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcSecondaryProdRadiusK0s;       //! Histo
+//  TH1F        *fHistAsMcSecondaryProdRadiusK0sMI;       //! Histo
+  TH1F        *fHistAsMcSecondaryProdRadiusLambda;       //! Histo
+//  TH1F        *fHistAsMcSecondaryProdRadiusLambdaMI;       //! Histo
+  TH1F        *fHistAsMcSecondaryProdRadiusAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcSecondaryProdRadiusAntiLambdaMI;       //! Histo
+
+  TH2F        *fHistAsMcSecondaryProdRadiusXvsYK0s;       //! Histo
+//  TH2F        *fHistAsMcSecondaryProdRadiusXvsYK0sMI;       //! Histo
+  TH2F        *fHistAsMcSecondaryProdRadiusXvsYLambda;       //! Histo
+//  TH2F        *fHistAsMcSecondaryProdRadiusXvsYLambdaMI;       //! Histo
+  TH2F        *fHistAsMcSecondaryProdRadiusXvsYAntiLambda;       //! Histo
+//  TH2F        *fHistAsMcSecondaryProdRadiusXvsYAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcSecondaryMotherPdgCodeK0s;       //! Histo
+//  TH1F        *fHistAsMcSecondaryMotherPdgCodeK0sMI;       //! Histo
+  TH1F        *fHistAsMcSecondaryMotherPdgCodeLambda;       //! Histo
+//  TH1F        *fHistAsMcSecondaryMotherPdgCodeLambdaMI;       //! Histo
+  TH1F        *fHistAsMcSecondaryMotherPdgCodeAntiLambda;       //! Histo
+//  TH1F        *fHistAsMcSecondaryMotherPdgCodeAntiLambdaMI;       //! Histo
+
+  TH1F        *fHistAsMcSecondaryPtLambdaFromSigma;       //! Histo
+//  TH1F        *fHistAsMcSecondaryPtLambdaFromSigmaMI;       //! Histo
+  TH1F        *fHistAsMcSecondaryPtAntiLambdaFromSigma;       //! Histo
+//  TH1F        *fHistAsMcSecondaryPtAntiLambdaFromSigmaMI;       //! Histo
+
+  AliAnalysisTaskPerformanceStrange(const AliAnalysisTaskPerformanceStrange&); 
+  AliAnalysisTaskPerformanceStrange& operator=(const AliAnalysisTaskPerformanceStrange&); 
+
+  ClassDef(AliAnalysisTaskPerformanceStrange, 1); 
+};
+
+#endif
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AnalysisStrange.C b/PWG2/SPECTRA/LambdaK0PbPb/AnalysisStrange.C
new file mode 100644 (file)
index 0000000..815194b
--- /dev/null
@@ -0,0 +1,155 @@
+const char *anatype = "ESD";
+
+void AnalysisStrange()
+{
+// Analysis using ESD data
+// Automatically generated analysis steering macro executed in grid subjobs
+
+   TStopwatch timer;
+   timer.Start();
+
+// Set temporary merging directory to current one
+   gSystem->Setenv("TMPDIR", gSystem->pwd());
+
+// Reset existing include path and add current directory first in the search
+   gSystem->SetIncludePath("-I.");
+// load base root libraries
+   gSystem->Load("libTree");
+   gSystem->Load("libGeom");
+   gSystem->Load("libVMC");
+   gSystem->Load("libPhysics");
+
+   gSystem->Load("libMinuit");
+
+// Load analysis framework libraries
+   if (!SetupPar("STEERBase")) return;
+   if (!SetupPar("ESD")) return;
+   if (!SetupPar("AOD")) return;
+   if (!SetupPar("ANALYSIS")) return;
+   if (!SetupPar("ANALYSISalice")) return;
+   if (!SetupPar("CORRFW")) return;
+
+// Compile other par packages
+   if (!SetupPar("OADB.par")) return;
+// include path
+   TString intPath = gInterpreter->GetIncludePath();
+   TObjArray *listpaths = intPath.Tokenize(" ");
+   TIter nextpath(listpaths);
+   TObjString *pname;
+   while ((pname=(TObjString*)nextpath())) {
+      TString current = pname->GetName();
+      if (current.Contains("AliRoot") || current.Contains("ALICE_ROOT")) continue;
+      gSystem->AddIncludePath(current);
+   }
+   if (listpaths) delete listpaths;
+   gROOT->ProcessLine(".include $ALICE_ROOT/include");
+   printf("Include path: %s\n", gSystem->GetIncludePath());
+
+// Add aditional AliRoot libraries
+
+// analysis source to be compiled at runtime (if any)
+   gROOT->ProcessLine(".L AliAnalysisTaskPerformanceStrange.cxx+g");
+
+// connect to AliEn and make the chain
+   if (!TGrid::Connect("alien://")) return;
+// read the analysis manager from file
+   TFile *file = TFile::Open("analysis.root");
+   if (!file) return;
+   TIter nextkey(file->GetListOfKeys());
+   AliAnalysisManager *mgr = 0;
+   TKey *key;
+   while ((key=(TKey*)nextkey())) {
+      if (!strcmp(key->GetClassName(), "AliAnalysisManager"))
+         mgr = (AliAnalysisManager*)file->Get(key->GetName());
+   };
+   if (!mgr) {
+      ::Error("AnalysisStrange", "No analysis manager found in file analysis.root");
+      return;
+   }
+
+   mgr->PrintStatus();
+   AliLog::SetGlobalLogLevel(AliLog::kError);
+   TChain *chain = CreateChain("wn.xml", anatype);
+
+   mgr->StartAnalysis("localfile", chain);
+   timer.Stop();
+   timer.Print();
+}
+
+//________________________________________________________________________________
+TChain* CreateChain(const char *xmlfile, const char *type="ESD")
+{
+// Create a chain using url's from xml file
+   TString filename;
+   Int_t run = 0;
+   TString treename = type;
+   treename.ToLower();
+   treename += "Tree";
+   printf("***************************************\n");
+   printf("    Getting chain of trees %s\n", treename.Data());
+   printf("***************************************\n");
+   TAlienCollection *coll = TAlienCollection::Open(xmlfile);
+   if (!coll) {
+      ::Error("CreateChain", "Cannot create an AliEn collection from %s", xmlfile);
+      return NULL;
+   }
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   TChain *chain = new TChain(treename);
+   coll->Reset();
+   while (coll->Next()) {
+      filename = coll->GetTURL();
+      if (mgr) {
+         Int_t nrun = AliAnalysisManager::GetRunFromAlienPath(filename);
+         if (nrun && nrun != run) {
+            printf("### Run number detected from chain: %d\n", nrun);
+            mgr->SetRunFromPath(nrun);
+            run = nrun;
+         }
+      }
+      chain->Add(filename);
+   }
+   if (!chain->GetNtrees()) {
+      ::Error("CreateChain", "No tree found from collection %s", xmlfile);
+      return NULL;
+   }
+   return chain;
+}
+
+//________________________________________________________________________________
+Bool_t SetupPar(const char *package) {
+// Compile the package and set it up.
+   TString pkgdir = package;
+   pkgdir.ReplaceAll(".par","");
+   gSystem->Exec(TString::Format("tar xvzf %s.par", pkgdir.Data()));
+   TString cdir = gSystem->WorkingDirectory();
+   gSystem->ChangeDirectory(pkgdir);
+   // Check for BUILD.sh and execute
+   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+      printf("*******************************\n");
+      printf("*** Building PAR archive    ***\n");
+      printf("*******************************\n");
+      if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+         ::Error("SetupPar", "Cannot build par archive %s", pkgdir.Data());
+         gSystem->ChangeDirectory(cdir);
+         return kFALSE;
+      }
+   } else {
+      ::Error("SetupPar","Cannot access PROOF-INF/BUILD.sh for package %s", pkgdir.Data());
+      gSystem->ChangeDirectory(cdir);
+      return kFALSE;
+   }
+   // Check for SETUP.C and execute
+   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+      printf("*******************************\n");
+      printf("***    Setup PAR archive    ***\n");
+      printf("*******************************\n");
+      gROOT->Macro("PROOF-INF/SETUP.C");
+   } else {
+      ::Error("SetupPar","Cannot access PROOF-INF/SETUP.C for package %s", pkgdir.Data());
+      gSystem->ChangeDirectory(cdir);
+      return kFALSE;
+   }
+   // Restore original workdir
+   gSystem->ChangeDirectory(cdir);
+   return kTRUE;
+}
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/CreateAlienHandler.C b/PWG2/SPECTRA/LambdaK0PbPb/CreateAlienHandler.C
new file mode 100644 (file)
index 0000000..89c0c12
--- /dev/null
@@ -0,0 +1,96 @@
+
+AliAnalysisGrid* CreateAlienHandler(const char * runlist, TList * listCode, const char * mode, Bool_t isMC)
+{
+  // FIXME: different paths
+  // Check if user has a valid token, otherwise make one. This has limitations.
+  // One can always follow the standard procedure of calling alien-token-init then
+  //   source /tmp/gclient_env_$UID in the current shell.
+  if (!AliAnalysisGrid::CreateToken()) return NULL;
+  AliAnalysisAlien *plugin = new AliAnalysisAlien(); 
+  plugin->SetOverwriteMode();
+  // Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+  plugin->SetRunMode(mode);
+  // Set versions of used packages
+  // FIXME: PAR FILES OPTIONAL?
+  plugin->SetAPIVersion("V1.1x");
+  plugin->SetROOTVersion("v5-27-06b");
+  //  plugin->SetAliROOTVersion("v4-21-09-AN");
+// PAR files: I'm using a modified ANALYSISalice package, so I need to load par files for everything
+   plugin->EnablePackage("STEERBase");
+   plugin->EnablePackage("ESD");
+   plugin->EnablePackage("AOD");
+   plugin->EnablePackage("CORRFW");
+   plugin->EnablePackage("ANALYSIS");
+   plugin->EnablePackage("ANALYSISalice");
+   plugin->EnablePackage("OADB");
+
+
+  // Declare input data to be processed.
+  // Method 1: Create automatically XML collections using alien 'find' command.
+  // Define production directory LFN
+   if(!isMC) 
+     plugin->SetGridDataDir("/alice/data/2010/LHC10h/");
+   else
+     plugin->SetGridDataDir("/alice/sim/LHC10h8/");
+
+   // Set data search pattern
+   if(!isMC)
+     plugin->SetDataPattern("*/pass1/*AliESDs.root");
+   else
+     plugin->SetDataPattern("AliESDs.root");
+
+   plugin->SetNrunsPerMaster(1); // One run per master job
+
+   // ...then add run numbers to be considered   
+   plugin->AddRunNumber(runlist);
+
+
+   plugin->SetGridWorkingDir("LambdaK0/");
+   plugin->SetGridOutputDir("out");
+
+   // plugin->SetDefaultOutputs(kFALSE);
+   // plugin->SetOutputFiles("EventStat_temp.root lambdak0.root");
+   // plugin->SetOutputArchive("log_archive.zip:stdout,stderr");
+
+  // Declare the analysis source files names separated by blancs. To be compiled runtime
+  // using ACLiC on the worker nodes.
+
+   // Load additional stuff
+   TString cxxToLoad, hToLoad;
+   TIterator * iter = listToLoad->MakeIterator();
+   TObjString * str = 0;
+   Bool_t first = kTRUE;
+   while (str = (TObjString*) iter->Next()) {
+     cxxToLoad = cxxToLoad + (first ? "" : " ") + str->String();
+     str->String().ReplaceAll("cxx","h");
+     hToLoad = hToLoad + (first ? "" : " ") + str->String();
+     first = kFALSE;
+   }
+   plugin->SetAnalysisSource(cxxToLoad.Data());
+   plugin->SetAdditionalLibs((hToLoad + " " + cxxToLoad).Data());
+
+
+
+  // Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+  plugin->SetAnalysisMacro("AnalysisStrange.C");
+  // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+  //plugin->SetSplitMaxInputFileNumber(100);
+  plugin->SetSplitMaxInputFileNumber(50);
+  // Optionally set number of failed jobs that will trigger killing waiting sub-jobs.
+  plugin->SetMaxInitFailed(5);
+  // Optionally resubmit threshold.
+  plugin->SetMasterResubmitThreshold(90);
+  // Optionally set time to live (default 30000 sec)
+  plugin->SetTTL(90000);
+  // Optionally set input format (default xml-single)
+  plugin->SetInputFormat("xml-single");
+  // Optionally modify the name of the generated JDL (default analysis.jdl)
+  plugin->SetJDLName("TaskStrange.jdl");
+  // Optionally modify job price (default 1)
+  plugin->SetPrice(1);      
+  // Optionally modify split mode (default 'se')    
+  plugin->SetSplitMode("se");
+  
+  return plugin;
+}
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/run.C b/PWG2/SPECTRA/LambdaK0PbPb/run.C
new file mode 100644 (file)
index 0000000..415d63e
--- /dev/null
@@ -0,0 +1,275 @@
+// #include <iostream>
+// #include "AliAnalysisManager.h"
+// #include "AliESDInputHandler.h"
+// #include "AliMCEventHandler.h"
+// #include "AliAnalysisGrid.h"
+// #include "AliCentralitySelectionTask.h"
+// #include "AliAnalysisCentralitySelector.h"
+// #include "AliAnalysisTaskPerformanceStrange.h"
+// #include "TString.h"
+// #include "TChain.h"
+// #include "TAlienCollection.h"
+// #include <fstream>
+// #include "TObjString.h"
+// #include "TIterator.h"
+// #include "TGrid.h"
+// #include "TROOT.h"
+
+// #include "CreateAlienHandler.C"
+// #include 
+
+using namespace std;
+
+enum { kMyRunModeLocal = 0, kMyRunModeCAF, kMyRunModeGRID};
+
+TList * listToLoad = new TList(); // Additional classes to be loaded, see InitAndLoadLibs
+
+TChain * GetAnalysisChain(const char * incollection);
+void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug=0) ;
+
+void run(const char * data, const char * passOrPath, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kFALSE, Int_t runMode = 0, Bool_t isMC = 0, const char* option = "",TString customSuffix = "", Int_t workers = -1, const char * gridMode="full")
+{
+  // runMode:
+  //
+  // 0 local 
+  // 1 proof
+  // 2 grid
+
+  if (nev < 0)
+    nev = 1234567890;
+  InitAndLoadLibs(runMode,workers,debug);
+
+  // Create the analysis manager
+  AliAnalysisManager * mgr = new AliAnalysisManager;
+
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler;
+  // Do I need any of this? 
+  // esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks AliESDTZERO ALIESDACORDE MuonTracks TrdTracks");
+  mgr->SetInputEventHandler(esdH);
+
+  if(isMC) {
+    AliMCEventHandler* handler = new AliMCEventHandler;
+    handler->SetPreReadMode(AliMCEventHandler::kLmPreRead);
+    mgr->SetMCtruthEventHandler(handler);
+  }
+
+
+  // If we are running on grid, we need the alien handler
+  if (runMode == kMyRunModeGRID) {
+    // Create and configure the alien handler plugin
+    TGrid::Connect("alien://");// Why do I need this? Without a get a bus error...
+    //    gROOT->LoadMacro("CreateAlienHandler.C");
+    AliAnalysisGrid *alienHandler = CreateAlienHandler(data, listToLoad, gridMode, isMC);  
+    if (!alienHandler) {
+      cout << "Cannot create alien handler" << endl;    
+      exit(1);
+    }
+    mgr->SetGridHandler(alienHandler);  
+  }
+  
+
+  // Physics selection
+  gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  AliPhysicsSelectionTask * physicsSelectionTask = AddTaskPhysicsSelection(isMC);
+
+  // Centrality
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+  AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
+
+  AliAnalysisCentralitySelector * centrSelector = new AliAnalysisCentralitySelector();
+  centrSelector->SetIsMC(isMC);
+  centrSelector->SetCentralityEstimator("V0M"); // Todo: add parameter to macro?
+
+  // Parse option strings
+  TString optionStr(option);
+  
+  // remove SAVE option if set
+  Bool_t doSave = kFALSE;
+
+  if (optionStr.Contains("SAVE"))
+    {
+      optionStr = optionStr(0,optionStr.Index("SAVE")) + optionStr(optionStr.Index("SAVE")+4, optionStr.Length());
+      doSave = kTRUE;
+    }
+  TString pathsuffix = "";
+  // Not used, but may be useful
+  Bool_t useMCKinematics = isMC;
+  if (optionStr.Contains("NOMCKIN")) {
+    cout << ">>>> Ignoring MC kinematics" << endl;
+    useMCKinematics=kFALSE;
+    pathsuffix+="_NOMCKIN";
+  }
+  
+  gROOT->ProcessLine(".L AddTaskLambdaK0PbPb.C");
+  Int_t nbin = 0; // will contain the number of centrality bins
+  AliAnalysisTaskPerformanceStrange ** task = AddTaskLambdaK0PbPb("lambdak0.root", centrSelector, nbin); // FIXME also pass cuts, centrality bin type selection(5,10% percentiles, ranges...)
+  // configure task
+  //  else if (iAODanalysis) task->SetAnalysisType("AOD");
+  // FIXME: add options to macro
+  // FIXME: put physics selection inside the task
+  cout << nbin << endl;
+  
+  for(Int_t ibin = 0; ibin < nbin; ibin++){
+    cout << "ibin " << ibin << "  "<< endl;//task[ibin] << endl;
+    
+    task[ibin]->SetAnalysisType("ESD");
+    cout << "1" << endl;
+    task[ibin]->SetAnalysisMC(isMC); // 0 or 1
+    cout << "2" << endl;
+    task[ibin]->SetCollidingSystems(1); // 0 =pp, 1=AA
+    cout << "3" << endl;
+    task[ibin]->SetAnalysisCut("no");
+    cout << "4" << endl;
+    task[ibin]->SetUsePID("withoutPID"); // withPID or withoutPID
+    cout << "5" << endl;
+  }
+
+  // Init and run the analy
+  if (!mgr->InitAnalysis()) return;
+
+  mgr->PrintStatus();
+  
+  if (runMode == kMyRunModeLocal ) {
+    // If running in local mode, create chain of ESD files
+    cout << "RUNNING LOCAL, CHAIN" << endl;    
+    TChain * chain = GetAnalysisChain(data);
+    //    chain->Print();
+    mgr->StartAnalysis("local",chain,nev);
+  } else if (runMode == kMyRunModeCAF) {
+    mgr->StartAnalysis("proof",TString(passOrPath)+TString(data)+"#esdTree",nev);
+  } else if (runMode == kMyRunModeGRID) {
+    mgr->StartAnalysis("grid");
+  } else {
+    cout << "ERROR: unknown run mode" << endl;        
+  }
+
+  pathsuffix += customSuffix;
+
+  if (doSave) MoveOutput(data, pathsuffix.Data());
+
+  
+}
+
+void MoveOutput(const char * data, const char * suffix = ""){
+
+  TString path("output/");
+  path = path + TString(data).Tokenize("/")->Last()->GetName() + suffix;
+  
+  TString fileName = "lambdak0.root";
+  gSystem->mkdir(path, kTRUE);
+  gSystem->Rename(fileName, path + "/" + fileName);
+  for(Int_t ibin = 0; ibin < 20; ibin++){
+    TString fileBin = fileName;
+    fileBin.ReplaceAll(".root",Form("_%2.2d.root",ibin));
+    gSystem->Rename(fileBin, path + "/" + fileBin);    
+  }
+  
+  gSystem->Rename("event_stat.root", path + "/event_stat.root");      
+  gSystem->Rename("EventStat_temp.root", path + "/EventStat_temp.root");      
+  Printf(">>>>> Moved files to %s", path.Data());
+}  
+
+
+
+TChain * GetAnalysisChain(const char * incollection){
+  // Builds a chain of esd files
+  // incollection can be
+  // - a single root file
+  // - an xml collection of files on alien
+  // - a ASCII containing a list of local root files
+
+  TChain* analysisChain = 0;
+  // chain
+  analysisChain = new TChain("esdTree");
+  if (TString(incollection).Contains(".root")){
+    analysisChain->Add(incollection);
+  }
+  else if (TString(incollection).Contains("xml")){
+    TGrid::Connect("alien://");
+    TGridCollection * coll = TAlienCollection::Open (incollection);
+    while(coll->Next()){
+      analysisChain->Add(TString("alien://")+coll->GetLFN());
+    }
+  } else {
+    ifstream file_collect(incollection);
+    TString line;
+    while (line.ReadLine(file_collect) ) {
+      analysisChain->Add(line.Data());
+    }
+  }
+  analysisChain->GetListOfFiles()->Print();
+
+  return analysisChain;
+}
+
+
+void InitAndLoadLibs(Int_t runMode, Int_t workers,Bool_t debug) {
+  // Loads libs and par files + custom task and classes (the order is important)
+  listToLoad->Add(new TObjString("$ALICE_ROOT/STEER/AliCentrality.cxx")); // FIXME: why do I have to load it?!?
+  listToLoad->Add(new TObjString("AliAnalysisCentralitySelector.cxx"));
+  listToLoad->Add(new TObjString("AliAnalysisTaskPerformanceStrange.cxx"));
+
+  if (runMode == kMyRunModeCAF)
+    {
+      cout << "Init in CAF mode" << endl;
+    
+      gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+      TProof * p = TProof::Open("alice-caf.cern.ch", workers>0 ? Form("workers=%d",workers) : "");
+      //TProof * p = TProof::Open("skaf.saske.sk", workers>0 ? Form("workers=%d",workers) : "");    
+      p->Exec("TObject *o = gEnv->GetTable()->FindObject(\"Proof.UseMergers\"); gEnv->GetTable()->Remove(o);", kTRUE); // avoid submerging
+      //      gProof->EnablePackage("VO_ALICE@AliRoot::v4-21-11-AN");
+
+      // Enable the needed package
+      // FIXME: what if I don't want to use par files?
+      gSystem->AddIncludePath("-I${ALICE_ROOT}/include/");
+      gSystem->AddIncludePath("-I${ALICE_ROOT}/STEER/");
+      gProof->UploadPackage("$ALICE_ROOT/obj/STEERBase");
+      gProof->EnablePackage("$ALICE_ROOT/obj/STEERBase");
+      gProof->UploadPackage("$ALICE_ROOT/obj/ESD");
+      gProof->EnablePackage("$ALICE_ROOT/obj/ESD");
+      gProof->UploadPackage("$ALICE_ROOT/obj/AOD");
+      gProof->EnablePackage("$ALICE_ROOT/obj/AOD");
+      gProof->UploadPackage("$ALICE_ROOT/obj/ANALYSIS");
+      gProof->EnablePackage("$ALICE_ROOT/obj/ANALYSIS");
+      gProof->UploadPackage("$ALICE_ROOT/obj/ANALYSISalice");
+      gProof->EnablePackage("$ALICE_ROOT/obj/ANALYSISalice");
+      gProof->UploadPackage("$ALICE_ROOT/obj/CORRFW");
+      gProof->EnablePackage("$ALICE_ROOT/obj/CORRFW");
+      gProof->UploadPackage("~/Desktop/OADB");//FIXME
+      gProof->EnablePackage("~/Desktop/OADB");//FIXME
+      
+    }
+  else
+    {
+      cout << "Init in Local or Grid mode" << endl;
+      gSystem->Load("libCore.so");  
+      gSystem->Load("libTree.so");
+      gSystem->Load("libGeom.so");
+      gSystem->Load("libVMC.so");
+      gSystem->Load("libPhysics.so");
+      gSystem->Load("libSTEERBase");
+      gSystem->Load("libESD");
+      gSystem->Load("libAOD");
+      gSystem->Load("libANALYSIS");
+      gSystem->Load("libANALYSISalice");   
+
+      // gSystem->Load("libVMC");
+      gROOT->ProcessLine(".include $ALICE_ROOT/include");
+      gROOT->ProcessLine(".include $ALICE_ROOT/STEER");
+    }
+  // Load helper classes
+  TIterator * iter = listToLoad->MakeIterator();
+  TObjString * name = 0;
+  while ((name = (TObjString *)iter->Next())) {
+    gSystem->ExpandPathName(name->String());
+    cout << name->String().Data() << endl;
+    if (runMode == kMyRunModeCAF) {
+      gProof->Load(name->String()+(debug?"++g":"+"));   
+    } else {
+      gROOT->LoadMacro(name->String()+(debug?"++g":"+"));   
+    }
+  }
+
+}
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/run.sh b/PWG2/SPECTRA/LambdaK0PbPb/run.sh
new file mode 100755 (executable)
index 0000000..dfedcd0
--- /dev/null
@@ -0,0 +1,119 @@
+#!/bin/bash
+
+run=137161
+pass=pass1_5plus
+mc=0
+mode="full"
+nev=1234566789
+workers=26
+ROPT=""
+listfile=""
+offset=0
+debug=kTRUE
+option="SAVE" #FIXME:set option
+suffix=""
+
+give_help() {
+
+cat <<ENDOFGUIDE
+This scripts runs the the physics selection and centrality on a specified run/dataset
+
+Available options:
+ Mode control, at least one of the following options should be used
+  -r <mode>                    Run the task
+                               Modes (compulsory):
+                                  0 local
+                                  1 caf
+                                  2 grid    
+  -d <run or dataset>          Run number(s) (grid) or dataset (caf) or file name (local)
+                               Local filename can be an xml collection of files on alies, 
+                               a single esd, or a text file with an ESD filename per line
+  -m                           Set if runnning on MC
+  -t <rootopt>                 Options passed to root
+  -l <list.txt>                Process sequentially all runs/dataset listed in the file 
+                               (one entry per line). If you use this option you don't 
+                               need the -d. in the case you are running on CAF, they 
+                               must have the same path
+ Grid only options
+  -g <gridmode>                Plugin Mode [default=$mode]
+  -p <recopass>                Reconstruction pass [default=$pass]       
+ CAF only options
+  -p <path>                    Data set path
+  -n <nev>                     Number of events
+  -w <workers>                 Number of workers [default=$workers]
+ENDOFGUIDE
+
+}
+
+while getopts "r:hd:mg:p:n:w:t:l:" opt; do
+  case $opt in
+    r)
+      runMode=$OPTARG
+      ;;
+    d)
+      run=$OPTARG
+      ;;
+    l)
+      listfile=$OPTARG
+      ;;
+    t)
+      ROPT=$OPTARG
+      ;;
+    n)
+      nev=$OPTARG
+      ;;
+    w)
+      workers=$OPTARG
+      ;;
+    m)
+      mc=kTRUE
+      ;;
+    g) 
+      mode=$OPTARG
+      ;;
+    p)
+      pass=$OPTARG
+      ;;
+    h)
+      give_help
+      exit 1
+      ;;
+    \?)
+      echo "Invalid option: -$OPTARG" >&2
+      give_help
+      exit 1
+      ;;
+    :)
+      echo "Option -$OPTARG requires an argument." >&2
+      give_help
+      exit 1
+      ;;
+  esac
+done
+
+
+runlist=$run
+if [ "$listfile" != "" ]
+    then
+    runlist=""
+    while read line
+    do
+       runlist="$runlist $line"
+    done < $listfile
+    
+fi
+
+echo "Run list: $runlist"
+
+
+if [ "$runMode" = "2" ]
+then
+    echo root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,$option,$suffix,$workers,\"$mode\"\)
+    root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,\"$option\",\"$suffix\",$workers,\"$mode\"\)
+else
+    for run in $runlist 
+    do
+       echo root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,$option,$suffix,$workers,\"$mode\"\)
+       root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,\"$option\",\"$suffix\",$workers,\"$mode\"\)
+    done
+fi