diff --git a/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py b/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py
index 73b62951b9b981b712098d4a90a8a78822296822..978d1b7ce229588db4fb88609a96acf6948f81d6 100755
--- a/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py
+++ b/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py
@@ -290,13 +290,6 @@ def PrUpgradeChecking(defTracks={},
         seq.Members += [trassociator]
         GaudiSequencer("MCLinksTrSeq").Members += [seq]
 
-    from Configurables import LHCb__Converters__RecVertex__v2__fromVectorLHCbRecVertex as FromVectorLHCbRecVertex
-    vertexConverter = FromVectorLHCbRecVertex("VertexConverter")
-    vertexConverter.InputVerticesName = "Rec/Vertex/Vector/Primary"
-    vertexConverter.InputTracksName = defTracksConverted["Velo"]["Location"]
-    vertexConverter.OutputVerticesName = "Rec/Vertex/Primary"
-    GaudiSequencer("MCLinksTrSeq").Members += [vertexConverter]
-
     from Configurables import PrTrackChecker, PrUTHitChecker
     from Configurables import LoKi__Hybrid__MCTool
     MCHybridFactory = LoKi__Hybrid__MCTool("MCHybridFactory")
@@ -468,9 +461,7 @@ def PrUpgradeChecking(defTracks={},
         from Configurables import PrimaryVertexChecker
         PVCheck = PrimaryVertexChecker("PVChecker")
         if not PVCheck.isPropertySet('inputVerticesName'):
-            PVCheck.inputVerticesName = "Rec/Vertex/Primary"
-        if not PVCheck.isPropertySet('matchByTracks'):
-            PVCheck.matchByTracks = False
+            PVCheck.inputVerticesName = "Rec/Vertex/Vector/Primary"
         if not PVCheck.isPropertySet('nTracksToBeRecble'):
             PVCheck.nTracksToBeRecble = 4
         if not PVCheck.isPropertySet('inputTracksName'):
diff --git a/Tr/PatChecker/CMakeLists.txt b/Tr/PatChecker/CMakeLists.txt
index 257b46b2ce2cd4ea38d93bf1265de44bafa9d9c0..812a343f44cc7f2654a69319d77864f6c6f984ed 100644
--- a/Tr/PatChecker/CMakeLists.txt
+++ b/Tr/PatChecker/CMakeLists.txt
@@ -13,7 +13,8 @@
 ################################################################################
 gaudi_subdir(PatChecker v3r16p1)
 
-gaudi_depends_on_subdirs(Event/DigiEvent
+gaudi_depends_on_subdirs(Det/VPDet
+	                 Event/DigiEvent
                          Event/LinkerEvent
                          Event/MCEvent
                          Event/PhysEvent
@@ -32,5 +33,5 @@ include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS})
 gaudi_add_module(PatChecker
                  src/*.cpp
                  INCLUDE_DIRS AIDA GSL Event/DigiEvent Kernel/MCInterfaces Tf/TfKernel
-                 LINK_LIBRARIES GSL LinkerEvent MCEvent PhysEvent GaudiAlgLib GaudiKernel)
+                 LINK_LIBRARIES GSL LinkerEvent MCEvent PhysEvent GaudiAlgLib GaudiKernel VPDetLib)
 
diff --git a/Tr/PatChecker/src/PrimaryVertexChecker.cpp b/Tr/PatChecker/src/PrimaryVertexChecker.cpp
index f15f5b97b344e0a0aa96120bec09240eddd0bc0d..0e3b53b9a19e6203cf94ef225960c7524821b728 100644
--- a/Tr/PatChecker/src/PrimaryVertexChecker.cpp
+++ b/Tr/PatChecker/src/PrimaryVertexChecker.cpp
@@ -8,217 +8,250 @@
 * granted to it by virtue of its status as an Intergovernmental Organization  *
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
-// Include files
-// from Gaudi
+
 #include "AIDA/IHistogram1D.h"
+
 #include "DetDesc/Condition.h"
+#include "DetDesc/ConditionAccessorHolder.h"
+
+#include "Event/MCHeader.h"
+#include "Event/MCParticle.h"
+#include "Event/MCProperty.h"
+#include "Event/MCTrackInfo.h"
+#include "Event/MCVertex.h"
+
 #include "Event/RecVertex.h"
+#include "Event/RecVertex_v2.h"
 #include "Event/State.h"
 #include "Event/Track.h"
+#include "Event/Track_v2.h"
+
+#include "GaudiAlg/Consumer.h"
+#include "GaudiAlg/GaudiTupleAlg.h"
 #include "GaudiAlg/Tuples.h"
+
 #include "GaudiKernel/PhysicalConstants.h"
 #include "GaudiKernel/SystemOfUnits.h"
 #include "GaudiUtils/HistoStats.h"
-#include "Linker/AllLinks.h"
-#include "Linker/LinkerWithKey.h"
-#include <Event/MCTrackInfo.h>
-#include <Linker/LinkedTo.h>
-// local
-#include "PrimaryVertexChecker.h"
+#include "Kernel/STLExtensions.h"
 
-bool sortmlt( MCPVInfo first, MCPVInfo second ) { return first.nRecTracks > second.nRecTracks; }
+#include <string>
+#include <tuple>
+#include <vector>
 
 //-----------------------------------------------------------------------------
 // Implementation file for class : PrimaryVertexChecker
 //-----------------------------------------------------------------------------
 
-// Declaration of the Algorithm Factory
-DECLARE_COMPONENT( PrimaryVertexChecker )
+namespace {
+
+  struct MCPVInfo {
+    LHCb::MCVertex*                pMCPV;             // pointer to MC PV
+    int                            nRecTracks;        // number of reconstructed tracks from this MCPV
+    int                            nRecBackTracks;    // number of reconstructed backward tracks
+    int                            indexRecPVInfo;    // index to reconstructed PVInfo (-1 if not reco)
+    int                            nCorrectTracks;    // correct tracks belonging to reconstructed PV
+    int                            multClosestMCPV;   // multiplicity of closest reconstructable MCPV
+    double                         distToClosestMCPV; // distance to closest reconstructible MCPV
+    bool                           decayCharm;        // type of mother particle
+    bool                           decayBeauty;
+    bool                           decayStrange;
+    std::vector<LHCb::MCParticle*> m_mcPartInMCPV;
+    std::vector<LHCb::Track*>      m_recTracksInMCPV;
+  };
+
+  struct RecPVInfo {
+  public:
+    int                               nTracks; // number of tracks
+    double                            chi2;
+    double                            nDoF;
+    int                               mother;
+    Gaudi::XYZPoint                   position;      // position
+    Gaudi::XYZPoint                   positionSigma; // position sigmas
+    int                               indexMCPVInfo; // index to MCPVInfo
+    const LHCb::Event::v2::RecVertex* pRECPV;
+  };
+
+  inline const std::string beamSpotCond = "/dd/Conditions/Online/Velo/MotionSystem";
+
+  struct Beamline_t {
+    double X = std::numeric_limits<double>::signaling_NaN();
+    double Y = std::numeric_limits<double>::signaling_NaN();
+    Beamline_t( Condition const& c )
+        : X{( c.param<double>( "ResolPosRC" ) + c.param<double>( "ResolPosLA" ) ) / 2}
+        , Y{c.param<double>( "ResolPosY" )} {}
+  };
+
+  bool sortmlt( MCPVInfo const& first, MCPVInfo const& second ) { return first.nRecTracks > second.nRecTracks; }
+
+  enum class recoAs {
+    all,
+    isolated,
+    close,
+    ntracks_low,
+    ntracks_high,
+    z_low,
+    z_middle,
+    z_high,
+    beauty,
+    charm,
+    strange,
+    other,
+    first,
+    second,
+    third,
+    fourth,
+    fifth
+  };
+
+  constexpr auto All = std::array{
+      recoAs::all,      recoAs::isolated, recoAs::close,  recoAs::ntracks_low, recoAs::ntracks_high, recoAs::z_low,
+      recoAs::z_middle, recoAs::z_high,   recoAs::beauty, recoAs::charm,       recoAs::strange,      recoAs::other,
+      recoAs::first,    recoAs::second,   recoAs::third,  recoAs::fourth,      recoAs::fifth};
+  constexpr auto Part  = std::array{recoAs::all, recoAs::beauty, recoAs::charm, recoAs::strange, recoAs::other};
+  constexpr auto Basic = std::array{recoAs::all,          recoAs::isolated, recoAs::close,    recoAs::ntracks_low,
+                                    recoAs::ntracks_high, recoAs::z_low,    recoAs::z_middle, recoAs::z_high,
+                                    recoAs::beauty,       recoAs::charm,    recoAs::strange,  recoAs::other};
+
+  constexpr int size_basic  = static_cast<int>( recoAs::other ) + 1;
+  constexpr int size_recoAs = static_cast<int>( recoAs::fifth ) + 1;
+  constexpr int size_multi  = size_recoAs - size_basic;
+  constexpr int begin_multi = static_cast<int>( recoAs::first );
+
+} // namespace
+
+class PrimaryVertexChecker
+    : public Gaudi::Functional::Consumer<void( const LHCb::Tracks&, LHCb::Event::v2::RecVertices const&,
+                                               const LHCb::MCVertices&, const LHCb::MCParticles&, const LHCb::MCHeader&,
+                                               const LHCb::MCProperty&, const Beamline_t& ),
+                                         LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, Beamline_t>> {
+public:
+  /// Standard constructor
+  PrimaryVertexChecker( const std::string& name, ISvcLocator* pSvcLocator )
+      : Consumer{name,
+                 pSvcLocator,
+                 {
+                     KeyValue{"inputTracksName", LHCb::TrackLocation::Default},
+                     KeyValue{"inputVerticesName", LHCb::Event::v2::RecVertexLocation::Primary},
+                     KeyValue{"MCVertexInput", LHCb::MCVertexLocation::Default},
+                     KeyValue{"MCParticleInput", LHCb::MCParticleLocation::Default},
+                     KeyValue{"MCHeaderLocation", LHCb::MCHeaderLocation::Default},
+                     KeyValue{"MCPropertyInput", LHCb::MCPropertyLocation::TrackInfo},
+                     KeyValue{"BeamSpotLocation", "AlgorithmSpecific-" + name + "-beamspot"},
+                 }} {}
+
+  StatusCode initialize() override; ///< Algorithm initialization
+  void       operator()( const LHCb::Tracks& tracks, LHCb::Event::v2::RecVertices const& vertices,
+                   const LHCb::MCVertices& mcvtx, const LHCb::MCParticles& mcps, const LHCb::MCHeader& mcheader,
+                   const LHCb::MCProperty& flags,
+                   const Beamline_t& ) const override; ///< Algorithm execution
+  StatusCode finalize() override;                            ///< Algorithm finalization
+
+private:
+  // bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
+
+  Gaudi::Property<int>    m_nTracksToBeRecble{this, "nTracksToBeRecble", 4};  // min number of tracks in PV
+  Gaudi::Property<bool>   m_produceHistogram{this, "produceHistogram", true}; // producing histograms (light version)
+  Gaudi::Property<bool>   m_produceNtuple{this, "produceNtuple", false};      // producing NTuples (full version)
+  Gaudi::Property<bool>   m_requireVelo{this, "RequireVelo", true};           // requiring VELO for tracks
+  Gaudi::Property<double> m_dzIsolated{this, "dzIsolated", 10.0 * Gaudi::Units::mm}; // split close/isolated PVs
+  Gaudi::Property<int>    m_nTracksToPrint{this, "nTracksToPrint", 10};              // split low/high multiplicity PVs
+  Gaudi::Property<double> m_zToPrint{this, "zToPrint", 50.0 * Gaudi::Units::mm};     // split in z
+
+  mutable Gaudi::Accumulators::Counter<>                            m_nevt{this, "nEvents"};
+  mutable std::map<recoAs, Gaudi::Accumulators::BinomialCounter<>>  m_false;     // False PVs vs Reco PVs
+  mutable std::map<recoAs, Gaudi::Accumulators::BinomialCounter<>>  m_eff;       // MC reconstructible vs Reconstructed
+  mutable std::map<recoAs, Gaudi::Accumulators::BinomialCounter<>>  m_mcpv;      // MC vs MC reconstrucible
+  mutable std::map<recoAs, Gaudi::Accumulators::AveragingCounter<>> m_av_mcp;    // average mc particles in MCPV
+  mutable std::map<recoAs, Gaudi::Accumulators::AveragingCounter<>> m_av_tracks; // average tracks in RecoPV
+
+  std::vector<MCPVInfo>::iterator closestMCPV( std::vector<MCPVInfo>& rblemcpv, const MCPVInfo& mc ) const;
+  // std::vector<MCPVInfo>::iterator& itmc ) const;
+  std::vector<MCPVInfo>::iterator closestMCPV( std::vector<MCPVInfo>& rblemcpv, const RecPVInfo& rec ) const;
+
+  void collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx, std::vector<LHCb::MCParticle*>& allprods ) const;
+  void printRat( const std::string name, const Gaudi::Accumulators::BinomialCounter<>& m_mcpv,
+                 const Gaudi::Accumulators::BinomialCounter<>& m_eff,
+                 const Gaudi::Accumulators::BinomialCounter<>& m_false );
+  void printAvTracks( const std::string name, const Gaudi::Accumulators::AveragingCounter<>& m_av_tracks,
+                      const Gaudi::Accumulators::AveragingCounter<>& m_av_mcp );
+  void printRes( std::string mes, double x, double y, double z );
+  std::string toString( const recoAs& n ) const;
+  bool        checkCondition( const MCPVInfo& MCPV, const recoAs& n ) const;
+
+  void   match_mc_vertex_by_distance( int ipv, std::vector<RecPVInfo>& rinfo, std::vector<MCPVInfo>& mcpvvec ) const;
+  int    count_velo_tracks( const std::vector<LHCb::Event::v2::WeightedTrack>& tracksPV );
+  void   count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec, const LHCb::MCProperty& flags ) const;
+  double check_histogram( const AIDA::IHistogram1D* h, bool var );
+};
 
-//=============================================================================
-// Standard constructor, initializes variables
-//=============================================================================
-PrimaryVertexChecker::PrimaryVertexChecker( const std::string& name, ISvcLocator* pSvcLocator )
-    : GaudiTupleAlg( name, pSvcLocator ) {
-  declareProperty( "nTracksToBeRecble", m_nTracksToBeRecble = 5 );
-  declareProperty( "produceNtuple", m_produceNtuple = false );
-  declareProperty( "RequireVelo", m_requireVelo = true );
-  declareProperty( "produceHistogram", m_produceHistogram = true );
-  declareProperty( "dzIsolated", m_dzIsolated = 10.0 * Gaudi::Units::mm );
-  declareProperty( "matchByTracks", m_matchByTracks = true );
-  declareProperty( "inputTracksName", m_inputTracksName = LHCb::TrackLocation::Default );
-  declareProperty( "inputVerticesName", m_inputVerticesName = LHCb::RecVertexLocation::Primary );
-}
+// Declaration of the Algorithm Factory
+DECLARE_COMPONENT_WITH_ID( PrimaryVertexChecker, "PrimaryVertexChecker" )
 
-//=============================================================================
-// Initialization
-//=============================================================================
 StatusCode PrimaryVertexChecker::initialize() {
-  StatusCode sc = GaudiTupleAlg::initialize(); // Must be executed first
-  if ( sc.isFailure() ) debug() << "==> Initialize" << endmsg;
-
-  m_nevt                 = 0;
-  m_nMCPV                = 0;
-  m_nRecMCPV             = 0;
-  m_nMCPV_isol           = 0;
-  m_nRecMCPV_isol        = 0;
-  m_nMCPV_close          = 0;
-  m_nRecMCPV_close       = 0;
-  m_nFalsePV             = 0;
-  m_nFalsePV_real        = 0;
-  m_nMCPV_1mult          = 0;
-  m_nRecMCPV_1mult       = 0;
-  m_nMCPV_isol_1mult     = 0;
-  m_nRecMCPV_isol_1mult  = 0;
-  m_nMCPV_close_1mult    = 0;
-  m_nRecMCPV_close_1mult = 0;
-  m_nMCPV_2mult          = 0;
-  m_nRecMCPV_2mult       = 0;
-  m_nMCPV_isol_2mult     = 0;
-  m_nRecMCPV_isol_2mult  = 0;
-  m_nMCPV_close_2mult    = 0;
-  m_nRecMCPV_close_2mult = 0;
-  m_nMCPV_3mult          = 0;
-  m_nRecMCPV_3mult       = 0;
-  m_nMCPV_isol_3mult     = 0;
-  m_nRecMCPV_isol_3mult  = 0;
-  m_nMCPV_close_3mult    = 0;
-  m_nRecMCPV_close_3mult = 0;
-  m_nBFalse              = 0;
-  m_nRecBFalse           = 0;
-
-  return StatusCode::SUCCESS;
+  return Consumer::initialize().andThen(
+      [&] { addConditionDerivation<Beamline_t( Condition const& )>( {beamSpotCond}, inputLocation<Beamline_t>() ); } );
 }
 
 //=============================================================================
 // Main execution
 //=============================================================================
-StatusCode PrimaryVertexChecker::execute() {
-  debug() << "==> Execute" << endmsg;
+void PrimaryVertexChecker::operator()( const LHCb::Tracks& tracks, LHCb::Event::v2::RecVertices const& recoVtx,
+                                       const LHCb::MCVertices& mcvtx, const LHCb::MCParticles& mcps,
+                                       const LHCb::MCHeader& mcheader, const LHCb::MCProperty& flags,
+                                       const Beamline_t& beamline ) const {
 
+  debug() << "==> Execute" << endmsg;
   // Event
   m_nevt++;
 
-  // Get input
-  std::vector<LHCb::Track*>     vecOfTracks;
-  std::vector<LHCb::RecVertex*> vecOfVertices;
-  std::string                   trackLoc;
-  bool                          tracks_ok   = getInputTracks( vecOfTracks, trackLoc );
-  bool                          vertices_ok = getInputVertices( vecOfVertices );
-
-  if ( !tracks_ok && m_matchByTracks ) {
-    return StatusCode::SUCCESS; // return SUCCESSS anyway
-  }
-
-  if ( !vertices_ok ) { return StatusCode::SUCCESS; }
-
-  std::string m_beamSpotCond = "/dd/Conditions/Online/Velo/MotionSystem";
-  //  Condition *myCond =  get<Condition>(detSvc(), m_beamSpotCond );
-  const double xRC = 0.0; // myCond -> paramAsDouble ( "ResolPosRC" ) ;
-  const double xLA = 0.0; // myCond -> paramAsDouble ( "ResolPosLA" ) ;
-  const double Y   = 0.0; // myCond -> paramAsDouble ( "ResolPosY"  ) ;
-
-  double m_beamSpotX = ( xRC + xLA ) / 2;
-  double m_beamSpotY = Y;
-
-  if ( debugLevel() )
-    debug() << trackLoc << " # tracks: " << vecOfTracks.size() << "  # vertices: " << vecOfVertices.size() << endmsg;
+  if ( msgLevel( MSG::DEBUG ) )
+    debug() << " # tracks: " << tracks.size() << "  # vertices: " << recoVtx.size() << endmsg;
 
   // Fill reconstucted PV info
-  std::vector<RecPVInfo>                  recpvvec;
-  std::vector<LHCb::RecVertex*>::iterator itRecV;
-  for ( itRecV = vecOfVertices.begin(); vecOfVertices.end() != itRecV; itRecV++ ) {
-    LHCb::RecVertex* pv;
-    pv = *itRecV;
+  std::vector<RecPVInfo> recpvvec;
+  recpvvec.reserve( recoVtx.size() );
+
+  for ( const auto& pv : recoVtx ) {
     RecPVInfo recinfo;
-    recinfo.pRECPV            = pv;
-    recinfo.position          = pv->position();
-    Gaudi::SymMatrix3x3 covPV = pv->covMatrix();
+    recinfo.pRECPV            = ( &pv );
+    recinfo.position          = pv.position();
+    Gaudi::SymMatrix3x3 covPV = pv.covMatrix();
     double              sigx  = sqrt( covPV( 0, 0 ) );
     double              sigy  = sqrt( covPV( 1, 1 ) );
     double              sigz  = sqrt( covPV( 2, 2 ) );
     Gaudi::XYZPoint     a3d( sigx, sigy, sigz );
     recinfo.positionSigma = a3d;
-    recinfo.nTracks       = pv->tracks().size();
-    double minRD          = 99999.;
-    double maxRD          = -99999.;
-    double chi2           = pv->chi2();
-    double nDoF           = pv->nDoF();
-
-    std::string                               container = "container";
-    typedef const SmartRefVector<LHCb::Track> PVTRACKS;
-    PVTRACKS&                                 tracksIn = pv->tracks();
-    PVTRACKS::const_iterator                  itin;
-
-    int    mother    = 0;
-    int    velo      = 0;
-    int    lg        = 0;
-    double d0        = 0;
-    double mind0     = 99999.0;
-    double maxd0     = -99999.0;
-    double trackChi2 = 0.0;
-    int    tr        = 0;
-
-    for ( itin = tracksIn.begin(); itin != tracksIn.end(); itin++ ) {
-      tr++;
-      const LHCb::Track* ptr = *itin;
-
-      if ( ptr->type() == LHCb::Track::Types::Long ) { lg++; }
-      if ( ptr->type() == LHCb::Track::Types::Velo ) { velo++; }
-
-      double           x        = ptr->position().x();
-      double           y        = ptr->position().y();
-      double           r2       = x * x + y * y;
-      Gaudi::XYZVector udir     = ptr->firstState().slopes().Unit();
-      Gaudi::XYZVector distance = ptr->firstState().position() - pv->position();
-      Gaudi::XYZVector d0i      = udir.Cross( distance.Cross( udir ) );
-
-      d0 = d0 + d0i.Mag2();
-      if ( d0i.Mag2() > maxd0 ) { maxd0 = d0i.Mag2(); }
-      if ( d0i.Mag2() < mind0 ) { mind0 = d0i.Mag2(); }
-
-      if ( r2 > maxRD ) { maxRD = r2; }
-      if ( r2 < minRD ) { minRD = r2; }
-
-      double trChi2 = ptr->chi2();
-      trackChi2     = trackChi2 + trChi2;
-    }
-
-    if ( m_matchByTracks ) { mother = check_mother_particle( pv, trackLoc ); }
+    recinfo.nTracks       = pv.tracks().size();
+    recinfo.chi2          = pv.chi2();
+    recinfo.nDoF          = pv.nDoF();
+    recinfo.indexMCPVInfo = -1;
 
-    recinfo.minTrackRD    = minRD;
-    recinfo.maxTrackRD    = maxRD;
+    int mother            = 0;
     recinfo.mother        = mother;
-    recinfo.chi2          = chi2;
-    recinfo.nDoF          = nDoF;
-    recinfo.d0            = d0;
-    recinfo.d0nTr         = (double)d0 / (double)tr;
-    recinfo.chi2nTr       = (double)trackChi2 / (double)tr;
-    recinfo.mind0         = mind0;
-    recinfo.maxd0         = maxd0;
-    recinfo.nVeloTracks   = velo;
-    recinfo.nLongTracks   = lg;
     recinfo.indexMCPVInfo = -1;
     recpvvec.push_back( recinfo );
   }
 
   // Fill MC PV info
   std::vector<MCPVInfo> mcpvvec;
-  LHCb::MCVertices*     mcvtx = get<LHCb::MCVertices>( LHCb::MCVertexLocation::Default );
-  for ( LHCb::MCVertices::const_iterator itMCV = mcvtx->begin(); mcvtx->end() != itMCV; itMCV++ ) {
-    const LHCb::MCParticle* motherPart = ( *itMCV )->mother();
+  mcpvvec.reserve( mcvtx.size() );
+
+  for ( const auto& mcpv : mcvtx ) {
+    const LHCb::MCParticle* motherPart = mcpv->mother();
     if ( 0 == motherPart ) {
-      if ( ( *itMCV )->type() == LHCb::MCVertex::MCVertexType::ppCollision ) {
+      if ( mcpv->type() == LHCb::MCVertex::MCVertexType::ppCollision ) {
         MCPVInfo mcprimvert;
-        mcprimvert.pMCPV             = *itMCV;
+        mcprimvert.pMCPV             = mcpv;
         mcprimvert.nRecTracks        = 0;
         mcprimvert.nRecBackTracks    = 0;
         mcprimvert.indexRecPVInfo    = -1;
         mcprimvert.nCorrectTracks    = 0;
         mcprimvert.multClosestMCPV   = 0;
         mcprimvert.distToClosestMCPV = 999999.;
-        mcprimvert.decayBeauty       = 0;
-        mcprimvert.decayCharm        = 0;
+        mcprimvert.decayBeauty       = false;
+        mcprimvert.decayCharm        = false;
+        mcprimvert.decayStrange      = false;
         mcprimvert.m_mcPartInMCPV.clear();
         mcprimvert.m_recTracksInMCPV.clear();
         mcpvvec.push_back( mcprimvert );
@@ -226,356 +259,252 @@ StatusCode PrimaryVertexChecker::execute() {
     }
   }
 
-  const LHCb::MCParticle::Container* mcps = get<LHCb::MCParticle::Container>( LHCb::MCParticleLocation::Default );
-
-  int sMCP = 0;
-  int sTr  = vecOfTracks.size();
-  for ( LHCb::MCParticle::Container::const_iterator imc = mcps->begin(); mcps->end() != imc; ++imc ) {
-    const LHCb::MCParticle* dec = *imc;
-    double                  z   = dec->originVertex()->position().z();
-    if ( dec->particleID().threeCharge() != 0 && ( z < 400.0 && z > -600.0 ) ) { sMCP++; }
-  }
-
-  if ( m_matchByTracks ) {
-
-    count_reconstructed_tracks( mcpvvec, vecOfTracks, trackLoc );
-
-  } else {
-
-    count_reconstructible_mc_particles( mcpvvec );
-  }
+  count_reconstructible_mc_particles( mcpvvec, flags );
 
   std::sort( mcpvvec.begin(), mcpvvec.end(), sortmlt );
 
   std::vector<MCPVInfo> rblemcpv;
+  std::vector<MCPVInfo> notrblemcpv;
   std::vector<MCPVInfo> not_rble_but_visible;
   std::vector<MCPVInfo> not_rble;
   int                   nmrc = 0;
 
   std::vector<MCPVInfo>::iterator itmc;
   for ( itmc = mcpvvec.begin(); mcpvvec.end() != itmc; itmc++ ) {
-    rblemcpv.push_back( *itmc );
-    if ( itmc->nRecTracks < m_nTracksToBeRecble ) { nmrc++; }
+    if ( itmc->nRecTracks >= m_nTracksToBeRecble ) {
+      rblemcpv.push_back( *itmc );
+    } else {
+      notrblemcpv.push_back( *itmc );
+    }
     if ( itmc->nRecTracks < m_nTracksToBeRecble && itmc->nRecTracks > 1 ) { not_rble_but_visible.push_back( *itmc ); }
     if ( itmc->nRecTracks < m_nTracksToBeRecble && itmc->nRecTracks < 2 ) { not_rble.push_back( *itmc ); }
   }
 
-  // match MC and REC PVs
-  if ( m_matchByTracks ) {
-
-    for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) { match_mc_vertex_by_tracks( ipv, recpvvec, rblemcpv ); }
+  for ( int ipv = 0; ipv < static_cast<int>( recpvvec.size() ); ipv++ ) {
+    match_mc_vertex_by_distance( ipv, recpvvec, rblemcpv );
+  }
 
-  } else {
+  rblemcpv.insert( rblemcpv.end(), notrblemcpv.begin(), notrblemcpv.end() );
 
-    for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) { match_mc_vertex_by_distance( ipv, recpvvec, rblemcpv ); }
+  debug() << endmsg << " MC vertices " << endmsg;
+  debug() << " ===================================" << endmsg;
+  for ( const auto& [ipv, rmcpv] : LHCb::range::enumerate( rblemcpv ) ) {
+    std::string     ff   = " ";
+    LHCb::MCVertex* mcpv = rmcpv.pMCPV;
+    if ( rmcpv.indexRecPVInfo < 0 ) ff = "  NOTRECO";
+    debug() << format( " %3d %3d  xyz ( %7.4f %7.4f %8.3f )   nrec = %4d", ipv, rmcpv.indexRecPVInfo,
+                       mcpv->position().x(), mcpv->position().y(), mcpv->position().z(), rmcpv.nRecTracks )
+            << ff << endmsg;
   }
 
-  if ( debugLevel() ) {
-    debug() << endmsg << " MC vertices " << endmsg;
-    debug() << " ===================================" << endmsg;
-    for ( int imcpv = 0; imcpv < (int)rblemcpv.size(); imcpv++ ) {
-      std::string     ff   = " ";
-      LHCb::MCVertex* mcpv = rblemcpv[imcpv].pMCPV;
-      if ( rblemcpv[imcpv].indexRecPVInfo < 0 ) ff = "  NOTRECO";
-      debug() << format( " %3d %3d  xyz ( %7.4f %7.4f %8.3f )   nrec = %4d", imcpv, rblemcpv[imcpv].indexRecPVInfo,
-                         mcpv->position().x(), mcpv->position().y(), mcpv->position().z(), rblemcpv[imcpv].nRecTracks )
-              << ff << endmsg;
-    }
-    debug() << " -----------------------------------" << endmsg << endmsg;
-
-    debug() << endmsg << " REC vertices " << endmsg;
-    debug() << " ===================================" << endmsg;
-    for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) {
-      std::string ff = " ";
-      if ( recpvvec[ipv].indexMCPVInfo < 0 ) ff = "  FALSE  ";
-      debug()
-          << format(
-                 " %3d %3d  xyz ( %7.4f %7.4f %8.3f )  ntra = %4d   sigxyz ( %7.4f %7.4f %8.4f )   chi2/NDF = %7.2f",
-                 ipv, recpvvec[ipv].indexMCPVInfo, recpvvec[ipv].position.x(), recpvvec[ipv].position.y(),
-                 recpvvec[ipv].position.z(), recpvvec[ipv].nTracks, recpvvec[ipv].nVeloTracks,
-                 recpvvec[ipv].positionSigma.x(), recpvvec[ipv].positionSigma.y(), recpvvec[ipv].positionSigma.z(),
-                 recpvvec[ipv].pRECPV->chi2PerDoF() )
-          << ff << endmsg;
-    }
-    debug() << " -----------------------------------" << endmsg;
+  debug() << " -----------------------------------" << endmsg << endmsg;
+
+  debug() << endmsg << " REC vertices " << endmsg;
+  debug() << " ===================================" << endmsg;
+  for ( const auto& [ipv, recpv] : LHCb::range::enumerate( recpvvec ) ) {
+    std::string ff = " ";
+    if ( recpvvec[ipv].indexMCPVInfo < 0 ) ff = "  FALSE  ";
+    debug() << format(
+                   " %3d %3d  xyz ( %7.4f %7.4f %8.3f )  ntra = %4d   sigxyz ( %7.4f %7.4f %8.4f )   chi2/NDF = %7.2f",
+                   ipv, recpv.indexMCPVInfo, recpv.position.x(), recpv.position.y(), recpv.position.z(), recpv.nTracks,
+                   recpv.positionSigma.x(), recpv.positionSigma.y(), recpv.positionSigma.z(),
+                   recpv.pRECPV->chi2PerDoF() )
+            << ff << endmsg;
   }
+  debug() << " -----------------------------------" << endmsg;
 
   // find nr of false PV
-
-  int nFalsePV      = 0;
   int nFalsePV_real = 0;
-  for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) {
-    int    fake       = 0;
-    double x          = recpvvec[ipv].position.x();
-    double y          = recpvvec[ipv].position.y();
-    double z          = recpvvec[ipv].position.z();
-    double r          = std::sqrt( x * x + y * y );
-    double errx       = recpvvec[ipv].positionSigma.x();
-    double erry       = recpvvec[ipv].positionSigma.y();
-    double errz       = recpvvec[ipv].positionSigma.z();
-    double errr       = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) );
-    double minRDTrack = recpvvec[ipv].minTrackRD;
-    double maxRDTrack = recpvvec[ipv].maxTrackRD;
-    int    mother     = recpvvec[ipv].mother;
-    double velo       = recpvvec[ipv].nVeloTracks;
-    double lg         = recpvvec[ipv].nLongTracks;
-    double d0         = recpvvec[ipv].d0;
-    double d0nTr      = recpvvec[ipv].d0nTr;
-    double chi2nTr    = recpvvec[ipv].chi2nTr;
-    double mind0      = recpvvec[ipv].mind0;
-    double maxd0      = recpvvec[ipv].maxd0;
-    double chi2       = recpvvec[ipv].chi2;
-    double nDoF       = recpvvec[ipv].nDoF;
-
-    if ( recpvvec[ipv].indexMCPVInfo < 0 ) {
-      nFalsePV++;
-      fake           = 1;
+  for ( const auto& [ipv, recpv] : LHCb::range::enumerate( recpvvec ) ) {
+    int    fake   = 0;
+    double x      = recpv.position.x();
+    double y      = recpv.position.y();
+    double z      = recpv.position.z();
+    double r      = std::sqrt( x * x + y * y );
+    double errx   = recpv.positionSigma.x();
+    double erry   = recpv.positionSigma.y();
+    double errz   = recpv.positionSigma.z();
+    double errr   = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) );
+    int    mother = recpv.mother;
+    double chi2   = recpv.chi2;
+    double nDoF   = recpv.nDoF;
+
+    if ( recpv.indexMCPVInfo < 0 ) {
+      fake        = 1;
+      auto   cmc  = closestMCPV( rblemcpv, recpv );
+      double dist = 999999.;
+
+      if ( cmc != rblemcpv.end() ) {
+        dist = ( cmc->pMCPV->position() - recpv.pRECPV->position() ).R();
+
+        for ( const auto& n : Part ) m_false[n] += checkCondition( *cmc, n );
+
+        if ( dist > m_dzIsolated.value() ) {
+          m_false[recoAs::isolated] += true;
+        } else {
+          m_false[recoAs::close] += true;
+        }
+        if ( recpv.nTracks >= m_nTracksToPrint.value() ) {
+          m_false[recoAs::ntracks_high] += true;
+        } else {
+          m_false[recoAs::ntracks_low] += true;
+        }
+        if ( recpv.position.z() < -m_zToPrint.value() ) {
+          m_false[recoAs::z_low] += true;
+        } else if ( recpv.position.z() < m_zToPrint.value() ) {
+          m_false[recoAs::z_middle] += true;
+        } else {
+          m_false[recoAs::z_high] += true;
+        }
+
+        int idx = std::distance( rblemcpv.begin(), cmc );
+        if ( idx < size_multi ) { m_false[All[begin_multi + idx]] += true; }
+      }
+
       bool vis_found = false;
       for ( unsigned int imc = 0; imc < not_rble_but_visible.size(); imc++ ) {
         if ( not_rble_but_visible[imc].indexRecPVInfo > -1 ) continue;
-        double dist = fabs( mcpvvec[imc].pMCPV->position().z() - recpvvec[ipv].position.z() );
-        if ( dist < 5.0 * recpvvec[ipv].positionSigma.z() ) {
+        double dist = fabs( mcpvvec[imc].pMCPV->position().z() - recpv.position.z() );
+        if ( dist < 5.0 * recpv.positionSigma.z() ) {
           vis_found                                = true;
           not_rble_but_visible[imc].indexRecPVInfo = 10;
           break;
         }
       } // imc
+
       if ( !vis_found ) nFalsePV_real++;
     }
     if ( m_produceNtuple ) {
       Tuple myTuple2 = nTuple( 102, "PV_nTuple2", CLID_ColumnWiseTuple );
-      myTuple2->column( "fake", double( fake ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "r", double( r ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "x", double( x ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "y", double( y ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "z", double( z ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "errr", double( errr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "errz", double( errz ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "errx", double( errx ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "erry", double( erry ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "minRDTrack", double( minRDTrack ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "maxRDTrack", double( maxRDTrack ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "mother", double( mother ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "velo", double( velo ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "long", double( lg ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "d0", double( d0 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "d0nTr", double( d0nTr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "chi2nTr", double( chi2nTr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "mind0", double( mind0 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "maxd0", double( maxd0 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "chi2", double( chi2 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->column( "nDoF", double( nDoF ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple2->write().ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
+      myTuple2->column( "fake", double( fake ) ).ignore();
+      myTuple2->column( "r", double( r ) ).ignore();
+      myTuple2->column( "x", double( x ) ).ignore();
+      myTuple2->column( "y", double( y ) ).ignore();
+      myTuple2->column( "z", double( z ) ).ignore();
+      myTuple2->column( "errr", double( errr ) ).ignore();
+      myTuple2->column( "errz", double( errz ) ).ignore();
+      myTuple2->column( "errx", double( errx ) ).ignore();
+      myTuple2->column( "erry", double( erry ) ).ignore();
+      myTuple2->column( "mother", double( mother ) ).ignore();
+      myTuple2->column( "chi2", double( chi2 ) ).ignore();
+      myTuple2->column( "nDoF", double( nDoF ) ).ignore();
+      myTuple2->write().ignore();
     }
   }
   // Fill distance to closest recble MC PV and its multiplicity
-  std::vector<MCPVInfo>::iterator itmcl;
-  for ( itmcl = rblemcpv.begin(); rblemcpv.end() != itmcl; itmcl++ ) {
-    std::vector<MCPVInfo>::iterator cmc  = closestMCPV( rblemcpv, itmcl );
+  for ( auto& mcpv : rblemcpv ) {
+    std::vector<MCPVInfo>::iterator cmc  = closestMCPV( rblemcpv, mcpv );
     double                          dist = 999999.;
     int                             mult = 0;
     if ( cmc != rblemcpv.end() ) {
-      dist = ( cmc->pMCPV->position() - itmcl->pMCPV->position() ).R();
+      dist = ( cmc->pMCPV->position() - mcpv.pMCPV->position() ).R();
       mult = cmc->nRecTracks;
     }
-    itmcl->distToClosestMCPV = dist;
-    itmcl->multClosestMCPV   = mult;
-  }
-
-  // count non-reconstructible close and isolated PVs
-  int nmrc_isol  = 0;
-  int nmrc_close = 0;
-
-  // Counters
-  int nMCPV                = rblemcpv.size() - nmrc;
-  int nRecMCPV             = 0;
-  int nMCPV_isol           = 0;
-  int nRecMCPV_isol        = 0;
-  int nMCPV_close          = 0;
-  int nRecMCPV_close       = 0;
-  int nMCPV_1mult          = 0;
-  int nRecMCPV_1mult       = 0;
-  int nMCPV_isol_1mult     = 0;
-  int nRecMCPV_isol_1mult  = 0;
-  int nMCPV_close_1mult    = 0;
-  int nRecMCPV_close_1mult = 0;
-  int nRecMCPV_wrong_1mult = 0;
-  int nMCPV_2mult          = 0;
-  int nRecMCPV_2mult       = 0;
-  int nMCPV_isol_2mult     = 0;
-  int nRecMCPV_isol_2mult  = 0;
-  int nMCPV_close_2mult    = 0;
-  int nRecMCPV_close_2mult = 0;
-  int nRecMCPV_wrong_2mult = 0;
-  int nMCPV_3mult          = 0;
-  int nRecMCPV_3mult       = 0;
-  int nMCPV_isol_3mult     = 0;
-  int nRecMCPV_isol_3mult  = 0;
-  int nMCPV_close_3mult    = 0;
-  int nRecMCPV_close_3mult = 0;
-  int nRecMCPV_wrong_3mult = 0;
-
-  for ( itmc = rblemcpv.begin(); rblemcpv.end() != itmc; itmc++ ) {
-    if ( itmc->distToClosestMCPV > m_dzIsolated ) nMCPV_isol++;
-    if ( itmc->distToClosestMCPV > m_dzIsolated && itmc->nRecTracks < m_nTracksToBeRecble ) nmrc_isol++;
-    if ( itmc->distToClosestMCPV < m_dzIsolated ) nMCPV_close++;
-    if ( itmc->distToClosestMCPV < m_dzIsolated && itmc->nRecTracks < m_nTracksToBeRecble ) nmrc_close++;
-
-    if ( itmc->indexRecPVInfo > -1 ) {
-      nRecMCPV++;
-      if ( itmc->distToClosestMCPV > m_dzIsolated ) nRecMCPV_isol++;
-      if ( itmc->distToClosestMCPV < m_dzIsolated ) nRecMCPV_close++;
-    }
-  }
-
-  // rblemcpv is already sorted
-
-  // highest mult
-  if ( rblemcpv.size() > 0 ) {
-    nMCPV_1mult++;
-    double dist = rblemcpv[0].distToClosestMCPV;
-    if ( dist > m_dzIsolated ) nMCPV_isol_1mult++;
-    if ( dist < m_dzIsolated ) nMCPV_close_1mult++;
-    if ( rblemcpv[0].indexRecPVInfo > -1 ) {
-      nRecMCPV_1mult++;
-      if ( dist > m_dzIsolated ) nRecMCPV_isol_1mult++;
-      if ( dist < m_dzIsolated ) nRecMCPV_close_1mult++;
-    } else {
-      nRecMCPV_wrong_1mult++;
-    }
+    mcpv.distToClosestMCPV = dist;
+    mcpv.multClosestMCPV   = mult;
   }
 
-  // second highest
-  if ( rblemcpv.size() > 1 ) {
-    nMCPV_2mult++;
-    double dist = rblemcpv[1].distToClosestMCPV;
-    if ( dist > m_dzIsolated ) nMCPV_isol_2mult++;
-    if ( dist < m_dzIsolated ) nMCPV_close_2mult++;
-    if ( rblemcpv[1].indexRecPVInfo > -1 ) {
-      nRecMCPV_2mult++;
-      if ( dist > m_dzIsolated ) nRecMCPV_isol_2mult++;
-      if ( dist < m_dzIsolated ) nRecMCPV_close_2mult++;
-    } else {
-      nRecMCPV_wrong_2mult++;
-    }
-  }
-  // third highest
-  if ( rblemcpv.size() > 2 ) {
-    nMCPV_3mult++;
-    double dist = rblemcpv[2].distToClosestMCPV;
-    if ( dist > m_dzIsolated ) nMCPV_isol_3mult++;
-    if ( dist < m_dzIsolated ) nMCPV_close_3mult++;
-    if ( rblemcpv[2].indexRecPVInfo > -1 ) {
-      nRecMCPV_3mult++;
-      if ( dist > m_dzIsolated ) nRecMCPV_isol_3mult++;
-      if ( dist < m_dzIsolated ) nRecMCPV_close_3mult++;
-    } else {
-      nRecMCPV_wrong_3mult++;
+  for ( const auto& itmc : rblemcpv ) {
+    for ( const auto& n : Basic ) {
+      bool cut = false;
+      cut      = checkCondition( itmc, n );
+      if ( cut ) {
+        if ( itmc.nRecTracks < m_nTracksToBeRecble ) {
+          m_mcpv[n] += false;
+        } else {
+          m_mcpv[n] += true;
+          m_av_mcp[n] += itmc.nRecTracks;
+          if ( itmc.indexRecPVInfo < 0 ) {
+            m_eff[n] += false;
+          } else {
+            m_eff[n] += true;
+            m_false[n] += false;
+          }
+        }
+        if ( itmc.indexRecPVInfo > -1 ) { m_av_tracks[n] += recpvvec[itmc.indexRecPVInfo].nTracks; }
+      }
     }
   }
 
-  nMCPV_isol  = nMCPV_isol - nmrc_isol;
-  nMCPV_close = nMCPV_close - nmrc_close;
-
-  m_nMCPV += nMCPV;
-  m_nRecMCPV += nRecMCPV;
-  m_nMCPV_isol += nMCPV_isol;
-  m_nRecMCPV_isol += nRecMCPV_isol;
-  m_nMCPV_close += nMCPV_close;
-  m_nRecMCPV_close += nRecMCPV_close;
-  m_nFalsePV += nFalsePV;
-  m_nFalsePV_real += nFalsePV_real;
-  m_nMCPV_1mult += nMCPV_1mult;
-  m_nRecMCPV_1mult += nRecMCPV_1mult;
-  m_nMCPV_isol_1mult += nMCPV_isol_1mult;
-  m_nRecMCPV_isol_1mult += nRecMCPV_isol_1mult;
-  m_nMCPV_close_1mult += nMCPV_close_1mult;
-  m_nRecMCPV_close_1mult += nRecMCPV_close_1mult;
-  m_nRecMCPV_wrong_1mult += nRecMCPV_wrong_1mult;
-  m_nMCPV_2mult += nMCPV_2mult;
-  m_nRecMCPV_2mult += nRecMCPV_2mult;
-  m_nMCPV_isol_2mult += nMCPV_isol_2mult;
-  m_nRecMCPV_isol_2mult += nRecMCPV_isol_2mult;
-  m_nMCPV_close_2mult += nMCPV_close_2mult;
-  m_nRecMCPV_wrong_2mult += nRecMCPV_wrong_2mult;
-  m_nRecMCPV_close_2mult += nRecMCPV_close_2mult;
-  m_nMCPV_3mult += nMCPV_3mult;
-  m_nRecMCPV_3mult += nRecMCPV_3mult;
-  m_nMCPV_isol_3mult += nMCPV_isol_3mult;
-  m_nRecMCPV_isol_3mult += nRecMCPV_isol_3mult;
-  m_nMCPV_close_3mult += nMCPV_close_3mult;
-  m_nRecMCPV_close_3mult += nRecMCPV_close_3mult;
-  m_nRecMCPV_wrong_3mult += nRecMCPV_wrong_3mult;
-
+  int mcpv = 0;
   int high = 0;
 
-  for ( itmc = rblemcpv.begin(); rblemcpv.end() != itmc; itmc++ ) {
-    double x                 = -99999.;
-    double y                 = -99999.;
-    double dx                = -99999.;
-    double dy                = -99999.;
-    double dz                = -99999.;
-    double r                 = -99999.;
-    double zMC               = -99999.;
-    double yMC               = -99999.;
-    double xMC               = -99999.;
-    double rMC               = -99999.;
-    double z                 = -99999.;
-    double errx              = -99999.;
-    double erry              = -99999.;
-    double errz              = -99999.;
-    double errr              = -99999.;
-    double minRDTrack        = 99999.;
-    double maxRDTrack        = -99999.;
-    double chi2              = -999999.;
-    double nDoF              = -999999.;
-    int    indRec            = itmc->indexRecPVInfo;
-    int    reconstructed     = 0;
-    int    ntracks_pvrec     = 0;
-    int    nvelotracks_pvrec = 0;
-    int    ntracks_pvmc      = 0;
-    int    dtrcks            = 0;
-    int    pvevt             = 0;
-    int    mother            = 0;
-    int    assoctrks         = 0;
-    int    nassoctrks        = 0;
-
-    zMC = itmc->pMCPV->position().z();
-    yMC = itmc->pMCPV->position().y();
-    xMC = itmc->pMCPV->position().x();
-    rMC = std::sqrt( ( xMC - m_beamSpotX ) * ( xMC - m_beamSpotX ) + ( yMC - m_beamSpotY ) * ( yMC - m_beamSpotY ) );
+  for ( auto& itmc : rblemcpv ) {
+    double x             = -99999.;
+    double y             = -99999.;
+    double dx            = -99999.;
+    double dy            = -99999.;
+    double dz            = -99999.;
+    double r             = -99999.;
+    double zMC           = -99999.;
+    double yMC           = -99999.;
+    double xMC           = -99999.;
+    double rMC           = -99999.;
+    double z             = -99999.;
+    double errx          = -99999.;
+    double erry          = -99999.;
+    double errz          = -99999.;
+    double errr          = -99999.;
+    double chi2          = -999999.;
+    double nDoF          = -999999.;
+    int    indRec        = itmc.indexRecPVInfo;
+    int    reconstructed = 0;
+    int    ntracks_pvrec = 0;
+    int    ntracks_pvmc  = 0;
+    int    dtrcks        = 0;
+    int    pvevt         = 0;
+    int    mother        = 0;
+    int    assoctrks     = 0;
+    int    nassoctrks    = 0;
+
+    zMC = itmc.pMCPV->position().z();
+    yMC = itmc.pMCPV->position().y();
+    xMC = itmc.pMCPV->position().x();
+    rMC = std::sqrt( ( xMC - beamline.X ) * ( xMC - beamline.X ) + ( yMC - beamline.Y ) * ( yMC - beamline.Y ) );
+
+    if ( mcpv < size_multi ) {
+      if ( itmc.nRecTracks < m_nTracksToBeRecble ) {
+        m_mcpv[All[begin_multi + mcpv]] += false;
+      } else {
+        m_mcpv[All[begin_multi + mcpv]] += true;
+        m_av_mcp[All[begin_multi + mcpv]] += itmc.nRecTracks;
+        if ( itmc.indexRecPVInfo < 0 ) {
+          m_eff[All[begin_multi + mcpv]] += false;
+        } else {
+          m_eff[All[begin_multi + mcpv]] += true;
+          m_false[All[begin_multi + mcpv]] += false;
+        }
+      }
+    }
 
     if ( indRec > -1 ) {
       high++;
       pvevt++;
       reconstructed = 1;
-      dx            = recpvvec[indRec].position.x() - itmc->pMCPV->position().x();
-      dy            = recpvvec[indRec].position.y() - itmc->pMCPV->position().y();
-      dz            = recpvvec[indRec].position.z() - itmc->pMCPV->position().z();
+      dx            = recpvvec[indRec].position.x() - itmc.pMCPV->position().x();
+      dy            = recpvvec[indRec].position.y() - itmc.pMCPV->position().y();
+      dz            = recpvvec[indRec].position.z() - itmc.pMCPV->position().z();
       x             = recpvvec[indRec].position.x();
       y             = recpvvec[indRec].position.y();
       z             = recpvvec[indRec].position.z();
-      // zMC = itmc->pMCPV->position().z();
-      r    = std::sqrt( ( x - m_beamSpotX ) * ( x - m_beamSpotX ) + ( y - m_beamSpotY ) * ( y - m_beamSpotY ) );
-      errx = recpvvec[indRec].positionSigma.x();
-      erry = recpvvec[indRec].positionSigma.y();
-      errz = recpvvec[indRec].positionSigma.z();
-      errr = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) );
-      ntracks_pvrec     = recpvvec[indRec].nTracks;
-      nvelotracks_pvrec = recpvvec[indRec].nVeloTracks;
-      ntracks_pvmc      = itmc->pMCPV->products().size();
-      dtrcks            = ntracks_pvmc - ntracks_pvrec;
-      minRDTrack        = recpvvec[indRec].minTrackRD;
-      maxRDTrack        = recpvvec[indRec].maxTrackRD;
-      mother            = recpvvec[indRec].mother;
-      chi2              = recpvvec[indRec].chi2;
-      nDoF              = recpvvec[indRec].nDoF;
-
+      r             = std::sqrt( ( x - beamline.X ) * ( x - beamline.X ) + ( y - beamline.Y ) * ( y - beamline.Y ) );
+      errx          = recpvvec[indRec].positionSigma.x();
+      erry          = recpvvec[indRec].positionSigma.y();
+      errz          = recpvvec[indRec].positionSigma.z();
+      errr          = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) );
+      ntracks_pvrec = recpvvec[indRec].nTracks;
+      ntracks_pvmc  = itmc.pMCPV->products().size();
+      dtrcks        = ntracks_pvmc - ntracks_pvrec;
+      mother        = recpvvec[indRec].mother;
+      chi2          = recpvvec[indRec].chi2;
+      nDoF          = recpvvec[indRec].nDoF;
+
+      if ( mcpv < size_multi ) { m_av_tracks[All[begin_multi + mcpv]] += recpvvec[indRec].nTracks; }
       // Filling histograms
       if ( m_produceHistogram ) {
-        plot( itmc->pMCPV->position().x(), 1001, "xmc", -0.25, 0.25, 50 );
-        plot( itmc->pMCPV->position().y(), 1002, "ymc", -0.25, 0.25, 50 );
-        plot( itmc->pMCPV->position().z(), 1003, "zmc", -20, 20, 50 );
+        plot( itmc.pMCPV->position().x(), 1001, "xmc", -0.25, 0.25, 50 );
+        plot( itmc.pMCPV->position().y(), 1002, "ymc", -0.25, 0.25, 50 );
+        plot( itmc.pMCPV->position().z(), 1003, "zmc", -20, 20, 50 );
         plot( recpvvec[indRec].position.x(), 1011, "xrd", -0.25, 0.25, 50 );
         plot( recpvvec[indRec].position.y(), 1012, "yrd", -0.25, 0.25, 50 );
         plot( recpvvec[indRec].position.z(), 1013, "zrd", -20, 20, 50 );
@@ -585,119 +514,109 @@ StatusCode PrimaryVertexChecker::execute() {
         plot( dx / errx, 1031, "pullx", -5., 5., 50 );
         plot( dy / erry, 1032, "pully", -5., 5., 50 );
         plot( dz / errz, 1033, "pullz", -5., 5., 50 );
+        if ( itmc.nRecTracks < 10 ) {
+          plot( dx, 1101, "dx", -0.25, 0.25, 50 );
+          plot( dy, 1102, "dy", -0.25, 0.25, 50 );
+          plot( dz, 1103, "dz", -0.5, 0.5, 50 );
+        } else if ( itmc.nRecTracks >= 10 && itmc.nRecTracks < 30 ) {
+          plot( dx, 1111, "dx", -0.25, 0.25, 50 );
+          plot( dy, 1112, "dy", -0.25, 0.25, 50 );
+          plot( dz, 1113, "dz", -0.5, 0.5, 50 );
+        } else {
+          plot( dx, 1121, "dx", -0.25, 0.25, 50 );
+          plot( dy, 1122, "dy", -0.25, 0.25, 50 );
+          plot( dz, 1123, "dz", -0.5, 0.5, 50 );
+        }
+
+        if ( itmc.pMCPV->position().z() < -m_zToPrint ) {
+          plot( dx, 1201, "dx", -0.25, 0.25, 50 );
+          plot( dy, 1202, "dy", -0.25, 0.25, 50 );
+          plot( dz, 1203, "dz", -0.5, 0.5, 50 );
+        } else if ( itmc.pMCPV->position().z() < m_zToPrint ) {
+          plot( dx, 1211, "dx", -0.25, 0.25, 50 );
+          plot( dy, 1212, "dy", -0.25, 0.25, 50 );
+          plot( dz, 1213, "dz", -0.5, 0.5, 50 );
+        } else {
+          plot( dx, 1221, "dx", -0.25, 0.25, 50 );
+          plot( dy, 1222, "dy", -0.25, 0.25, 50 );
+          plot( dz, 1223, "dz", -0.5, 0.5, 50 );
+        }
+
         plot( double( ntracks_pvrec ), 1041, "ntracks_pvrec", 0., 150., 50 );
         plot( double( dtrcks ), 1042, "mcrdtracks", 0., 150., 50 );
-        plot( double( nvelotracks_pvrec ), 1043, "nvelotracks_pvrec", 0., 150., 50 );
         if ( pvevt == 1 ) {
           plot( double( recpvvec.size() ), 1051, "nPVperEvt", -0.5, 5.5, 6 );
           for ( int ipvrec = 0; ipvrec < (int)recpvvec.size(); ipvrec++ ) {
             assoctrks = assoctrks + recpvvec[ipvrec].nTracks;
           }
-          nassoctrks = vecOfTracks.size() - assoctrks;
+          nassoctrks = tracks.size() - assoctrks;
           plot( double( nassoctrks ), 1052, "nassoctrks", 0., 150., 50 );
         }
       }
     }
+    mcpv++;
 
     int    isolated     = 0;
-    double dist_closest = itmc->distToClosestMCPV;
-    if ( dist_closest > m_dzIsolated ) { isolated = 1; }
+    double dist_closest = itmc.distToClosestMCPV;
+    if ( dist_closest > m_dzIsolated.value() ) { isolated = 1; }
 
     // Filling ntuple
     if ( m_produceNtuple ) {
       Tuple myTuple = nTuple( 101, "PV_nTuple", CLID_ColumnWiseTuple );
-      myTuple->column( "reco", double( reconstructed ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "isol", double( isolated ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "ntracks", double( ntracks_pvrec ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "nvelotracks", double( nvelotracks_pvrec ) )
-          .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "nrectrmc", double( itmc->nRecTracks ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "dzclose", dist_closest ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "nmcpv", double( rblemcpv.size() ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "nmcallpv", double( mcpvvec.size() ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "nrecpv", double( recpvvec.size() ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "decayCharm", double( itmc->decayCharm ) )
-          .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "decayBeauty", double( itmc->decayBeauty ) )
-          .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "multi", double( high ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "dx", dx ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "dy", dy ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "dz", dz ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "x", x ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "y", y ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "r", r ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "zMC", zMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "yMC", yMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "xMC", xMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "rMC", rMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "z", z ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "xBeam", m_beamSpotX ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "yBeam", m_beamSpotY ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "errx", errx ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "erry", erry ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "errz", errz ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "errr", errr ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "minRDTrack", minRDTrack ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "maxRDTrack", maxRDTrack ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "mother", double( mother ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "evtnr", double( m_nevt ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "chi2", double( chi2 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "nDoF", double( nDoF ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "size_tracks", double( sTr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "size_mcp", double( sMCP ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->column( "mcpvrec", double( nmrc ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
-      myTuple->write().ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
+      myTuple->column( "reco", double( reconstructed ) ).ignore();
+      myTuple->column( "isol", double( isolated ) ).ignore();
+      myTuple->column( "ntracks", double( ntracks_pvrec ) ).ignore();
+      myTuple->column( "nrectrmc", double( itmc.nRecTracks ) ).ignore();
+      myTuple->column( "dzclose", dist_closest ).ignore();
+      myTuple->column( "nmcpv", double( rblemcpv.size() ) ).ignore();
+      myTuple->column( "mtruemcpv", double( mcheader.numOfPrimaryVertices() ) ).ignore();
+      myTuple->column( "nmcallpv", double( mcpvvec.size() ) ).ignore();
+      myTuple->column( "nrecpv", double( recpvvec.size() ) ).ignore();
+      myTuple->column( "decayCharm", double( itmc.decayCharm ) ).ignore();
+      myTuple->column( "decayBeauty", double( itmc.decayBeauty ) ).ignore();
+      myTuple->column( "decayStrange", double( itmc.decayStrange ) ).ignore();
+      myTuple->column( "multirec", double( high ) ).ignore();
+      myTuple->column( "multimc", double( mcpv ) ).ignore();
+      myTuple->column( "dx", dx ).ignore();
+      myTuple->column( "dy", dy ).ignore();
+      myTuple->column( "dz", dz ).ignore();
+      myTuple->column( "x", x ).ignore();
+      myTuple->column( "y", y ).ignore();
+      myTuple->column( "r", r ).ignore();
+      myTuple->column( "zMC", zMC ).ignore();
+      myTuple->column( "yMC", yMC ).ignore();
+      myTuple->column( "xMC", xMC ).ignore();
+      myTuple->column( "rMC", rMC ).ignore();
+      myTuple->column( "z", z ).ignore();
+      myTuple->column( "xBeam", beamline.X ).ignore();
+      myTuple->column( "yBeam", beamline.Y ).ignore();
+      myTuple->column( "errx", errx ).ignore();
+      myTuple->column( "erry", erry ).ignore();
+      myTuple->column( "errz", errz ).ignore();
+      myTuple->column( "errr", errr ).ignore();
+      myTuple->column( "mother", double( mother ) ).ignore();
+      myTuple->column( "evtnr", double( m_nevt.nEntries() ) ).ignore();
+      myTuple->column( "chi2", double( chi2 ) ).ignore();
+      myTuple->column( "nDoF", double( nDoF ) ).ignore();
+      myTuple->column( "size_tracks", double( tracks.size() ) ).ignore();
+      myTuple->column( "size_mcp", double( mcps.size() ) ).ignore();
+      myTuple->column( "mcpvrec", double( nmrc ) ).ignore();
+      myTuple->write().ignore();
     }
   }
 
-  return StatusCode::SUCCESS;
-}
-
-void PrimaryVertexChecker::match_mc_vertex_by_tracks( int ipv, std::vector<RecPVInfo>& rinfo,
-                                                      std::vector<MCPVInfo>& mcpvvec ) {
-
-  LHCb::RecVertex*                          invtx = rinfo[ipv].pRECPV;
-  typedef const SmartRefVector<LHCb::Track> PVTRACKS;
-  PVTRACKS&                                 tracksIn = invtx->tracks();
-  PVTRACKS::const_iterator                  itin;
-
-  int    indexmc  = -1;
-  double ratiomax = 0.;
-  for ( int imc = 0; imc < (int)mcpvvec.size(); imc++ ) {
-    if ( mcpvvec[imc].indexRecPVInfo > -1 ) continue;
-    int ntrin = 0;
-    for ( std::vector<LHCb::Track*>::iterator itr = mcpvvec[imc].m_recTracksInMCPV.begin();
-          itr != mcpvvec[imc].m_recTracksInMCPV.end(); itr++ ) {
-      for ( itin = tracksIn.begin(); itin != tracksIn.end(); itin++ ) {
-        const LHCb::Track* ptr = *itin;
-        if ( ptr == *itr ) {
-          ntrin++;
-          break;
-        }
-      }
-    }
-    double ratio = 1. * ntrin / tracksIn.size();
-    if ( ratio > ratiomax ) {
-      ratiomax = ratio;
-      indexmc  = imc;
-    }
-  } // imc
-  if ( ratiomax > 0.05 ) {
-    rinfo[ipv].indexMCPVInfo        = indexmc;
-    mcpvvec[indexmc].indexRecPVInfo = ipv;
-  }
+  return;
 }
 
 void PrimaryVertexChecker::match_mc_vertex_by_distance( int ipv, std::vector<RecPVInfo>& rinfo,
-                                                        std::vector<MCPVInfo>& mcpvvec ) {
+                                                        std::vector<MCPVInfo>& mcpvvec ) const {
 
   double mindist = 999999.;
   int    indexmc = -1;
 
-  for ( int imc = 0; imc < (int)mcpvvec.size(); imc++ ) {
-    if ( mcpvvec[imc].indexRecPVInfo > -1 ) continue;
-    double dist = fabs( mcpvvec[imc].pMCPV->position().z() - rinfo[ipv].position.z() );
+  for ( const auto& [imc, mcpv] : LHCb::range::enumerate( mcpvvec ) ) {
+    if ( mcpv.indexRecPVInfo > -1 ) continue;
+    double dist = fabs( mcpv.pMCPV->position().z() - rinfo[ipv].position.z() );
     if ( dist < mindist ) {
       mindist = dist;
       indexmc = imc;
@@ -711,16 +630,16 @@ void PrimaryVertexChecker::match_mc_vertex_by_distance( int ipv, std::vector<Rec
   }
 }
 
-std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<MCPVInfo>&           rblemcpv,
-                                                                   std::vector<MCPVInfo>::iterator& itmc ) {
+std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<MCPVInfo>& rblemcpv,
+                                                                   const MCPVInfo&        mc ) const {
 
-  std::vector<MCPVInfo>::iterator itret   = rblemcpv.end();
-  double                          mindist = 999999.;
+  auto   itret   = rblemcpv.end();
+  double mindist = 999999.;
   if ( rblemcpv.size() < 2 ) return itret;
   std::vector<MCPVInfo>::iterator it;
   for ( it = rblemcpv.begin(); it != rblemcpv.end(); it++ ) {
-    if ( it->pMCPV != itmc->pMCPV ) {
-      double dist = ( it->pMCPV->position() - itmc->pMCPV->position() ).R();
+    if ( it->pMCPV != mc.pMCPV ) {
+      double dist = ( it->pMCPV->position() - mc.pMCPV->position() ).R();
       if ( dist < mindist ) {
         mindist = dist;
         itret   = it;
@@ -730,20 +649,20 @@ std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<M
   return itret;
 }
 
-void PrimaryVertexChecker::collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx,
-                                             std::vector<LHCb::MCParticle*>& allprods ) {
-
-  SmartRefVector<LHCb::MCParticle>           daughters = mcvtx->products();
-  SmartRefVector<LHCb::MCParticle>::iterator idau;
-  for ( idau = daughters.begin(); idau != daughters.end(); idau++ ) {
-    double dv2 = ( mcpv->position() - ( *idau )->originVertex()->position() ).Mag2();
-    if ( dv2 > ( 100. * Gaudi::Units::mm ) * ( 100. * Gaudi::Units::mm ) ) continue;
-    LHCb::MCParticle* pmcp = *idau;
-    allprods.push_back( pmcp );
-    SmartRefVector<LHCb::MCVertex>           decays = ( *idau )->endVertices();
-    SmartRefVector<LHCb::MCVertex>::iterator ivtx;
-    for ( ivtx = decays.begin(); ivtx != decays.end(); ivtx++ ) { collectProductss( mcpv, *ivtx, allprods ); }
+std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<MCPVInfo>& rblemcpv,
+                                                                   const RecPVInfo&       rec ) const {
+  auto   itret   = rblemcpv.end();
+  double mindist = 999999.;
+  if ( rblemcpv.size() < 2 ) return itret;
+  std::vector<MCPVInfo>::iterator it;
+  for ( it = rblemcpv.begin(); it != rblemcpv.end(); it++ ) {
+    double dist = ( it->pMCPV->position() - rec.pRECPV->position() ).R();
+    if ( dist < mindist ) {
+      mindist = dist;
+      itret   = it;
+    }
   }
+  return itret;
 }
 
 //=============================================================================
@@ -752,290 +671,204 @@ void PrimaryVertexChecker::collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVerte
 StatusCode PrimaryVertexChecker::finalize() {
 
   debug() << "==> Finalize" << endmsg;
+  info() << "     ************************************ " << endmsg;
+  info() << " MC PV is reconstructible if at least " << m_nTracksToBeRecble.value() << "  tracks are reconstructed"
+         << endmsg;
+  info() << " MC PV is isolated if dz to closest reconstructible MC PV >  " << m_dzIsolated.value() << " mm" << endmsg;
+  info() << " MC efficiency split by tracks with threshold: (" << m_nTracksToBeRecble.value() << ","
+         << m_nTracksToPrint.value() << "), >= " << m_nTracksToPrint.value() << endmsg;
+  info() << " MC efficiency split by z position: <-" << m_zToPrint.value() << ", (-" << m_zToPrint.value() << ","
+         << m_zToPrint.value() << "), >" << m_zToPrint.value() << endmsg;
+  std::string ff = "by distance";
+  info() << " REC and MC vertices matched:  " << ff << endmsg;
 
-  info() << " ============================================" << endmsg;
-  info() << " Efficiencies for reconstructible MC vertices: " << endmsg;
-  info() << " ============================================" << endmsg;
+  info() << " " << endmsg;
+  for ( const auto& n : All ) { printRat( toString( n ), m_mcpv[n], m_eff[n], m_false[n] ); }
+  info() << " " << endmsg;
+  for ( const auto& n : All ) { printAvTracks( toString( n ), m_av_tracks[n], m_av_mcp[n] ); }
   info() << " " << endmsg;
 
-  info() << " MC PV is reconstructible if at least " << m_nTracksToBeRecble << "  tracks are reconstructed" << endmsg;
-  info() << " MC PV is isolated if dz to closest reconstructible MC PV >  " << m_dzIsolated << " mm" << endmsg;
-  std::string ff = "by counting tracks";
-  if ( !m_matchByTracks ) ff = "by dz distance";
-  info() << " REC and MC vertices matched:  " << ff << endmsg;
+  printRes( "1_res_all", check_histogram( histo( HistoID( 1021 ) ), true ),
+            check_histogram( histo( HistoID( 1022 ) ), true ), check_histogram( histo( HistoID( 1023 ) ), true ) );
+
+  printRes( "2_res_ntracks<10", check_histogram( histo( HistoID( 1101 ) ), true ),
+            check_histogram( histo( HistoID( 1102 ) ), true ), check_histogram( histo( HistoID( 1103 ) ), true ) );
+  printRes( "3_res_ntracks(10,30)", check_histogram( histo( HistoID( 1111 ) ), true ),
+            check_histogram( histo( HistoID( 1112 ) ), true ), check_histogram( histo( HistoID( 1113 ) ), true ) );
+  printRes( "4_res_ntracks>30", check_histogram( histo( HistoID( 1121 ) ), true ),
+            check_histogram( histo( HistoID( 1122 ) ), true ), check_histogram( histo( HistoID( 1123 ) ), true ) );
+
+  printRes( "5_res_z<-50", check_histogram( histo( HistoID( 1201 ) ), true ),
+            check_histogram( histo( HistoID( 1202 ) ), true ), check_histogram( histo( HistoID( 1203 ) ), true ) );
+  printRes( "6_res_z(-50,50)", check_histogram( histo( HistoID( 1211 ) ), true ),
+            check_histogram( histo( HistoID( 1212 ) ), true ), check_histogram( histo( HistoID( 1213 ) ), true ) );
+  printRes( "7_res_z>50", check_histogram( histo( HistoID( 1221 ) ), true ),
+            check_histogram( histo( HistoID( 1222 ) ), true ), check_histogram( histo( HistoID( 1223 ) ), true ) );
 
   info() << " " << endmsg;
 
-  printRat( "All", m_nRecMCPV, m_nMCPV );
-  printRat( "Isolated", m_nRecMCPV_isol, m_nMCPV_isol );
-  printRat( "Close", m_nRecMCPV_close, m_nMCPV_close );
-  printRat( "False rate", m_nFalsePV, m_nRecMCPV + m_nFalsePV );
-
-  if ( debugLevel() ) { printRat( "Real false rate", m_nFalsePV_real, m_nRecMCPV + m_nFalsePV_real ); }
-
-  info() << endmsg;
-  printRat( "False PV as B", m_nRecBFalse, m_nBFalse );
-  info() << endmsg;
-
-  info() << "      --------------------------------------------" << endmsg;
-  info() << "           Substatistics: " << endmsg;
-  info() << "      --------------------------------------------" << endmsg;
-  info() << "      1st PV (highest multiplicity): " << endmsg;
-  printRat( "All", m_nRecMCPV_1mult, m_nMCPV_1mult );
-  printRat( "Isolated", m_nRecMCPV_isol_1mult, m_nMCPV_isol_1mult );
-  printRat( "Close", m_nRecMCPV_close_1mult, m_nMCPV_close_1mult );
-
-  info() << "      ---------------------------------------" << endmsg;
-  info() << "      2nd PV: " << endmsg;
-  printRat( "All", m_nRecMCPV_2mult, m_nMCPV_2mult );
-  printRat( "Isolated", m_nRecMCPV_isol_2mult, m_nMCPV_isol_2mult );
-  printRat( "Close", m_nRecMCPV_close_2mult, m_nMCPV_close_2mult );
-
-  info() << "      ---------------------------------------" << endmsg;
-  info() << "      3rd PV: " << endmsg;
-  printRat( "All", m_nRecMCPV_3mult, m_nMCPV_3mult );
-  printRat( "Isolated", m_nRecMCPV_isol_3mult, m_nMCPV_isol_3mult );
-  printRat( "Close", m_nRecMCPV_close_3mult, m_nMCPV_close_3mult );
+  printRes( "1_pull_width_all", check_histogram( histo( HistoID( 1031 ) ), true ),
+            check_histogram( histo( HistoID( 1032 ) ), true ), check_histogram( histo( HistoID( 1033 ) ), true ) );
+  printRes( "1_pull_mean_all", check_histogram( histo( HistoID( 1031 ) ), false ),
+            check_histogram( histo( HistoID( 1032 ) ), false ), check_histogram( histo( HistoID( 1033 ) ), false ) );
 
   info() << " " << endmsg;
-  if ( debugLevel() ) {
-    info() << " * Real false rate means: no visible MC PV within 5 sigma of REC PV."
-           << " Visible MC PV: 2 tracks reconstructed" << endmsg;
-    info() << " " << endmsg;
-  }
-  const AIDA::IHistogram1D* dx    = histo( HistoID( 1021 ) );
-  const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
-  const AIDA::IHistogram1D* dy    = histo( HistoID( 1022 ) );
-  const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
-  const AIDA::IHistogram1D* dz    = histo( HistoID( 1023 ) );
-  const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
-  if ( dx ) {
-    info() << "      ---------------------------------------" << endmsg;
-    info() << "dx:    "
-           << format( "mean =  %5.3f +/- %5.3f, RMS =  %5.3f +/- %5.3f", dx->mean(),
-                      Gaudi::Utils::HistoStats::meanErr( dx ), dx->rms(), Gaudi::Utils::HistoStats::rmsErr( dx ) )
-           << endmsg;
-  }
-  if ( dy ) {
-    info() << "dy:    "
-           << format( "mean =  %5.3f +/- %5.3f, RMS =  %5.3f +/- %5.3f", dy->mean(),
-                      Gaudi::Utils::HistoStats::meanErr( dy ), dy->rms(), Gaudi::Utils::HistoStats::rmsErr( dy ) )
-           << endmsg;
-  }
-  if ( dz ) {
-    info() << "dz:    "
-           << format( "mean =  %5.3f +/- %5.3f, RMS =  %5.3f +/- %5.3f", dz->mean(),
-                      Gaudi::Utils::HistoStats::meanErr( dz ), dz->rms(), Gaudi::Utils::HistoStats::rmsErr( dz ) )
-           << endmsg;
-  }
-  info() << "      ---------------------------------------" << endmsg;
-  if ( pullx ) {
-    info() << "pullx: "
-           << format( "mean =  %5.3f +/- %5.3f, RMS =  %5.3f +/- %5.3f", pullx->mean(),
-                      Gaudi::Utils::HistoStats::meanErr( pullx ), pullx->rms(),
-                      Gaudi::Utils::HistoStats::rmsErr( pullx ) )
-           << endmsg;
-  }
-  if ( pully ) {
-    info() << "pully: "
-           << format( "mean =  %5.3f +/- %5.3f, RMS =  %5.3f +/- %5.3f", pully->mean(),
-                      Gaudi::Utils::HistoStats::meanErr( pully ), pully->rms(),
-                      Gaudi::Utils::HistoStats::rmsErr( pully ) )
-           << endmsg;
-  }
-  if ( pullz ) {
-    info() << "pullz: "
-           << format( "mean =  %5.3f +/- %5.3f, RMS =  %5.3f +/- %5.3f", pullz->mean(),
-                      Gaudi::Utils::HistoStats::meanErr( pullz ), pullz->rms(),
-                      Gaudi::Utils::HistoStats::rmsErr( pullz ) )
-           << endmsg;
-  }
-  info() << " ============================================" << endmsg;
 
   return GaudiTupleAlg::finalize(); // Must be called after all other actions
 }
 
-void PrimaryVertexChecker::printRat( std::string mes, int a, int b ) {
+void PrimaryVertexChecker::printRat( const std::string name, const Gaudi::Accumulators::BinomialCounter<>& m_mcpv,
+                                     const Gaudi::Accumulators::BinomialCounter<>& m_eff,
+                                     const Gaudi::Accumulators::BinomialCounter<>& m_false ) {
+  unsigned int len  = 25;
+  std::string  pmes = name;
+  if ( pmes.size() < len ) pmes.resize( len, ' ' );
+  pmes += " : ";
+  info() << pmes
+         << format( "%8d from %8d (%8d-%-8d) [ %5.2f %], false %4d from reco. %8d (%8d+%-4d) [ %5.2f %] ",
+                    m_eff.nTrueEntries(), m_eff.nEntries(), m_mcpv.nEntries(), m_mcpv.nFalseEntries(),
+                    m_eff.efficiency() * 100, m_false.nTrueEntries(), m_false.nEntries(), m_false.nFalseEntries(),
+                    m_false.nTrueEntries(), m_false.efficiency() * 100 )
+         << endmsg;
+}
 
-  double rat = 0.;
-  if ( b > 0 ) rat = 1.0 * a / b;
+void PrimaryVertexChecker::printAvTracks( const std::string                              name,
+                                          const Gaudi::Accumulators::AveragingCounter<>& m_av_tracks,
+                                          const Gaudi::Accumulators::AveragingCounter<>& m_av_mcp ) {
+  unsigned int len  = 25;
+  std::string  pmes = name;
+  if ( pmes.size() < len ) pmes.resize( len, ' ' );
+  pmes += " : ";
+  info() << pmes << format( "av. PV tracks: %6.2f [MC: %6.2f]", m_av_tracks.mean(), m_av_mcp.mean() ) << endmsg;
+}
 
-  // reformat message
-  unsigned int len  = 20;
+void PrimaryVertexChecker::printRes( std::string mes, double x, double y, double z ) {
+  unsigned int len  = 25;
   std::string  pmes = mes;
   while ( pmes.length() < len ) { pmes += " "; }
   pmes += " : ";
-
-  info() << pmes << format( " %6.3f ( %7d / %8d )", rat, a, b ) << endmsg;
+  info() << pmes << format( " x: %+5.3f, y: %+5.3f, z: %+5.3f", x, y, z ) << endmsg;
 }
 
-void PrimaryVertexChecker::count_reconstructed_tracks( std::vector<MCPVInfo>&     mcpvvec,
-                                                       std::vector<LHCb::Track*>& vecOfTracks, std::string trackLoc ) {
-
-  LinkedTo<LHCb::MCParticle> trackMClink( eventSvc(), msgSvc(), trackLoc );
-
-  // Find # of reconstructed tracks of every MC PV
-  std::vector<MCPVInfo>::iterator itinfomc;
-  for ( itinfomc = mcpvvec.begin(); mcpvvec.end() != itinfomc; itinfomc++ ) {
-    LHCb::MCVertex* avtx = itinfomc->pMCPV;
-
-    SmartRefVector<LHCb::MCParticle> parts = avtx->products();
-    std::vector<LHCb::MCParticle*>   allproducts;
-    collectProductss( avtx, avtx, allproducts );
+double PrimaryVertexChecker::check_histogram( const AIDA::IHistogram1D* h, bool rms ) {
+  return h ? ( rms ? h->rms() : h->mean() ) : 0.;
+}
 
-    LHCb::Track*                             recTrack = 0;
-    std::vector<LHCb::MCParticle*>::iterator imcp;
-    for ( imcp = allproducts.begin(); allproducts.end() != imcp; imcp++ ) {
-      LHCb::MCParticle* pmcp   = *imcp;
-      int               isReco = 0;
-      for ( std::vector<LHCb::Track*>::iterator itrec = vecOfTracks.begin(); itrec != vecOfTracks.end(); itrec++ ) {
-        LHCb::MCParticle* partmc = trackMClink.first( ( *itrec )->key() );
-        if ( partmc && partmc == pmcp ) {
-          recTrack = ( *itrec );
-          if ( !m_requireVelo || recTrack->hasVelo() ) {
-            isReco = 1;
-            break;
-          }
-        }
-      }
-      if ( pmcp->particleID().threeCharge() != 0 && isReco ) {
-        double dv2 = ( avtx->position() - pmcp->originVertex()->position() ).Mag2();
-        if ( dv2 < 0.0001 ) {
-          itinfomc->m_mcPartInMCPV.push_back( pmcp );
-          itinfomc->m_recTracksInMCPV.push_back( recTrack );
-        }
-      }
-    }
-    itinfomc->nRecTracks = itinfomc->m_mcPartInMCPV.size();
-  }
+int PrimaryVertexChecker::count_velo_tracks( const std::vector<LHCb::Event::v2::WeightedTrack>& tracksPV ) {
+  return std::count_if( tracksPV.begin(), tracksPV.end(), []( const auto& t ) { return t.track->hasVelo(); } );
 }
 
-int PrimaryVertexChecker::count_velo_tracks( LHCb::RecVertex* RecVtx ) {
-  SmartRefVector<LHCb::Track> vtx_tracks  = RecVtx->tracks();
-  int                         nVeloTracks = 0;
-  for ( unsigned int it = 0; it < vtx_tracks.size(); it++ ) {
-    const LHCb::Track* ptr = vtx_tracks[it];
-    if ( ptr->hasVelo() ) nVeloTracks++;
+void PrimaryVertexChecker::collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx,
+                                             std::vector<LHCb::MCParticle*>& allprods ) const {
+
+  SmartRefVector<LHCb::MCParticle>           daughters = mcvtx->products();
+  SmartRefVector<LHCb::MCParticle>::iterator idau;
+  for ( idau = daughters.begin(); idau != daughters.end(); idau++ ) {
+    double dv2 = ( mcpv->position() - ( *idau )->originVertex()->position() ).Mag2();
+    if ( dv2 > ( 100. * Gaudi::Units::mm ) * ( 100. * Gaudi::Units::mm ) ) continue;
+    LHCb::MCParticle* pmcp = *idau;
+    allprods.push_back( pmcp );
+    SmartRefVector<LHCb::MCVertex>           decays = ( *idau )->endVertices();
+    SmartRefVector<LHCb::MCVertex>::iterator ivtx;
+    for ( ivtx = decays.begin(); ivtx != decays.end(); ivtx++ ) { collectProductss( mcpv, *ivtx, allprods ); }
   }
-  return nVeloTracks;
 }
 
-void PrimaryVertexChecker::count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec ) {
+void PrimaryVertexChecker::count_reconstructible_mc_particles( std::vector<MCPVInfo>&  mcpvvec,
+                                                               const LHCb::MCProperty& flags ) const {
 
-  const auto trInfo = MCTrackInfo{*get<LHCb::MCProperty>( LHCb::MCPropertyLocation::TrackInfo )};
-  for ( auto itinfomc = mcpvvec.begin(); mcpvvec.end() != itinfomc; itinfomc++ ) {
-    LHCb::MCVertex*                  avtx = itinfomc->pMCPV;
+  const MCTrackInfo trInfo = {flags};
+  for ( auto& itinfomc : mcpvvec ) {
+    LHCb::MCVertex*                  avtx = itinfomc.pMCPV;
     std::vector<LHCb::MCParticle*>   mcPartInMCPV;
     SmartRefVector<LHCb::MCParticle> parts = avtx->products();
     std::vector<LHCb::MCParticle*>   allproducts;
     collectProductss( avtx, avtx, allproducts );
 
-    for ( LHCb::MCParticle* pmcp : allproducts ) {
-      if ( pmcp->particleID().threeCharge() != 0 && ( !m_requireVelo || trInfo.hasVelo( pmcp ) ) ) {
-
-        if ( pmcp->particleID().hasBottom() ) { itinfomc->decayBeauty = 1.0; }
-        if ( pmcp->particleID().hasCharm() ) { itinfomc->decayCharm = 1.0; }
-        if ( pmcp->particleID().threeCharge() != 0 && ( !m_requireVelo || trInfo.hasVelo( pmcp ) ) ) {
-
-          double dv2 = ( avtx->position() - pmcp->originVertex()->position() ).Mag2();
-          if ( dv2 < 0.0000001 && pmcp->p() > 100. * Gaudi::Units::MeV ) { mcPartInMCPV.push_back( pmcp ); }
-        }
+    for ( const auto pmcp : allproducts ) {
+      if ( pmcp->particleID().isMeson() || pmcp->particleID().isBaryon() ) {
+        if ( pmcp->particleID().hasBottom() ) { itinfomc.decayBeauty = true; }
+        if ( pmcp->particleID().hasCharm() ) { itinfomc.decayCharm = true; }
+        if ( pmcp->particleID().hasStrange() ) { itinfomc.decayStrange = true; }
       }
-      itinfomc->nRecTracks = mcPartInMCPV.size();
-    }
-  }
-}
-int PrimaryVertexChecker::check_mother_particle( LHCb::RecVertex* pv, std::string trackLoc ) {
-
-  int         motherPart = 0;
-  const auto& tracksIn   = pv->tracks();
-
-  LinkedTo<LHCb::MCParticle, LHCb::Track> link( eventSvc(), msgSvc(), trackLoc );
-  for ( const LHCb::Track* ptr : tracksIn ) {
-
-    LHCb::MCParticle* part = nullptr;
-    part                   = link.first( ptr->key() );
-    if ( part != nullptr ) {
-      const LHCb::MCParticle* part1 = part;
-      int                     i     = 0;
-      while ( part1 != nullptr ) {
-        // std::cout<<"X"<<std::endl;
-        const LHCb::MCParticle* mother = part1->mother();
-        if ( mother != nullptr ) {
-          i     = i + 1;
-          part1 = mother;
-          /// std::cout<<i<<": "<<part1->particleID().pid()<<std::endl;
-        } else {
-          break;
-        }
+      if ( ( pmcp->particleID().isMeson() || pmcp->particleID().isBaryon() ) &&
+           ( !m_requireVelo || trInfo.hasVelo( pmcp ) ) ) {
+        double dv2 = ( avtx->position() - pmcp->originVertex()->position() ).Mag2();
+        if ( dv2 < 0.0000001 && pmcp->p() > 100. * Gaudi::Units::MeV ) { mcPartInMCPV.push_back( pmcp ); }
       }
-
-      if ( part1->particleID().hasBottom() == true ) { motherPart = 2.0; }
-      if ( part1->particleID().hasCharm() == true ) { motherPart = 1.0; }
+      itinfomc.nRecTracks = mcPartInMCPV.size();
     }
   }
-  return motherPart;
 }
 
-bool PrimaryVertexChecker::getInputTracks( std::vector<LHCb::Track*>& vecOfTracks, std::string& trackLoc ) {
-
-  std::string tracksName = m_inputTracksName;
-
-  if ( m_inputTracksName == "Offline" ) {
-    tracksName = LHCb::TrackLocation::Default;
-  } else if ( m_inputTracksName == "3D" ) {
-    tracksName = LHCb::TrackLocation::Velo;
-  } else if ( m_inputTracksName == "2D" ) {
-    tracksName = LHCb::TrackLocation::RZVelo;
-  }
-
-  if ( tracksName == "none" ) {
-    debug() << " Tracks not specified " << tracksName << endmsg;
+bool PrimaryVertexChecker::checkCondition( const MCPVInfo& MCPV, const recoAs& n ) const {
+  switch ( n ) {
+  case recoAs::all:
+    return true;
+  case recoAs::isolated:
+    return ( MCPV.distToClosestMCPV > m_dzIsolated );
+  case recoAs::close:
+    return ( MCPV.distToClosestMCPV <= m_dzIsolated );
+  case recoAs::ntracks_low:
+    return ( MCPV.nRecTracks >= m_nTracksToBeRecble && MCPV.nRecTracks < m_nTracksToPrint );
+  case recoAs::ntracks_high:
+    return ( MCPV.nRecTracks >= m_nTracksToPrint );
+  case recoAs::z_low:
+    return ( MCPV.pMCPV->position().z() < -m_zToPrint );
+  case recoAs::z_middle:
+    return ( MCPV.pMCPV->position().z() >= -m_zToPrint && MCPV.pMCPV->position().z() < m_zToPrint );
+  case recoAs::z_high:
+    return ( MCPV.pMCPV->position().z() >= m_zToPrint );
+  case recoAs::beauty:
+    return ( MCPV.decayBeauty );
+  case recoAs::charm:
+    return ( MCPV.decayCharm );
+  case recoAs::strange:
+    return ( MCPV.decayStrange );
+  case recoAs::other:
+    return ( !( MCPV.decayBeauty ) && !( MCPV.decayCharm ) && !( MCPV.decayStrange ) );
+  default:
     return false;
   }
-
-  trackLoc = tracksName;
-
-  LHCb::Tracks* usedTracks;
-
-  usedTracks = getOrCreate<LHCb::Tracks, LHCb::Tracks>( tracksName );
-  if ( usedTracks->size() == 0 ) return false;
-
-  std::vector<LHCb::Track*>::const_iterator itT;
-  for ( itT = usedTracks->begin(); usedTracks->end() != itT; itT++ ) {
-    LHCb::Track* ptr = ( *itT );
-    vecOfTracks.push_back( ptr );
-  }
-
-  return true;
 }
 
-bool PrimaryVertexChecker::getInputVertices( std::vector<LHCb::RecVertex*>& vecOfVertices ) {
-
-  std::string verticesName = m_inputVerticesName;
-
-  if ( m_inputVerticesName == "Offline" ) {
-    verticesName = LHCb::RecVertexLocation::Primary;
-  } else if ( m_inputVerticesName == "3D" ) {
-    verticesName = LHCb::RecVertexLocation::Velo3D;
-  } else if ( m_inputVerticesName == "2D" ) {
-    verticesName = LHCb::RecVertexLocation::Velo2D;
-  }
-
-  if ( verticesName == "none" ) {
-    debug() << " Vertices not specified " << verticesName << endmsg;
-    return false;
+std::string PrimaryVertexChecker::toString( const recoAs& n ) const {
+  switch ( n ) {
+  case recoAs::all:
+    return format( "0%d all", int( recoAs::all ) );
+  case recoAs::isolated:
+    return format( "0%d isolated", int( recoAs::isolated ) );
+  case recoAs::close:
+    return format( "0%d close", int( recoAs::close ) );
+  case recoAs::ntracks_low:
+    return format( "0%d ntracks<%d", int( recoAs::ntracks_low ), int( m_nTracksToPrint ) );
+  case recoAs::ntracks_high:
+    return format( "0%d ntracks>=%d", int( recoAs::ntracks_high ), int( m_nTracksToPrint ) );
+  case recoAs::z_low:
+    return format( "0%d z<%2.1f", int( recoAs::z_low ), -m_zToPrint );
+  case recoAs::z_middle:
+    return format( "0%d z in (%2.1f, %2.1f)", int( recoAs::z_middle ), -m_zToPrint, +m_zToPrint );
+  case recoAs::z_high:
+    return format( "0%d z >=%2.1f", int( recoAs::z_high ), +m_zToPrint );
+  case recoAs::beauty:
+    return format( "0%d decayBeauty", int( recoAs::beauty ) );
+  case recoAs::charm:
+    return format( "0%d decayCharm", int( recoAs::charm ) );
+  case recoAs::strange:
+    return format( "%d decayStrange", int( recoAs::strange ) );
+  case recoAs::other:
+    return format( "%d other", int( recoAs::other ) );
+  case recoAs::first:
+    return format( "%d 1MCPV", int( recoAs::first ) );
+  case recoAs::second:
+    return format( "%d 2MCPV", int( recoAs::second ) );
+  case recoAs::third:
+    return format( "%d 3MCPV", int( recoAs::third ) );
+  case recoAs::fourth:
+    return format( "%d 4MCPV", int( recoAs::fourth ) );
+  case recoAs::fifth:
+    return format( "%d 5MCPV", int( recoAs::fifth ) );
+  default:
+    return "not defined";
   }
-
-  LHCb::RecVertices* recoVertices;
-
-  recoVertices = getOrCreate<LHCb::RecVertices, LHCb::RecVertices>( verticesName );
-
-  std::vector<LHCb::RecVertex*>::const_iterator itVer;
-  for ( itVer = recoVertices->begin(); recoVertices->end() != itVer; itVer++ ) {
-    LHCb::RecVertex* ppv = ( *itVer );
-    vecOfVertices.push_back( ppv );
-  }
-
-  return true;
 }
diff --git a/Tr/PatChecker/src/PrimaryVertexChecker.h b/Tr/PatChecker/src/PrimaryVertexChecker.h
deleted file mode 100644
index b5d62bd502903e1df6e3fd205b15253964dd814b..0000000000000000000000000000000000000000
--- a/Tr/PatChecker/src/PrimaryVertexChecker.h
+++ /dev/null
@@ -1,111 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration      *
-*                                                                             *
-* This software is distributed under the terms of the GNU General Public      *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
-*                                                                             *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization  *
-* or submit itself to any jurisdiction.                                       *
-\*****************************************************************************/
-#ifndef PRIMARYVERTEXCHECKER_H
-#define PRIMARYVERTEXCHECKER_H 1
-
-#include "Event/MCVertex.h"
-#include "Event/RecVertex.h"
-#include "Event/Track.h"
-#include "MCInterfaces/IForcedBDecayTool.h"
-
-#include "GaudiAlg/GaudiTupleAlg.h"
-
-typedef struct {
-  LHCb::MCVertex*                pMCPV;             // pointer to MC PV
-  int                            nRecTracks;        // number of reconstructed tracks from this MCPV
-  int                            nRecBackTracks;    // number of reconstructed backward tracks
-  int                            indexRecPVInfo;    // index to reconstructed PVInfo (-1 if not reco)
-  int                            nCorrectTracks;    // correct tracks belonging to reconstructed PV
-  int                            multClosestMCPV;   // multiplicity of closest reconstructable MCPV
-  double                         distToClosestMCPV; // distance to closest reconstructible MCPV
-  int                            decayCharm;        // type of mother particle
-  int                            decayBeauty;
-  std::vector<LHCb::MCParticle*> m_mcPartInMCPV;
-  std::vector<LHCb::Track*>      m_recTracksInMCPV;
-} MCPVInfo;
-
-typedef struct {
-public:
-  int              nTracks;     // number of tracks
-  int              nVeloTracks; // number of velo tracks in a vertex
-  int              nLongTracks;
-  double           minTrackRD; //
-  double           maxTrackRD; //
-  double           chi2;
-  double           nDoF;
-  double           d0;
-  double           d0nTr;
-  double           chi2nTr;
-  double           mind0;
-  double           maxd0;
-  int              mother;
-  Gaudi::XYZPoint  position;      // position
-  Gaudi::XYZPoint  positionSigma; // position sigmas
-  int              indexMCPVInfo; // index to MCPVInfo
-  LHCb::RecVertex* pRECPV;        // pointer to REC PV
-} RecPVInfo;
-
-class PrimaryVertexChecker : public GaudiTupleAlg {
-public:
-  /// Standard constructor
-  PrimaryVertexChecker( const std::string& name, ISvcLocator* pSvcLocator );
-
-  StatusCode initialize() override; ///< Algorithm initialization
-  StatusCode execute() override;    ///< Algorithm execution
-  StatusCode finalize() override;   ///< Algorithm finalization
-
-protected:
-  bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
-
-private:
-  int         m_nTracksToBeRecble;
-  bool        m_produceHistogram;
-  bool        m_produceNtuple;
-  bool        m_requireVelo;
-  double      m_dzIsolated;
-  bool        m_matchByTracks;
-  std::string m_inputTracksName;
-  std::string m_inputVerticesName;
-
-  int m_nevt;
-  int m_nMCPV, m_nRecMCPV;
-  int m_nMCPV_isol, m_nRecMCPV_isol;
-  int m_nMCPV_close, m_nRecMCPV_close;
-  int m_nFalsePV, m_nFalsePV_real;
-  int m_nMCPV_1mult, m_nRecMCPV_1mult;
-  int m_nMCPV_isol_1mult, m_nRecMCPV_isol_1mult;
-  int m_nMCPV_close_1mult, m_nRecMCPV_close_1mult;
-  int m_nRecMCPV_wrong_1mult;
-  int m_nMCPV_2mult, m_nRecMCPV_2mult;
-  int m_nMCPV_isol_2mult, m_nRecMCPV_isol_2mult;
-  int m_nMCPV_close_2mult, m_nRecMCPV_close_2mult;
-  int m_nRecMCPV_wrong_2mult;
-  int m_nMCPV_3mult, m_nRecMCPV_3mult;
-  int m_nMCPV_isol_3mult, m_nRecMCPV_isol_3mult;
-  int m_nMCPV_close_3mult, m_nRecMCPV_close_3mult;
-  int m_nRecMCPV_wrong_3mult;
-  int m_nBFalse, m_nRecBFalse;
-
-  std::vector<MCPVInfo>::iterator closestMCPV( std::vector<MCPVInfo>& rblemcpv, std::vector<MCPVInfo>::iterator& itmc );
-
-  void collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx, std::vector<LHCb::MCParticle*>& allprods );
-  void printRat( std::string mes, int a, int b );
-  void match_mc_vertex_by_tracks( int ipv, std::vector<RecPVInfo>& rinfo, std::vector<MCPVInfo>& mcpvvec );
-  void match_mc_vertex_by_distance( int ipv, std::vector<RecPVInfo>& rinfo, std::vector<MCPVInfo>& mcpvvec );
-  void count_reconstructed_tracks( std::vector<MCPVInfo>& mcpvvec, std::vector<LHCb::Track*>& vecOfTracks,
-                                   std::string trackLoc );
-  int  count_velo_tracks( LHCb::RecVertex* RecVtx );
-  void count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec );
-  bool getInputTracks( std::vector<LHCb::Track*>& vecOfTracks, std::string& trackLoc );
-  bool getInputVertices( std::vector<LHCb::RecVertex*>& vecOfVertices );
-  int  check_mother_particle( LHCb::RecVertex* RecVtx, std::string trackLoc );
-};
-#endif // PRIMARYVERTEXCHECKER_H