diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/CMakeLists.txt b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6ac8bd25ed028d7dcb73160f75b8eaa473ee37d5
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/CMakeLists.txt
@@ -0,0 +1,41 @@
+################################################################################
+# Package: InDetZVTOPVxFinder
+################################################################################
+
+# Declare the package name:
+atlas_subdir( InDetZVTOPVxFinder )
+
+# declare the package's dependencies:
+atlas_depends_on_subdirs( PUBLIC
+                          Control/AthenaBaseComps
+                          GaudiKernel
+                          InnerDetector/InDetRecTools/InDetRecToolInterfaces
+                          Tracking/TrkEvent/VxVertex
+                          Tracking/TrkEvent/TrkTrack
+                          Reconstruction/Particle
+                          Tracking/TrkEvent/TrkParticleBase
+                          Tracking/TrkEvent/TrkParameters
+                          Event/xAOD/xAODTracking
+                          PRIVATE
+                          External/AtlasCLHEP
+                          Control/DataModel
+                          InnerDetector/InDetConditions/InDetBeamSpotService
+                          Tracking/TrkExtrapolation/TrkExInterfaces
+                          Tracking/TrkDetDescr/TrkSurfaces
+                          Tracking/TrkVertexFitter/TrkVertexFitterInterfaces
+                          Tracking/TrkEvent/VxSecVertex
+                          Event/EventPrimitives
+                          Tracking/TrkVertexFitter/TrkVxEdmCnv )
+
+# External dependencies
+find_package( HepMC )
+
+# Component(s) in the package:
+atlas_add_component( InDetZVTOPVxFinder
+                     src/*.cxx
+                     src/components/*.cxx
+                     INCLUDE_DIRS ${HEPMC_INCLUDE_DIRS}
+                     LINK_LIBRARIES ${HEPMC_LIBRARIES} AthenaBaseComps GaudiKernel InDetRecToolInterfaces VxVertex Particle TrkParticleBase TrkParameters xAODTracking DataModel InDetBeamSpotService TrkExInterfaces TrkSurfaces TrkVertexFitterInterfaces VxSecVertex EventPrimitives TrkVxEdmCnvLib )
+
+# Install files from the package:
+atlas_install_headers( InDetZVTOPVxFinder )
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_AmbiguitySolver.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_AmbiguitySolver.h
new file mode 100755
index 0000000000000000000000000000000000000000..7471f1f60e22d8b331007a6ce7b50b8c527f14a2
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_AmbiguitySolver.h
@@ -0,0 +1,43 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// IZVTOP_AmbiguitySolver.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef IINDETZVTOP_AMBIGUITYSOLVER_H
+#define IINDETZVTOP_AMBIGUITYSOLVER_H
+
+#include "GaudiKernel/IAlgTool.h"
+
+//xAOD includes --David S.
+#include "xAODTracking/Vertex.h"
+
+//namespace Trk --David S.
+//{
+// class VxCandidate;
+//}
+
+namespace InDet 
+{
+
+
+  static const InterfaceID IID_IZVTOP_AmbiguitySolver("InDet::IZVTOP_AmbiguitySolver", 1, 0);
+
+  class IZVTOP_AmbiguitySolver : virtual public IAlgTool {
+  public:
+    static const InterfaceID& interfaceID( ) ;
+
+    // enter declaration of your interface-defining member functions here
+    //virtual  std::vector<Trk::VxCandidate*> solveAmbiguities(std::vector<Trk::VxCandidate*> VxContainer) = 0; --David S.
+    virtual std::vector< xAOD::Vertex* > solveAmbiguities(std::vector< xAOD::Vertex* > VertexContainer) = 0;
+  };
+
+  inline const InterfaceID& InDet::IZVTOP_AmbiguitySolver::interfaceID()
+    { 
+      return IID_IZVTOP_AmbiguitySolver; 
+    }
+
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h
new file mode 100755
index 0000000000000000000000000000000000000000..015ce332649afa9b58b508a169039eda4534423a
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h
@@ -0,0 +1,35 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// IZVTOP_SimpleVtxProbCalc.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef IINDETZVTOP_SIMPLEVTXPROBCALC_H
+#define IINDETZVTOP_SIMPLEVTXPROBCALC_H
+
+#include "GaudiKernel/IAlgTool.h"
+
+
+namespace InDet 
+{
+
+
+  static const InterfaceID IID_IZVTOP_SimpleVtxProbCalc("InDet::IZVTOP_SimpleVtxProbCalc", 1, 0);
+
+  class IZVTOP_SimpleVtxProbCalc : virtual public IAlgTool {
+  public:
+    static const InterfaceID& interfaceID( ) ;
+
+    // enter declaration of your interface-defining member functions here
+    virtual double calcVtxProb(double & trk_func_sum, double & trk_func_sum2) = 0;
+  };
+
+  inline const InterfaceID& InDet::IZVTOP_SimpleVtxProbCalc::interfaceID()
+    { 
+      return IID_IZVTOP_SimpleVtxProbCalc; 
+    }
+
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_SpatialPointFinder.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_SpatialPointFinder.h
new file mode 100755
index 0000000000000000000000000000000000000000..86fab39fa4e89ce03656c0198a6fd5b7f082d649
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_SpatialPointFinder.h
@@ -0,0 +1,52 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// IZVTOP_SpatialPointFinder.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef IINDETZVTOP_SPATIALPOINTFINDER_H
+#define IINDETZVTOP_SPATIALPOINTFINDER_H
+
+#include "GaudiKernel/IAlgTool.h"
+#include "TrkParameters/TrackParameters.h"
+
+namespace Rec {
+  class TrackParticle;
+}
+namespace Trk {
+  class Track;
+  class TrackParticleBase;
+  class RecVertex;
+  class Vertex;
+  //class TrackParameters;
+}
+namespace InDet 
+{
+
+
+  static const InterfaceID IID_IZVTOP_SpatialPointFinder("InDet::IZVTOP_SpatialPointFinder", 1, 0);
+
+  class IZVTOP_SpatialPointFinder : virtual public IAlgTool {
+  public:
+    static const InterfaceID& interfaceID( ) ;
+
+    // declaration of interface-defining member functions
+    virtual Trk::Vertex* findSpatialPoint(const Trk::Track* trk_1, const Trk::Track* trk_2) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Rec::TrackParticle* trk_1, const Rec::TrackParticle* trk_2) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Trk::TrackParticleBase* trk_1, const Trk::TrackParticleBase* trk_2) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::Track* trk_1) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Rec::TrackParticle* trk_1) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParticleBase* trk_1) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Trk::TrackParameters* perigee_1, const Trk::TrackParameters* perigee_2) = 0;
+    virtual Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParameters* perigee_1) = 0;
+  };
+
+  inline const InterfaceID& InDet::IZVTOP_SpatialPointFinder::interfaceID()
+    { 
+      return IID_IZVTOP_SpatialPointFinder; 
+    }
+
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h
new file mode 100755
index 0000000000000000000000000000000000000000..e28bd3f4b01d5e55f114af93b7ff91b4f84e17a2
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h
@@ -0,0 +1,53 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// IZVTOP_TrkProbTubeCalc.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef IINDETZVTOP_TRKPROBTUBECALC_H
+#define IINDETZVTOP_TRKPROBTUBECALC_H
+
+#include "GaudiKernel/IAlgTool.h"
+#include "TrkParameters/TrackParameters.h"
+
+namespace Trk
+{
+  class Track;
+  class RecVertex;
+  class Vertex;
+  class TrackParticleBase;
+}
+
+namespace Rec {
+  class TrackParticle;
+}
+namespace InDet 
+{
+
+
+  static const InterfaceID IID_IZVTOP_TrkProbTubeCalc("InDet::IZVTOP_TrkProbTubeCalc", 1, 0);
+
+  class IZVTOP_TrkProbTubeCalc : virtual public IAlgTool {
+  public:
+    static const InterfaceID& interfaceID( ) ;
+
+    // declaration of interface-defining member functions
+    //trk
+    virtual double calcProbTube(const Trk::Track& trk, Trk::Vertex& vec) = 0;
+    virtual double calcProbTube(const Rec::TrackParticle& trk, Trk::Vertex& vec) = 0;
+    virtual double calcProbTube(const Trk::TrackParticleBase& trk, Trk::Vertex& vec) = 0;
+    virtual double calcProbTube(const Trk::Perigee* trk, Trk::Vertex& vec) = 0;
+    //beam spot
+    virtual double calcProbTube(const Trk::RecVertex& vtx, Trk::Vertex& vec) = 0;
+
+  };
+
+  inline const InterfaceID& InDet::IZVTOP_TrkProbTubeCalc::interfaceID()
+    { 
+      return IID_IZVTOP_TrkProbTubeCalc; 
+    }
+
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_AmbiguitySolver.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_AmbiguitySolver.h
new file mode 100755
index 0000000000000000000000000000000000000000..4db855fdf749834cae4a454a7cb111852b89e481
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_AmbiguitySolver.h
@@ -0,0 +1,61 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_AmbiguitySolver.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_AMBIGUITYSOLVER_H
+#define INDETZVTOP_AMBIGUITYSOLVER_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "InDetZVTOPVxFinder/IZVTOP_AmbiguitySolver.h"
+namespace Trk
+{
+  //class VxCandidate; --David S.
+  class Track;
+}
+namespace InDet 
+{
+
+  class IZVTOP_TrkProbTubeCalc;
+  class IZVTOP_SimpleVtxProbCalc;
+
+  /** @class ZVTOP_AmbiguitySolver 
+
+      This is for the Doxygen-Documentation.  
+      Please delete these lines and fill in information about
+      the Algorithm!
+      Please precede every member function declaration with a
+      short Doxygen comment stating the purpose of this function.
+      
+      @author  Tatjana Lenz <tatjana.lenz@cern.ch>
+  */  
+
+  class ZVTOP_AmbiguitySolver : virtual public IZVTOP_AmbiguitySolver, public AthAlgTool
+    {
+    public:
+      ZVTOP_AmbiguitySolver(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_AmbiguitySolver ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+      //std::vector<Trk::VxCandidate*> solveAmbiguities(std::vector<Trk::VxCandidate*> VxContainer); --David S.
+      std::vector< xAOD::Vertex* > solveAmbiguities(std::vector< xAOD::Vertex* > VertexContainer);
+      
+    private:
+
+      ToolHandle <InDet::IZVTOP_TrkProbTubeCalc>	m_TrkProbTubeCalc;
+      ToolHandle <InDet::IZVTOP_SimpleVtxProbCalc>	m_VtxProbCalc;
+
+      
+    }; 
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SecVtxTool.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SecVtxTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..ad2f420fa67f30cd850d5e56d6a1fdc44f5d5d0a
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SecVtxTool.h
@@ -0,0 +1,100 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_Tool.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_SECVTXTOOL_H
+#define INDETZVTOP_SECVTXTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "InDetRecToolInterfaces/ISecVertexInJetFinder.h"
+
+/* Forward declarations */
+namespace Trk
+{
+  class Track;
+  class VxSecVertexInfo;
+  class RecVertex;
+  class TrackParticleBase;
+  class IVertexFitter;
+}
+namespace InDet 
+{
+    class ISecVertexInJetFinder;
+    class IZVTOP_SpatialPointFinder;
+    class IZVTOP_TrkProbTubeCalc;
+    class IZVTOP_SimpleVtxProbCalc;
+    class IZVTOP_AmbiguitySolver;
+    class ZVTOP_TrkPartBaseVertexState;
+    class ZVTOP_VertexState;
+
+  /** @class ZVTOP_Tool 
+
+      ---Topological Vertex Finder Tool ---
+      This tool reconstructs a set of topological vertices each 
+      associated with an independent subset of the charged tracks.
+      Vertices are reconstructed by associating tracks with 3D spatial regions
+      according to the vertex probability function which is based on the trajectories
+      and position resolution of the tracks.
+      Docu: "A Topological Vertex Reconstruction Algorithm for Hadronic Jets" 
+            by David J. Jackson, SLAC-PUB-7215, December 1996 
+
+      Following properties can be set/changed via jobOptions:
+ 
+      ---Input: TrackCollection
+      - \c TrackParticleName: The name of the StoreGate input container from
+      which the TrackParticle are read. The default is "???", the container
+      from the legacy converters/ambiguity processor.
+      
+      ---Output: VxContainer
+      - \c VxCollectionOutputName: The name of the StoreGate container where the fit
+      results are put. default is "VxCollection".
+
+      ---Vertex Fitter: Billoir FullVertexFitter
+      - \c FitRoutine: The routine which should be used for the fitting. The
+      default is "FullVertexFitter".
+      
+
+      
+      @author  Tatjana Lenz <tatjana.lenz@NOSPAMcern.ch>
+  */  
+  class ZVTOP_SecVtxTool : virtual public ISecVertexInJetFinder, public AthAlgTool
+    {
+    public:
+      ZVTOP_SecVtxTool(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_SecVtxTool ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+ 
+     virtual const Trk::VxSecVertexInfo* findSecVertex(const Trk::RecVertex & primaryVertex,const TLorentzVector & jetMomentum,const std::vector<const Trk::TrackParticleBase*> & inputTracks) const;
+
+     //Added purely to satisfy new inheritance in ISecVertexInJetFinder, not yet implemented --David S.
+     virtual const Trk::VxSecVertexInfo* findSecVertex(const xAOD::Vertex & primaryVertex, const TLorentzVector & jetMomentum,const std::vector<const xAOD::IParticle*> & inputTracks) const;
+
+     virtual const Trk::VxSecVertexInfo* findSecVertex(const Trk::RecVertex & primaryVertex,const TLorentzVector & jetMomentum,const std::vector<const Trk::Track*> & inputTracks) const;
+
+    private:
+      
+      ToolHandle <Trk::IVertexFitter>					m_iVertexFitter;
+      ToolHandle <InDet::IZVTOP_SpatialPointFinder>	m_iSpatialPointFinder;
+      ToolHandle <InDet::IZVTOP_TrkProbTubeCalc>	m_iTrkProbTubeCalc;
+      ToolHandle <InDet::IZVTOP_SimpleVtxProbCalc>	m_iVtxProbCalc;
+      ToolHandle <InDet::IZVTOP_AmbiguitySolver>	m_iAmbiguitySolver;
+      double										m_resolvingCut;
+      double										m_resolvingStep;
+      double										m_vtxProb_cut;
+      double										m_trk_chi2_cut;
+    }; 
+} // end of namespace
+
+#endif 
+ 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SimpleVtxProbCalc.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SimpleVtxProbCalc.h
new file mode 100755
index 0000000000000000000000000000000000000000..ee4d80f2c91321e54bfc0578e1dbcc3583342b9a
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SimpleVtxProbCalc.h
@@ -0,0 +1,45 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_SimpleVtxProbCalc.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_SIMPLEVTXPROBCALC_H
+#define INDETZVTOP_SIMPLEVTXPROBCALC_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h"
+
+
+namespace InDet 
+{
+
+  /** @class ZVTOP_SimpleVtxProbCalc 
+      
+      @author  Tatjana Lenz <tatjana.lenz@cern.ch>
+  */  
+
+  class ZVTOP_SimpleVtxProbCalc : virtual public IZVTOP_SimpleVtxProbCalc, public AthAlgTool
+    {
+    public:
+      ZVTOP_SimpleVtxProbCalc(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_SimpleVtxProbCalc ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+      
+      double calcVtxProb(double & trk_func_sum, double & trk_func_sum2);
+    private:
+
+     
+      
+    }; 
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SlowSpatialPointFinder.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SlowSpatialPointFinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..6e1fbfa88f290fb3e7090f32734a7bf16f10604e
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SlowSpatialPointFinder.h
@@ -0,0 +1,78 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_SlowSpatialPointFinder.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_SLOWSPATIALPOINTFINDER_H
+#define INDETZVTOP_SLOWSPATIALPOINTFINDER_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SpatialPointFinder.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "TrkParameters/TrackParameters.h"
+namespace Trk
+{
+  class Track;
+  class RecVertex;
+  class TrackParticleBase;
+  class IVertexLinearizedTrackFactory;
+  class Vertex;
+}
+
+namespace Rec 
+{
+  class TrackParticle;
+}
+
+namespace InDet 
+{
+
+  /** @class ZVTOP_SlowSpatialPointFinder 
+
+
+      HelperClass for Topological Vertex Finder -ZVTOP-
+      Calculates the spatial point, which define the region of interest,
+      the region where we later look for the maxima of the vertex probability
+      function. Since at least two of the tracks (or RecVertex and Trk::Track) 
+      must contribute to a maximum in the vertex probability function, 
+      we calculate for each track pair the spatial point.
+      The found spatial point is rejected if chi² of this point is larger as chi²_cut.
+      
+      @author  Tatjana Lenz <tatjana.lenz@cern.ch>
+  */  
+
+  class ZVTOP_SlowSpatialPointFinder : virtual public IZVTOP_SpatialPointFinder, public AthAlgTool
+    {
+    public:
+      ZVTOP_SlowSpatialPointFinder(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_SlowSpatialPointFinder ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+
+      Trk::Vertex* findSpatialPoint(const Trk::Track* trk_1, const Trk::Track* trk_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::Track* trk_1);
+      Trk::Vertex* findSpatialPoint(const Rec::TrackParticle* trk_1, const Rec::TrackParticle* trk_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Rec::TrackParticle* trk_1);
+      Trk::Vertex* findSpatialPoint(const Trk::TrackParticleBase* trk_1, const Trk::TrackParticleBase* trk_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParticleBase* trk_1);
+      Trk::Vertex* findSpatialPoint(const Trk::TrackParameters* perigee_1, const Trk::TrackParameters* perigee_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParameters* perigee_1);
+
+    private:
+      
+    	ToolHandle< Trk::IVertexLinearizedTrackFactory > m_linFactory;
+      double				m_chi2;
+      
+    }; 
+} // end of namespace
+
+#endif 
+ 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SpatialPointFinder.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SpatialPointFinder.h
new file mode 100755
index 0000000000000000000000000000000000000000..1c8f2f955cbbadcd380f544bf5f530c62f21946b
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_SpatialPointFinder.h
@@ -0,0 +1,72 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_SpatialPointFinder.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_SPATIALPOINTFINDER_H
+#define INDETZVTOP_SPATIALPOINTFINDER_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SpatialPointFinder.h"
+
+namespace Trk
+{
+  class Track;
+  class RecVertex;
+  class TrackParticleBase;
+  class Vertex;
+}
+
+namespace Rec {
+  class TrackParticle;
+}
+
+namespace InDet 
+{
+
+  /** @class ZVTOP_SpatialPointFinder 
+
+
+      HelperClass for Topological Vertex Finder -ZVTOP-
+      Calculates the spatial point, which define the region of interest,
+      the region where we later look for the maxima of the vertex probability
+      function. Since at least two of the tracks (or RecVertex and Trk::Track) 
+      must contribute to a maximum in the vertex probability function, 
+      we calculate for each track pair the spatial point.
+      The found spatial point is rejected if chi² of this point is larger as chi²_cut.
+      
+      @author  Tatjana Lenz <tatjana.lenz@cern.ch>
+  */  
+
+  class ZVTOP_SpatialPointFinder : virtual public IZVTOP_SpatialPointFinder, public AthAlgTool
+    {
+    public:
+      ZVTOP_SpatialPointFinder(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_SpatialPointFinder ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+
+      Trk::Vertex* findSpatialPoint(const Trk::Track* trk_1, const Trk::Track* trk_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::Track* trk_1);
+      Trk::Vertex* findSpatialPoint(const Rec::TrackParticle* trk_1, const Rec::TrackParticle* trk_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Rec::TrackParticle* trk_1);
+      Trk::Vertex* findSpatialPoint(const Trk::TrackParticleBase* trk_1, const Trk::TrackParticleBase* trk_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParticleBase* trk_1);
+      Trk::Vertex* findSpatialPoint(const Trk::TrackParameters* perigee_1, const Trk::TrackParameters* perigee_2);
+      Trk::Vertex* findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParameters* perigee_1);
+
+    private:
+      double				m_chi2;
+      
+    }; 
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_Tool.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_Tool.h
new file mode 100755
index 0000000000000000000000000000000000000000..8b125b6ba1ffe6dfaca822f8e3d8a7e60fa86f3f
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_Tool.h
@@ -0,0 +1,95 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_Tool.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_TOOL_H
+#define INDETZVTOP_TOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "InDetRecToolInterfaces/IVertexFinder.h"
+#include "xAODTracking/Vertex.h"
+#include "xAODTracking/TrackParticle.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/TrackParticleContainer.h"
+
+class IBeamCondSvc;
+
+namespace Trk
+{
+  class IVertexFitter;
+}
+
+namespace InDet 
+{
+
+  class IZVTOP_SpatialPointFinder;
+  class IZVTOP_TrkProbTubeCalc;
+  class IZVTOP_SimpleVtxProbCalc;
+
+  /** @class ZVTOP_Tool 
+
+      ---Topological Vertex Finder Tool ---
+      This tool reconstructs a set of topological vertices each 
+      associated with an independent subset of the charged tracks.
+      Vertices are reconstructed by associating tracks with 3D spatial regions
+      according to the vertex probability function which is based on the trajectories
+      and position resolution of the tracks.
+      Docu: "A Topological Vertex Reconstruction Algorithm for Hadronic Jets" 
+            by David J. Jackson, SLAC-PUB-7215, December 1996 
+
+      Following properties can be set/changed via jobOptions:
+ 
+      ---Input: TrackCollection
+      - \c TracksName: The name of the StoreGate input container from
+      which the tracks are read. The default is "Tracks", the container
+      from the legacy converters/ambiguity processor.
+      
+      ---Output: VxContainer
+      - \c VxCollectionOutputName: The name of the StoreGate container where the fit
+      results are put. default is "VxCollection".
+
+      ---Vertex Fitter: Billoir FullVertexFitter
+      - \c FitRoutine: The routine which should be used for the fitting. The
+      default is "FullVertexFitter".
+      
+
+      
+      @author  Tatjana Lenz <tatjana.lenz@cern.ch>
+  */  
+  class ZVTOP_Tool : virtual public IVertexFinder, public AthAlgTool
+    {
+    public:
+      ZVTOP_Tool(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_Tool ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+      std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> findVertex(const TrackCollection* trackTES);
+      std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> findVertex(const Trk::TrackParticleBaseCollection* trackTES);
+      std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> findVertex(const xAOD::TrackParticleContainer* trackParticles);
+    private:
+      
+      ToolHandle <Trk::IVertexFitter>			m_iVertexFitter;
+      ToolHandle <InDet::IZVTOP_SpatialPointFinder>	m_iSpatialPointFinder;
+      ToolHandle <InDet::IZVTOP_TrkProbTubeCalc>	m_iTrkProbTubeCalc;
+      ToolHandle <InDet::IZVTOP_SimpleVtxProbCalc>	m_iVtxProbCalc;
+      
+      double						m_resolvingCut;
+      double						m_resolvingStep;
+      double						m_vtxProb_cut;
+      double						m_trk_chi2_cut;
+      ServiceHandle<IBeamCondSvc>                       m_iBeamCondSvc; //!< pointer to the beam condition service
+    };
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkPartBaseVertexState.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkPartBaseVertexState.h
new file mode 100644
index 0000000000000000000000000000000000000000..945764981b6863abd8d61311e0d8f4bce13dda80
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkPartBaseVertexState.h
@@ -0,0 +1,81 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef IINDETZVTOP_TRKPARTBASEVERTEXSTATE_H
+#define IINDETZVTOP_TRKPARTBASEVERTEXSTATE_H 
+
+#include "VxVertex/Vertex.h"
+#include "TrkParticleBase/TrackParticleBase.h"
+
+/**
+ *  ZVTOP_Tool helper class
+ *  to store the full vertex/tracks & vtx_probability
+ * information
+ *  @author:      Tatjana Lenz <Tatjana.Lenz@cern.ch>
+ */
+
+
+namespace InDet
+{
+
+ class ZVTOP_TrkPartBaseVertexState
+ {
+public:
+
+    /**Default Constructor */
+    ZVTOP_TrkPartBaseVertexState();
+    ZVTOP_TrkPartBaseVertexState(double& vtx_prob, Trk::Vertex& vertex, bool& beam_spot, std::vector<const Trk::TrackParticleBase*>& tracks);
+
+   /**virtual destructor of a class*/
+   virtual ~ZVTOP_TrkPartBaseVertexState();
+
+    /**Copy Constructor */
+    ZVTOP_TrkPartBaseVertexState (const ZVTOP_TrkPartBaseVertexState& vs); //copy-constructor
+    //ZVTOP_VertexState &operator=  (const ZVTOP_VertexState &);  //assignement operator
+    
+    double& vtx_prob(void);
+
+    Trk::Vertex& vertex(void);
+
+    bool& beam_spot(void);
+
+    std::vector<const Trk::TrackParticleBase*>&  tracks(void);
+ 
+    inline void set_beam_spot(bool);
+
+protected:
+    //vertex probability
+    double m_vtx_prob;
+    //associated vertex 
+    Trk::Vertex m_vertex;
+    bool m_beam_spot;
+    // associated tracks
+    std::vector<const Trk::TrackParticleBase*> m_tracks; 
+};//end of class definitions
+
+inline double&  ZVTOP_TrkPartBaseVertexState::vtx_prob(void)
+  {
+    return m_vtx_prob;
+  }
+
+inline Trk::Vertex& ZVTOP_TrkPartBaseVertexState::vertex(void)
+  {
+    return m_vertex;
+  }
+
+inline std::vector<const Trk::TrackParticleBase*>& ZVTOP_TrkPartBaseVertexState::tracks(void)
+  {
+    return m_tracks;
+  }
+inline bool& ZVTOP_TrkPartBaseVertexState::beam_spot(void)
+  {
+    return m_beam_spot;
+  }
+inline void ZVTOP_TrkPartBaseVertexState::set_beam_spot(bool beam_spot)
+  {
+    m_beam_spot = beam_spot;
+  }
+}//end of namespace definitions
+
+#endif
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkPartVertexState.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkPartVertexState.h
new file mode 100644
index 0000000000000000000000000000000000000000..077b83a810ccabb6bc4be979be041e3b2590a521
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkPartVertexState.h
@@ -0,0 +1,80 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef IINDETZVTOP_TRKPARTVERTEXSTATE_H
+#define IINDETZVTOP_TRKPARTVERTEXSTATE_H 
+
+#include "VxVertex/Vertex.h"
+#include "Particle/TrackParticle.h"
+/**
+ *  ZVTOP_Tool helper class
+ *  to store the full vertex/tracks & vtx_probability
+ * information
+ *  @author:      Tatjana Lenz <Tatjana.Lenz@cern.ch>
+ */
+
+
+namespace InDet
+{
+
+ class ZVTOP_TrkPartVertexState
+ {
+public:
+
+    /**Default Constructor */
+    ZVTOP_TrkPartVertexState();
+    ZVTOP_TrkPartVertexState(double& vtx_prob, Trk::Vertex& vertex, bool& IP, std::vector<const Rec::TrackParticle*>& tracks);
+
+   /**virtual destructor of a class*/
+   virtual ~ZVTOP_TrkPartVertexState();
+
+    /**Copy Constructor */
+    ZVTOP_TrkPartVertexState (const ZVTOP_TrkPartVertexState& vs); //copy-constructor
+    //ZVTOP_VertexState &operator=  (const ZVTOP_VertexState &);  //assignement operator
+    
+    double& vtx_prob(void);
+
+    Trk::Vertex& vertex(void);
+
+    bool& IP(void);
+
+    std::vector<const Rec::TrackParticle*>&  tracks(void);
+ 
+    inline void set_IP(bool);
+
+protected:
+    //vertex probability
+    double m_vtx_prob;
+    //associated vertex 
+    Trk::Vertex m_vertex;
+    bool m_IP;
+    // associated tracks
+    std::vector<const Rec::TrackParticle*> m_tracks; 
+};//end of class definitions
+
+inline double&  ZVTOP_TrkPartVertexState::vtx_prob(void)
+  {
+    return m_vtx_prob;
+  }
+
+inline Trk::Vertex& ZVTOP_TrkPartVertexState::vertex(void)
+  {
+    return m_vertex;
+  }
+
+inline std::vector<const Rec::TrackParticle*>& ZVTOP_TrkPartVertexState::tracks(void)
+  {
+    return m_tracks;
+  }
+inline bool& ZVTOP_TrkPartVertexState::IP(void)
+  {
+    return m_IP;
+  }
+inline void ZVTOP_TrkPartVertexState::set_IP(bool IP)
+  {
+    m_IP = IP;
+  }
+}//end of namespace definitions
+
+#endif
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkProbTubeCalc.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkProbTubeCalc.h
new file mode 100755
index 0000000000000000000000000000000000000000..c6910ed5ad8235a943dc7a6e91ea93c502d5b2e9
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_TrkProbTubeCalc.h
@@ -0,0 +1,59 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_TrkProbTubeCalc.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef INDETZVTOP_TRKPROBTUBECALC_H
+#define INDETZVTOP_TRKPROBTUBECALC_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h"
+#include "TrkTrack/Track.h"
+#include "VxVertex/RecVertex.h"
+#include"Particle/TrackParticle.h"
+namespace Trk
+{
+  class IExtrapolator;
+}
+namespace InDet 
+{
+
+  /** @class ZVTOP_TrkProbTubeCalc 
+
+      ZVTOP helper class, calculates the Gaussian probability
+      tube for the track trajectory.
+      
+      @author  Tatjana Lenz <tatjana.lenz@cern.ch>
+  */  
+
+  class ZVTOP_TrkProbTubeCalc : virtual public IZVTOP_TrkProbTubeCalc, public AthAlgTool
+    {
+    public:
+      ZVTOP_TrkProbTubeCalc(const std::string&,const std::string&,const IInterface*);
+
+       /** default destructor */
+      virtual ~ZVTOP_TrkProbTubeCalc ();
+      
+       /** standard Athena-Algorithm method */
+      virtual StatusCode initialize();
+       /** standard Athena-Algorithm method */
+      virtual StatusCode finalize  ();
+      
+      virtual double calcProbTube(const Trk::Track& trk, Trk::Vertex& vec);
+      virtual double calcProbTube(const Rec::TrackParticle& trk, Trk::Vertex& vec);
+      virtual double calcProbTube(const Trk::TrackParticleBase& trk, Trk::Vertex& vec);
+      virtual double calcProbTube(const Trk::Perigee* trk, Trk::Vertex& vec);
+      virtual double calcProbTube(const Trk::RecVertex& vtx, Trk::Vertex& vec);
+    private:
+      
+      /** member variables for algorithm properties: */
+     ToolHandle< Trk::IExtrapolator > m_extrapolator;
+      
+    }; 
+} // end of namespace
+
+#endif 
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_VertexState.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_VertexState.h
new file mode 100755
index 0000000000000000000000000000000000000000..96d0d74425c2f7a4b711f3be1cf3de6dbc1f8830
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/InDetZVTOPVxFinder/ZVTOP_VertexState.h
@@ -0,0 +1,80 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef IINDETZVTOP_VERTEXSTATE_H
+#define IINDETZVTOP_VERTEXSTATE_H 
+
+#include "VxVertex/Vertex.h"
+#include "TrkTrack/Track.h"
+/**
+ *  ZVTOP_Tool helper class
+ *  to store the full vertex/tracks & vtx_probability
+ * information
+ *  @author:      Tatjana Lenz <Tatjana.Lenz@cern.ch>
+ */
+
+
+namespace InDet
+{
+
+ class ZVTOP_VertexState
+ {
+public:
+
+    /**Default Constructor */
+    ZVTOP_VertexState();
+    ZVTOP_VertexState(double& vtx_prob, Trk::Vertex& vertex, bool& beam_spot, std::vector<const Trk::Track*>& tracks);
+
+   /**virtual destructor of a class*/
+   virtual ~ZVTOP_VertexState();
+
+    /**Copy Constructor */
+    ZVTOP_VertexState (const ZVTOP_VertexState& vs); //copy-constructor
+    //ZVTOP_VertexState &operator=  (const ZVTOP_VertexState &);  //assignement operator
+    
+    double& vtx_prob(void);
+
+    Trk::Vertex& vertex(void);
+
+    bool& beam_spot(void);
+
+    std::vector<const Trk::Track*>&  tracks(void);
+ 
+    inline void set_beam_spot(bool);
+
+protected:
+    //vertex probability
+    double m_vtx_prob;
+    //associated vertex 
+    Trk::Vertex m_vertex;
+    bool m_beam_spot;
+    // associated tracks
+    std::vector<const Trk::Track*> m_tracks; 
+};//end of class definitions
+
+inline double&  ZVTOP_VertexState::vtx_prob(void)
+  {
+    return m_vtx_prob;
+  }
+
+inline Trk::Vertex& ZVTOP_VertexState::vertex(void)
+  {
+    return m_vertex;
+  }
+
+inline std::vector<const Trk::Track*>& ZVTOP_VertexState::tracks(void)
+  {
+    return m_tracks;
+  }
+inline bool& ZVTOP_VertexState::beam_spot(void)
+  {
+    return m_beam_spot;
+  }
+inline void ZVTOP_VertexState::set_beam_spot(bool beam_spot)
+  {
+    m_beam_spot = beam_spot;
+  }
+}//end of namespace definitions
+
+#endif
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/cmt/requirements b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/cmt/requirements
new file mode 100755
index 0000000000000000000000000000000000000000..3c2277331ca874983c71503d287a45d9b3e478d8
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/cmt/requirements
@@ -0,0 +1,35 @@
+package InDetZVTOPVxFinder
+
+author Tatjana Lenz <tatjana.lenz@cern.ch>
+
+public
+
+use AthenaBaseComps           AthenaBaseComps-*              Control
+use AtlasPolicy		      AtlasPolicy-* 
+use GaudiInterface	      GaudiInterface-*		     External
+use InDetRecToolInterfaces    InDetRecToolInterfaces-*	     InnerDetector/InDetRecTools
+use VxVertex		      VxVertex-*		     Tracking/TrkEvent
+use TrkTrack		      TrkTrack-*		     Tracking/TrkEvent
+use Particle		      Particle-*		     Reconstruction
+use TrkParticleBase	      TrkParticleBase-*        	     Tracking/TrkEvent
+use TrkParameters	      TrkParameters-*                Tracking/TrkEvent
+use xAODTracking              xAODTracking-*                  Event/xAOD
+private
+
+apply_pattern component_library
+library InDetZVTOPVxFinder *.cxx components/*.cxx
+
+private
+use AtlasCLHEP				AtlasCLHEP-*					External
+use DataModel				DataModel-*					Control
+use InDetBeamSpotService		InDetBeamSpotService-*				InnerDetector/InDetConditions
+use TrkExInterfaces			TrkExInterfaces-*				Tracking/TrkExtrapolation
+use TrkSurfaces				TrkSurfaces-*					Tracking/TrkDetDescr
+use TrkVertexFitterInterfaces	        TrkVertexFitterInterfaces-*  	                Tracking/TrkVertexFitter
+use VxSecVertex				VxSecVertex-*					Tracking/TrkEvent
+use EventPrimitives                     EventPrimitives-*                               Event
+#use TrkVxEdmCnv                         TrkVxEdmCnv-*                                   Tracking/TrkVertexFitter
+
+##macro cppdebugflags '$(cppdebugflags_s)'
+##macro_remove componentshr_linkopts "-Wl,-s"
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/doc/mainpage.h b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/doc/mainpage.h
new file mode 100755
index 0000000000000000000000000000000000000000..bd6f4197ccdb31324a0e9b99b6607a7e582fc7ac
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/doc/mainpage.h
@@ -0,0 +1,29 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+@mainpage InDetZVTOPVxFinder Package
+
+Contains:
+- ZVTOP_AmbiguitySolver.cxx (dummy) <br>
+- ZVTOP_SimpleVtxProbCalc.cxx (helper class, calculates vtx probability at a given point) <br>
+- ZVTOP_SpatialPointFinder.cxx (helper class, seed finder for vertex candidates) <br>
+- ZVTOP_Tool.cxx (main tool) <br>
+- ZVTOP_TrkProbTubeCalc.cxx (helper class, calculates track probabilities at a given point) <br>
+- ZVTOP_VertexState.cxx (helper class) <br>
+
+Documentation: David J. Jackson, Nucl.Instrum.Meth.A388: 247-253, 1997
+
+@author Tatjana Lenz <tatjana.lenz@cern.ch>
+
+@section InDetZVTOPVxFinderIntro Introduction
+
+This AlgTool reconstructs a set of vertices and associated subset of the charged tracks, the association accords to a vertex probability function which is based on the trajectories and position resolution of the tracks. The track probability function, \f$ f_i \f$ is defined as a Gaussian tube around a track: the width is given by a track covariance (= error) matrix. The vertex probability function is \f$ \sum f_i - \frac{ \sum f_i^2}{\sum f_i} \f$. The job is to search for local maxima of the vertex probability function. So you have to calculate the vertex probability function at euch space point in detector. To reduce the number of spatial points, all possible combinations of two track pairs are build. For each pair the spatial point is calculated (see ZVTOP_SpatialPointFinder), all spatial points are candidates for a vertex. To reduce the number of vertex candidates the vertex candidate cluster are build. The tracks associated with each cluster are fitted, the vertex fitter can be set/changed via jobOptions, default Billoir FullFitter.
+On the fitted VxCandidate some quality cuts can be applied. The result is stored in a VxContainer. 
+
+@section InDetZVTOPVxFinderReq Requirements
+
+@include requirements
+
+*/
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_AmbiguitySolver.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_AmbiguitySolver.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..66293b0f3e564e9a806652c71e4d01a4828333cb
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_AmbiguitySolver.cxx
@@ -0,0 +1,149 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_AmbiguitySolver.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_AmbiguitySolver.h"
+#include "InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h"
+#include "VxVertex/VxCandidate.h"
+#include "TrkTrack/Track.h"
+#include "TrkParameters/TrackParameters.h"
+#include "VxVertex/VxTrackAtVertex.h"
+
+//================ Constructor =================================================
+
+InDet::ZVTOP_AmbiguitySolver::ZVTOP_AmbiguitySolver(const std::string& t,
+			  const std::string& n,
+			  const IInterface*  p )
+  :
+  AthAlgTool(t,n,p),
+  m_TrkProbTubeCalc("InDet::ZVTOP_TrkProbTubeCalc"),
+  m_VtxProbCalc("InDet::ZVTOP_VtxProbCalc")
+{
+  declareInterface<IZVTOP_AmbiguitySolver>(this);
+
+  //  template for property decalration
+  declareProperty("TrkProbTubeCalcTool",m_TrkProbTubeCalc);
+  declareProperty("VtxProbCalcTool",m_VtxProbCalc);
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_AmbiguitySolver::~ZVTOP_AmbiguitySolver()
+{}
+
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_AmbiguitySolver::initialize()
+{
+  
+  StatusCode sc = AlgTool::initialize();
+
+  if (sc.isFailure()) return sc;
+  //Gaussian Probability Tube for the Track Trajectory
+  if ( m_TrkProbTubeCalc.retrieve().isFailure() ) {
+      msg(MSG::FATAL) << "Failed to retrieve tool " << m_TrkProbTubeCalc<< endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_TrkProbTubeCalc << endreq;
+
+  //Vertex Probability Function
+  if ( m_VtxProbCalc.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_VtxProbCalc<< endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_VtxProbCalc << endreq;
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//std::vector<Trk::VxCandidate*> InDet::ZVTOP_AmbiguitySolver::solveAmbiguities(std::vector<Trk::VxCandidate*> VxContainer) --David S.
+std::vector<xAOD::Vertex*> InDet::ZVTOP_AmbiguitySolver::solveAmbiguities(std::vector<xAOD::Vertex*> VertexContainer)
+{
+  //std::vector<Trk::VxCandidate*> newVxContainer; --David S.
+  std::vector<xAOD::Vertex*> newVertexContainer;
+  //std::vector<Trk::VxCandidate*>::iterator itr = VxContainer.begin(); --David S.
+  std::vector<xAOD::Vertex*>::iterator itr = VertexContainer.begin();
+  //std::map<double, Trk::VxCandidate*> vxprob_map; --David S.
+  std::map<double, xAOD::Vertex*> vxprob_map;
+  double sum_TrkProbTube = 0.;
+  double sum2_TrkProbTube = 0.;
+  //for (; itr != VxContainer.end(); ++itr){ --David S.
+  for (; itr != VertexContainer.end(); ++itr){
+    //Amg::Vector3D pos = (*itr)->recVertex().position(); --David S.
+    Amg::Vector3D pos = (*itr)->position();
+    Trk::Vertex vtx_pos(pos);
+    //std::vector<Trk::VxTrackAtVertex*> trksAtVtx = *((*itr)->vxTrackAtVertex()); --David S.
+    std::vector<Trk::VxTrackAtVertex> trksAtVtx = (*itr)->vxTrackAtVertex();
+    //std::vector<Trk::VxTrackAtVertex*>::iterator trk_itr = trksAtVtx.begin(); --David S.
+    std::vector<Trk::VxTrackAtVertex>::iterator trk_itr = trksAtVtx.begin();
+    for (; trk_itr != trksAtVtx.end(); ++trk_itr ){
+      const Trk::Perigee*	perigee = dynamic_cast<const Trk::Perigee*>((*trk_itr).initialPerigee());
+      double TrkProbTube = m_TrkProbTubeCalc->calcProbTube(perigee, vtx_pos);
+      sum_TrkProbTube  += TrkProbTube;
+      sum2_TrkProbTube += pow(TrkProbTube,2);
+    }
+    double vtx_prob = m_VtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+    //vxprob_map.insert(std::map<double, Trk::VxCandidate*>::value_type(vtx_prob,(*itr))); --David S.
+    vxprob_map.insert(std::map<double, xAOD::Vertex*>::value_type(vtx_prob,(*itr)));
+  }
+  
+  //std::map<double, Trk::VxCandidate*>::reverse_iterator  mapItr = vxprob_map.rbegin(); --David S.
+  std::map<double, xAOD::Vertex*>::reverse_iterator mapItr = vxprob_map.rbegin();
+  //std::map<double, Trk::VxCandidate*> initial_perigee_pos_map; --David S.
+  std::map<double, xAOD::Vertex*> initial_perigee_pos_map;
+  for(; mapItr!=vxprob_map.rend(); ++mapItr){
+    //Trk::VxCandidate* vtx = (*mapItr).second; --David S.
+    xAOD::Vertex* vtx = (*mapItr).second;
+    //std::vector< Trk::VxTrackAtVertex*> trkAtVtx = *(vtx->vxTrackAtVertex()); --David S.
+    std::vector<Trk::VxTrackAtVertex> trkAtVtx = vtx->vxTrackAtVertex();
+    //std::vector< Trk::VxTrackAtVertex*>::iterator trkAtVtx_Iter = trkAtVtx.begin(); --David S.
+    std::vector<Trk::VxTrackAtVertex>::iterator trkAtVtx_Iter = trkAtVtx.begin();
+    for (; trkAtVtx_Iter!= trkAtVtx.end(); ++trkAtVtx_Iter){
+      const Trk::Perigee*	perigee = dynamic_cast<const Trk::Perigee*>((*trkAtVtx_Iter).initialPerigee());
+      Amg::Vector3D pos = perigee->position();
+      //initial_perigee_pos_map.insert(std::map<double, Trk::VxCandidate*>::value_type(pos.mag(),vtx)); --David S.
+      initial_perigee_pos_map.insert(std::map<double, xAOD::Vertex*>::value_type(pos.mag(),vtx));
+    }
+  }
+  mapItr = vxprob_map.rbegin();
+  for (; mapItr!= vxprob_map.rend(); ++mapItr){
+    //Trk::VxCandidate* vtx = (*mapItr).second; --David S.
+    xAOD::Vertex* vtx = (*mapItr).second;
+    //std::vector< Trk::VxTrackAtVertex*> trkAtVtx = *(vtx->vxTrackAtVertex()); --David S.
+    std::vector<Trk::VxTrackAtVertex> trkAtVtx = vtx->vxTrackAtVertex();
+    //std::vector< Trk::VxTrackAtVertex*>::iterator trkAtVtx_Iter = trkAtVtx.begin(); --David S.
+    std::vector<Trk::VxTrackAtVertex>::iterator trkAtVtx_Iter = trkAtVtx.begin();
+    for (; trkAtVtx_Iter!= trkAtVtx.end(); ++trkAtVtx_Iter){
+      const Trk::Perigee*	perigee = dynamic_cast<const Trk::Perigee*>((*trkAtVtx_Iter).initialPerigee());
+      double pos_mag = perigee->position().mag();
+      //std::map<double, Trk::VxCandidate*>::iterator pos_magItr = initial_perigee_pos_map.find(pos_mag); --David S.
+      std::map<double, xAOD::Vertex*>::iterator pos_magItr = initial_perigee_pos_map.find(pos_mag);
+      if (pos_magItr != initial_perigee_pos_map.end()) {
+	if (vtx != (*pos_magItr).second){
+	  trkAtVtx.erase(trkAtVtx_Iter);
+	  if (trkAtVtx_Iter == trkAtVtx.end()) break;
+	}
+      }
+    }
+    if (trkAtVtx.size() > 1) { 
+      //newVxContainer.push_back(vtx); --David S.
+      newVertexContainer.push_back(vtx);
+    }
+  }
+  //return newVxContainer; --David S.
+  return newVertexContainer;
+}
+//================ Finalisation =================================================
+
+StatusCode InDet::ZVTOP_AmbiguitySolver::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+
+//============================================================================================
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SecVtxTool.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SecVtxTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..03807a7f35ac6afe5cbade61b672b163f7ae98d3
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SecVtxTool.cxx
@@ -0,0 +1,689 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_Tool.cxx, (c) ATLAS Detector software
+// begin   : 30-10-2006
+// authors : Tatjana Lenz
+// email   : tatjana.lenz@cern.ch
+// changes :
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_SecVtxTool.h"
+#include "TrkVertexFitterInterfaces/IVertexFitter.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SpatialPointFinder.h"
+#include "InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h"
+#include "InDetZVTOPVxFinder/IZVTOP_AmbiguitySolver.h"
+#include "InDetZVTOPVxFinder/ZVTOP_TrkPartBaseVertexState.h"
+#include "InDetZVTOPVxFinder/ZVTOP_VertexState.h"
+#include "VxVertex/Vertex.h"
+#include "VxVertex/VxTrackAtVertex.h"
+#include "TrkParticleBase/TrackParticleBaseCollection.h"
+#include "TrkTrack/Track.h"
+#include "VxSecVertex/VxSecVertexInfo.h"
+#include <map>
+#include <utility>
+//================ Constructor =================================================
+
+InDet::ZVTOP_SecVtxTool::ZVTOP_SecVtxTool(const std::string& t,
+			  const std::string& n,
+			  const IInterface*  p )
+  :
+  AthAlgTool(t,n,p),
+  m_iVertexFitter("Trk::FullVertexFitter"),
+  m_iSpatialPointFinder("InDet::SpatialPointFinder"),
+  m_iTrkProbTubeCalc("InDet::TrkProbTubeCalc"),
+  m_iVtxProbCalc("InDet::VtxProbCalc"),
+	m_iAmbiguitySolver("InDet::ZVTOP_AmbiguitySolver")
+{
+  //  template for property declaration
+  declareInterface<ISecVertexInJetFinder>(this);
+  declareProperty("VertexFitterTool",m_iVertexFitter);
+  declareProperty("SpatialPointFinderTool",m_iSpatialPointFinder);
+  declareProperty("TrkProbTubeCalcTool",m_iTrkProbTubeCalc);
+  declareProperty("VtxProbCalcTool",m_iVtxProbCalc);
+  declareProperty("AmbiguitySolverTool",m_iAmbiguitySolver);
+  declareProperty("ResolvingCutValue", m_resolvingCut = 0.6);
+  declareProperty("ResolvingStepValue", m_resolvingStep = 0.3);
+  declareProperty("VtxProbCut", m_vtxProb_cut = 0.00001);
+  declareProperty("MaxChi2PerTrack", m_trk_chi2_cut = 5.);
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_SecVtxTool::~ZVTOP_SecVtxTool(){}
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_SecVtxTool::initialize()
+{
+  StatusCode sc = AlgTool::initialize();
+
+  msg (MSG::INFO) << name() << " initialize()" << endreq;
+  if (sc.isFailure()) return sc;
+
+  /* Retrieve Tools*/
+  //SpatialPointFinder
+  if ( m_iSpatialPointFinder.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iSpatialPointFinder << endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iSpatialPointFinder << endreq;
+
+  //Gaussian Probability Tube for the Track Trajectory
+  if ( m_iTrkProbTubeCalc.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iTrkProbTubeCalc<< endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iTrkProbTubeCalc << endreq;
+
+  //Vertex Probability Function
+  if ( m_iVtxProbCalc.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iVtxProbCalc<< endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iVtxProbCalc << endreq;
+
+  //VxFitter
+  if ( m_iVertexFitter.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iVertexFitter << endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iVertexFitter << endreq;
+
+  if ( m_iAmbiguitySolver.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iAmbiguitySolver << endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iAmbiguitySolver << endreq;
+
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================================
+
+StatusCode InDet::ZVTOP_SecVtxTool::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+
+//===========================================================================================
+// TrackParticle
+//===========================================================================================
+const Trk::VxSecVertexInfo* InDet::ZVTOP_SecVtxTool::findSecVertex(const Trk::RecVertex & primaryVertex,
+								const TLorentzVector & jetMomentum,
+								const std::vector<const Trk::TrackParticleBase*> & inputTracks) const
+{
+  //std::vector<Trk::VxCandidate*> secVertices; --David S.
+  std::vector<xAOD::Vertex*> secVertices;
+  //std::vector<Trk::VxCandidate*> new_secVertices(0); --David S.
+  std::vector<xAOD::Vertex*> new_secVertices(0);
+  //if (msgLvl(MSG::DEBUG)) msg () << "jet momentum Et s= "<<jetMomentum.et() <<  endreq; --David S.
+  if (msgLvl(MSG::DEBUG)) msg () << "jet momentum Et s= "<<jetMomentum.Et() <<  endreq;
+  if (inputTracks.size() != 0) {
+    //some variables
+    typedef std::vector<const Trk::TrackParticleBase*>::const_iterator TrkPartVecIter;
+    std::vector<const Trk::TrackParticleBase*> trkColl; // the collection of tracks, which are assigned to one spatial point
+    std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> vtxState_org;// vector of tracks and vertices, tracks can be assigned to more than one vertex, these ambiguities will be resolved later
+    std::multimap<double, InDet::ZVTOP_TrkPartBaseVertexState*> vtxProb_map;
+    //---------------------------------------------------------------------------//  
+    //-------step1: find spatial points & calculate vtx probabilities-------//
+    //--------------------------------------------------------------------------//
+    std::vector< std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> > vertex_objects; //vector of vtx candidate & track collection pairs
+    if (msgLvl(MSG::DEBUG)) msg () << "step ONE search for the spatial points, number of tracks = "<<inputTracks.size() <<  endreq;
+    for (TrkPartVecIter itr_1  = inputTracks.begin(); itr_1 != inputTracks.end()-1; itr_1++)
+  {
+    double sum_TrkProbTube = 0.;
+    double sum2_TrkProbTube = 0.;
+    double vtxProb = 0.;
+    // trk-trk combination
+    for (TrkPartVecIter itr_2 = itr_1+1; itr_2 != inputTracks.end(); itr_2++)
+    {
+      Trk::Vertex* spatialPoint;
+      spatialPoint = m_iSpatialPointFinder->findSpatialPoint((*itr_1),(*itr_2));
+      if (spatialPoint != 0) 
+      {
+	double TrkProbTube_1 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+	double TrkProbTube_2 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_2),*spatialPoint);
+	sum_TrkProbTube = TrkProbTube_1 + TrkProbTube_2;
+	sum2_TrkProbTube = pow(TrkProbTube_1,2.) + pow(TrkProbTube_2,2.);
+	if (sum_TrkProbTube != 0. )
+	  { 
+            vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+            trkColl.push_back((*itr_1));
+            trkColl.push_back((*itr_2));
+            bool IP = false;
+            if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_TrkPartBaseVertexState(vtxProb, *spatialPoint, IP, trkColl ));
+            trkColl.clear();
+	  }
+      }
+      delete spatialPoint;
+    }
+    //trk-IP combination
+    Trk::Vertex* spatialPoint = 0;
+    spatialPoint = m_iSpatialPointFinder->findSpatialPoint(primaryVertex,(*itr_1));
+    if (spatialPoint != 0) 
+    {
+      double IPprobTube = m_iTrkProbTubeCalc->calcProbTube(primaryVertex,*spatialPoint);
+      double TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+      sum_TrkProbTube = IPprobTube + TrkProbTube;
+      sum2_TrkProbTube = pow(TrkProbTube,2.) + pow(IPprobTube,2.);
+      if (sum_TrkProbTube != 0. ) 
+      {
+         vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+         trkColl.push_back((*itr_1));
+         bool IP = true;
+         if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_TrkPartBaseVertexState(vtxProb, *spatialPoint, IP, trkColl ));
+         trkColl.clear();
+      }
+    }
+    delete spatialPoint;
+  }
+  if (msgLvl(MSG::DEBUG)) msg () << " number of spatial points = "<<vtxState_org.size()<<endreq;
+  //reduce the spatial point collection
+  typedef std::vector<InDet::ZVTOP_TrkPartBaseVertexState*>::iterator Sp_Iter;
+  std::vector< InDet::ZVTOP_TrkPartBaseVertexState*> vtxState;
+  for (Sp_Iter itr1 = vtxState_org.begin(); itr1 != vtxState_org.end(); itr1++)
+  {
+
+    if (vtxState.size() == 0) vtxState.push_back(*itr1);
+    else 
+      {
+	Trk::Vertex vtx_new = (*itr1)->vertex();
+	bool can_be_resolved = false;
+	for (Sp_Iter itr2 = vtxState.begin(); itr2 != vtxState.end(); itr2++)
+	  {
+	    Trk::Vertex vtx_in = (*itr2)->vertex();
+            double r_diff = sqrt(pow(vtx_new.position().x(),2.) + pow(vtx_new.position().y(),2.)) - sqrt(pow(vtx_in.position().x(),2.) + pow(vtx_in.position().y(),2.));
+            double z_diff = vtx_new.position().z() - vtx_in.position().z();
+            if (fabs(r_diff) > 0.001 && fabs(z_diff) > 0.001 && r_diff != 0. && z_diff != 0.) can_be_resolved = true;
+            else 
+	      {
+                for (unsigned int trk_itr = 0; trk_itr < (*itr1)->tracks().size(); trk_itr++) (*itr2)->tracks().push_back((*itr1)->tracks()[trk_itr]); 
+                can_be_resolved = false;
+                break;
+	      }
+	  }
+	if (can_be_resolved)  vtxState.push_back(*itr1);
+      }
+  }
+  for (Sp_Iter itr_sort = vtxState.begin();  itr_sort!= vtxState.end(); itr_sort++)  vtxProb_map.insert(std::multimap<double, InDet::ZVTOP_TrkPartBaseVertexState*>::value_type((*itr_sort)->vtx_prob(), *itr_sort));
+
+
+//-------------------------------------------------------------------------------//
+//----step2: vertex clustering, aim: to cluster all vertices together, ----//
+//----which can not be resolved. The found vertex seeds and ---------//
+//----associated tracks are taken as an input for the vertex fit. --------//
+//------------------------------------------------------------------------------//
+
+  if (vtxState.size() != 0 ){
+    if (msgLvl(MSG::DEBUG)) msg () << " step TWO vertex clustering, number of reduced vertices = "<<vtxState.size()<< endreq;
+    //sort the vtxState collection, starting with a highest vtx probability
+    std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> vtxState_sorted;
+    std::multimap<double, InDet::ZVTOP_TrkPartBaseVertexState*>::reverse_iterator s= vtxProb_map.rbegin();
+    for (;  s!=vtxProb_map.rend();  s++)  {
+      InDet::ZVTOP_TrkPartBaseVertexState* vtx_state = (*s).second;
+      vtxState_sorted.push_back(vtx_state);
+    }
+    //store first element in the vertex_objects ---->vector<vector<VertexState>>
+    std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> vtx_coll;
+    vtx_coll.push_back(*(vtxState_sorted.begin()));
+    vertex_objects.push_back(vtx_coll);//store first seed
+    std::vector< InDet::ZVTOP_TrkPartBaseVertexState*>  vtxState_new(vtxState_sorted); //copy of the vtxState_sorted
+    vtxState_new.erase(vtxState_new.begin()); // without first element
+    //loop over vtxState_sorted
+    for (Sp_Iter org_itr= vtxState_sorted.begin(); org_itr != vtxState_sorted.end(); org_itr++)
+    {
+      Trk::Vertex seed = (*org_itr)->vertex();
+      //found where the seed is already stored
+      bool vtx_is_stored = false;
+      unsigned int vertex_objectsPosition = 0; 
+      for (unsigned int vertex_objectsItr = 0; vertex_objectsItr!= vertex_objects.size(); vertex_objectsItr++)
+        {
+          //stored vertices in the line vertex_objectsItr
+          std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> stored_vertices = vertex_objects[vertex_objectsItr];
+           //---->inner loop
+          for (unsigned int stored_vtx_itr = 0; stored_vtx_itr!= stored_vertices.size(); stored_vtx_itr++)
+          {
+	    Trk::Vertex storedVtx = stored_vertices[stored_vtx_itr]->vertex();
+	    double diff = ((*org_itr)->vertex().position().x() - storedVtx.position().x()) + ((*org_itr)->vertex().position().y() - storedVtx.position().y()) + ((*org_itr)->vertex().position().z() - storedVtx.position().z());
+	    if (fabs(diff) < 0.0001)
+              {
+		vertex_objectsPosition = vertex_objectsItr;
+		vtx_is_stored = true;
+		break; // break inner loop if found
+              }
+	  }
+	  if (vtx_is_stored == true) break; // break outer loop if found
+	}
+      if (!vtx_is_stored) {
+	//if not stored
+	std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> new_vtx_coll;
+	new_vtx_coll.push_back(*org_itr);
+	vertex_objects.push_back(new_vtx_coll);
+	vertex_objectsPosition = vertex_objects.size() - 1;
+      }
+      for (Sp_Iter new_itr = vtxState_new.begin(); new_itr!= vtxState_new.end(); new_itr++)
+	{
+	  Trk::Vertex cand = (*new_itr)->vertex();
+	  //calculate the vtx prob function for this seed
+	  std::multimap<double, Trk::Vertex > vtx_prob_coll;
+	  vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*new_itr)->vtx_prob(),(*new_itr)->vertex()));
+	  vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*org_itr)->vtx_prob(),(*org_itr)->vertex()));
+	  for (double alpha = m_resolvingStep; alpha <1.; alpha += m_resolvingStep)
+	    {
+	      Trk::Vertex SP_line = Trk::Vertex((*org_itr)->vertex().position() + alpha*((*new_itr)->vertex().position() - (*org_itr)->vertex().position()));
+	      double sum_TrkProbTube = 0.;
+	      double sum2_TrkProbTube = 0.;
+	      for (TrkPartVecIter trk_itr = inputTracks.begin(); trk_itr != inputTracks.end(); trk_itr++)
+		{
+		  double TrkProbTube = 0.;
+		  TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*trk_itr),SP_line);
+		  sum_TrkProbTube += TrkProbTube;
+		  sum2_TrkProbTube += pow(TrkProbTube,2.);
+		}
+	      double IPprobTube = m_iTrkProbTubeCalc->calcProbTube(primaryVertex,SP_line);
+	      sum_TrkProbTube += IPprobTube;
+	      sum2_TrkProbTube += pow(IPprobTube,2.);
+	      double vtx_prob = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+	      vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type(vtx_prob, SP_line));
+	    }
+	  double vtx_prob_ratio = (*vtx_prob_coll.begin()).first/(*new_itr)->vtx_prob();
+	  if (vtx_prob_ratio >= m_resolvingCut)  {
+	    //check if the vertices can be resolved
+	    vertex_objects[vertex_objectsPosition].push_back(*new_itr);
+	    vtxState_new.erase(new_itr);
+	    if (new_itr != vtxState_new.end()) --new_itr;
+	    if (new_itr == vtxState_new.end()) break;
+	  }
+	}
+    }//loop over orig reduced coll
+  }//if spatial point collection size > 0
+
+  
+   //--------------------------------------------------------------------------------//
+  //--------------------step3: call vertex fitter----------------------------------//
+  //--------------------------------------------------------------------------------//
+
+  if (msgLvl(MSG::DEBUG)) msg () << " step THREE, vertex fit, vertex_objects size = "<<vertex_objects.size()<< endreq;
+  //std::vector<Trk::VxCandidate* > theVxContainer; --David S.
+  std::vector<xAOD::Vertex*> theVertexContainer;
+  //prepare trk coll --> remove double tracks
+  std::vector<std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> >::iterator vtx_itr = vertex_objects.begin();
+  for (; vtx_itr!= vertex_objects.end(); vtx_itr++)
+    {
+      bool with_IP = false;
+      std::vector <const Trk::TrackParticleBase*> trk_coll(0); //--->trk_coll for the fit
+      std::vector<InDet::ZVTOP_TrkPartBaseVertexState*>::iterator itr = (*vtx_itr).begin();
+      trk_coll.push_back((*itr)->tracks()[0]);
+      for (; itr != (*vtx_itr).end(); itr++)
+	{ 
+	  for ( std::vector <const Trk::TrackParticleBase*>::iterator vs_itr = ((*itr)->tracks()).begin(); vs_itr!= ((*itr)->tracks()).end(); vs_itr++)
+	    {
+	      if (msgLvl(MSG::DEBUG)) msg () <<"old trk coll size = "<<trk_coll.size()<<", new track = "<<(*vs_itr)<<endreq;
+	      bool found = false;
+	      for (std::vector <const Trk::TrackParticleBase*>::iterator trk_itr = trk_coll.begin();  trk_itr!= trk_coll.end(); trk_itr++)  {
+		if (*vs_itr == *trk_itr)  found = true;
+	      }
+	      if (!found)  trk_coll.push_back(*vs_itr);
+	    }
+	  if ((*itr)->beam_spot()) with_IP = true;
+	  if (msgLvl(MSG::DEBUG)) msg () <<"new track collection size = "<<trk_coll.size()<<endreq;
+	}
+      //call the fitter
+      //Trk::VxCandidate * myVxCandidate(0); --David S.
+      xAOD::Vertex * myxAODVertex(0);
+      //const Amg::Vector3D p(0.,0.,0.); --David S.
+      const Amg::Vector3D startingPoint(0.,0.,0.);
+      //Trk::Vertex startingPoint(p); --David S.
+      //if (!with_IP) myVxCandidate = m_iVertexFitter->fit(trk_coll,startingPoint); --David S.
+      if (!with_IP) myxAODVertex = m_iVertexFitter->fit(trk_coll,startingPoint);
+      bool bad_chi2 = true;
+      //if (myVxCandidate) --David S.
+      if (myxAODVertex)
+	{
+	  while (bad_chi2)
+	    {
+	      //if (msgLvl(MSG::DEBUG)) msg () <<"rec Vertex = "<<myVxCandidate->recVertex().position()<<endreq; --David S.
+	      if (msgLvl(MSG::DEBUG)) msg () << "Vertex pos = " << myxAODVertex->position() << endreq;
+	      //std::vector< Trk::VxTrackAtVertex*> trkAtVtx = *(myVxCandidate->vxTrackAtVertex()); --David S.
+	      std::vector<Trk::VxTrackAtVertex> trkAtVtx = myxAODVertex->vxTrackAtVertex();
+	      //std::vector< Trk::VxTrackAtVertex*>::iterator trkAtVtx_Iter = trkAtVtx.begin(); --David S.
+	      std::vector<Trk::VxTrackAtVertex>::iterator trkAtVtx_Iter = trkAtVtx.begin();
+	      std::vector< const Trk::TrackParticleBase*>::iterator trk_Iter = trk_coll.begin();
+	      double largest_chi2 = 0.;
+	      std::vector< const Trk::TrackParticleBase*>::iterator index;
+	      for (; trkAtVtx_Iter!= trkAtVtx.end(); trkAtVtx_Iter++)
+		{
+		  double chi2 = (*trkAtVtx_Iter).trackQuality().chiSquared();
+		  if (chi2 > largest_chi2) {
+		    largest_chi2 = chi2;
+		    index = trk_Iter;
+		  }
+		  trk_Iter++;
+		}
+	      if (largest_chi2 > m_trk_chi2_cut)
+		{
+		  if (trk_coll.size() < 3) break;
+		  else {
+		    trk_coll.erase(index);
+		    if (trk_coll.size() >= 2) {
+		      //if (myVxCandidate!=0) { delete myVxCandidate; myVxCandidate=0; } --David S.
+		      if (myxAODVertex!=0) { delete myxAODVertex; myxAODVertex=0; }
+		      //myVxCandidate = m_iVertexFitter->fit(trk_coll, startingPoint); --David S.
+		      myxAODVertex = m_iVertexFitter->fit(trk_coll, startingPoint);
+		    }
+		    //if (myVxCandidate == 0) break; --David S.
+		    if (myxAODVertex == 0) break;
+		  }
+		} else bad_chi2 = false;
+	    }
+	}
+      //if (myVxCandidate && bad_chi2 == false) secVertices.push_back(myVxCandidate); --David S.
+      if (myxAODVertex && bad_chi2 == false) secVertices.push_back(myxAODVertex);
+    }
+  new_secVertices = m_iAmbiguitySolver->solveAmbiguities(secVertices);
+  if (msgLvl(MSG::DEBUG)) msg () <<"vertex container size = "<<secVertices.size()<<endreq;
+  for (Sp_Iter iter = vtxState_org.begin(); iter != vtxState_org.end(); iter++) delete *iter;
+  
+  } else {
+    if (msgLvl(MSG::DEBUG)) msg () <<"No tracks in this event, provide to the next event" <<  endreq;
+  }
+  
+  return new Trk::VxSecVertexInfo(secVertices);
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////
+//for xAOD::IParticle --David S.
+//Added purely to satisfy new inheritance in ISecVertexInJetFinder, not yet implemented
+/////////////////////////////////////////////////////////////////////////////////////////
+const Trk::VxSecVertexInfo* InDet::ZVTOP_SecVtxTool::findSecVertex(const xAOD::Vertex & primaryVertex, const TLorentzVector & jetMomentum,const std::vector<const xAOD::IParticle*> & inputTracks) const {
+
+  if(msgLvl(MSG::DEBUG)){
+    msg(MSG::DEBUG) << "No ZVTOP_SecVtxTool implementation for xAOD::IParticle" << endreq;
+    msg(MSG::DEBUG) << "Returning null VxSecVertexInfo*" << endreq;
+  }
+
+  Trk::VxSecVertexInfo* returnInfo(0);
+  return returnInfo;
+
+}
+
+////////////////////////////////////////////////////////////
+//for Trk::Track
+///////////////////////////////////////////////////////////
+const Trk::VxSecVertexInfo* InDet::ZVTOP_SecVtxTool::findSecVertex(const Trk::RecVertex & primaryVertex,const TLorentzVector & jetMomentum, const std::vector<const Trk::Track*> & inputTracks) const
+{
+  //std::vector<Trk::VxCandidate*> secVertices; --David S.
+  std::vector<xAOD::Vertex*> secVertices;
+  //if (msgLvl(MSG::DEBUG)) msg () << "jet momentum Et = "<<jetMomentum.et() <<  endreq; --David S.
+  if (msgLvl(MSG::DEBUG)) msg () << "jet momentum Et = "<<jetMomentum.Et() <<  endreq;
+  if (inputTracks.size() != 0) {
+    //some variables
+    typedef std::vector<const Trk::Track*>::const_iterator TrackDataVecIter;
+    std::vector<const Trk::Track*> trkColl(0); // the collection of tracks, which are assigned to one spatial point
+    std::vector<InDet::ZVTOP_VertexState*> vtxState_org(0);// vector of tracks and vertices, tracks can be assigned to more than one vertex, these ambiguities will be resolved later
+    std::multimap<double, InDet::ZVTOP_VertexState*> vtxProb_map;
+    //---------------------------------------------------------------------------//  
+    //-------step1: find spatial points & calculate vtx probabilities-------//
+    //--------------------------------------------------------------------------//
+
+    std::vector< std::vector<InDet::ZVTOP_VertexState*> > vertex_objects(0); //vector of vtx candidate & track collection pairs
+    
+    if (msgLvl(MSG::DEBUG)) msg () << "step ONE search for the spatial points, number of tracks = "<<inputTracks.size() <<  endreq;
+    for (TrackDataVecIter itr_1  = inputTracks.begin(); itr_1 != inputTracks.end()-1; itr_1++)
+      {
+	double sum_TrkProbTube = 0.;
+	double sum2_TrkProbTube = 0.;
+	double vtxProb = 0.;
+	// trk-trk combination
+	for (TrackDataVecIter itr_2 = itr_1+1; itr_2 != inputTracks.end(); itr_2++)
+	  {
+	    Trk::Vertex* spatialPoint;
+	    spatialPoint = m_iSpatialPointFinder->findSpatialPoint((*itr_1),(*itr_2));
+	    if (spatialPoint != 0) 
+	      {
+		double TrkProbTube_1 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+		double TrkProbTube_2 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_2),*spatialPoint);
+		sum_TrkProbTube = TrkProbTube_1 + TrkProbTube_2;
+		sum2_TrkProbTube = pow(TrkProbTube_1,2.) + pow(TrkProbTube_2,2.);
+		if (sum_TrkProbTube != 0. )
+		  { 
+		    vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+		    trkColl.push_back((*itr_1));
+		    trkColl.push_back((*itr_2));
+		    bool IP = false;
+		    if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_VertexState(vtxProb, *spatialPoint, IP, trkColl ));
+		    trkColl.clear();
+		  }
+	      }
+	    delete spatialPoint;
+	  }
+	//trk-IP combination
+	Trk::Vertex* spatialPoint = 0;
+	spatialPoint = m_iSpatialPointFinder->findSpatialPoint(primaryVertex,(*itr_1));
+	if (spatialPoint != 0) 
+	  {
+	    double BeamProbTube = m_iTrkProbTubeCalc->calcProbTube(primaryVertex,*spatialPoint);
+	    double TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+	    sum_TrkProbTube = BeamProbTube + TrkProbTube;
+	    sum2_TrkProbTube = pow(TrkProbTube,2.) + pow(BeamProbTube,2.);
+	    if (sum_TrkProbTube != 0. ) 
+	      {
+		vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+		trkColl.push_back((*itr_1));
+		bool IP = true;
+		if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_VertexState(vtxProb, *spatialPoint, IP, trkColl ));
+		trkColl.clear();
+	      }
+	  }
+	delete spatialPoint;
+      }
+    if (msgLvl(MSG::DEBUG)) msg () << " number of spatial points = "<<vtxState_org.size()<<endreq;
+    //reduce the spatial point collection
+    typedef std::vector<InDet::ZVTOP_VertexState*>::iterator Sp_Iter;
+    std::vector< InDet::ZVTOP_VertexState*> vtxState;
+    for (Sp_Iter itr1 = vtxState_org.begin(); itr1 != vtxState_org.end(); itr1++)
+      {
+	
+	if (vtxState.size() == 0) vtxState.push_back(*itr1);
+	else {
+	  Trk::Vertex vtx_new = (*itr1)->vertex();
+	  bool can_be_resolved = false;
+	  for (Sp_Iter itr2 = vtxState.begin(); itr2 != vtxState.end(); itr2++)
+	    {
+	      Trk::Vertex vtx_in = (*itr2)->vertex();
+	      double r_diff = sqrt(pow(vtx_new.position().x(),2.) + pow(vtx_new.position().y(),2.)) - sqrt(pow(vtx_in.position().x(),2.) + pow(vtx_in.position().y(),2.));
+	      double z_diff = vtx_new.position().z() - vtx_in.position().z();
+	      if (fabs(r_diff) > 0.001 && fabs(z_diff) > 0.001 && r_diff != 0. && z_diff != 0.) can_be_resolved = true;
+	      else 
+		{
+		  for (unsigned int trk_itr = 0; trk_itr < (*itr1)->tracks().size(); trk_itr++) (*itr2)->tracks().push_back((*itr1)->tracks()[trk_itr]); 
+		  can_be_resolved = false;
+		  break;
+		}
+	    }
+	  if (can_be_resolved)  vtxState.push_back(*itr1);
+	}
+      }
+    for (Sp_Iter itr_sort = vtxState.begin();  itr_sort!= vtxState.end(); itr_sort++)  vtxProb_map.insert(std::multimap<double, InDet::ZVTOP_VertexState*>::value_type((*itr_sort)->vtx_prob(), *itr_sort));
+
+
+    //-------------------------------------------------------------------------------//
+    //----step2: vertex clustering, aim: to cluster all vertices together, ----//
+    //----which can not be resolved. The found vertex seeds and ---------//
+    //----associated tracks are taken as an input for the vertex fit. --------//
+    //------------------------------------------------------------------------------//
+
+
+    if (vtxState.size() != 0 ){
+      if (msgLvl(MSG::DEBUG)) msg () << " step TWO vertex clustering, number of reduced vertices = "<<vtxState.size()<< endreq;
+      //sort the vtxState collection, starting with a highest vtx probability
+      std::vector<InDet::ZVTOP_VertexState*> vtxState_sorted;
+      unsigned int IP_counter = 0;
+      std::multimap<double, InDet::ZVTOP_VertexState*>::reverse_iterator s= vtxProb_map.rbegin();
+      for (;  s!=vtxProb_map.rend();  s++)  {
+	InDet::ZVTOP_VertexState* vtx_state = (*s).second;
+	if (vtx_state->beam_spot()) IP_counter += 1;
+	if (IP_counter > 1) vtx_state->set_beam_spot(false);
+	vtxState_sorted.push_back(vtx_state);
+      }
+      //store first element in the vertex_objects ---->vector<vector<VertexState>>
+      std::vector<InDet::ZVTOP_VertexState*> vtx_coll;
+      vtx_coll.push_back(*(vtxState_sorted.begin()));
+      vertex_objects.push_back(vtx_coll);//store first seed
+      std::vector< InDet::ZVTOP_VertexState*>  vtxState_new(vtxState_sorted); //copy of the vtxState_sorted
+      vtxState_new.erase(vtxState_new.begin()); // without first element
+      //loop over vtxState_sorted
+      for (Sp_Iter org_itr= vtxState_sorted.begin(); org_itr != vtxState_sorted.end(); org_itr++)
+	{
+	  Trk::Vertex seed = (*org_itr)->vertex();
+	  //found where the seed is already stored
+	  bool vtx_is_stored = false;
+	  unsigned int vertex_objectsPosition = 0; 
+	  for (unsigned int vertex_objectsItr = 0; vertex_objectsItr!= vertex_objects.size(); vertex_objectsItr++)
+	    {
+	      //stored vertices in the line vertex_objectsItr
+	      std::vector<InDet::ZVTOP_VertexState*> stored_vertices = vertex_objects[vertex_objectsItr];
+	      //---->inner loop
+	      for (unsigned int stored_vtx_itr = 0; stored_vtx_itr!= stored_vertices.size(); stored_vtx_itr++)
+		{
+		  Trk::Vertex storedVtx = stored_vertices[stored_vtx_itr]->vertex();
+		  double diff = ((*org_itr)->vertex().position().x() - storedVtx.position().x()) + ((*org_itr)->vertex().position().y() - storedVtx.position().y()) + ((*org_itr)->vertex().position().z() - storedVtx.position().z());
+		  if (fabs(diff) < 0.0001)
+		    {
+		      vertex_objectsPosition = vertex_objectsItr;
+		      vtx_is_stored = true;
+		      break; // break inner loop if found
+		    }
+		}
+	      if (vtx_is_stored == true) break; // break outer loop if found
+	    }
+	  if (!vtx_is_stored) {
+	    //if not stored
+	    std::vector<InDet::ZVTOP_VertexState*> new_vtx_coll;
+	    new_vtx_coll.push_back(*org_itr);
+	    vertex_objects.push_back(new_vtx_coll);
+	    vertex_objectsPosition = vertex_objects.size() - 1;
+	  }
+	  for (Sp_Iter new_itr = vtxState_new.begin(); new_itr!= vtxState_new.end(); new_itr++)
+	    {
+	      Trk::Vertex cand = (*new_itr)->vertex();
+	      //calculate the vtx prob function for this seed
+	      std::multimap<double, Trk::Vertex > vtx_prob_coll;
+	      vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*new_itr)->vtx_prob(),(*new_itr)->vertex()));
+	      vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*org_itr)->vtx_prob(),(*org_itr)->vertex()));
+	      for (double alpha = m_resolvingStep; alpha <1.; alpha += m_resolvingStep)
+		{
+		  Trk::Vertex SP_line = Trk::Vertex((*org_itr)->vertex().position() + alpha*((*new_itr)->vertex().position() - (*org_itr)->vertex().position()));
+                  double sum_TrkProbTube = 0.;
+                  double sum2_TrkProbTube = 0.;
+                  for (TrackDataVecIter trk_itr = inputTracks.begin(); trk_itr != inputTracks.end(); trk_itr++)
+		    {
+	              double TrkProbTube = 0.;
+                      TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*trk_itr),SP_line);
+		      sum_TrkProbTube += TrkProbTube;
+		      sum2_TrkProbTube += pow(TrkProbTube,2.);
+		    }
+		  double BeamProbTube = m_iTrkProbTubeCalc->calcProbTube(primaryVertex,SP_line);
+		  sum_TrkProbTube += BeamProbTube;
+		  sum2_TrkProbTube += pow(BeamProbTube,2.);
+		  double vtx_prob = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+		  vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type(vtx_prob, SP_line));
+		}
+              double vtx_prob_ratio = (*vtx_prob_coll.begin()).first/(*new_itr)->vtx_prob();
+	      if (vtx_prob_ratio >= m_resolvingCut)  {
+		//check if the vertices can be resolved
+		vertex_objects[vertex_objectsPosition].push_back(*new_itr);
+		vtxState_new.erase(new_itr);
+		if (new_itr != vtxState_new.end()) --new_itr;
+		if (new_itr == vtxState_new.end()) break;
+	      }
+	    }
+	}//loop over orig reduced coll
+    }//if spatial point collection size > 0
+
+    
+    //--------------------------------------------------------------------------------//
+    //--------------------step3: call vertex fitter----------------------------------//
+    //--------------------------------------------------------------------------------//
+    
+    if (msgLvl(MSG::DEBUG)) msg () << " step THREE, vertex fit, vertex_objects size = "<<vertex_objects.size()<< endreq;
+    //std::vector<Trk::VxCandidate* > theVxContainer; --David S.
+    std::vector<xAOD::Vertex*> theVertexContainer;
+    //prepare trk coll --> remove double tracks
+    std::vector<std::vector<InDet::ZVTOP_VertexState*> >::iterator vtx_itr = vertex_objects.begin();
+    for (; vtx_itr!= vertex_objects.end(); vtx_itr++)
+      {
+	bool with_IP = false;
+	std::vector <const Trk::Track*> trk_coll(0); //--->trk_coll for the fit
+	std::vector<InDet::ZVTOP_VertexState*>::iterator itr = (*vtx_itr).begin();
+	trk_coll.push_back((*itr)->tracks()[0]);
+	for (; itr != (*vtx_itr).end(); itr++)
+	  { 
+	    for ( std::vector <const Trk::Track*>::iterator vs_itr = ((*itr)->tracks()).begin(); vs_itr!= ((*itr)->tracks()).end(); vs_itr++)
+	      {
+		bool found = false;
+		for (std::vector <const Trk::Track*>::iterator trk_itr = trk_coll.begin();  trk_itr!= trk_coll.end(); trk_itr++)  {
+		  if (*vs_itr == *trk_itr)  found = true;
+		}
+		if (!found)  trk_coll.push_back(*vs_itr);
+	      }
+	    if ((*itr)->beam_spot()) with_IP = true;
+	  }
+	//call the fitter
+	//Trk::VxCandidate * myVxCandidate(0); --David S.
+	xAOD::Vertex * myxAODVertex(0);
+	//if (!with_IP) myVxCandidate = m_iVertexFitter->fit(trk_coll); --David S.
+	if (!with_IP) myxAODVertex = m_iVertexFitter->fit(trk_coll);
+	bool bad_chi2 = true;
+	//if (myVxCandidate) { --David S.
+	if (myxAODVertex) {
+          while (bad_chi2)
+	    {
+	      //std::vector< Trk::VxTrackAtVertex*> trkAtVtx = *(myVxCandidate->vxTrackAtVertex()); --David S.
+	      std::vector<Trk::VxTrackAtVertex> trkAtVtx = myxAODVertex->vxTrackAtVertex();
+	      //std::vector< Trk::VxTrackAtVertex*>::iterator trkAtVtx_Iter = trkAtVtx.begin(); --David S.
+	      std::vector<Trk::VxTrackAtVertex>::iterator trkAtVtx_Iter = trkAtVtx.begin();
+	      std::vector< const Trk::Track*>::iterator trk_Iter = trk_coll.begin();
+	      double largest_chi2 = 0.;
+	      std::vector< const Trk::Track*>::iterator index;
+	      for (; trkAtVtx_Iter!= trkAtVtx.end(); trkAtVtx_Iter++)
+		{
+		  double chi2 = (*trkAtVtx_Iter).trackQuality().chiSquared();
+		  if (chi2 > largest_chi2) {
+		    largest_chi2 = chi2;
+		    index = trk_Iter;
+		  }
+		  trk_Iter++;
+		}
+	      if (largest_chi2 > m_trk_chi2_cut)
+		{
+		  if (trk_coll.size() < 3) break;
+		  else {
+		    trk_coll.erase(index);
+		    if (trk_coll.size() >= 2) {
+		      //if (myVxCandidate!=0) { delete myVxCandidate; myVxCandidate=0; } --David S.
+		      if (myxAODVertex!=0) { delete myxAODVertex; myxAODVertex=0; }
+		      //myVxCandidate = m_iVertexFitter->fit(trk_coll); --David S.
+		      myxAODVertex = m_iVertexFitter->fit(trk_coll);
+		    }
+		    //if (myVxCandidate == 0) break; --David S.
+		    if (myxAODVertex == 0) break;
+		  }
+		} else bad_chi2 = false;
+	    }
+	}
+	//if (myVxCandidate && bad_chi2 == false) secVertices.push_back(myVxCandidate); --David S.
+	if (myxAODVertex && bad_chi2 == false) secVertices.push_back(myxAODVertex);
+      }
+    if (msgLvl(MSG::DEBUG)) msg () <<"vertex container size = "<<secVertices.size()<<endreq;
+    for (Sp_Iter iter = vtxState_org.begin(); iter != vtxState_org.end(); iter++) delete *iter;
+  } else {
+    if (msgLvl(MSG::DEBUG)) msg () <<"No tracks in this event, provide to the next event" <<  endreq;
+  }
+  
+  return new Trk::VxSecVertexInfo(secVertices);
+  
+
+}
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SimpleVtxProbCalc.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SimpleVtxProbCalc.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..1d34c95e7814c848eabd3915e7cc6484b4fd1bd8
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SimpleVtxProbCalc.cxx
@@ -0,0 +1,59 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_SimpleVtxProbCalc.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_SimpleVtxProbCalc.h"
+
+//================ Constructor =================================================
+
+InDet::ZVTOP_SimpleVtxProbCalc::ZVTOP_SimpleVtxProbCalc(const std::string& t,
+			  const std::string& n,
+			  const IInterface*  p )
+  :
+  AthAlgTool(t,n,p)
+{
+  declareInterface<IZVTOP_SimpleVtxProbCalc>(this);
+
+  //  template for property decalration
+  //declareProperty("PropertyName", m_propertyName);
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_SimpleVtxProbCalc::~ZVTOP_SimpleVtxProbCalc()
+{}
+
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_SimpleVtxProbCalc::initialize()
+{
+  
+  StatusCode sc = AlgTool::initialize();
+
+  if (sc.isFailure()) return sc;
+
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================
+
+StatusCode InDet::ZVTOP_SimpleVtxProbCalc::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+
+//============================================================================================
+double InDet::ZVTOP_SimpleVtxProbCalc::calcVtxProb(double & trk_func_sum, double & trk_func_sum2)
+{
+  double vtx_prob = 0.;
+  if (trk_func_sum != 0.) vtx_prob = trk_func_sum - trk_func_sum2/trk_func_sum;
+  return vtx_prob;
+}
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SlowSpatialPointFinder.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SlowSpatialPointFinder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..be39abd1ff88354830a04339a03f302107e64780
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SlowSpatialPointFinder.cxx
@@ -0,0 +1,256 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_SpatialPointFinder.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_SlowSpatialPointFinder.h"
+#include "TrkParameters/TrackParameters.h"
+#include "VxVertex/Vertex.h"
+#include "Particle/TrackParticle.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include "VxVertex/LinearizedTrack.h"
+#include "TrkVertexFitterInterfaces/IVertexLinearizedTrackFactory.h"
+
+//================ Constructor =================================================
+
+InDet::ZVTOP_SlowSpatialPointFinder::ZVTOP_SlowSpatialPointFinder(const std::string& t, const std::string& n, const IInterface*  p )
+  :
+  AthAlgTool(t,n,p)
+{
+  declareInterface<IZVTOP_SpatialPointFinder>(this);
+  //  template for property declaration
+  declareProperty("Chi2_cut_value", m_chi2);
+  declareProperty ( "LinearizedTrackFactory", m_linFactory );
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_SlowSpatialPointFinder::~ZVTOP_SlowSpatialPointFinder()
+{}
+
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_SlowSpatialPointFinder::initialize()
+{
+  
+  StatusCode sc = AlgTool::initialize();
+
+  if (sc.isFailure()) return sc;
+  // get tool svc
+  if(m_linFactory.retrieve().isFailure())
+    {
+      msg (MSG::ERROR) <<"Could not find ToolSvc."<<endreq;
+      return sc;
+    } else msg (MSG::INFO) << "Retrieved tool " << m_linFactory << endreq;
+  
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================
+
+StatusCode InDet::ZVTOP_SlowSpatialPointFinder::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::Track* trk_1, const Trk::Track* trk_2)
+{
+  const Trk::TrackParameters* perigee_1(dynamic_cast<const Trk::TrackParameters*>(trk_1->perigeeParameters()));
+  const Trk::TrackParameters* perigee_2(dynamic_cast<const Trk::TrackParameters*>(trk_2->perigeeParameters()));
+  if (!perigee_1 | !perigee_2) {
+    if (msgLvl(MSG::VERBOSE)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+		Trk::Vertex* vertex = findSpatialPoint(perigee_1, perigee_2);
+		return vertex;
+  }//if measured perigee
+}
+
+//============================================================================================
+
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Trk::Track* trk_1)
+{
+  const Trk::TrackParameters *perigee_1(dynamic_cast<const Trk::TrackParameters*>(trk_1->perigeeParameters()));
+  if (!perigee_1) {
+    if (msgLvl(MSG::VERBOSE)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    //we need Trk::Vertex
+    Trk::Vertex* vertex = findSpatialPoint(vtx, perigee_1);
+    return vertex;
+  }
+}
+
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Rec::TrackParticle* trk_1, const Rec::TrackParticle* trk_2)
+{
+  const Trk::TrackParameters* perigee_1(trk_1->measuredPerigee());
+  const Trk::TrackParameters* perigee_2(trk_2->measuredPerigee());
+  if (!perigee_1 | !perigee_2) {
+    if (msgLvl(MSG::VERBOSE)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    //we need Trk::Vertex
+    Trk::Vertex* vertex = findSpatialPoint(perigee_1, perigee_2);
+    return vertex;
+  }//if measured perigee
+}
+
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Rec::TrackParticle* trk_1)
+{
+  const Trk::TrackParameters* perigee_1(trk_1->measuredPerigee());
+  if (!perigee_1) {
+    if (msgLvl(MSG::VERBOSE)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex* vertex = findSpatialPoint(vtx, perigee_1);
+    return vertex;
+  }
+}
+
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::TrackParticleBase* trk_1, const Trk::TrackParticleBase* trk_2)
+{
+  const Trk::TrackParameters* perigee_1 = dynamic_cast<const Trk::TrackParameters*>(&(trk_1)->definingParameters());
+  const Trk::TrackParameters* perigee_2 = dynamic_cast<const Trk::TrackParameters*>(&(trk_2)->definingParameters());
+  if (!perigee_1 | !perigee_2) {
+    if (msgLvl(MSG::VERBOSE)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    //we need Trk::Vertex
+    Trk::Vertex* vertex = findSpatialPoint(perigee_1, perigee_2);
+    return vertex;
+  }//if measured perigee
+}
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParticleBase* trk_1)
+{
+  const Trk::TrackParameters* perigee_1 = dynamic_cast<const Trk::TrackParameters*>(&(trk_1)->definingParameters());
+  if (!perigee_1) {
+    if (msgLvl(MSG::VERBOSE)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex* vertex = findSpatialPoint(vtx, perigee_1);
+    return vertex;
+  }
+}
+
+//=============================================================================================
+
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::TrackParameters* perigee_1, const Trk::TrackParameters* perigee_2)
+{
+  Amg::Vector3D spatialPoint;
+  //we need Trk::Vertex
+  double chi2(0.);
+  const Amg::Vector3D lin_point(0.,0.,0.); // we start our linearization always at ATLAS origin!
+  Trk::Vertex linPoint(lin_point);
+  for (unsigned int i = 0; i < 3; i++){
+    Trk::LinearizedTrack* linTrack1 = m_linFactory->linearizedTrack (perigee_1, linPoint.position() );
+    Trk::LinearizedTrack* linTrack2 = m_linFactory->linearizedTrack (perigee_2, linPoint.position() );
+    if (linTrack1 && linTrack2)
+      {
+	Amg::Vector3D locXpVec_1 = linTrack1->expectedPositionAtPCA();
+	locXpVec_1[0] = locXpVec_1[0] - linPoint.position()[0];
+	locXpVec_1[1] = locXpVec_1[1] - linPoint.position()[1];
+	locXpVec_1[2] = locXpVec_1[2] - linPoint.position()[2];
+	
+	Amg::Vector3D locXpVec_2 = linTrack2->expectedPositionAtPCA();
+	locXpVec_2[0] = locXpVec_2[0] - linPoint.position()[0];
+	locXpVec_2[1] = locXpVec_2[1] - linPoint.position()[1];
+	locXpVec_2[2] = locXpVec_2[2] - linPoint.position()[2];
+	
+	AmgMatrix(2,2) billoirCovMat_1 = linTrack1->expectedCovarianceAtPCA().block<2,2>(0,0);
+	AmgMatrix(2,2) billoirCovMat_2 = linTrack2->expectedCovarianceAtPCA().block<2,2>(0,0);
+	AmgMatrix(2,2) billoirWeightMat_1 = billoirCovMat_1.inverse().eval();
+	AmgMatrix(2,2) billoirWeightMat_2 = billoirCovMat_2.inverse().eval();
+	AmgMatrix(3,3) DtWD_sum;
+	Amg::Vector3D DtWDx_sum;
+	AmgMatrix(3,3) DtWD_1;
+	AmgMatrix(3,3) DtWD_2;
+	// D matrix for d0 and z0
+	AmgMatrix(2,3) D_1 = linTrack1->positionJacobian().block<2,3>(0,0);
+	AmgMatrix(2,3) D_2 = linTrack2->positionJacobian().block<2,3>(0,0);
+	// Calculate DtWD and DtWD*x and sum them
+	DtWD_1 = D_1.transpose()*billoirWeightMat_1*D_1;
+	DtWD_2 = D_2.transpose()*billoirWeightMat_2*D_2;
+	DtWD_sum = DtWD_1 + DtWD_2;
+	DtWDx_sum = DtWD_1 * locXpVec_1 + DtWD_2* locXpVec_2;
+	AmgMatrix(3,3) cov_spatialPoint = DtWD_sum.inverse();
+	spatialPoint = cov_spatialPoint * DtWDx_sum;
+	Amg::Vector3D XpVec_sp_1 = locXpVec_1 - spatialPoint;
+	Amg::Vector3D XpVec_sp_2 = locXpVec_2 - spatialPoint;
+	chi2 =  XpVec_sp_1.dot(DtWD_1*XpVec_sp_1) + XpVec_sp_2.dot(DtWD_2*XpVec_sp_2);
+	spatialPoint[0] += linPoint.position()[0];
+	spatialPoint[1] += linPoint.position()[1];
+	spatialPoint[2] += linPoint.position()[2];
+	if (msgLvl(MSG::VERBOSE)) msg() <<"found spatial point = ("<<spatialPoint[0]<<", "<<spatialPoint[1]<<", "<<spatialPoint[2]<<")"<< endreq;
+	if (msgLvl(MSG::VERBOSE)) msg() <<"chi2 = "<<chi2<<endreq;
+	linPoint = Trk::Vertex(spatialPoint);
+      }	
+    delete linTrack1; linTrack1=0; delete linTrack2; linTrack2=0;
+  }// two iterations
+  if (chi2 <= m_chi2) return new Trk::Vertex(spatialPoint);
+  else {
+    if (msgLvl(MSG::VERBOSE)) msg() <<"found spatial point candidate doesn't pass chi2_cut" << endreq;
+    return 0;
+  }
+}
+
+//========================================================================================================
+Trk::Vertex* InDet::ZVTOP_SlowSpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParameters* perigee_1)
+{
+  Amg::Vector3D spatialPoint;
+  //we need Trk::Vertex
+  double chi2(0.);
+  const Amg::Vector3D lin_point(0.,0.,0.); // we start our linearization always at ATLAS origin!
+  Trk::Vertex linPoint(lin_point);
+  for (unsigned int i = 0; i < 3; i++){
+    Trk::LinearizedTrack* linTrack1 = m_linFactory->linearizedTrack (perigee_1, linPoint.position() );
+    if (linTrack1)
+      {
+	Amg::Vector3D locXpVec_1(0.,0.,0.);
+	AmgMatrix(3,3) DtWD_1;
+	AmgMatrix(3,3) vtx_weight;
+	locXpVec_1[0]= locXpVec_1[0] - linPoint.position()[0];
+	locXpVec_1[1]= locXpVec_1[1] - linPoint.position()[1];
+	locXpVec_1[2]= locXpVec_1[2] - linPoint.position()[2];
+
+	AmgMatrix(3,3) DtWD_sum;
+	Amg::Vector3D DtWDx_sum;
+	AmgMatrix(2,2) billoirCovMat_1 = linTrack1->expectedCovarianceAtPCA().block<2,2>(0,0);
+	AmgMatrix(2,2) billoirWeightMat_1 = billoirCovMat_1.inverse().eval();
+	// D matrix for d0 and z0
+	AmgMatrix(2,3) D_1 = linTrack1->positionJacobian().block<2,3>(0,0);
+	// Calculate DtWD and DtWD*x and sum them
+	DtWD_1 = D_1.transpose()*billoirWeightMat_1*D_1 ;
+	vtx_weight = vtx.covariancePosition().inverse();
+	DtWD_sum = DtWD_1 + vtx_weight;
+	Amg::Vector3D vtx_pos(3,0); vtx_pos[0] = vtx.position().x(); vtx_pos[1] = vtx.position().y(); vtx_pos[2] = vtx.position().z();
+	DtWDx_sum =  DtWD_1 * locXpVec_1 + vtx_weight*vtx_pos;
+	AmgMatrix(3,3) cov_spatialPoint = DtWD_sum.inverse();
+	spatialPoint = cov_spatialPoint * DtWDx_sum;
+	Amg::Vector3D XpVec_sp_1 = locXpVec_1 - spatialPoint;
+	Amg::Vector3D sp_vtx = spatialPoint - vtx_pos;
+	chi2 =  XpVec_sp_1.dot(DtWD_1*XpVec_sp_1) + sp_vtx.dot(vtx_weight*sp_vtx);
+	spatialPoint[0] += linPoint.position()[0];
+	spatialPoint[1] += linPoint.position()[1];
+	spatialPoint[2] += linPoint.position()[2];
+	if (msgLvl(MSG::VERBOSE)) msg() <<"found spatial point = ("<<spatialPoint[0]<<", "<<spatialPoint[1]<<", "<<spatialPoint[2]<<")"<< endreq;
+	if (msgLvl(MSG::VERBOSE)) msg() <<"chi2 = "<<chi2<<endreq;
+	linPoint = Trk::Vertex (spatialPoint);
+      }
+    delete linTrack1; linTrack1=0; 
+  }
+  if (chi2 <= m_chi2) return new Trk::Vertex(spatialPoint);
+  else {
+    if (msgLvl(MSG::VERBOSE)) msg() <<"found spatial point candidate doesn't pass chi2_cut" << endreq;
+    return 0;
+  }
+}
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SpatialPointFinder.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SpatialPointFinder.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..50cee085697a3fc0ead2f12ca1e7cde54f02804a
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_SpatialPointFinder.cxx
@@ -0,0 +1,250 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_SpatialPointFinder.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_SpatialPointFinder.h"
+#include "TrkParameters/TrackParameters.h"
+#include "VxVertex/Vertex.h"
+#include "Particle/TrackParticle.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include "TrkParticleBase/TrackParticleBase.h"
+//================ Constructor =================================================
+
+InDet::ZVTOP_SpatialPointFinder::ZVTOP_SpatialPointFinder(const std::string& t,
+			  const std::string& n,
+			  const IInterface*  p )
+  :
+  AthAlgTool(t,n,p)
+{
+  declareInterface<IZVTOP_SpatialPointFinder>(this);
+
+  //  template for property declaration
+  declareProperty("Chi2_cut_value", m_chi2);
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_SpatialPointFinder::~ZVTOP_SpatialPointFinder()
+{}
+
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_SpatialPointFinder::initialize()
+{
+  
+  StatusCode sc = AlgTool::initialize();
+
+  if (sc.isFailure()) return sc;
+
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================
+
+StatusCode InDet::ZVTOP_SpatialPointFinder::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::Track* trk_1, const Trk::Track* trk_2)
+{
+  const Trk::TrackParameters* perigee_1(dynamic_cast<const Trk::TrackParameters*>(trk_1->perigeeParameters()));
+  const Trk::TrackParameters* perigee_2(dynamic_cast<const Trk::TrackParameters*>(trk_2->perigeeParameters()));
+  if (!perigee_1 | !perigee_2) {
+    if (msgLvl(MSG::DEBUG)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex * vertex = findSpatialPoint(perigee_1,perigee_2);
+    return vertex;
+  }//if measured perigee
+}
+
+//===============================================================================================
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Trk::Track* trk_1)
+{
+  const Trk::TrackParameters *perigee_1(dynamic_cast<const Trk::TrackParameters*>(trk_1->perigeeParameters()));
+  if (!perigee_1) {
+    if (msgLvl(MSG::DEBUG)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex * vertex = findSpatialPoint(vtx,perigee_1);
+    return vertex;
+  }
+}
+
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Rec::TrackParticle* trk_1, const Rec::TrackParticle* trk_2)
+{
+  const Trk::TrackParameters* perigee_1(trk_1->measuredPerigee());
+  const Trk::TrackParameters* perigee_2(trk_2->measuredPerigee());
+  if (!perigee_1 | !perigee_2) {
+    if (msgLvl(MSG::DEBUG)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+		Trk::Vertex * vertex = findSpatialPoint(perigee_1,perigee_2);
+		return vertex;
+  }//if measured perigee
+}
+
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Rec::TrackParticle* trk_1)
+{
+  const Trk::TrackParameters* perigee_1(trk_1->measuredPerigee());
+  if (!perigee_1) {
+    if (msgLvl(MSG::DEBUG)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex * vertex = findSpatialPoint(vtx,perigee_1);
+    return vertex;
+  }
+}
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::TrackParticleBase* trk_1, const Trk::TrackParticleBase* trk_2)
+{
+  const Trk::TrackParameters* perigee_1 = dynamic_cast<const Trk::TrackParameters*>(&trk_1->definingParameters());
+  const Trk::TrackParameters* perigee_2 = dynamic_cast<const Trk::TrackParameters*>(&trk_2->definingParameters());
+  if (!perigee_1 | !perigee_2) {
+    if (msgLvl(MSG::DEBUG)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex * vertex = findSpatialPoint(perigee_1,perigee_2);
+    return vertex;
+  }//if measured perigee
+}
+
+//============================================================================================
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParticleBase* trk_1)
+{
+  const Trk::TrackParameters* perigee_1 = dynamic_cast<const Trk::TrackParameters*>(&trk_1->definingParameters());
+  if (!perigee_1) {
+    if (msgLvl(MSG::DEBUG)) msg() << "Dynamic cast to MeasuredPerigee failed. Skipping this pair" << endreq;
+    return 0;
+  } else {
+    Trk::Vertex * vertex = findSpatialPoint(vtx,perigee_1);
+    return vertex;
+  }
+}
+
+//===============================================================================================
+
+
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::TrackParameters* param_1, const Trk::TrackParameters* param_2)
+{
+  Amg::Vector3D spatialPoint;
+  const Trk::Perigee* perigee_1 =dynamic_cast<const Trk::Perigee*>(param_1);
+  const Trk::Perigee* perigee_2 =dynamic_cast<const Trk::Perigee*>(param_2);
+
+  double cot_theta_1 = 1./tan(perigee_1->parameters()[Trk::theta]);
+  double sphi_1 = sin(perigee_1->parameters()[Trk::phi]);
+  double cphi_1 = cos(perigee_1->parameters()[Trk::phi]);
+  Amg::Vector3D locXpVec_1;
+  locXpVec_1[0]= -perigee_1->parameters()[Trk::d0]*sphi_1;
+  locXpVec_1[1]= perigee_1->parameters()[Trk::d0]*cphi_1;
+  locXpVec_1[2]= perigee_1->parameters()[Trk::z0];
+  
+  double cot_theta_2 = 1./tan(perigee_2->parameters()[Trk::theta]);
+  double sphi_2 = sin(perigee_2->parameters()[Trk::phi]);
+  double cphi_2 = cos(perigee_2->parameters()[Trk::phi]);
+  Amg::Vector3D locXpVec_2;
+  locXpVec_2[0]= -perigee_2->parameters()[Trk::d0]*sphi_2;
+  locXpVec_2[1]= perigee_2->parameters()[Trk::d0]*cphi_2;
+  locXpVec_2[2]= perigee_2->parameters()[Trk::z0];
+  
+  AmgMatrix(3,3) DtWD_sum;
+  Amg::Vector3D DtWDx_sum;
+  double chi2(0.);
+  AmgMatrix(3,3) DtWD_1;
+  AmgMatrix(3,3) DtWD_2;
+  AmgMatrix(2,2) cov_mat_1;
+  AmgMatrix(2,2) cov_mat_2;
+  cov_mat_1 =(* perigee_1->covariance()).block<2,2>(0,0);
+  cov_mat_2 =(* perigee_2->covariance()).block<2,2>(0,0);
+
+  // and the weightmatrix
+  AmgMatrix(2,2) weight_mat_1 = cov_mat_1.inverse().eval();
+  AmgMatrix(2,2) weight_mat_2 = cov_mat_2.inverse().eval();
+
+  // calculate the D matrix
+  AmgMatrix(2,3) D_1;
+  AmgMatrix(2,3) D_2;
+  D_1(0,0) = -sphi_1              ; D_1(0,1) = cphi_1              ; D_1(0,2) = 0.;
+  D_1(1,0) = -cphi_1*cot_theta_1  ; D_1(1,1) = -sphi_1*cot_theta_1 ; D_1(1,2) = 1.;
+  D_2(0,0) = -sphi_2              ; D_2(0,1) = cphi_2              ; D_2(0,2) = 0.;
+  D_2(1,0) = -cphi_2*cot_theta_2  ; D_2(1,1) = -sphi_2*cot_theta_2 ; D_2(1,2) = 1.;
+  // Calculate DtWD and DtWD*x and sum them
+  DtWD_1 = D_1.transpose()*weight_mat_1*D_1;
+  DtWD_2 = D_2.transpose()*weight_mat_2*D_2;
+  DtWD_sum = DtWD_1 + DtWD_2;
+  DtWDx_sum = DtWD_1*locXpVec_1 + DtWD_2*locXpVec_2;
+  AmgMatrix(3,3) cov_spatialPoint = DtWD_sum.inverse();
+  spatialPoint = cov_spatialPoint * DtWDx_sum;
+  Amg::Vector3D XpVec_sp_1 = locXpVec_1 - spatialPoint;
+  Amg::Vector3D XpVec_sp_2 = locXpVec_2 - spatialPoint;
+  chi2 =  XpVec_sp_1.dot(DtWD_1*XpVec_sp_1) + XpVec_sp_2.dot(DtWD_2*XpVec_sp_2);
+  if (chi2 <= m_chi2) 
+    { 
+      if (msgLvl(MSG::DEBUG)) msg() <<"found spatial point = ("<<spatialPoint[0]<<", "<<spatialPoint[1]<<", "<<spatialPoint[2]<<")"<< endreq;
+      return new Trk::Vertex(spatialPoint);
+    } else 
+    {
+      if (msgLvl(MSG::DEBUG)) msg() <<"found spatial point candidate doesn't pass chi2_cut" << endreq;
+      return 0;
+    }
+     	
+}
+
+
+//============================================================================================
+
+Trk::Vertex* InDet::ZVTOP_SpatialPointFinder::findSpatialPoint(const Trk::RecVertex vtx, const Trk::TrackParameters* param_1)
+{
+  Amg::Vector3D spatialPoint;
+  const Trk::Perigee* perigee_1 =dynamic_cast<const Trk::Perigee*>(param_1);
+  double cot_theta_1 = 1./tan(perigee_1->parameters()[Trk::theta]);
+  double sphi_1 = sin(perigee_1->parameters()[Trk::phi]);
+  double cphi_1 = cos(perigee_1->parameters()[Trk::phi]);
+  Amg::Vector3D locXpVec_1;
+  AmgMatrix(3,3) DtWD_1;
+  AmgMatrix(3,3) vtx_weight;
+  double chi2(0.);
+  locXpVec_1[0]= -perigee_1->parameters()[Trk::d0]*sphi_1;
+  locXpVec_1[1]= perigee_1->parameters()[Trk::d0]*cphi_1;
+  locXpVec_1[2]= perigee_1->parameters()[Trk::z0];
+
+  AmgMatrix(3,3) DtWD_sum;
+  Amg::Vector3D DtWDx_sum;
+  AmgMatrix(2,2) cov_mat_1 = (*perigee_1->covariance()).block<2,2>(0,0);
+  // and the weightmatrix
+  AmgMatrix(2,2) weight_mat_1 = cov_mat_1.inverse();
+  AmgMatrix(2,3) D_1;
+  D_1(0,0) = -sphi_1              ; D_1(0,1) = cphi_1              ; D_1(0,2) = 0.;
+  D_1(1,0) = -cphi_1*cot_theta_1  ; D_1(1,1) = -sphi_1*cot_theta_1 ; D_1(1,2) = 1.;
+  DtWD_1 = D_1.transpose()*weight_mat_1*D_1;
+  vtx_weight = vtx.covariancePosition().inverse();
+  DtWD_sum = DtWD_1 + vtx_weight;
+  Amg::Vector3D vtx_pos; vtx_pos[0] = vtx.position().x(); vtx_pos[1] = vtx.position().y(); vtx_pos[2] = vtx.position().z();
+  DtWDx_sum = DtWD_1*locXpVec_1 + vtx_weight*vtx_pos;
+  AmgMatrix(3,3) cov_spatialPoint = DtWD_sum.inverse();
+  spatialPoint = cov_spatialPoint * DtWDx_sum;
+  Amg::Vector3D XpVec_sp_1 = locXpVec_1 - spatialPoint;
+  Amg::Vector3D sp_vtx = spatialPoint - vtx_pos;
+  chi2 =  XpVec_sp_1.dot(DtWD_1*XpVec_sp_1) + sp_vtx.dot(vtx_weight*sp_vtx);
+  if (chi2 <= m_chi2) 
+    { 
+      if (msgLvl(MSG::DEBUG)) msg() <<"found spatial point = ("<<spatialPoint[0]<<", "<<spatialPoint[1]<<", "<<spatialPoint[2]<<")"<< endreq;
+      return new Trk::Vertex(spatialPoint);
+    } else 
+    {
+      if (msgLvl(MSG::DEBUG)) msg() <<"found spatial point candidate doesn't pass chi2_cut" << endreq;
+      return 0;
+    }
+
+}
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_Tool.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_Tool.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..c1c280d04884900480fcdedcfb7f89756c3eea9e
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_Tool.cxx
@@ -0,0 +1,853 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_Tool.cxx, (c) ATLAS Detector software
+// begin   : 30-10-2006
+// authors : Tatjana Lenz
+// email   : tatjana.lenz@cern.ch
+// changes :
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_Tool.h"
+#include "TrkVertexFitterInterfaces/IVertexFitter.h"
+#include "InDetZVTOPVxFinder/ZVTOP_SpatialPointFinder.h"
+#include "InDetZVTOPVxFinder/ZVTOP_TrkPartBaseVertexState.h"
+#include "InDetZVTOPVxFinder/ZVTOP_VertexState.h"
+#include "InDetZVTOPVxFinder/IZVTOP_TrkProbTubeCalc.h"
+#include "InDetZVTOPVxFinder/IZVTOP_SimpleVtxProbCalc.h"
+#include "VxVertex/Vertex.h"
+#include "VxVertex/VxTrackAtVertex.h"
+#include "TrkTrack/LinkToTrack.h"
+#include "InDetBeamSpotService/IBeamCondSvc.h"
+#include "DataModel/DataVector.h"
+#include "TrkParticleBase/LinkToTrackParticleBase.h"
+#include "TrkParticleBase/TrackParticleBaseCollection.h"
+#include "TrkTrack/TrackCollection.h"
+#include "VxVertex/VxContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include <map>
+#include <utility>
+//================ Constructor =================================================
+
+InDet::ZVTOP_Tool::ZVTOP_Tool(const std::string& t,
+			  const std::string& n,
+			  const IInterface*  p )
+  :
+  AthAlgTool(t,n,p),
+  m_iVertexFitter("Trk::FullVertexFitter"),
+  m_iSpatialPointFinder("InDet::SpatialPointFinder"),
+  m_iTrkProbTubeCalc("InDet::TrkProbTubeCalc"),
+  m_iVtxProbCalc("InDet::VtxProbCalc"),
+  m_iBeamCondSvc("BeamCondSvc",n)
+{
+  //  template for property declaration
+  declareInterface<IVertexFinder>(this);
+  declareProperty("VertexFitterTool",m_iVertexFitter);
+  declareProperty("SpatialPointFinderTool",m_iSpatialPointFinder);
+  declareProperty("TrkProbTubeCalcTool",m_iTrkProbTubeCalc);
+  declareProperty("VtxProbCalcTool",m_iVtxProbCalc);
+  declareProperty("ResolvingCutValue", m_resolvingCut = 0.6);
+  declareProperty("ResolvingStepValue", m_resolvingStep = 0.3);
+  declareProperty("VtxProbCut", m_vtxProb_cut = 0.00001);
+  declareProperty("MaxChi2PerTrack", m_trk_chi2_cut = 5.);
+  declareProperty("BeamPositionSvc", m_iBeamCondSvc);
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_Tool::~ZVTOP_Tool()
+{}
+
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_Tool::initialize()
+{
+  StatusCode sc = AlgTool::initialize();
+
+  msg (MSG::INFO) << name() << " initialize()" << endreq;
+  if (sc.isFailure()) return sc;
+
+  /* Retrieve Tools*/
+  //SpatialPointFinder
+  if ( m_iSpatialPointFinder.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iSpatialPointFinder << endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iSpatialPointFinder << endreq;
+
+  //Gaussian Probability Tube for the Track Trajectory
+  if ( m_iTrkProbTubeCalc.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iTrkProbTubeCalc<< endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iTrkProbTubeCalc << endreq;
+
+  //Vertex Probability Function
+  if ( m_iVtxProbCalc.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iVtxProbCalc<< endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iVtxProbCalc << endreq;
+
+  //Beam Spot
+  sc = m_iBeamCondSvc.retrieve();
+  if (sc.isFailure())
+  {
+    msg (MSG::ERROR) << "Could not find BeamCondSvc." << endreq;
+    return sc;
+  }
+  //VxFitter
+  if ( m_iVertexFitter.retrieve().isFailure() ) {
+      msg (MSG::FATAL) << "Failed to retrieve tool " << m_iVertexFitter << endreq;
+      return StatusCode::FAILURE;
+  } else msg (MSG::INFO) << "Retrieved tool " << m_iVertexFitter << endreq;
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================
+
+StatusCode InDet::ZVTOP_Tool::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+
+//============================================================================================
+//VxContainer* InDet::ZVTOP_Tool::findVertex(const TrackCollection* trackTES) --David S.
+std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> InDet::ZVTOP_Tool::findVertex(const TrackCollection* trackTES)
+{
+//VxContainer* newVxContainer = new VxContainer; --David S.
+xAOD::VertexContainer* newVertexContainer = new xAOD::VertexContainer;
+xAOD::VertexAuxContainer* newVertexAuxContainer = new xAOD::VertexAuxContainer;
+newVertexContainer->setStore( newVertexAuxContainer );
+if (trackTES->size() != 0) {
+  //some variables
+  typedef DataVector<Trk::Track>::const_iterator TrackDataVecIter;
+  const Trk::RecVertex beam_spot = m_iBeamCondSvc->beamVtx();
+
+  // new Fitter EDM takes xAOD::Vertex instead of Trk::RecVertex as beam_spot constraint, so create one out of beamVtx() --David S.
+  xAOD::Vertex theconstraint = xAOD::Vertex(); // Default constructor creates a private store
+  theconstraint.setPosition( beam_spot.position() );
+  theconstraint.setCovariancePosition( beam_spot.covariancePosition() );
+  theconstraint.setFitQuality( beam_spot.fitQuality().chiSquared(), beam_spot.fitQuality().doubleNumberDoF() );
+
+  std::vector<const Trk::Track*> trkColl(0); // the collection of tracks, which are assigned to one spatial point
+  std::vector<InDet::ZVTOP_VertexState*> vtxState_org(0);// vector of tracks and vertices, tracks can be assigned to more than one vertex, these ambiguities will be resolved later
+  std::multimap<double, InDet::ZVTOP_VertexState*> vtxProb_map;
+  //---------------------------------------------------------------------------//  
+  //-------step1: find spatial points & calculate vtx probabilities-------//
+  //--------------------------------------------------------------------------//
+
+  std::vector< std::vector<InDet::ZVTOP_VertexState*> > vertex_objects(0); //vector of vtx candidate & track collection pairs
+  
+  if (msgLvl(MSG::DEBUG)) msg() << "step ONE search for the spatial points, number of tracks = "<<(*trackTES).size() <<  endreq;
+  for (TrackDataVecIter itr_1  = (*trackTES).begin(); itr_1 != (*trackTES).end()-1; itr_1++)
+  {
+    double sum_TrkProbTube = 0.;
+    double sum2_TrkProbTube = 0.;
+    double vtxProb = 0.;
+    // trk-trk combination
+    for (TrackDataVecIter itr_2 = itr_1+1; itr_2 != (*trackTES).end(); itr_2++)
+    {
+      Trk::Vertex* spatialPoint;
+      spatialPoint = m_iSpatialPointFinder->findSpatialPoint((*itr_1),(*itr_2));
+      if (spatialPoint != 0) 
+      {
+         double TrkProbTube_1 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+         double TrkProbTube_2 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_2),*spatialPoint);
+         sum_TrkProbTube = TrkProbTube_1 + TrkProbTube_2;
+         sum2_TrkProbTube = pow(TrkProbTube_1,2.) + pow(TrkProbTube_2,2.);
+         if (sum_TrkProbTube != 0. )
+         { 
+            vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+            trkColl.push_back((*itr_1));
+            trkColl.push_back((*itr_2));
+            bool beam_spot = false;
+            if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_VertexState(vtxProb, *spatialPoint, beam_spot, trkColl ));
+            trkColl.clear();
+         }
+      }
+      delete spatialPoint;
+    }
+    //trk-IP combination
+    Trk::Vertex* spatialPoint = 0;
+    spatialPoint = m_iSpatialPointFinder->findSpatialPoint(beam_spot,(*itr_1));
+    if (spatialPoint != 0) 
+    {
+      double BeamProbTube = m_iTrkProbTubeCalc->calcProbTube(beam_spot,*spatialPoint);
+      double TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+      sum_TrkProbTube = BeamProbTube + TrkProbTube;
+      sum2_TrkProbTube = pow(TrkProbTube,2.) + pow(BeamProbTube,2.);
+      if (sum_TrkProbTube != 0. ) 
+      {
+         vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+         trkColl.push_back((*itr_1));
+         bool beam_spot = true;
+         if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_VertexState(vtxProb, *spatialPoint, beam_spot, trkColl ));
+         trkColl.clear();
+      }
+    }
+    delete spatialPoint;
+  }
+  if (msgLvl(MSG::DEBUG)) msg() << " number of spatial points = "<<vtxState_org.size()<<endreq;
+  //reduce the spatial point collection
+  typedef std::vector<InDet::ZVTOP_VertexState*>::iterator Sp_Iter;
+  std::vector< InDet::ZVTOP_VertexState*> vtxState;
+  for (Sp_Iter itr1 = vtxState_org.begin(); itr1 != vtxState_org.end(); itr1++)
+  {
+
+     if (vtxState.size() == 0) vtxState.push_back(*itr1);
+     else {
+         Trk::Vertex vtx_new = (*itr1)->vertex();
+         bool can_be_resolved = false;
+         for (Sp_Iter itr2 = vtxState.begin(); itr2 != vtxState.end(); itr2++)
+        {
+           Trk::Vertex vtx_in = (*itr2)->vertex();
+            double r_diff = sqrt(pow(vtx_new.position().x(),2.) + pow(vtx_new.position().y(),2.)) - sqrt(pow(vtx_in.position().x(),2.) + pow(vtx_in.position().y(),2.));
+            double z_diff = vtx_new.position().z() - vtx_in.position().z();
+            if (fabs(r_diff) > 0.001 && fabs(z_diff) > 0.001 && r_diff != 0. && z_diff != 0.) can_be_resolved = true;
+            else {
+                for (unsigned int trk_itr = 0; trk_itr < (*itr1)->tracks().size(); trk_itr++) (*itr2)->tracks().push_back((*itr1)->tracks()[trk_itr]); 
+                can_be_resolved = false;
+                break;
+            }
+         }
+         if (can_be_resolved)  vtxState.push_back(*itr1);
+      }
+   }
+   for (Sp_Iter itr_sort = vtxState.begin();  itr_sort!= vtxState.end(); itr_sort++)  vtxProb_map.insert(std::multimap<double, InDet::ZVTOP_VertexState*>::value_type((*itr_sort)->vtx_prob(), *itr_sort));
+
+
+//-------------------------------------------------------------------------------//
+//----step2: vertex clustering, aim: to cluster all vertices together, ----//
+//----which can not be resolved. The found vertex seeds and ---------//
+//----associated tracks are taken as an input for the vertex fit. --------//
+//------------------------------------------------------------------------------//
+
+
+  if (vtxState.size() != 0 ){
+    if (msgLvl(MSG::DEBUG)) msg() << " step TWO vertex clustering, number of reduced vertices = "<<vtxState.size()<< endreq;
+    //sort the vtxState collection, starting with a highest vtx probability
+    std::vector<InDet::ZVTOP_VertexState*> vtxState_sorted;
+    unsigned int IP_counter = 0;
+    std::multimap<double, InDet::ZVTOP_VertexState*>::reverse_iterator s= vtxProb_map.rbegin();
+    for (;  s!=vtxProb_map.rend();  s++)  {
+           InDet::ZVTOP_VertexState* vtx_state = (*s).second;
+           if (vtx_state->beam_spot()) IP_counter += 1;
+           if (IP_counter > 1) vtx_state->set_beam_spot(false);
+           vtxState_sorted.push_back(vtx_state);
+     }
+    //store first element in the vertex_objects ---->vector<vector<VertexState>>
+    std::vector<InDet::ZVTOP_VertexState*> vtx_coll;
+    vtx_coll.push_back(*(vtxState_sorted.begin()));
+    vertex_objects.push_back(vtx_coll);//store first seed
+    std::vector< InDet::ZVTOP_VertexState*>  vtxState_new(vtxState_sorted); //copy of the vtxState_sorted
+    vtxState_new.erase(vtxState_new.begin()); // without first element
+    //loop over vtxState_sorted
+    for (Sp_Iter org_itr= vtxState_sorted.begin(); org_itr != vtxState_sorted.end(); org_itr++)
+    {
+       Trk::Vertex seed = (*org_itr)->vertex();
+        //found where the seed is already stored
+        bool vtx_is_stored = false;
+        unsigned int vertex_objectsPosition = 0; 
+        for (unsigned int vertex_objectsItr = 0; vertex_objectsItr!= vertex_objects.size(); vertex_objectsItr++)
+        {
+          //stored vertices in the line vertex_objectsItr
+          std::vector<InDet::ZVTOP_VertexState*> stored_vertices = vertex_objects[vertex_objectsItr];
+           //---->inner loop
+          for (unsigned int stored_vtx_itr = 0; stored_vtx_itr!= stored_vertices.size(); stored_vtx_itr++)
+          {
+              Trk::Vertex storedVtx = stored_vertices[stored_vtx_itr]->vertex();
+              double diff = ((*org_itr)->vertex().position().x() - storedVtx.position().x()) + ((*org_itr)->vertex().position().y() - storedVtx.position().y()) + ((*org_itr)->vertex().position().z() - storedVtx.position().z());
+              if (fabs(diff) < 0.0001)
+              {
+                 vertex_objectsPosition = vertex_objectsItr;
+                 vtx_is_stored = true;
+                 break; // break inner loop if found
+              }
+            }
+            if (vtx_is_stored == true) break; // break outer loop if found
+         }
+         if (!vtx_is_stored) {
+                 //if not stored
+                std::vector<InDet::ZVTOP_VertexState*> new_vtx_coll;
+                new_vtx_coll.push_back(*org_itr);
+                vertex_objects.push_back(new_vtx_coll);
+                vertex_objectsPosition = vertex_objects.size() - 1;
+         }
+         for (Sp_Iter new_itr = vtxState_new.begin(); new_itr!= vtxState_new.end(); new_itr++)
+         {
+             Trk::Vertex cand = (*new_itr)->vertex();
+             //calculate the vtx prob function for this seed
+             std::multimap<double, Trk::Vertex > vtx_prob_coll;
+             vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*new_itr)->vtx_prob(),(*new_itr)->vertex()));
+             vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*org_itr)->vtx_prob(),(*org_itr)->vertex()));
+             for (double alpha = m_resolvingStep; alpha <1.; alpha += m_resolvingStep)
+             {
+                 Trk::Vertex SP_line = Trk::Vertex((*org_itr)->vertex().position() + alpha*((*new_itr)->vertex().position() - (*org_itr)->vertex().position()));
+                  double sum_TrkProbTube = 0.;
+                  double sum2_TrkProbTube = 0.;
+                  for (TrackDataVecIter trk_itr = (*trackTES).begin(); trk_itr != (*trackTES).end(); trk_itr++)
+                  {
+	              double TrkProbTube = 0.;
+                      TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*trk_itr),SP_line);
+                       sum_TrkProbTube += TrkProbTube;
+                       sum2_TrkProbTube += pow(TrkProbTube,2.);
+                   }
+                   double BeamProbTube = m_iTrkProbTubeCalc->calcProbTube(beam_spot,SP_line);
+                   sum_TrkProbTube += BeamProbTube;
+                   sum2_TrkProbTube += pow(BeamProbTube,2.);
+                   double vtx_prob = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+                    vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type(vtx_prob, SP_line));
+               }
+              double vtx_prob_ratio = (*vtx_prob_coll.begin()).first/(*new_itr)->vtx_prob();
+                if (vtx_prob_ratio >= m_resolvingCut)  {
+                       //check if the vertices can be resolved
+                       vertex_objects[vertex_objectsPosition].push_back(*new_itr);
+                       vtxState_new.erase(new_itr);
+                       if (new_itr != vtxState_new.end()) --new_itr;
+                       if (new_itr == vtxState_new.end()) break;
+                }
+             }
+      }//loop over orig reduced coll
+   }//if spatial point collection size > 0
+
+
+   //--------------------------------------------------------------------------------//
+  //--------------------step3: call vertex fitter----------------------------------//
+  //--------------------------------------------------------------------------------//
+
+    if (msgLvl(MSG::DEBUG)) msg() << " step THREE, vertex fit, vertex_objects size = "<<vertex_objects.size()<< endreq;
+    //std::vector<Trk::VxCandidate* > theVxContainer; --David S.
+    std::vector< xAOD::Vertex* > theVertexContainer;
+    //prepare trk coll --> remove double tracks
+    std::vector<std::vector<InDet::ZVTOP_VertexState*> >::iterator vtx_itr = vertex_objects.begin();
+    for (; vtx_itr!= vertex_objects.end(); vtx_itr++)
+    {
+      bool with_beam_spot = false;
+      std::vector <const Trk::Track*> trk_coll(0); //--->trk_coll for the fit
+      std::vector<InDet::ZVTOP_VertexState*>::iterator itr = (*vtx_itr).begin();
+      trk_coll.push_back((*itr)->tracks()[0]);
+      for (; itr != (*vtx_itr).end(); itr++)
+      { 
+        for ( std::vector <const Trk::Track*>::iterator vs_itr = ((*itr)->tracks()).begin(); vs_itr!= ((*itr)->tracks()).end(); vs_itr++)
+        {
+          bool found = false;
+	  for (std::vector <const Trk::Track*>::iterator trk_itr = trk_coll.begin();  trk_itr!= trk_coll.end(); trk_itr++)  {
+             if (*vs_itr == *trk_itr)  found = true;
+          }
+          if (!found)  trk_coll.push_back(*vs_itr);
+        }
+       if ((*itr)->beam_spot()) with_beam_spot = true;
+       }
+       //call the fitter
+       //Trk::VxCandidate * myVxCandidate(0); --David S.
+       xAOD::Vertex * myxAODVertex(0);
+       //if (with_beam_spot) myVxCandidate = m_iVertexFitter->fit(trk_coll,beam_spot); --David S.
+       if (with_beam_spot) myxAODVertex = m_iVertexFitter->fit(trk_coll,theconstraint);
+       //else myVxCandidate = m_iVertexFitter->fit(trk_coll); --David S.
+       else myxAODVertex = m_iVertexFitter->fit(trk_coll);
+       bool bad_chi2 = true;
+       //if (myVxCandidate) { --David S.
+       if (myxAODVertex) {
+          while (bad_chi2)
+          {
+             //std::vector< Trk::VxTrackAtVertex*> trkAtVtx = *(myVxCandidate->vxTrackAtVertex()); --David S.
+             std::vector< Trk::VxTrackAtVertex > trkAtVtx = myxAODVertex->vxTrackAtVertex();
+             //std::vector< Trk::VxTrackAtVertex*>::iterator trkAtVtx_Iter = trkAtVtx.begin(); --David S.
+             std::vector< Trk::VxTrackAtVertex >::iterator trkAtVtx_Iter = trkAtVtx.begin();
+             std::vector< const Trk::Track*>::iterator trk_Iter = trk_coll.begin();
+             double largest_chi2 = 0.;
+             std::vector< const Trk::Track*>::iterator index;
+             for (; trkAtVtx_Iter!= trkAtVtx.end(); trkAtVtx_Iter++)
+             {
+               double chi2 = (*trkAtVtx_Iter).trackQuality().chiSquared();
+               if (chi2 > largest_chi2) {
+                   largest_chi2 = chi2;
+                   index = trk_Iter;
+               }
+               trk_Iter++;
+             }
+             if (largest_chi2 > m_trk_chi2_cut)
+             {
+               if (trk_coll.size() < 3) break;
+               else {
+                  trk_coll.erase(index);
+                  if (trk_coll.size() >= 2) {
+                    //if (myVxCandidate!=0) { delete myVxCandidate; myVxCandidate=0; } --David S.
+                    if (myxAODVertex!=0) { delete myxAODVertex; myxAODVertex=0; }
+                    //if (with_beam_spot) myVxCandidate = m_iVertexFitter->fit(trk_coll, beam_spot); --David S.
+                    if (with_beam_spot) myxAODVertex = m_iVertexFitter->fit(trk_coll, theconstraint);
+                    //else myVxCandidate = m_iVertexFitter->fit(trk_coll); --David S.
+                    else myxAODVertex = m_iVertexFitter->fit(trk_coll);
+                  }
+                  //if (myVxCandidate == 0) break; --David S.
+                  if (myxAODVertex == 0) break;
+               }
+             } else bad_chi2 = false;
+           }
+      }
+      //if (myVxCandidate && bad_chi2 == false && with_beam_spot) newVxContainer->push_back(myVxCandidate); --David S.
+      if (myxAODVertex && bad_chi2 == false && with_beam_spot) newVertexContainer->push_back(myxAODVertex);
+      //if (myVxCandidate && bad_chi2 == false && !with_beam_spot) theVxContainer.push_back(myVxCandidate); --David S.
+      if (myxAODVertex && bad_chi2 == false && !with_beam_spot) theVertexContainer.push_back(myxAODVertex);
+    }
+    //if (msgLvl(MSG::DEBUG)) msg() <<"vertex container size = "<<theVxContainer.size()<<endreq; --David S.
+    if (msgLvl(MSG::DEBUG)) msg() << "vertex container size = " << theVertexContainer.size() << endreq;
+
+   //ambiguity solving for two track vertices
+   //typedef std::vector<Trk::VxCandidate* >::iterator theVxContainerIter; --David S.
+   typedef std::vector< xAOD::Vertex* >::iterator theVertexContainerIter;
+   //theVxContainerIter iter = theVxContainer.begin(); --David S.
+   theVertexContainerIter iter = theVertexContainer.begin();
+   //theVxContainerIter iter_end = theVxContainer.end(); --David S.
+   theVertexContainerIter iter_end = theVertexContainer.end();
+   //std::vector<Trk::VxCandidate*> twoTrkContainer(0); --David S.
+   std::vector< xAOD::Vertex* > twoTrkContainer(0);
+   for (; iter!= iter_end; iter++)
+   {
+     //Trk::VxCandidate* theVxCandidate(*iter); --David S.
+     xAOD::Vertex* thexAODVertex(*iter);
+     //if (theVxCandidate->vxTrackAtVertex()->size() == 2 ) twoTrkContainer.push_back(theVxCandidate); --David S.
+     if (thexAODVertex->vxTrackAtVertex().size() == 2 ) twoTrkContainer.push_back(thexAODVertex);
+     //else newVxContainer->push_back(theVxCandidate); --David S.
+     else newVertexContainer->push_back(thexAODVertex);
+   }
+   //std::vector<Trk::VxCandidate*>::iterator twoTrk_itr = twoTrkContainer.begin(); --David S.
+   std::vector< xAOD::Vertex* >::iterator twoTrk_itr = twoTrkContainer.begin();
+   for (; twoTrk_itr!= twoTrkContainer.end(); twoTrk_itr++)
+   {
+     // TODO: make sure that links are actually set in VxTrackAtVertices! --David S.
+     //Trk::ITrackLink* trklink1 = (*(*twoTrk_itr)->vxTrackAtVertex())[0]->trackOrParticleLink(); --David S.
+     Trk::ITrackLink* trklink1 = (*twoTrk_itr)->vxTrackAtVertex()[0].trackOrParticleLink();
+     Trk::LinkToTrack* linkToTrack1 = dynamic_cast<Trk::LinkToTrack*>(trklink1);
+     const Trk::Track* first_trk = linkToTrack1->cachedElement();
+     //double first_trk_chi2 = (*(*twoTrk_itr)->vxTrackAtVertex())[0]->trackQuality().chiSquared(); --David S.
+     double first_trk_chi2 = (*twoTrk_itr)->vxTrackAtVertex()[0].trackQuality().chiSquared();
+     //Trk::ITrackLink* trklink2 = (*(*twoTrk_itr)->vxTrackAtVertex())[1]->trackOrParticleLink(); --David S.
+     Trk::ITrackLink* trklink2 = (*twoTrk_itr)->vxTrackAtVertex()[1].trackOrParticleLink();
+     Trk::LinkToTrack* linkToTrack2 = dynamic_cast<Trk::LinkToTrack*>(trklink2);
+     const Trk::Track* second_trk = linkToTrack2->cachedElement();
+     //double second_trk_chi2 = (*(*twoTrk_itr)->vxTrackAtVertex())[1]->trackQuality().chiSquared(); --David S.
+     double second_trk_chi2 = (*twoTrk_itr)->vxTrackAtVertex()[1].trackQuality().chiSquared();
+     
+     //VxContainer::iterator Citer = newVxContainer->begin(); --David S.
+     xAOD::VertexContainer::iterator Citer = newVertexContainer->begin();
+     //VxContainer::iterator Citer_end = newVxContainer->end(); --David S.
+     xAOD::VertexContainer::iterator Citer_end = newVertexContainer->end();
+     bool first_found = false;
+     bool second_found = false;
+     ////loop over all VxCandidates --David S.
+     //loop over all Vertices
+     for (; Citer != Citer_end; Citer++)
+     {
+       //unsigned int size = (*Citer)->vxTrackAtVertex()->size(); --David S.
+       unsigned int size = (*Citer)->vxTrackAtVertex().size();
+       for (unsigned int counter = 0; counter < size; counter++)
+       {  
+         // TODO: make sure that links are actually set in VxTrackAtVertices! --David S.
+         //Trk::ITrackLink* trklink = (*(*Citer)->vxTrackAtVertex())[counter]->trackOrParticleLink(); --David S.
+         Trk::ITrackLink* trklink = (*Citer)->vxTrackAtVertex()[counter].trackOrParticleLink();
+         Trk::LinkToTrack* linkToTrack = dynamic_cast<Trk::LinkToTrack*>(trklink);
+         const Trk::Track* trk= linkToTrack->cachedElement();
+         //double trk_chi2 = (*(*Citer)->vxTrackAtVertex())[counter]->trackQuality().chiSquared(); --David S.
+         double trk_chi2 = (*Citer)->vxTrackAtVertex()[counter].trackQuality().chiSquared();
+         if (trk == first_trk && trk_chi2 < first_trk_chi2 ) first_found = true;
+         if (trk == second_trk && trk_chi2 < second_trk_chi2) second_found = true;
+        }
+     }
+     //if (first_found==false && second_found==false) newVxContainer->push_back(*twoTrk_itr); --David S.
+     if (first_found==false && second_found==false) newVertexContainer->push_back(*twoTrk_itr);
+     else delete *twoTrk_itr;
+   //}//end loop over two track Vx Candidates --David S.
+   }//end loop over two track Vertices
+   //if (msgLvl(MSG::DEBUG)) msg() <<"new VxContainer size = "<<newVxContainer->size()<<endreq; --David S.
+   if (msgLvl(MSG::DEBUG)) msg() << "new VertexContainer size = " << newVertexContainer->size() << endreq;
+   for (Sp_Iter iter = vtxState_org.begin(); iter != vtxState_org.end(); iter++) delete *iter;
+} else {
+    if (msgLvl(MSG::DEBUG)) msg() <<"No tracks in this event, provide to the next event" <<  endreq;
+  }
+
+//return newVxContainer; --David S.
+return std::make_pair(newVertexContainer, newVertexAuxContainer);
+}
+//////////////////////////////////////////////////
+///run on AOD
+/////////////////////////////////////////////////
+//VxContainer* InDet::ZVTOP_Tool::findVertex(const Trk::TrackParticleBaseCollection* trackTES) --David S.
+std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> InDet::ZVTOP_Tool::findVertex(const Trk::TrackParticleBaseCollection* trackTES)
+{
+//VxContainer* newVxContainer = new VxContainer; --David S.
+xAOD::VertexContainer* newVertexContainer = new xAOD::VertexContainer;
+xAOD::VertexAuxContainer* newVertexAuxContainer = new xAOD::VertexAuxContainer;
+newVertexContainer->setStore( newVertexAuxContainer );
+if (trackTES->size() != 0) {
+  //some variables
+  typedef Trk::TrackParticleBaseCollection::const_iterator TrackDataVecIter;
+  const Trk::RecVertex beam_spot = m_iBeamCondSvc->beamVtx();
+
+  // new Fitter EDM takes xAOD::Vertex instead of Trk::RecVertex as beam_spot constraint, so create one out of beamVtx() --David S.
+  xAOD::Vertex theconstraint = xAOD::Vertex(); // Default constructor creates a private store
+  theconstraint.setPosition( beam_spot.position() );
+  theconstraint.setCovariancePosition( beam_spot.covariancePosition() );
+  theconstraint.setFitQuality( beam_spot.fitQuality().chiSquared(), beam_spot.fitQuality().doubleNumberDoF() );
+
+  std::vector<const Trk::TrackParticleBase*> trkColl(0); // the collection of tracks, which are assigned to one spatial point
+  std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> vtxState_org(0);// vector of tracks and vertices, tracks can be assigned to more than one vertex, these ambiguities will be resolved later
+  std::multimap<double, InDet::ZVTOP_TrkPartBaseVertexState*> vtxProb_map;
+  //---------------------------------------------------------------------------//  
+  //-------step1: find spatial points & calculate vtx probabilities-------//
+  //--------------------------------------------------------------------------//
+
+  std::vector< std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> > vertex_objects(0); //vector of vtx candidate & track collection pairs
+  
+  if (msgLvl(MSG::DEBUG)) msg() << "step ONE search for the spatial points, number of tracks = "<<(*trackTES).size() <<  endreq;
+  for (TrackDataVecIter itr_1  = (*trackTES).begin(); itr_1 != (*trackTES).end()-1; itr_1++)
+  {
+    double sum_TrkProbTube = 0.;
+    double sum2_TrkProbTube = 0.;
+    double vtxProb = 0.;
+    // trk-trk combination
+    for (TrackDataVecIter itr_2 = itr_1+1; itr_2 != (*trackTES).end(); itr_2++)
+    {
+      Trk::Vertex* spatialPoint;
+      spatialPoint = m_iSpatialPointFinder->findSpatialPoint((*itr_1),(*itr_2));
+      if (spatialPoint != 0) 
+      {
+         double TrkProbTube_1 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+         double TrkProbTube_2 = m_iTrkProbTubeCalc->calcProbTube(*(*itr_2),*spatialPoint);
+         sum_TrkProbTube = TrkProbTube_1 + TrkProbTube_2;
+         sum2_TrkProbTube = pow(TrkProbTube_1,2.) + pow(TrkProbTube_2,2.);
+         if (sum_TrkProbTube != 0. )
+         { 
+            vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+            trkColl.push_back((*itr_1));
+            trkColl.push_back((*itr_2));
+            bool beam_spot = false;
+            if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_TrkPartBaseVertexState(vtxProb, *spatialPoint, beam_spot, trkColl ));
+            trkColl.clear();
+         }
+      }
+      delete spatialPoint;
+    }
+    //trk-IP combination
+    Trk::Vertex* spatialPoint = 0;
+    spatialPoint = m_iSpatialPointFinder->findSpatialPoint(beam_spot,(*itr_1));
+    if (spatialPoint != 0) 
+    {
+      double BeamProbTube = m_iTrkProbTubeCalc->calcProbTube(beam_spot,*spatialPoint);
+      double TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*itr_1),*spatialPoint);
+      sum_TrkProbTube = BeamProbTube + TrkProbTube;
+      sum2_TrkProbTube = pow(TrkProbTube,2.) + pow(BeamProbTube,2.);
+      if (sum_TrkProbTube != 0. ) 
+      {
+         vtxProb = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+         trkColl.push_back((*itr_1));
+         bool beam_spot = true;
+         if (vtxProb > m_vtxProb_cut) vtxState_org.push_back(new InDet::ZVTOP_TrkPartBaseVertexState(vtxProb, *spatialPoint, beam_spot, trkColl ));
+         trkColl.clear();
+      }
+    }
+    delete spatialPoint;
+  }
+  if (msgLvl(MSG::DEBUG)) msg() << " number of spatial points = "<<vtxState_org.size()<<endreq;
+  //reduce the spatial point collection
+  typedef std::vector<InDet::ZVTOP_TrkPartBaseVertexState*>::iterator Sp_Iter;
+  std::vector< InDet::ZVTOP_TrkPartBaseVertexState*> vtxState;
+  for (Sp_Iter itr1 = vtxState_org.begin(); itr1 != vtxState_org.end(); itr1++)
+  {
+
+     if (vtxState.size() == 0) vtxState.push_back(*itr1);
+     else {
+         Trk::Vertex vtx_new = (*itr1)->vertex();
+         bool can_be_resolved = false;
+         for (Sp_Iter itr2 = vtxState.begin(); itr2 != vtxState.end(); itr2++)
+        {
+           Trk::Vertex vtx_in = (*itr2)->vertex();
+            double r_diff = sqrt(pow(vtx_new.position().x(),2.) + pow(vtx_new.position().y(),2.)) - sqrt(pow(vtx_in.position().x(),2.) + pow(vtx_in.position().y(),2.));
+            double z_diff = vtx_new.position().z() - vtx_in.position().z();
+            if (fabs(r_diff) > 0.001 && fabs(z_diff) > 0.001 && r_diff != 0. && z_diff != 0.) can_be_resolved = true;
+            else {
+                for (unsigned int trk_itr = 0; trk_itr < (*itr1)->tracks().size(); trk_itr++) (*itr2)->tracks().push_back((*itr1)->tracks()[trk_itr]); 
+                can_be_resolved = false;
+                break;
+            }
+         }
+         if (can_be_resolved)  vtxState.push_back(*itr1);
+      }
+   }
+   for (Sp_Iter itr_sort = vtxState.begin();  itr_sort!= vtxState.end(); itr_sort++)  vtxProb_map.insert(std::multimap<double, InDet::ZVTOP_TrkPartBaseVertexState*>::value_type((*itr_sort)->vtx_prob(), *itr_sort));
+
+
+//-------------------------------------------------------------------------------//
+//----step2: vertex clustering, aim: to cluster all vertices together, ----//
+//----which can not be resolved. The found vertex seeds and ---------//
+//----associated tracks are taken as an input for the vertex fit. --------//
+//------------------------------------------------------------------------------//
+
+
+  if (vtxState.size() != 0 ){
+    if (msgLvl(MSG::DEBUG)) msg() << " step TWO vertex clustering, number of reduced vertices = "<<vtxState.size()<< endreq;
+    //sort the vtxState collection, starting with a highest vtx probability
+    std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> vtxState_sorted;
+    unsigned int IP_counter = 0;
+    std::multimap<double, InDet::ZVTOP_TrkPartBaseVertexState*>::reverse_iterator s= vtxProb_map.rbegin();
+    for (;  s!=vtxProb_map.rend();  s++)  {
+           InDet::ZVTOP_TrkPartBaseVertexState* vtx_state = (*s).second;
+           if (vtx_state->beam_spot()) IP_counter += 1;
+           if (IP_counter > 1) vtx_state->set_beam_spot(false);
+           vtxState_sorted.push_back(vtx_state);
+     }
+    //store first element in the vertex_objects ---->vector<vector<VertexState>>
+    std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> vtx_coll;
+    vtx_coll.push_back(*(vtxState_sorted.begin()));
+    vertex_objects.push_back(vtx_coll);//store first seed
+    std::vector< InDet::ZVTOP_TrkPartBaseVertexState*>  vtxState_new(vtxState_sorted); //copy of the vtxState_sorted
+    vtxState_new.erase(vtxState_new.begin()); // without first element
+    //loop over vtxState_sorted
+    for (Sp_Iter org_itr= vtxState_sorted.begin(); org_itr != vtxState_sorted.end(); org_itr++)
+    {
+       Trk::Vertex seed = (*org_itr)->vertex();
+        //found where the seed is already stored
+        bool vtx_is_stored = false;
+        unsigned int vertex_objectsPosition = 0; 
+        for (unsigned int vertex_objectsItr = 0; vertex_objectsItr!= vertex_objects.size(); vertex_objectsItr++)
+        {
+          //stored vertices in the line vertex_objectsItr
+          std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> stored_vertices = vertex_objects[vertex_objectsItr];
+           //---->inner loop
+          for (unsigned int stored_vtx_itr = 0; stored_vtx_itr!= stored_vertices.size(); stored_vtx_itr++)
+          {
+              Trk::Vertex storedVtx = stored_vertices[stored_vtx_itr]->vertex();
+              double diff = ((*org_itr)->vertex().position().x() - storedVtx.position().x()) + ((*org_itr)->vertex().position().y() - storedVtx.position().y()) + ((*org_itr)->vertex().position().z() - storedVtx.position().z());
+              if (fabs(diff) < 0.0001)
+              {
+                 vertex_objectsPosition = vertex_objectsItr;
+                 vtx_is_stored = true;
+                 break; // break inner loop if found
+              }
+            }
+            if (vtx_is_stored == true) break; // break outer loop if found
+         }
+         if (!vtx_is_stored) {
+                 //if not stored
+                std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> new_vtx_coll;
+                new_vtx_coll.push_back(*org_itr);
+                vertex_objects.push_back(new_vtx_coll);
+                vertex_objectsPosition = vertex_objects.size() - 1;
+         }
+         for (Sp_Iter new_itr = vtxState_new.begin(); new_itr!= vtxState_new.end(); new_itr++)
+         {
+             Trk::Vertex cand = (*new_itr)->vertex();
+             //calculate the vtx prob function for this seed
+             std::multimap<double, Trk::Vertex > vtx_prob_coll;
+             vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*new_itr)->vtx_prob(),(*new_itr)->vertex()));
+             vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type((*org_itr)->vtx_prob(),(*org_itr)->vertex()));
+             for (double alpha = m_resolvingStep; alpha <1.; alpha += m_resolvingStep)
+             {
+                 Trk::Vertex SP_line = Trk::Vertex((*org_itr)->vertex().position() + alpha*((*new_itr)->vertex().position() - (*org_itr)->vertex().position()));
+                  double sum_TrkProbTube = 0.;
+                  double sum2_TrkProbTube = 0.;
+                  for (TrackDataVecIter trk_itr = (*trackTES).begin(); trk_itr != (*trackTES).end(); trk_itr++)
+                  {
+	              double TrkProbTube = 0.;
+                      TrkProbTube = m_iTrkProbTubeCalc->calcProbTube(*(*trk_itr),SP_line);
+                       sum_TrkProbTube += TrkProbTube;
+                       sum2_TrkProbTube += pow(TrkProbTube,2.);
+                   }
+                   double BeamProbTube = m_iTrkProbTubeCalc->calcProbTube(beam_spot,SP_line);
+                   sum_TrkProbTube += BeamProbTube;
+                   sum2_TrkProbTube += pow(BeamProbTube,2.);
+                   double vtx_prob = m_iVtxProbCalc->calcVtxProb(sum_TrkProbTube, sum2_TrkProbTube);
+                    vtx_prob_coll.insert(std::multimap<double, Trk::Vertex >::value_type(vtx_prob, SP_line));
+               }
+              double vtx_prob_ratio = (*vtx_prob_coll.begin()).first/(*new_itr)->vtx_prob();
+                if (vtx_prob_ratio >= m_resolvingCut)  {
+                       //check if the vertices can be resolved
+                       vertex_objects[vertex_objectsPosition].push_back(*new_itr);
+                       vtxState_new.erase(new_itr);
+                       if (new_itr != vtxState_new.end()) --new_itr;
+                       if (new_itr == vtxState_new.end()) break;
+                } 
+             }
+      }//loop over orig reduced coll
+   }//if spatial point collection size > 0
+
+
+   //--------------------------------------------------------------------------------//
+  //--------------------step3: call vertex fitter----------------------------------//
+  //--------------------------------------------------------------------------------//
+
+    if (msgLvl(MSG::DEBUG)) msg() << " step THREE, vertex fit, vertex_objects size = "<<vertex_objects.size()<< endreq;
+    //std::vector<Trk::VxCandidate* > theVxContainer; --David S.
+    std::vector< xAOD::Vertex* > theVertexContainer;
+    //prepare trk coll --> remove double tracks
+    std::vector<std::vector<InDet::ZVTOP_TrkPartBaseVertexState*> >::iterator vtx_itr = vertex_objects.begin();
+    for (; vtx_itr!= vertex_objects.end(); vtx_itr++)
+    {
+      bool with_beam_spot = false;
+      std::vector <const Trk::TrackParticleBase*> trk_coll(0); //--->trk_coll for the fit
+      std::vector<InDet::ZVTOP_TrkPartBaseVertexState*>::iterator itr = (*vtx_itr).begin();
+      trk_coll.push_back((*itr)->tracks()[0]);
+      for (; itr != (*vtx_itr).end(); itr++)
+      { 
+        for ( std::vector <const Trk::TrackParticleBase*>::iterator vs_itr = ((*itr)->tracks()).begin(); vs_itr!= ((*itr)->tracks()).end(); vs_itr++)
+        {
+          bool found = false;
+	  for (std::vector <const Trk::TrackParticleBase*>::iterator trk_itr = trk_coll.begin();  trk_itr!= trk_coll.end(); trk_itr++)  {
+             if (*vs_itr == *trk_itr)  found = true;
+          }
+          if (!found)  trk_coll.push_back(*vs_itr);
+        }
+       if ((*itr)->beam_spot()) with_beam_spot = true;
+       }
+       //call the fitter
+       //Trk::VxCandidate * myVxCandidate(0); --David S.
+       xAOD::Vertex * myxAODVertex(0);
+       //const Amg::Vector3D p(0.,0.,0.); --David S.
+       const Amg::Vector3D startingPoint(0.,0.,0.);
+       //Trk::Vertex startingPoint(p); --David S.
+       //if (with_beam_spot) myVxCandidate = m_iVertexFitter->fit(trk_coll,beam_spot); --David S.
+       if (with_beam_spot) myxAODVertex = m_iVertexFitter->fit(trk_coll,theconstraint);
+       //else myVxCandidate = m_iVertexFitter->fit(trk_coll,startingPoint); --David S.
+       else myxAODVertex = m_iVertexFitter->fit(trk_coll,startingPoint);
+       bool bad_chi2 = true;
+       //if (myVxCandidate) { --David S.
+       if (myxAODVertex) {
+          while (bad_chi2)
+          {
+             //std::vector< Trk::VxTrackAtVertex*> trkAtVtx = *(myVxCandidate->vxTrackAtVertex()); --David S.
+             std::vector< Trk::VxTrackAtVertex > trkAtVtx = myxAODVertex->vxTrackAtVertex();
+             //std::vector< Trk::VxTrackAtVertex*>::iterator trkAtVtx_Iter = trkAtVtx.begin(); --David S.
+             std::vector< Trk::VxTrackAtVertex >::iterator trkAtVtx_Iter = trkAtVtx.begin();
+             std::vector< const Trk::TrackParticleBase*>::iterator trk_Iter = trk_coll.begin();
+             double largest_chi2 = 0.;
+             std::vector< const Trk::TrackParticleBase*>::iterator index;
+             for (; trkAtVtx_Iter!= trkAtVtx.end(); trkAtVtx_Iter++)
+             {
+               double chi2 = (*trkAtVtx_Iter).trackQuality().chiSquared();
+
+               if (chi2 > largest_chi2) {
+                   largest_chi2 = chi2;
+                   index = trk_Iter;
+               }
+               trk_Iter++;
+             }
+             if (largest_chi2 > m_trk_chi2_cut)
+             {
+               if (trk_coll.size() < 3) break;
+               else {
+                  trk_coll.erase(index);
+                  if (trk_coll.size() >= 2) {
+                    //if (myVxCandidate!=0) { delete myVxCandidate; myVxCandidate=0; } --David S.
+                    if (myxAODVertex!=0) { delete myxAODVertex; myxAODVertex=0; }
+                    //if (with_beam_spot) myVxCandidate = m_iVertexFitter->fit(trk_coll, beam_spot); --David S.
+                    if (with_beam_spot) myxAODVertex = m_iVertexFitter->fit(trk_coll, theconstraint);
+                    //else myVxCandidate = m_iVertexFitter->fit(trk_coll,startingPoint); --David S.
+                    else myxAODVertex = m_iVertexFitter->fit(trk_coll,startingPoint);
+                  }
+                  //if (myVxCandidate == 0) break; --David S.
+                  if (myxAODVertex == 0) break;
+               }
+             } else bad_chi2 = false;
+           }
+      }
+      //if (myVxCandidate && bad_chi2 == false && with_beam_spot) newVxContainer->push_back(myVxCandidate); --David S.
+      if (myxAODVertex && bad_chi2 == false && with_beam_spot) newVertexContainer->push_back(myxAODVertex);
+      //if (myVxCandidate && bad_chi2 == false && !with_beam_spot) theVxContainer.push_back(myVxCandidate); --David S.
+      if (myxAODVertex && bad_chi2 == false && !with_beam_spot) theVertexContainer.push_back(myxAODVertex);
+    }
+    //if (msgLvl(MSG::DEBUG)) msg() <<"vertex container size = "<<theVxContainer.size()<<endreq; --David S.
+    if (msgLvl(MSG::DEBUG)) msg() << "vertex container size = " << theVertexContainer.size() << endreq;
+
+   //ambiguity solving for two track vertices
+   //typedef std::vector<Trk::VxCandidate* >::iterator theVxContainerIter; --David S.
+   typedef std::vector< xAOD::Vertex* >::iterator theVertexContainerIter;
+   //theVxContainerIter iter = theVxContainer.begin(); --David S.
+   theVertexContainerIter iter = theVertexContainer.begin();
+   //theVxContainerIter iter_end = theVxContainer.end(); --David S.
+   theVertexContainerIter iter_end = theVertexContainer.end();
+   //std::vector<Trk::VxCandidate*> twoTrkContainer(0); --David S.
+   std::vector< xAOD::Vertex* > twoTrkContainer(0);
+   for (; iter!= iter_end; iter++)
+   {
+     //Trk::VxCandidate* theVxCandidate(*iter); --David S.
+     xAOD::Vertex* thexAODVertex(*iter);
+     //if (theVxCandidate->vxTrackAtVertex()->size() == 2 ) twoTrkContainer.push_back(theVxCandidate); --David S.
+     if (thexAODVertex->vxTrackAtVertex().size() == 2 ) twoTrkContainer.push_back(thexAODVertex);
+     //else newVxContainer->push_back(theVxCandidate); --David S.
+     else newVertexContainer->push_back(thexAODVertex);
+   }
+   //std::vector<Trk::VxCandidate*>::iterator twoTrk_itr = twoTrkContainer.begin(); --David S.
+   std::vector< xAOD::Vertex* >::iterator twoTrk_itr = twoTrkContainer.begin();
+   for (; twoTrk_itr!= twoTrkContainer.end(); twoTrk_itr++)
+   {
+     // TODO: make sure that links are actually set in VxTrackAtVertices! --David S.
+     //Trk::ITrackLink* trklink1 = (*(*twoTrk_itr)->vxTrackAtVertex())[0]->trackOrParticleLink(); --David S.
+     Trk::ITrackLink* trklink1 = (*twoTrk_itr)->vxTrackAtVertex()[0].trackOrParticleLink();
+     Trk::LinkToTrackParticleBase* linkToTrack1 = dynamic_cast<Trk::LinkToTrackParticleBase*>(trklink1);
+     const Trk::TrackParticleBase* first_trk = linkToTrack1->cachedElement();
+     //double first_trk_chi2 = (*(*twoTrk_itr)->vxTrackAtVertex())[0]->trackQuality().chiSquared(); --David S.
+     double first_trk_chi2 = (*twoTrk_itr)->vxTrackAtVertex()[0].trackQuality().chiSquared();
+     //Trk::ITrackLink* trklink2 = (*(*twoTrk_itr)->vxTrackAtVertex())[1]->trackOrParticleLink(); --David S.
+     Trk::ITrackLink* trklink2 = (*twoTrk_itr)->vxTrackAtVertex()[1].trackOrParticleLink();
+     Trk::LinkToTrackParticleBase* linkToTrack2 = dynamic_cast<Trk::LinkToTrackParticleBase*>(trklink2);
+     const Trk::TrackParticleBase* second_trk = linkToTrack2->cachedElement();
+     //double second_trk_chi2 = (*(*twoTrk_itr)->vxTrackAtVertex())[1]->trackQuality().chiSquared(); --David S.
+     double second_trk_chi2 = (*twoTrk_itr)->vxTrackAtVertex()[1].trackQuality().chiSquared();
+     
+     //VxContainer::iterator Citer = newVxContainer->begin(); --David S.
+     xAOD::VertexContainer::iterator Citer = newVertexContainer->begin();
+     //VxContainer::iterator Citer_end = newVxContainer->end(); --David S.
+     xAOD::VertexContainer::iterator Citer_end = newVertexContainer->end();
+     bool first_found = false;
+     bool second_found = false;
+     ////loop over all VxCandidates --David S.
+     //loop over all Vertices
+     for (; Citer != Citer_end; Citer++)
+     {
+       //unsigned int size = (*Citer)->vxTrackAtVertex()->size(); --David S.
+       unsigned int size = (*Citer)->vxTrackAtVertex().size();
+       for (unsigned int counter = 0; counter < size; counter++)
+       {  
+         // TODO: make sure that links are actually set in VxTrackAtVertices! --David S.
+         //Trk::ITrackLink* trklink = (*(*Citer)->vxTrackAtVertex())[counter]->trackOrParticleLink(); --David S.
+         Trk::ITrackLink* trklink = (*Citer)->vxTrackAtVertex()[counter].trackOrParticleLink();
+         Trk::LinkToTrackParticleBase* linkToTrack = dynamic_cast<Trk::LinkToTrackParticleBase*>(trklink);
+         const Trk::TrackParticleBase* trk= linkToTrack->cachedElement();
+         //double trk_chi2 = (*(*Citer)->vxTrackAtVertex())[counter]->trackQuality().chiSquared(); --David S.
+         double trk_chi2 = (*Citer)->vxTrackAtVertex()[counter].trackQuality().chiSquared();
+         if (trk == first_trk && trk_chi2 < first_trk_chi2 ) first_found = true;
+         if (trk == second_trk && trk_chi2 < second_trk_chi2) second_found = true;
+        }
+     }
+     //if (first_found==false && second_found==false) newVxContainer->push_back(*twoTrk_itr); --David S.
+     if (first_found == false && second_found == false) newVertexContainer->push_back(*twoTrk_itr);
+     else delete *twoTrk_itr;
+   //}//end loop over two track Vx Candidates --David S.
+   }//end loop over two track Vertices
+   //if (msgLvl(MSG::DEBUG)) msg() <<"new VxContainer size = "<<newVxContainer->size()<<endreq; --David S.
+   if (msgLvl(MSG::DEBUG)) msg() << "new VertexContainer size = " << newVertexContainer->size() << endreq;
+   for (Sp_Iter iter = vtxState_org.begin(); iter != vtxState_org.end(); iter++) delete *iter;
+} else {
+    if (msgLvl(MSG::DEBUG)) msg() <<"No tracks in this event, provide to the next event" <<  endreq;
+  }
+
+//return newVxContainer; --David S.
+return std::make_pair(newVertexContainer, newVertexAuxContainer);
+}
+
+
+std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*>  InDet::ZVTOP_Tool::findVertex(const xAOD::TrackParticleContainer* trackParticles)
+{
+     if(msgLvl(MSG::DEBUG)){ 
+ 	msg(MSG::DEBUG) << "N="<<trackParticles->size()<<" xAOD::TrackParticles found" << endreq; 
+ 	msg(MSG::DEBUG) << "No ZVTOP_Tool implementation for xAOD::TrackParticle" << endreq; 
+     } 
+     xAOD::VertexContainer *xAODContainer(0);
+     xAOD::VertexAuxContainer *xAODAuxContainer(0);   
+     return std::make_pair(xAODContainer, xAODAuxContainer); 
+}
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkPartBaseVertexState.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkPartBaseVertexState.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ebd459a0bed68073bf78abec2ae765ff50c85638
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkPartBaseVertexState.cxx
@@ -0,0 +1,30 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+     ZVTOP_TrkPartBaseVertexState.cxx  -  Description
+                             -------------------
+    authors : Tatjana Lenz <tatjana.lenz@cern.ch>
+ ***************************************************************************/
+#include "InDetZVTOPVxFinder/ZVTOP_TrkPartBaseVertexState.h"
+#include "VxVertex/Vertex.h"
+#include "TrkParticleBase/TrackParticleBase.h"
+
+namespace InDet {
+
+  ZVTOP_TrkPartBaseVertexState::ZVTOP_TrkPartBaseVertexState():m_vtx_prob(),m_vertex(),m_beam_spot(), m_tracks(0){}
+
+  ZVTOP_TrkPartBaseVertexState::ZVTOP_TrkPartBaseVertexState(double& vtx_prob, Trk::Vertex& vertex, bool& beam_spot, std::vector<const Trk::TrackParticleBase*>& tracks): m_vtx_prob(vtx_prob), m_vertex(vertex), m_beam_spot(beam_spot), m_tracks(tracks){}
+
+  ZVTOP_TrkPartBaseVertexState::ZVTOP_TrkPartBaseVertexState(const  ZVTOP_TrkPartBaseVertexState& vs):
+  m_vtx_prob(vs.m_vtx_prob),
+  m_vertex(vs.m_vertex),
+  m_beam_spot(vs.m_beam_spot),
+  m_tracks(std::vector<const Trk::TrackParticleBase*>())
+  {
+         for (std::vector<const Trk::TrackParticleBase*>::const_iterator itr = vs.m_tracks.begin();
+         itr != vs.m_tracks.end(); ++itr) m_tracks.push_back(*itr); 
+  }
+  ZVTOP_TrkPartBaseVertexState::~ZVTOP_TrkPartBaseVertexState() {}
+} // end of namespace
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkPartVertexState.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkPartVertexState.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..dc606440fee7e35be433fad40ab55709e18a8e58
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkPartVertexState.cxx
@@ -0,0 +1,31 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+                           ZVTOP_TrkPartVertexState.cxx  -  Description
+                             -------------------
+    authors : Tatjana Lenz <tatjana.lenz@cern.ch>
+ ***************************************************************************/
+
+#include "InDetZVTOPVxFinder/ZVTOP_TrkPartVertexState.h"
+#include "VxVertex/Vertex.h"
+#include "Particle/TrackParticle.h"
+
+namespace InDet {
+
+  ZVTOP_TrkPartVertexState::ZVTOP_TrkPartVertexState():m_vtx_prob(),m_vertex(),m_IP(), m_tracks(0){}
+
+  ZVTOP_TrkPartVertexState::ZVTOP_TrkPartVertexState(double& vtx_prob, Trk::Vertex& vertex, bool& IP, std::vector<const Rec::TrackParticle*>& tracks): m_vtx_prob(vtx_prob), m_vertex(vertex), m_IP(IP), m_tracks(tracks){}
+
+  ZVTOP_TrkPartVertexState::ZVTOP_TrkPartVertexState(const  ZVTOP_TrkPartVertexState& vs):
+  m_vtx_prob(vs.m_vtx_prob),
+  m_vertex(vs.m_vertex),
+  m_IP(vs.m_IP),
+  m_tracks(std::vector<const Rec::TrackParticle*>())
+  {
+         for (std::vector<const Rec::TrackParticle*>::const_iterator itr = vs.m_tracks.begin();
+         itr != vs.m_tracks.end(); ++itr) m_tracks.push_back(*itr); 
+  }
+  ZVTOP_TrkPartVertexState::~ZVTOP_TrkPartVertexState() {}
+} // end of namespace
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkProbTubeCalc.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkProbTubeCalc.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..cbf6aa0f19da9037bd3559de5709752f9d46439e
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_TrkProbTubeCalc.cxx
@@ -0,0 +1,132 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ZVTOP_TrkProbTubeCalc.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "InDetZVTOPVxFinder/ZVTOP_TrkProbTubeCalc.h"
+#include "TrkSurfaces/PerigeeSurface.h"
+#include "TrkParameters/TrackParameters.h"
+#include "VxVertex/Vertex.h"
+#include "CLHEP/Matrix/SymMatrix.h"
+#include "CLHEP/Matrix/Vector.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include "TrkExInterfaces/IExtrapolator.h"
+
+//================ Constructor =================================================
+
+InDet::ZVTOP_TrkProbTubeCalc::ZVTOP_TrkProbTubeCalc(const std::string& t,
+			  const std::string& n,
+			  const IInterface*  p )
+  :
+  AthAlgTool(t,n,p)
+{
+  declareInterface<IZVTOP_TrkProbTubeCalc>(this);
+  //  template for property declaration
+  declareProperty("Extrapolator",		m_extrapolator);
+}
+
+//================ Destructor =================================================
+
+InDet::ZVTOP_TrkProbTubeCalc::~ZVTOP_TrkProbTubeCalc()
+{}
+
+
+//================ Initialisation =================================================
+
+StatusCode InDet::ZVTOP_TrkProbTubeCalc::initialize()
+{
+  
+  StatusCode sc = AlgTool::initialize();
+
+  if (sc.isFailure()) return sc;
+  
+  /* Get the right extrapolator tool from ToolSvc */
+  if ( m_extrapolator.retrieve().isFailure() ) {
+    msg (MSG::FATAL) << "Failed to retrieve tool " << m_extrapolator << endreq;
+    return StatusCode::FAILURE;
+  } else {
+    msg (MSG::INFO) << "Retrieved tool " << m_extrapolator << endreq;
+  }
+  msg (MSG::INFO) << "initialize() successful in " << name() << endreq;
+  return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================
+
+StatusCode InDet::ZVTOP_TrkProbTubeCalc::finalize()
+{
+  StatusCode sc = AlgTool::finalize();
+  return sc;
+}
+
+//============================================================================================
+double InDet::ZVTOP_TrkProbTubeCalc::calcProbTube(const Trk::Track& trk, Trk::Vertex& vec)
+{
+  double probTube = 0.;
+  //perigee surface 
+  const Trk::Perigee* trkPer(dynamic_cast<const Trk::Perigee*>(trk.perigeeParameters()));
+  probTube = calcProbTube(trkPer, vec);
+  return probTube;
+}
+
+//============================================================================================
+double InDet::ZVTOP_TrkProbTubeCalc::calcProbTube(const Trk::RecVertex& beam_spot, Trk::Vertex& vec)
+{
+  double probTube = 0.;
+  Amg::Vector3D diff; 
+  diff[0]= beam_spot.position()[0] - vec.position()[0];
+  diff[1]= beam_spot.position()[1] - vec.position()[1];
+  diff[2]= beam_spot.position()[2] - vec.position()[2];
+  AmgMatrix(3,3) beam_spot_weight = beam_spot.covariancePosition().inverse();
+  probTube = exp(-0.5*diff.transpose()*beam_spot_weight*diff);
+  return probTube;
+}
+//============================================================================================
+double InDet::ZVTOP_TrkProbTubeCalc::calcProbTube(const Rec::TrackParticle& trk, Trk::Vertex& vec)
+{
+  double probTube = 0.;
+  const Trk::Perigee* trkPer(trk.measuredPerigee());
+  probTube = calcProbTube(trkPer,vec);
+  return probTube;
+}
+
+double InDet::ZVTOP_TrkProbTubeCalc::calcProbTube(const Trk::TrackParticleBase& trk, Trk::Vertex& vec)
+{
+  double probTube = 0.;
+  //perigee surface 
+  const Trk::Perigee* trkPer = dynamic_cast<const Trk::Perigee*>(&(trk.definingParameters()));
+  probTube = calcProbTube(trkPer,vec);
+  return probTube;
+}
+
+double InDet::ZVTOP_TrkProbTubeCalc::calcProbTube(const Trk::Perigee* trkPer, Trk::Vertex& vec){
+
+  double probTube = 0.;
+  Amg::Vector3D lp =vec.position();
+  Trk::PerigeeSurface perigeeSurface(lp);
+  double f_phi0 = trkPer->parameters()[Trk::phi0];
+  double f_theta = trkPer->parameters()[Trk::theta];
+  double f_qOverP = trkPer->parameters()[Trk::qOverP];
+  Trk::Perigee * extrapolatedPerigee(const_cast<Trk::Perigee*>(dynamic_cast<const Trk::Perigee*>(m_extrapolator->extrapolateDirectly(*trkPer, perigeeSurface))));
+  if (extrapolatedPerigee) {
+    AmgVector(5) diff;
+    diff[0]= extrapolatedPerigee->parameters()[Trk::d0];
+    diff[1]= extrapolatedPerigee->parameters()[Trk::z0];
+    diff[2]= extrapolatedPerigee->parameters()[Trk::phi0] - f_phi0;
+    diff[3]= extrapolatedPerigee->parameters()[Trk::theta] - f_theta;
+    diff[4]= extrapolatedPerigee->parameters()[Trk::qOverP] - f_qOverP;
+    if (extrapolatedPerigee->covariance() != 0) {
+      AmgMatrix(5,5) exp_perigee_weight = (*extrapolatedPerigee->covariance()).inverse();
+      probTube = exp(-0.5*diff.transpose()*exp_perigee_weight*diff);
+    } else {
+      if (msgLvl(MSG::DEBUG)) msg()  << "extrapolateted perigee has NO information on the covariance matrix" << endreq;
+    }
+  } else {
+    if (msgLvl(MSG::DEBUG)) msg() << "Perigee was not extrapolated!" << endreq;
+  }
+  delete extrapolatedPerigee;
+  return probTube;
+}
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_VertexState.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_VertexState.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..353b6511d28e58e483d0289eaff2bc3f1aa8efe1
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/ZVTOP_VertexState.cxx
@@ -0,0 +1,31 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+                           ZVTOP_VertexState.cxx  -  Description
+                             -------------------
+    authors : Tatjana Lenz <tatjana.lenz@cern.ch>
+ ***************************************************************************/
+
+#include "InDetZVTOPVxFinder/ZVTOP_VertexState.h"
+#include "VxVertex/Vertex.h"
+#include "TrkTrack/Track.h"
+
+namespace InDet {
+
+  ZVTOP_VertexState::ZVTOP_VertexState():m_vtx_prob(),m_vertex(),m_beam_spot(), m_tracks(0){}
+
+  ZVTOP_VertexState::ZVTOP_VertexState(double& vtx_prob, Trk::Vertex& vertex, bool& beam_spot, std::vector<const Trk::Track*>& tracks): m_vtx_prob(vtx_prob), m_vertex(vertex), m_beam_spot(beam_spot), m_tracks(tracks){}
+
+  ZVTOP_VertexState::ZVTOP_VertexState(const  ZVTOP_VertexState& vs):
+  m_vtx_prob(vs.m_vtx_prob),
+  m_vertex(vs.m_vertex),
+  m_beam_spot(vs.m_beam_spot),
+  m_tracks(std::vector<const Trk::Track*>())
+  {
+         for (std::vector<const Trk::Track*>::const_iterator itr = vs.m_tracks.begin();
+         itr != vs.m_tracks.end(); ++itr) m_tracks.push_back(*itr); 
+  }
+  ZVTOP_VertexState::~ZVTOP_VertexState() {}
+} // end of namespace
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/components/InDetZVTOPVxFinder_entries.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/components/InDetZVTOPVxFinder_entries.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..ac1ecc1cd22e11088c74ff3b9d5f30e7a12f5b9d
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/components/InDetZVTOPVxFinder_entries.cxx
@@ -0,0 +1,28 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "InDetZVTOPVxFinder/ZVTOP_Tool.h"
+#include "InDetZVTOPVxFinder/ZVTOP_SecVtxTool.h"
+#include "InDetZVTOPVxFinder/ZVTOP_SpatialPointFinder.h"
+#include "InDetZVTOPVxFinder/ZVTOP_SlowSpatialPointFinder.h"
+#include "InDetZVTOPVxFinder/ZVTOP_TrkProbTubeCalc.h"
+#include "InDetZVTOPVxFinder/ZVTOP_SimpleVtxProbCalc.h"
+#include "InDetZVTOPVxFinder/ZVTOP_AmbiguitySolver.h"
+
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_Tool )
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_SecVtxTool )
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_SpatialPointFinder )
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_SlowSpatialPointFinder )
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_TrkProbTubeCalc )
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_SimpleVtxProbCalc )
+DECLARE_NAMESPACE_TOOL_FACTORY( InDet, ZVTOP_AmbiguitySolver )
+
+DECLARE_FACTORY_ENTRIES( InDetZVTOPVxFinder )
+{
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_Tool )
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_SecVtxTool )
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_SpatialPointFinder )
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_SlowSpatialPointFinder )
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_TrkProbTubeCalc )
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_SimpleVtxProbCalc )
+	DECLARE_NAMESPACE_TOOL( InDet, ZVTOP_AmbiguitySolver)
+}
+
diff --git a/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/components/InDetZVTOPVxFinder_load.cxx b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/components/InDetZVTOPVxFinder_load.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..cd14b7e1e729dafbc002c909880d36fde93dcd23
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetZVTOPVxFinder/src/components/InDetZVTOPVxFinder_load.cxx
@@ -0,0 +1,5 @@
+
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( InDetZVTOPVxFinder )
+