diff --git a/Phys/ParticleMatching/src/MassConstrainedVeloTrack.cpp b/Phys/ParticleMatching/src/MassConstrainedVeloTrack.cpp
index c3748f396e02b34969a596961cb1bb808cdc78a3..9a042d1a3a3c3f8dfdbdee71622470086d2c7a6b 100644
--- a/Phys/ParticleMatching/src/MassConstrainedVeloTrack.cpp
+++ b/Phys/ParticleMatching/src/MassConstrainedVeloTrack.cpp
@@ -18,7 +18,7 @@
 #include "TrackInterfaces/ITrackFitter.h"
 
 namespace {
-  using output_t = std::tuple<LHCb::Particles, LHCb::Particles, LHCb::ProtoParticles, LHCb::Tracks>;
+  using output_t = std::tuple<LHCb::Particles, LHCb::Particles, LHCb::ProtoParticles, LHCb::Tracks, LHCb::Vertices>;
 }
 // Given a decay chain, the Algorithm applies a mass constraint to the head of the chain
 // to infer the momentum of one of the children which is a Velo track.
@@ -35,7 +35,7 @@ public:
             pSvcLocator,
             { KeyValue{ "InputParticles", "" }, KeyValue{ "StandardGeometryTop", LHCb::standard_geometry_top } },
             { KeyValue{ "OutputParticles", "" }, KeyValue{ "OutputVeloParticles", "" }, KeyValue{ "VeloProtos", "" },
-              KeyValue{ "FittedVeloTracks", "" } } } {}
+              KeyValue{ "FittedVeloTracks", "" }, KeyValue{ "OutputVertices", "" } } } {}
 
   output_t operator()( const LHCb::Particle::Range&, const DetectorElement& ) const override;
 
@@ -58,7 +58,7 @@ output_t MassConstrainedVeloTrack::operator()( const LHCb::Particle::Range& comp
                                                const DetectorElement&       lhcb ) const {
 
   output_t ret{};
-  auto& [output, velo_output, pp_output, track_output] = ret;
+  auto& [output, velo_output, pp_output, track_output, vertices_output] = ret;
   m_nIn += composites.size();
 
   for ( const auto composite : composites ) {
@@ -123,8 +123,9 @@ output_t MassConstrainedVeloTrack::operator()( const LHCb::Particle::Range& comp
     auto velo_track_minus = std::make_unique<LHCb::Event::v1::Track>( *velo_particle_minus->proto()->track() );
 
     // function that handles fitting and output for both solutions
-    auto refit_track_and_add_outputs = [this, &track_output, &pp_output, &velo_output, &output, &track_pid, &lhcb](
-                                           auto velo_track, auto velo_particle, auto outp, auto p_probe, auto vtx ) {
+    auto refit_track_and_add_outputs = [this, &track_output, &pp_output, &velo_output, &output, &vertices_output,
+                                        &track_pid, &lhcb]( auto velo_track, auto velo_particle, auto outp,
+                                                            auto p_probe, auto vtx ) {
       for ( auto state : velo_track->states() ) {
         auto qop = velo_particle->charge() / p_probe;
         state->setQOverP( qop );
@@ -139,9 +140,10 @@ output_t MassConstrainedVeloTrack::operator()( const LHCb::Particle::Range& comp
       pp_output.insert( velo_proto.release() );
       outp->addToDaughters( velo_particle.get() );
       vtx->addToOutgoingParticles( velo_particle.get() );
-      outp->setEndVertex( vtx.release() );
+      outp->setEndVertex( vtx.get() );
       velo_output.insert( velo_particle.release() );
       output.insert( outp.release() );
+      vertices_output.insert( vtx.release() );
     };
     // Tests that these are the only relevant inefficiencies that the algorithm might introduce. I.e. the fraction of
     // failing fits is essentially 0.