Update EndVertex Definition for MCParticle

Event/MCEvent/include/Event/MCParticle.h

@@ -165,6 +165,8 @@ namespace LHCb {
     /// Retrieve (const)  Vector of pointers to decay vertices
     const SmartRefVector<LHCb::MCVertex>& endVertices() const;
 
+    const MCVertex* goodEndVertex() const;
+
     /// Update  Vector of pointers to decay vertices
     MCParticle& setEndVertices( SmartRefVector<LHCb::MCVertex> value );
 

Event/MCEvent/src/MCParticle.cpp

@@ -31,30 +31,57 @@ std::ostream& LHCb::MCParticle::fillStream( std::ostream& s ) const {
 }
 
 namespace LHCb {
+
   LinAlg::Vec3<SIMDWrapper::scalar::float_v> referencePoint( const MCParticle& p ) {
     auto const& pos = p.originVertex()->position();
     return LinAlg::Vec3<SIMDWrapper::scalar::float_v>{pos.x(), pos.y(), pos.z()};
   }
 
   /**
-   * @brief Get first end-vertex of MCParticle if it's from a decay or hadronic interaction
+   * @brief Get MCParticle's first end-vertex that likely destroyed the particle
+   *
+   * @return const LHCb::MCVertex*
+   * @note Choosing an endVertex is not trivial and might depend on the physics case.
+   * This implementation was presented in
+   * https://indico.cern.ch/event/1118300/contributions/4734729/attachments/2389881.
+   * It first tries to find the first vertex where the particle is actually destroyed.
+   * If no such vertex is found, it checks for an hadronic IA and returns this vertex or
+   * nullptr if nothing was found.
+   */
+  const MCVertex* MCParticle::goodEndVertex() const {
+    const auto& v            = endVertices();
+    const auto  destructiveV = std::find_if( v.begin(), v.end(), []( auto vtx ) {
+      switch ( vtx->type() ) {
+      case MCVertex::DecayVertex:
+      case MCVertex::OscillatedAndDecay:
+      case MCVertex::Annihilation:
+      case MCVertex::PairProduction:
+        return true;
+      default:
+        return false;
+      }
+    } );
+    if ( destructiveV != v.end() ) {
+      return destructiveV->data();
+    } else {
+      const auto hadronIAV =
+          std::find_if( v.begin(), v.end(), []( auto vtx ) { return vtx->type() == MCVertex::HadronicInteraction; } );
+      return hadronIAV != v.end() ? hadronIAV->data() : nullptr;
+    }
+  }
+
+  /**
+   * @brief Get position of a good endVertex for MCParticle
    *
    * @param p LHCb::MCParticle
    * @return LinAlg::Vec3<SIMDWrapper::scalar::float_v>
    * @note Returns LinAlg::Vec3<SIMDWrapper::scalar::float_v> filled with NaN if no end-vertex
    * found. This function is used by ThOr functors via ADL.
-   * Choosing an endVertex is not trivial and might depend on the physics case. This
-   * implementation is taken from
-   * https://gitlab.cern.ch/lhcb/Analysis/-/blob/v34r1/Phys/DecayTreeTupleMC/src/MCTupleToolKinematic.cpp#L112
    */
   LinAlg::Vec3<SIMDWrapper::scalar::float_v> endVertexPos( const MCParticle& p ) {
-    const auto& endVertices = p.endVertices();
-    const auto  mcV         = std::find_if( endVertices.begin(), endVertices.end(), []( auto vtx ) {
-      return vtx->type() == MCVertex::DecayVertex || vtx->type() == MCVertex::OscillatedAndDecay ||
-             vtx->type() == MCVertex::HadronicInteraction;
-    } );
-    if ( mcV != endVertices.end() ) {
-      const auto pos = ( *mcV )->position();
+    const auto* goodVtx = p.goodEndVertex();
+    if ( goodVtx ) {
+      const auto pos = goodVtx->position();
       return LinAlg::Vec3<SIMDWrapper::scalar::float_v>{pos.x(), pos.y(), pos.z()};
     } else {
       constexpr auto nan = std::numeric_limits<float>::quiet_NaN();