Merge branch 'SegmentFitVolVII' into 'main'

MuonRecoChainTester - Pointer based truth matching

See merge request atlas/athena!77743

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/python/MuonRecoChainTesterConfig.py

@@ -7,7 +7,7 @@ if __name__=="__main__":
     parser.set_defaults(nEvents = -1)
     parser.set_defaults(noMM=True)
     parser.set_defaults(noSTGC=True)
-    parser.set_defaults(outRootFile="HoughTransformTester.root")
+    parser.set_defaults(outRootFile="RecoChainTester.root")
     # parser.set_defaults(condTag="CONDBR2-BLKPA-2023-03")
     parser.set_defaults(inputFile=[
                                    "/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/MuonGeomRTT/R3SimHits.pool.root"

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/python/PatternTestConfig.py

@@ -165,6 +165,7 @@ def TrackTruthMatchCfg(flags):
     
     track_cols = ["MuonTracksR4", "MuonTracksFromHoughR4", "MuonSpectrometerTracks"]
     track_colstp = ["MuonSpectrometerTrackParticlesR4", "MuonSpectrometerTrackParticlesFromHoughR4", "MuonSpectrometerTrackParticles"]
+
     for trk in track_cols:
         result.merge(MuonDetailedTrackTruthMakerCfg(flags, name=f"MuonDetailedTruthTrkMaker{trk}",
                                                         TrackCollectionNames=[trk]))

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonRecoChainTester.cxx

@@ -4,12 +4,13 @@
 // Framework includes
 #include "MuonRecoChainTester.h"
 
-#include "MuonTesterTree/EventInfoBranch.h"
+#include "MuonTesterTree/MuonTesterTreeDict.h"
 #include "MuonTesterTree/TrackChi2Branch.h"
 #include "MuonPRDTest/SegmentVariables.h"
 #include "MuonPRDTest/ParticleVariables.h"
 #include "FourMomUtils/xAODP4Helpers.h"
 #include "AthContainers/ConstDataVector.h"
+#include "xAODTruth/xAODTruthHelpers.h"
 #include "StoreGate/ReadHandle.h"
 
 using namespace MuonVal;
@@ -17,20 +18,7 @@ using namespace MuonPRDTest;
 
 
 namespace {
-  template <class RefPartType, class SearchPartType>
-      const SearchPartType* findClosestParticle(const RefPartType* reference,
-                                                const DataVector<SearchPartType>& candidateContainer) {
-          const SearchPartType* best{nullptr};
-          for (const SearchPartType* candidate : candidateContainer) {
-              if (!best || xAOD::P4Helpers::deltaR2(reference, candidate) <
-                           xAOD::P4Helpers::deltaR2(reference, best)) {
-                  best = candidate;
-              }
-          }
-          return best;            
-      }
-  static const SG::Decorator<int> acc_truthMatched{"truthMatched"};
-
+    static const SG::Decorator<int> acc_truthMatched{"truthMatched"};
 }
 namespace MuonValR4{
     StatusCode MuonRecoChainTester::initialize() {
@@ -53,49 +41,36 @@ namespace MuonValR4{
 
         m_legacyTrks = std::make_shared<IParticleFourMomBranch>(m_tree, "LegacyMSTrks");
         m_legacyTrks->addVariable(std::make_shared<TrackChi2Branch>(*m_legacyTrks));
-        m_legacyTrks->addVariable<int>("truthMatched");
 
         m_TrksHoughR4 = std::make_shared<IParticleFourMomBranch>(m_tree, "HoughMSTrks");
         m_TrksHoughR4->addVariable(std::make_shared<TrackChi2Branch>(*m_TrksHoughR4));
-        m_TrksHoughR4->addVariable<int>("truthMatched");
 
         m_TrksSegmentR4 = std::make_shared<IParticleFourMomBranch>(m_tree, "MSTrksR4");
         m_TrksSegmentR4->addVariable(std::make_shared<TrackChi2Branch>(*m_TrksSegmentR4));
-        m_TrksSegmentR4->addVariable<int>("truthMatched");
 
         m_tree.addBranch(m_legacyTrks);
         m_tree.addBranch(m_TrksSegmentR4);
         m_tree.addBranch(m_TrksHoughR4);
         
         if (m_isMC) {
+            m_trkTruthLinks.emplace_back(m_legacyTrackKey, "truthParticleLink");
+            m_trkTruthLinks.emplace_back(m_TrackKeyHoughR4, "truthParticleLink");
+            m_trkTruthLinks.emplace_back(m_TrackKeyR4, "truthParticleLink");
+
             m_truthTrks = std::make_shared<IParticleFourMomBranch>(m_tree, "TruthMuons");
-            m_truthTrks->addVariable<int>("legacyMatched");
-            m_truthTrks->addVariable<int>("houghMatched");
-            m_truthTrks->addVariable<int>("r4Matched");
-            
+            BilateralLinkerBranch::connectCollections(m_legacyTrks, m_truthTrks, [](const xAOD::IParticle* trk){ 
+                                                        return xAOD::TruthHelpers::getTruthParticle(*trk); }, "truth", "LegacyMS");
+            BilateralLinkerBranch::connectCollections(m_TrksHoughR4, m_truthTrks, [](const xAOD::IParticle* trk){ 
+                                                            return xAOD::TruthHelpers::getTruthParticle(*trk); }, "truth", "HoughMS");
+            BilateralLinkerBranch::connectCollections(m_TrksSegmentR4, m_truthTrks, [](const xAOD::IParticle* trk){ 
+                                                                return xAOD::TruthHelpers::getTruthParticle(*trk); }, "truth", "MSTrksR4");
+        
             m_tree.addBranch(m_truthTrks);
         }
+        ATH_CHECK(m_trkTruthLinks.initialize());
         ATH_CHECK(m_tree.init(this));
         return StatusCode::SUCCESS;
     }
-    void MuonRecoChainTester::matchTrackToTruth(const xAOD::TruthParticle* truth,
-                                                const xAOD::TrackParticleContainer& tracks,
-                                                const SG::Decorator<int>& decRecoMatch,
-                                                const std::shared_ptr<MuonVal::IParticleFourMomBranch>& trkBranch) const {
-      constexpr double matchDR = 0.2;
-
-      decRecoMatch(*truth) = -1;
-      const xAOD::TrackParticle* closest = findClosestParticle(truth, tracks);
-
-      if (!closest || xAOD::P4Helpers::deltaR(truth, closest) > matchDR) {
-          return;
-      }
-      // Decorate the truth particle index to the reco particle
-      acc_truthMatched(*closest) = static_cast<int>(m_truthTrks->size());
-      trkBranch->push_back(closest);
-      /// Decorate the reco particle index to the truth particle
-      decRecoMatch(*truth) = trkBranch->find(closest); 
-    }
     void MuonRecoChainTester::fillBucketsPerStation(const MuonR4::SpacePointContainer& spContainer,
                                                     const StIdx station,
                                                     MuonVal::ScalarBranch<uint16_t>& outBranch) const{
@@ -107,11 +82,6 @@ namespace MuonValR4{
 
     StatusCode MuonRecoChainTester::execute() {
     
-      const SG::Decorator<int> acc_legacyMatched{"legacyMatched"};
-      const SG::Decorator<int> acc_houghMatched{"houghMatched"};
-      const SG::Decorator<int> acc_r4Matched{"r4Matched"};
-
-      
       const EventContext& ctx{Gaudi::Hive::currentContext()};
       SG::ReadHandle legacyTrks{m_legacyTrackKey, ctx};
       ATH_CHECK(legacyTrks.isPresent());
@@ -120,48 +90,27 @@ namespace MuonValR4{
       SG::ReadHandle trksR4{m_TrackKeyR4, ctx};
       ATH_CHECK(trksR4.isPresent());
       
+      ATH_MSG_DEBUG("Fill reconstructed tracks from "<<m_legacyTrackKey.fullKey());
       for (const xAOD::TrackParticle* trk : *legacyTrks) {
-          acc_truthMatched(*trk) = -1;
+        m_legacyTrks->push_back(trk);
       }
+      ATH_MSG_DEBUG("Fill reconstructed tracks from "<<m_TrackKeyR4.fullKey());
       for (const xAOD::TrackParticle* trk : *trksR4) {
-          acc_truthMatched(*trk) = -1;
+          m_TrksSegmentR4->push_back(trk);
       }
+      ATH_MSG_DEBUG("Fill reconstructed tracks from "<<m_TrackKeyHoughR4.fullKey());
       for (const xAOD::TrackParticle* trk : *trksFromHoughR4) {
-          acc_truthMatched(*trk) = -1;
-      }
-      
-      ConstDataVector<xAOD::TruthParticleContainer> truthParts{SG::VIEW_ELEMENTS};
+          m_TrksHoughR4->push_back(trk);
+      } 
+  
       if (!m_truthKey.empty()) {
           SG::ReadHandle readHandle{m_truthKey, ctx};
           ATH_CHECK(readHandle.isPresent());
+          ATH_MSG_DEBUG("Fill truth from "<<m_truthKey.fullKey());
           for (const xAOD::TruthParticle* truth : *readHandle) {
-              if (!truth->isMuon()) continue;
-              if (truth->status() != 1) continue;
-              truthParts.push_back(truth);
-          }
-      }
-      for (const xAOD::TruthParticle* truth : truthParts) {
-          matchTrackToTruth(truth, *legacyTrks, acc_legacyMatched, m_legacyTrks);
-          matchTrackToTruth(truth, *trksFromHoughR4, acc_houghMatched, m_TrksHoughR4);
-          matchTrackToTruth(truth, *trksR4, acc_r4Matched, m_TrksSegmentR4);
-          if (truth->eta() > 2. && acc_legacyMatched(*truth) != -1
-              && acc_r4Matched(*truth) == -1 ) {
-            ATH_MSG_VERBOSE("In event "<<ctx.eventID()<<" the new chain is inefficient. eta: "
-                          <<truth->eta()<<", pT: "<<truth->pt()
-                          <<" legacy: "<<acc_legacyMatched(*truth)<<", r4-hough: "<<
-                          acc_houghMatched(*truth)<<", r4 seg: "<<acc_r4Matched(*truth));
-          }
-          m_truthTrks->push_back(truth);
-      }
-      for (const xAOD::TrackParticle* trk : *legacyTrks) {
-          m_legacyTrks->push_back(trk);
-      }
-      for (const xAOD::TrackParticle* trk : *trksR4) {
-          m_TrksSegmentR4->push_back(trk);
+            m_truthTrks->push_back(truth);
+        }
       }
-      for (const xAOD::TrackParticle* trk : *trksFromHoughR4) {
-          m_TrksHoughR4->push_back(trk);
-      } 
       /** Fill the bucket summary counts */
       SG::ReadHandle spContainer{m_spacePointKey, ctx};
       ATH_CHECK(spContainer.isPresent());

MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MuonRecoChainTester.h

@@ -13,6 +13,7 @@
 #include "MuonTesterTree/MuonTesterTreeDict.h"
 
 #include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/ReadDecorHandleKeyArray.h"
 
 namespace MuonValR4{
 
@@ -26,17 +27,6 @@ namespace MuonValR4{
             virtual StatusCode finalize() override;
 
   private:
-        /** @brief Finds the closest track to the truth particle in terms of dR
-         *  @param truthPart: Pointer to the reference truth particle
-         *  @param tracks: MS-track collection where the closest track shall be searched
-         *  @param decRecoMatch: Decorator to decorate the track particle index in the output tree
-         *                       to the truth particle
-         *  @param trkBranch: Track branch collection to save the track particle  */
-        void matchTrackToTruth(const xAOD::TruthParticle* truthPart,
-                               const xAOD::TrackParticleContainer& tracks,
-                               const SG::Decorator<int>& decRecoMatch,
-                               const std::shared_ptr<MuonVal::IParticleFourMomBranch>& trkBranch) const;
-        
         using StIdx = Muon::MuonStationIndex::StIndex;
         /** @brief Counts how many buckets are in a particular station
          *  @param spContainer: Space point collection
@@ -86,7 +76,9 @@ namespace MuonValR4{
         SG::ReadHandleKey<xAOD::TrackParticleContainer> m_TrackKeyR4{this, "TrackKeyR4", "MuonSpectrometerTrackParticlesR4"};
 
         /** @brief Key to the truth particle collection */
-        SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthKey{this, "TruthKey", "TruthParticles"};
+        SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthKey{this, "TruthKey", "MuonTruthParticles"};
+        /** @brief Decoration dependency to the MS truth track links */
+        SG::ReadDecorHandleKeyArray<xAOD::TrackParticleContainer> m_trkTruthLinks{this, "TruthTrackLinks", {}};
         /** @brief Key to the space point container */
         SG::ReadHandleKey<MuonR4::SpacePointContainer> m_spacePointKey{this, "SpacePointContainer", "MuonSpacePoints"};
   

MuonSpectrometer/MuonValidation/MuonTesterTree/MuonTesterTree/IParticleFourMomBranch.h

@@ -124,8 +124,9 @@ namespace MuonVal {
     
         private:
             MuonTesterTree& m_parent;
-            std::string m_name;
+            std::string m_name{};
             bool m_init{false};
+            bool m_updated{false};
             VectorBranch<float> m_pt{m_parent.tree(), m_name +"_pt"};
             VectorBranch<float> m_eta{m_parent.tree(), m_name + "_eta"};
             VectorBranch<float> m_phi{m_parent.tree(), m_name + "_phi"};

MuonSpectrometer/MuonValidation/MuonTesterTree/MuonTesterTree/LinkerBranch.h

+/*
+  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef MUONTESTER_LINKERBRANCH_H
+#define MUONTESTER_LINKERBRANCH_H
+
+#include <MuonTesterTree/VectorBranch.h>
+#include <MuonTesterTree/IParticleFourMomBranch.h>
+
+namespace MuonVal{
+    /** @brief The linker branch is a helper class to establish linking indices
+     *         between dumped particle collections. Once the primary particle is
+     *         added to the collection, the linker branch calls the linker function
+     *         to find the related particle and to add it into it's forseen collection.
+     *         The position where the related particle is saved in the tree, is then 
+     *         saved by the linker branch */
+    class LinkerBranch : public VectorBranch<unsigned short>,
+                         virtual public IParticleDecorationBranch {
+    public:
+        /** @brief Abreviation of the pointer to the particle branch. */
+        using ParticleBranch_ptr  = std::shared_ptr<IParticleFourMomBranch>;
+
+        /** @brief Typedef of the linker function. The function takes as argument the primary
+         *         particle of interest and shall then return the particle to associate with */
+        using Linker_t = std::function<const xAOD::IParticle*(const xAOD::IParticle*)>;
+
+        /** @brief Standard constructor fo the LinkerBranch
+         *  @param parent: Particle branch collection where the linker branch is appended to. 
+         *                 The primary particle of interest is provided by this collection.
+         *  @param linkColl: Particle branch collection where the linked particle shall be appended
+         *  @param linker: Function to find the linked particle from the primary particle
+         *  @param altName: by default the name is <primaryColl>_linked<Secondary>. 
+         *                  Parse an alternative name */
+        LinkerBranch(IParticleFourMomBranch& parent,
+                     ParticleBranch_ptr linkColl,
+                     Linker_t  linker,
+                     const std::string& altName = "");
+        /** @brief Use the push back methods of the parent class */
+        using VectorBranch<unsigned short>::push_back;
+        /** @brief Interface methods to handle the particle */
+        void push_back(const xAOD::IParticle* p) override;
+        void push_back(const xAOD::IParticle& p) override;
+        void operator+=(const xAOD::IParticle* p) override;
+        void operator+=(const xAOD::IParticle& p) override;
+
+    private:
+        std::weak_ptr<IParticleFourMomBranch> m_linkColl;
+        Linker_t m_linkerFunc;
+    };
+
+    class BilateralLinkerBranch: public VectorBranch<unsigned short>,
+                                 virtual public IParticleDecorationBranch {
+        public:
+            using Linker_t = LinkerBranch::Linker_t;
+            using ParticleBranch_ptr = LinkerBranch::ParticleBranch_ptr;
+            /** @brief Use the push back methods of the parent class */
+            using VectorBranch<unsigned short>::push_back;
+            /** @brief Interface methods to handle the particle */
+            void push_back(const xAOD::IParticle* p) override;
+            void push_back(const xAOD::IParticle& p) override;
+            void operator+=(const xAOD::IParticle* p) override;
+            void operator+=(const xAOD::IParticle& p) override;
+        
+            bool fill(const EventContext& ctx) override;
+
+            static bool connectCollections(ParticleBranch_ptr primColl,
+                                           ParticleBranch_ptr secondColl,
+                                           Linker_t fromPrimToSec,
+                                           const std::string& altPrimName ="",
+                                           const std::string& altSecName = "");
+        private:
+            BilateralLinkerBranch(IParticleFourMomBranch& bilatColl,
+                                  ParticleBranch_ptr primColl,
+                                  Linker_t linker,
+                                  const std::string& altName);
+            const IParticleFourMomBranch& m_parent;
+            std::weak_ptr<IParticleFourMomBranch> m_linkColl;
+            Linker_t m_linkerFunc;
+    
+
+    };
+}
+#endif
\ No newline at end of file

MuonSpectrometer/MuonValidation/MuonTesterTree/MuonTesterTree/MuonTesterTreeDict.h

 /*
-Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
 */
 #ifndef MUONTESTER_MUONTESTERTREEDICT_H
 #define MUONTESTER_MUONTESTERTREEDICT_H
-
 #include <MuonTesterTree/ArrayBranch.h>
 #include <MuonTesterTree/AuxElementBranch.h>
 #include <MuonTesterTree/CoordTransformBranch.h>
 #include <MuonTesterTree/EventHashBranch.h>
+#include <MuonTesterTree/EventIDBranch.h>
 #include <MuonTesterTree/EventInfoBranch.h>
 #include <MuonTesterTree/FourVectorBranch.h>
 #include <MuonTesterTree/IMuonTesterBranch.h>
 #include <MuonTesterTree/IParticleFourMomBranch.h>
 #include <MuonTesterTree/IdentifierBranch.h>
+#include <MuonTesterTree/LinkerBranch.h>
 #include <MuonTesterTree/MatrixBranch.h>
 #include <MuonTesterTree/MuonTesterBranch.h>
 #include <MuonTesterTree/MuonTesterTree.h>
-#include <MuonTesterTree/MuonTesterTreeDict.h>
 #include <MuonTesterTree/ScalarBranch.h>
 #include <MuonTesterTree/SetBranch.h>
 #include <MuonTesterTree/ThreeVectorBranch.h>
+#include <MuonTesterTree/TrackChi2Branch.h>
 #include <MuonTesterTree/TwoVectorBranch.h>
 #include <MuonTesterTree/VectorBranch.h>
 #endif
+

MuonSpectrometer/MuonValidation/MuonTesterTree/Root/IParticleFourMomBranch.cxx

@@ -25,7 +25,7 @@ namespace MuonVal{
         return m_cached_particles;
     }
 
-    size_t IParticleFourMomBranch::size() const { return m_pt.size(); }
+    size_t IParticleFourMomBranch::size() const { return m_updated ? m_cached_particles.size() : 0; }
     std::string IParticleFourMomBranch::name() const { return m_name; }
     const TTree* IParticleFourMomBranch::tree() const { return m_pt.tree(); }
     TTree* IParticleFourMomBranch::tree() { return m_pt.tree(); }
@@ -33,9 +33,16 @@ namespace MuonVal{
     void IParticleFourMomBranch::operator+=(const xAOD::IParticle& p) {push_back(p); }
     void IParticleFourMomBranch::push_back(const xAOD::IParticle& p) { push_back(&p); }
     void IParticleFourMomBranch::push_back(const xAOD::IParticle* p) {
+        if (!m_updated) {
+            m_cached_particles.clear();
+            m_updated = true;
+        }
         /// Avoid that the particle is added twice to the tree
-        if (!p || find(p) < size())
+        if (!p || find(p) < size()){
+            if (p) ATH_MSG_VERBOSE("Rejected particle ("<<p<<") "<<p->pt()<<", "<<p->eta()<<", "<<p->phi()<<". Size: "<<size());
             return;
+        }
+
         m_cached_particles.push_back(p);
         m_pt.push_back(p->pt() * MeVtoGeV);
         m_eta.push_back(p->eta());
@@ -47,6 +54,8 @@ namespace MuonVal{
         } else if (p->type() == xAOD::Type::ObjectType::TrackParticle){
             q = static_cast<const xAOD::TrackParticle*>(p)->charge();
         }
+        ATH_MSG_VERBOSE("New particle ("<<p<<") "<<p->pt()<<", "<<p->eta()<<", "<<p->phi()<<", q: "<<q<<". Size: "<<size());
+
         m_q.push_back(q);
         for (const auto& var : m_variables) {
             var->push_back(p);
@@ -60,10 +69,15 @@ namespace MuonVal{
         if (!p) {
             return dummyIdx;
         }
-        return find([p](const xAOD::IParticle* cached) { return p == cached; });
+        return find([p](const xAOD::IParticle* cached) { 
+            return p == cached; 
+        });
     }
     size_t IParticleFourMomBranch::find(std::function<bool(const xAOD::IParticle*)> func) const {
-        size_t j = 0;
+        if (!m_updated) {
+            return dummyIdx;
+        }
+        size_t j{0};
         for (const xAOD::IParticle* p : m_cached_particles) {
             if (func(p)) {
                 return j;
@@ -73,6 +87,9 @@ namespace MuonVal{
         return dummyIdx;
     }
     bool IParticleFourMomBranch::fill(const EventContext& ctx) {
+        if (!m_updated) {
+            m_cached_particles.clear();
+        }
         if (!m_pt.fill(ctx) || !m_eta.fill(ctx) || !m_phi.fill(ctx) || !m_e.fill(ctx) || !m_q.fill(ctx)){
             return false;
         }
@@ -82,7 +99,7 @@ namespace MuonVal{
             }
         }
         m_init = true;
-        m_cached_particles.clear();
+        m_updated = false;
         return true;
     }
     bool IParticleFourMomBranch::initialized() const {

MuonSpectrometer/MuonValidation/MuonTesterTree/Root/LinkerBranch.cxx

+/*
+  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
+*/
+#include <MuonTesterTree/LinkerBranch.h>
+#include <format>
+namespace MuonVal{
+    LinkerBranch::LinkerBranch(IParticleFourMomBranch& parent,
+                               ParticleBranch_ptr linkColl,
+                               Linker_t  linker,
+                               const std::string& altName):
+        VectorBranch<unsigned short>{parent.tree(), 
+                                     std::format("{:}_{:}Link", parent.name(), altName.empty() ? linkColl->name() : altName)},
+        m_linkColl{linkColl},
+        m_linkerFunc{linker} {}
+
+    void LinkerBranch::operator+=(const xAOD::IParticle* p) {
+        push_back(p);
+    }
+    void LinkerBranch::operator+=(const xAOD::IParticle& p) {
+        push_back(p);
+    }
+    void LinkerBranch::push_back(const xAOD::IParticle& p) {
+        push_back(&p);
+    }
+    void LinkerBranch::push_back(const xAOD::IParticle* p) {
+        const xAOD::IParticle* related = m_linkerFunc(p);
+        if (related) {
+            ParticleBranch_ptr linkColl = m_linkColl.lock();
+            linkColl->push_back(related);
+            VectorBranch<unsigned short>::push_back(linkColl->find(related));
+        } else {
+            VectorBranch<unsigned short>::push_back(-1);
+        }
+    }
+
+    
+    bool BilateralLinkerBranch::connectCollections(ParticleBranch_ptr primColl,
+                                                   ParticleBranch_ptr secondColl,
+                                                   Linker_t fromPrimToSec,
+                                                   const std::string& altPrimName,
+                                                   const std::string& altSecName) {
+        return primColl->addVariable(std::make_unique<LinkerBranch>(*primColl, secondColl, fromPrimToSec, altPrimName)) &&
+               secondColl->addVariable(std::unique_ptr<IParticleDecorationBranch>{new BilateralLinkerBranch(*secondColl, primColl, fromPrimToSec, altSecName)});
+    }
+                     
+    BilateralLinkerBranch::BilateralLinkerBranch(IParticleFourMomBranch& parent,
+                                                 ParticleBranch_ptr primColl,
+                                                 Linker_t linker,
+                                                 const std::string& altName):
+        VectorBranch<unsigned short>{parent.tree(), 
+                                     std::format("{:}_{:}Link", parent.name(), altName.empty() ? primColl->name() : altName)},
+        m_parent{parent},
+        m_linkColl{primColl},
+        m_linkerFunc{linker} {
+        setDefault(-1);
+    }
+
+    void BilateralLinkerBranch::push_back(const xAOD::IParticle* /* p*/) {
+    }
+    void BilateralLinkerBranch::push_back(const xAOD::IParticle& p) {
+        push_back(&p);
+    }
+    void BilateralLinkerBranch::operator+=(const xAOD::IParticle* p){
+        push_back(p);
+    }
+    void BilateralLinkerBranch::operator+=(const xAOD::IParticle& p){
+        push_back(p);
+    }          
+    bool BilateralLinkerBranch::fill(const EventContext& ctx) {
+        ATH_MSG_VERBOSE("Fill "<<name()<<", size: "<<m_parent.size()<<", "<<m_parent.name());
+        if (m_parent.size()) {
+            /** Allocate the memory */
+            get(m_parent.size() -1);
+            ParticleBranch_ptr linkColl = m_linkColl.lock();
+            const std::vector<const xAOD::IParticle*>& linkeMe = linkColl->getCached();
+            for (std::size_t primToSec = 0 ; primToSec < linkeMe.size(); ++primToSec) {
+                const size_t linkIdx = m_parent.find(m_linkerFunc(linkeMe[primToSec]));
+                if (linkIdx < size()) {
+                    get(linkIdx) = primToSec;
+                }
+            }
+        }
+        return VectorBranch<unsigned short>::fill(ctx);
+    }
+}
\ No newline at end of file