diff --git a/Phys/DaVinciFilters/CMakeLists.txt b/Phys/DaVinciFilters/CMakeLists.txt
index 0979ea3f09c95deaf4ecb2271ccdc56fdc468749..5a39c90e9143554e1f4ec0355044c0af43d43297 100644
--- a/Phys/DaVinciFilters/CMakeLists.txt
+++ b/Phys/DaVinciFilters/CMakeLists.txt
@@ -17,6 +17,7 @@ gaudi_add_module(DaVinciFilters
     SOURCES
         src/Velo2LongB2jpsiKTrackEff_VeloTrackEffFilter.cpp
         src/Velo2Long_Ds2PhiPi_TrackEffFilter.cpp
+	src/Velo2Long_L2ProtonPi_TrackEffFilter.cpp
 	    src/KSLongVeloFilter.cpp
     LINK
         Boost::headers
diff --git a/Phys/DaVinciFilters/src/Velo2Long_L2ProtonPi_TrackEffFilter.cpp b/Phys/DaVinciFilters/src/Velo2Long_L2ProtonPi_TrackEffFilter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d5a2b21ebe4d9d7b6c5edf39b3244be271b21717
--- /dev/null
+++ b/Phys/DaVinciFilters/src/Velo2Long_L2ProtonPi_TrackEffFilter.cpp
@@ -0,0 +1,263 @@
+/*****************************************************************************\
+* (c) Copyright 2023 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.                                       *
+\*****************************************************************************/
+#include "Core/FloatComparison.h"
+#include "Event/Particle.h"
+#include "Event/RecVertex.h"
+#include "Kernel/IParticleDescendants.h"
+#include "LHCbAlgs/Transformer.h"
+#include "LHCbMath/GeomFun.h"
+#include "LHCbMath/Line.h"
+#include "LHCbMath/MatrixTransforms.h"
+#include "LHCbMath/MatrixUtils.h"
+
+#include "Kernel/IParticle2State.h"
+#include "Kernel/IParticleCombiner.h"
+
+#include "DetDesc/DetectorElement.h"
+#include "DetDesc/GenericConditionAccessorHolder.h"
+
+#include <functional>
+#include <string>
+
+namespace {
+  constexpr double protonMass  = 938.272 * Gaudi::Units::MeV;
+  constexpr double protonMass2 = protonMass * protonMass;
+  constexpr double LMass2      = 1115.683 * 1115.683 * Gaudi::Units::MeV * Gaudi::Units::MeV;
+
+  using in_pv           = LHCb::RecVertex::Range;
+  using in_parts        = LHCb::Particle::Range;
+  using output_t        = std::tuple<LHCb::Particle::Container,
+                              LHCb::Particle::Container,      // updated probe particles
+                              LHCb::Vertex::Container,        // refitted vertices
+                              LHCb::ProtoParticle::Container, // updated probe protos
+                              LHCb::Track::Container          // updated tracks
+                              >;
+  using filter_output_t = std::tuple<bool, LHCb::Particle::Container,
+                                     LHCb::Particle::Container,      // updated probe particles
+                                     LHCb::Vertex::Container,        // refitted vertices
+                                     LHCb::ProtoParticle::Container, // updated probe protos
+                                     LHCb::Track::Container>;
+  using base_t =
+      LHCb::Algorithm::MultiTransformerFilter<output_t( in_parts const&, in_pv const&, const DetectorElement& ),
+                                              LHCb::Algorithm::Traits::usesConditions<DetectorElement>>;
+
+  double impactParameterSquared( const Gaudi::XYZPoint& vertex, const Gaudi::XYZPoint& point,
+                                 const Gaudi::XYZVector& dir ) {
+    // Direction does not need to be normalised.
+    const auto distance = point - vertex;
+    return ( distance - dir * 1. / dir.Mag2() * ( distance.Dot( dir ) ) ).Mag2();
+  }
+
+  auto bpvAndMinImpactParameterSquared( const in_pv& pvs, const int nTracksMin, const Gaudi::XYZPoint& point,
+                                        const Gaudi::XYZVector& dir ) {
+    // Direction does not need to be normalised.
+    return std::accumulate( pvs.begin(), pvs.end(), std::pair<LHCb::RecVertex const*, double>{ nullptr, -1 },
+                            [&]( auto r, const auto& pv ) {
+                              if ( ( pv->nDoF() + 3 ) > 2 * nTracksMin ) {
+                                double ip2 = impactParameterSquared( pv->position(), point, dir );
+                                if ( !r.first || ip2 < r.second ) return std::pair{ pv, ip2 };
+                              }
+                              return r;
+                            } );
+  }
+
+  double transverseMomentumDir( const Gaudi::XYZVector& mom, const Gaudi::XYZVector& dir ) {
+    const double dmag2 = dir.Mag2();
+    if ( LHCb::essentiallyZero( dmag2 ) ) { return mom.R(); }
+    const Gaudi::XYZVector perp = mom - dir * ( mom.Dot( dir ) / dmag2 );
+    return perp.R();
+  }
+
+  Gaudi::XYZVector get_ip_vector( const Gaudi::XYZPoint& vertex_point, const Gaudi::XYZPoint& point,
+                                  const Gaudi::XYZVector& direction ) {
+    const Gaudi::Math::Line my_line( point, direction );
+    return Gaudi::Math::closestPoint( vertex_point, my_line ) - vertex_point;
+  }
+
+  double get_ip_perp( const Gaudi::XYZPoint& vertex_point, const Gaudi::XYZPoint& point,
+                      const Gaudi::XYZVector& direction, const Gaudi::XYZVector& momentum_vector_A,
+                      const Gaudi::XYZVector& momentum_vector_B ) {
+    const auto&      ipVector = get_ip_vector( vertex_point, point, direction );
+    Gaudi::XYZVector perp     = momentum_vector_A.Unit().Cross( momentum_vector_B.Unit() );
+    return perp.Dot( ipVector );
+  }
+
+} // namespace
+
+/*
+
+  using output_t        = std::tuple<LHCb::Particle::Container,
+                                  LHCb::Particle::Container,  // updated probe particles
+                                  LHCb::Vertex::Container, // refitted vertices
+                                  LHCb::ProtoParticle::Container, // updated probe protos
+                                  LHCb::Track::Container // updated tracks*/
+class Velo2Long_L2ProtonPi_TrackEffFilter final : public base_t {
+public:
+  Velo2Long_L2ProtonPi_TrackEffFilter( const std::string& name, ISvcLocator* pSvc )
+      : base_t( name, pSvc,
+                { KeyValue{ "InputParticle", "" }, KeyValue{ "InputPVs", "" },
+                  KeyValue{ "StandardGeometryTop", LHCb::standard_geometry_top } },
+                { KeyValue{ "OutputParticles", "" }, KeyValue{ "OutputProbeParticles", "" },
+                  KeyValue{ "OutputVertices", "" }, KeyValue{ "OutputProbePotos", "" },
+                  KeyValue{ "OutputProbeTracks", "" } } ) {}
+
+  filter_output_t operator()( const in_parts& parts, const in_pv& vertices,
+                              const DetectorElement& lhcb ) const override {
+
+    filter_output_t all_output;
+    auto& [event_accepted, out_heads, out_probe_particles, out_vertices, out_probe_protos, out_probe_tracks] =
+        all_output;
+
+    ++m_eventCount;
+    m_candidateCount += parts.size();
+    for ( const LHCb::Particle* orig_p : parts ) {
+      if ( orig_p->daughters().size() < 2 ) continue;
+
+      const LHCb::Particle* orig_probe_particle = ( orig_p->daughters() ).at( m_probe_index );
+      const LHCb::Particle* orig_tag_particle   = ( orig_p->daughters() ).at( m_tag_index );
+
+      if ( orig_probe_particle == nullptr || orig_tag_particle == nullptr ) continue;
+      if ( orig_p->endVertex() == nullptr ) continue;
+
+      auto my_probe_momentum               = orig_probe_particle->momentum();
+      auto mass_constrained_probe_momentum = orig_probe_particle->momentum();
+
+      const auto Etag   = orig_tag_particle->momentum().E();
+      const auto pprobe = std::sqrt( orig_probe_particle->momentum().Vect().Mag2() );
+      const auto ptagcostheta =
+          orig_tag_particle->momentum().Vect().Dot( orig_probe_particle->momentum().Vect() ) / pprobe;
+      const auto DeltaMSquared = LMass2 - protonMass2 - orig_tag_particle->momentum().M2();
+
+      auto my_new_momentum = 0.5 * ( DeltaMSquared ) / ( Etag - ptagcostheta );
+
+      auto Eprobe     = std::sqrt( 1. + ( protonMass / my_new_momentum ) * ( protonMass / my_new_momentum ) );
+      my_new_momentum = DeltaMSquared / ( 2.0 * ( Etag * Eprobe - ptagcostheta ) );
+      Eprobe          = std::sqrt( 1. + ( protonMass / my_new_momentum ) * ( protonMass / my_new_momentum ) );
+      my_new_momentum = DeltaMSquared / ( 2.0 * ( Etag * Eprobe - ptagcostheta ) );
+      Eprobe          = std::sqrt( 1. + ( protonMass / my_new_momentum ) * ( protonMass / my_new_momentum ) );
+      my_new_momentum = DeltaMSquared / ( 2.0 * ( Etag * Eprobe - ptagcostheta ) );
+
+      const auto scale_MConstr = my_new_momentum / pprobe;
+
+      mass_constrained_probe_momentum.SetPxPyPzE(
+          scale_MConstr * mass_constrained_probe_momentum.x(), scale_MConstr * mass_constrained_probe_momentum.y(),
+          scale_MConstr * mass_constrained_probe_momentum.z(),
+          std::sqrt( protonMass2 + scale_MConstr * scale_MConstr * mass_constrained_probe_momentum.P2() ) );
+
+      auto [bpv, ip2] = bpvAndMinImpactParameterSquared(
+          vertices, 3, orig_p->endVertex()->position(),
+          ( orig_p->momentum() - my_probe_momentum + mass_constrained_probe_momentum ).Vect() );
+      if ( bpv == nullptr ) {
+        ++m_noBPVfound;
+        continue;
+      }
+
+      const double ip_perp =
+          get_ip_perp( bpv->position(), orig_p->endVertex()->position(),
+                       ( orig_p->momentum() - my_probe_momentum + mass_constrained_probe_momentum ).Vect(),
+                       mass_constrained_probe_momentum.Vect(), orig_tag_particle->momentum().Vect() );
+
+      Gaudi::LorentzVector PV_constrained_probe_momentum( orig_p->momentum() );
+
+      const LHCb::VertexBase* vx         = orig_p->endVertex();
+      const Gaudi::XYZVector  dir        = vx->position() - bpv->position();
+      auto                    prb_p_PERP = transverseMomentumDir( orig_probe_particle->momentum().Vect(), dir );
+      auto                    tag_p_perp = transverseMomentumDir( orig_tag_particle->momentum().Vect(), dir );
+
+      if ( std::abs( prb_p_PERP ) < 1 ) continue;
+
+      // protect against corrupt data
+      if ( orig_probe_particle->proto() == nullptr || orig_probe_particle->proto()->track() == nullptr ) continue;
+
+      const auto scalePVConst = std::abs( ( tag_p_perp ) / ( prb_p_PERP ) );
+      PV_constrained_probe_momentum.SetPxPyPzE(
+          scalePVConst * my_probe_momentum.x(), scalePVConst * my_probe_momentum.y(),
+          scalePVConst * my_probe_momentum.z(),
+          std::sqrt( protonMass2 + scalePVConst * scalePVConst * my_probe_momentum.P2() ) );
+      Gaudi::LorentzVector PV_constrained_total_momentum =
+          PV_constrained_probe_momentum + orig_tag_particle->momentum();
+
+      if ( PV_constrained_total_momentum.M() > m_pvconstr_mass_min &&
+           PV_constrained_total_momentum.M() < m_pvconstr_mass_max &&
+           PV_constrained_probe_momentum.P() > m_probe_p_min && PV_constrained_probe_momentum.Pt() > m_probe_pt_min &&
+           PV_constrained_probe_momentum.Pt() < m_probe_pt_max && std::abs( ip_perp ) < ip_perp_max &&
+           std::sqrt( ip2 ) < ip_max ) {
+
+        auto   new_probe_track = std::make_unique<LHCb::Event::v1::Track>( *( orig_probe_particle->proto()->track() ) );
+        double qop             = -1. * (double)( orig_tag_particle->charge() ) / PV_constrained_probe_momentum.P();
+        new_probe_track->firstState().setQOverP( qop );
+
+        // assume a generous 30% uncertainty on the momentum
+        new_probe_track->firstState().setErrQOverP2( std::pow( 0.30 * qop, 2 ) );
+
+        auto new_probe_proto = std::unique_ptr<LHCb::ProtoParticle>( orig_probe_particle->proto()->clone() );
+        new_probe_proto->setTrack( new_probe_track.get() );
+
+        auto new_probe_particle = std::unique_ptr<LHCb::Particle>( orig_probe_particle->clone() );
+        new_probe_particle->setProto( new_probe_proto.get() );
+        m_p2s->state2Particle( new_probe_track->firstState(), *new_probe_particle ).ignore();
+
+        LHCb::Particle::ConstVector children;
+        children.reserve( 2 ); // Reserve space for 2 children
+        children.push_back( orig_tag_particle );
+        children.push_back( new_probe_particle.get() );
+
+        auto vtx      = std::unique_ptr<LHCb::Vertex>( orig_p->endVertex()->clone() );
+        auto new_part = std::make_unique<LHCb::Particle>( orig_p->particleID() );
+        auto combSc   = m_combiner->combine( children, *new_part, *vtx, *lhcb.geometry() );
+        if ( !( combSc.isSuccess() ) ) {
+          ++m_msg_fit_failure;
+          continue; // all new objects should be unique_ptrs and cleaned up
+        }
+
+        new_part->setPV( bpv );
+
+        out_heads.insert( new_part.release() );
+        out_vertices.insert( vtx.release() );
+        out_probe_particles.insert( new_probe_particle.release() );
+        out_probe_protos.insert( new_probe_proto.release() );
+        out_probe_tracks.insert( new_probe_track.release() );
+
+        ++m_newpartCounts;
+      }
+    }
+
+    event_accepted = !( out_heads.empty() );
+
+    return all_output;
+  };
+
+private:
+  // cuts
+  Gaudi::Property<double> m_probe_p_min{ this, "PVConstrainedProbePMin", 9000. * Gaudi::Units::MeV };
+  Gaudi::Property<double> m_probe_pt_min{ this, "PVConstrainedProbePtMin", 300. * Gaudi::Units::MeV };
+  Gaudi::Property<double> m_probe_pt_max{ this, "PVConstrainedProbePtMax", 100000. * Gaudi::Units::MeV };
+
+  Gaudi::Property<double> m_pvconstr_mass_min{ this, "PVConstrainedMassMin", 965 * Gaudi::Units::MeV };
+  Gaudi::Property<double> m_pvconstr_mass_max{ this, "PVConstrainedMassMax", 1225 * Gaudi::Units::MeV };
+
+  Gaudi::Property<double> ip_perp_max{ this, "IPperpendicularMax", 2 * Gaudi::Units::mm };
+  Gaudi::Property<double> ip_max{ this, "IPMax", 10 * Gaudi::Units::mm };
+
+  Gaudi::Property<unsigned int> m_probe_index{ this, "ProbeIndex", 1 };
+  Gaudi::Property<unsigned int> m_tag_index{ this, "TagIndex", 0 };
+
+  ToolHandle<IParticleCombiner> m_combiner{ this, "ParticleCombiner", "ParticleVertexFitter" };
+  ToolHandle<IParticle2State>   m_p2s{ this, "Particle2State", "Particle2State" };
+
+  mutable Gaudi::Accumulators::Counter<> m_eventCount{ this, "Events" };
+  mutable Gaudi::Accumulators::Counter<> m_candidateCount{ this, "Input Particles" };
+  mutable Gaudi::Accumulators::Counter<> m_newpartCounts{ this, "Accepted Particles" };
+  mutable Gaudi::Accumulators::Counter<> m_noBPVfound{ this, "Particles without BPV" };
+  mutable Gaudi::Accumulators::Counter<> m_msg_fit_failure{ this, "Failed vertex fit" };
+};
+
+DECLARE_COMPONENT( Velo2Long_L2ProtonPi_TrackEffFilter )