From 2059d1589e123ce57174cc5491dfe4da5e7f8255 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Wed, 27 Apr 2022 23:08:58 +0200
Subject: [PATCH 01/14] first working hit search algo

---
 Pr/PrPixel/CMakeLists.txt                   |   2 +
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 192 ++++++++++++++++++++
 2 files changed, 194 insertions(+)
 create mode 100644 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp

diff --git a/Pr/PrPixel/CMakeLists.txt b/Pr/PrPixel/CMakeLists.txt
index 7f56490f906..ce0a95d469b 100644
--- a/Pr/PrPixel/CMakeLists.txt
+++ b/Pr/PrPixel/CMakeLists.txt
@@ -20,6 +20,7 @@ gaudi_add_module(PrPixel
         src/VeloClusterTrackingSIMD.cpp
         src/VeloKalman.cpp
         src/PrVPHitsToVPLightClusters.cpp
+        src/VeloHeavyFlavourTracking.cpp
     LINK
         Gaudi::GaudiAlgLib
         LHCb::DAQEventLib
@@ -32,4 +33,5 @@ gaudi_add_module(PrPixel
         LHCb::VPKernel
         LHCb::VPDAQLib
         Rec::PrKernel
+        Rec::SelKernelLib
 )
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
new file mode 100644
index 00000000000..eb5f19e3886
--- /dev/null
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -0,0 +1,192 @@
+/*****************************************************************************\
+* (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.                                       *
+\*****************************************************************************/
+
+#include "DetDesc/GenericConditionAccessorHolder.h"
+#include "DetDesc/IConditionDerivationMgr.h"
+#include "Event/Particle.h"
+#include "Event/PrHits.h"
+#include "Event/PrimaryVertices.h"
+#include "Event/RawEvent.h"
+#include "Event/RelationTable.h"
+#include "GaudiAlg/Transformer.h"
+#include "Kernel/LHCbID.h"
+#include "LHCbMath/MatVec.h"
+#include "SelKernel/Utilities.h"
+#include "SelKernel/VertexRelation.h"
+#include "VPDet/DeVP.h"
+#include "VPDet/VPDetPaths.h"
+#include <vector>
+
+namespace LHCb::Pr::Velo {
+
+  namespace {
+
+    namespace VeloInfo {
+      int constexpr NSensors = 208;
+    } // namespace VeloInfo
+
+    namespace TrackInfo {
+      struct nInSensor : LHCb::Event::int_field {};
+      struct hasHitInFirstSensor : LHCb::Event::int_field {};
+      struct mCorr : LHCb::Event::float_field {};
+      struct cosAngle : LHCb::Event::float_field {};
+    } // namespace TrackInfo
+
+    // in/outputs
+    using Composites = LHCb::Particles;
+    using VeloHits   = LHCb::Pr::VP::Hits;
+    using PVs        = LHCb::Event::PV::PrimaryVertexContainer;
+    using Output     = LHCb::Event::RelationTable1D<Composites, TrackInfo::nInSensor, TrackInfo::hasHitInFirstSensor,
+                                                TrackInfo::mCorr, TrackInfo::cosAngle>;
+
+    // miscellaneous
+    using HitPos   = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>;
+    using SensorID = LHCb::Detector::VPChannelID::SensorID;
+
+    // to store relevant sensor information
+    struct SensorInfo {
+      float       z;
+      SensorID    id;
+      friend bool operator<( SensorInfo const& lhs, SensorInfo const& rhs ) {
+        return std::pair{lhs.z, lhs.id} < std::pair{rhs.z, rhs.id};
+      }
+    };
+    using SensorInfos = std::vector<SensorInfo>;
+
+    // to store relevant hit information
+    struct HitInfo {
+      float       d;
+      HitPos      position;
+      friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.d < rhs.d; }
+    };
+    using HitInfos = std::vector<HitInfo>;
+
+    // functions
+    template <typename Position3, typename Vector3>
+    HitInfos getHitsForFirstCrossedSensor( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV,
+                                           Vector3 const& tvec, float max_distance = 1. * Gaudi::Units::mm ) {
+      HitInfos poshits;
+      if ( !sensors.size() ) return poshits;
+      // get first sensor (in z) and its reference point
+      auto firstsensor = sensors.front();
+      auto ref         = PV + ( firstsensor.z - PV.Z() ) * tvec;
+      // look if this sensor has hits and if they are in search window
+      for ( auto hit : hits.scalar() ) {
+        auto id = LHCb::Detector::VPChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
+        if ( id.sensor() == firstsensor.id ) {
+          auto pos  = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
+          auto diff = ( pos.vec3() - ref ).mag();
+          if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3()} ); }
+        }
+      }
+      return poshits;
+    }
+
+  } // namespace
+
+  /**
+   * @class VeloHeavyFlavourTracking
+   * Search for hits of charged b/c hadrons based on primary and secondary vertices
+   *
+   * @author Maarten van Veghel (Nikhef, RUG)
+   */
+
+  class VeloHeavyFlavourTracking
+      : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const& ),
+                                              LHCb::DetDesc::usesConditions<DeVP>> {
+  public:
+    // standard constructor
+    VeloHeavyFlavourTracking( const std::string& name, ISvcLocator* pSvcLocator )
+        : Transformer( name, pSvcLocator,
+                       {KeyValue{"Composites", ""}, KeyValue{"PVs", ""}, KeyValue{"Hits", ""},
+                        KeyValue{"DeVP", LHCb::Det::VP::det_path}},
+                       {KeyValue{"OutputLocation", ""}} ) {}
+
+    // main execution
+    Output operator()( Composites const& composites, PVs const& pvs, VeloHits const& hits,
+                       DeVP const& velo ) const override {
+      // get z positions of sensors (should be in initialize TODO)
+      SensorInfos info_sensors( VeloInfo::NSensors );
+      auto        getSensorZs = [&info_sensors]( auto const& sensor ) {
+        info_sensors.push_back( {(float)sensor.z(), sensor.sensorNumber()} );
+      };
+      velo.runOnAllSensors( getSensorZs );
+      std::sort( info_sensors.begin(), info_sensors.end() );
+
+      // output
+      Output output_table( &composites );
+      output_table.reserve( composites.size() );
+
+      for ( auto const composite : composites ) {
+        // get SV and PV information
+        using Sel::Utils::endVertexPos;
+        auto const& endVertex = endVertexPos( *composite );
+        auto const& PV        = endVertexPos( Sel::getBestPV( *composite, pvs ) );
+        auto const  diff      = endVertex - PV;
+        auto const  tvec      = diff / diff.Z();
+
+        // search window crosses sensors?
+        int         nInSensor = 0;
+        SensorInfos crossedSensors;
+        auto        pos = PV;
+        for ( auto sinfo : info_sensors ) {
+          if ( sinfo.z < pos.Z() ) continue;
+          if ( sinfo.z > endVertex.Z() ) break;
+          // if in range look if search window intersects with sensor
+          pos = pos + ( sinfo.z - pos.Z() ) * tvec;
+          if ( velo.sensor( sinfo.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) {
+            nInSensor++;
+            crossedSensors.push_back( {pos.Z().cast(), sinfo.id} );
+          }
+        }
+
+        // look for hits in relevant sensors
+        // TODO maybe check all sensors in this plane, in case the search window is at the edge?
+        auto hits_first_sensor = getHitsForFirstCrossedSensor( crossedSensors, hits, PV, tvec );
+        std::sort( hits_first_sensor.begin(), hits_first_sensor.end() );
+        bool hasHitInFirstSensor = hits_first_sensor.size() > 0;
+
+        // determine mother particle direction based on closest hit in first crossed sensor
+        // or using SV-PV if not available
+        auto direction =
+            hasHitInFirstSensor ? ( hits_first_sensor.front().position - PV ).normalize() : diff.normalize();
+
+        // obtain corrected mass and opening angle
+        using Sel::Utils::mass2;
+        using Sel::Utils::threeMomentum;
+        using std::sqrt;
+        auto const mom      = threeMomentum( *composite );
+        auto const perp     = mom - direction * dot( mom, direction );
+        auto const pt       = perp.mag();
+        auto const mcorr    = pt + sqrt( mass2( *composite ) + pt * pt );
+        auto const cosangle = dot( mom, direction ) / mom.mag();
+
+        // statistics
+        m_nInSensor += nInSensor;
+        if ( nInSensor > 0 ) m_nHitCandidates += hits_first_sensor.size();
+
+        // set output
+        output_table.emplace_back<SIMDWrapper::InstructionSet::Scalar>().set( composite->key(), nInSensor,
+                                                                              hasHitInFirstSensor, mcorr, cosangle );
+      }
+
+      return output_table;
+    }
+
+  private:
+    // statistics monitoring
+    mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "number of sensors crossed"};
+    mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "number of hit candidates per crossed sensor"};
+  };
+
+} // namespace LHCb::Pr::Velo
+
+DECLARE_COMPONENT_WITH_ID( LHCb::Pr::Velo::VeloHeavyFlavourTracking, "VeloHeavyFlavourTracking" )
-- 
GitLab


From db358c253d48de604a367e6f9896668f9d259ec8 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Thu, 28 Apr 2022 16:46:35 +0200
Subject: [PATCH 02/14] searches whole module, not just sensor

---
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 104 ++++++++++++--------
 1 file changed, 63 insertions(+), 41 deletions(-)

diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index eb5f19e3886..81751211bed 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -48,8 +48,9 @@ namespace LHCb::Pr::Velo {
                                                 TrackInfo::mCorr, TrackInfo::cosAngle>;
 
     // miscellaneous
-    using HitPos   = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>;
-    using SensorID = LHCb::Detector::VPChannelID::SensorID;
+    using HitPos    = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>;
+    using SensorID  = LHCb::Detector::VPChannelID::SensorID;
+    using ChannelID = LHCb::Detector::VPChannelID;
 
     // to store relevant sensor information
     struct SensorInfo {
@@ -65,31 +66,68 @@ namespace LHCb::Pr::Velo {
     struct HitInfo {
       float       d;
       HitPos      position;
+      ChannelID   id;
       friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.d < rhs.d; }
     };
     using HitInfos = std::vector<HitInfo>;
 
     // functions
     template <typename Position3, typename Vector3>
-    HitInfos getHitsForFirstCrossedSensor( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV,
-                                           Vector3 const& tvec, float max_distance = 1. * Gaudi::Units::mm ) {
+    SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& endVertex,
+                                   Position3 const& PV, Vector3 const& tvec ) {
+      SensorInfos crossedSensors;
+      auto        pos = PV;
+      for ( auto sensor : sensors ) {
+        if ( sensor.z < pos.Z() ) continue;
+        if ( sensor.z > endVertex.Z() ) break;
+        // if in range look if search window intersects with sensor
+        pos = pos + ( sensor.z - pos.Z() ) * tvec;
+        if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) )
+          crossedSensors.push_back( {pos.Z().cast(), sensor.id} );
+      }
+      return crossedSensors;
+    }
+
+    template <typename Position3, typename Vector3>
+    HitInfos getHitsForFirstCrossedModule( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV,
+                                           Vector3 const& tvec, float max_distance ) {
       HitInfos poshits;
       if ( !sensors.size() ) return poshits;
-      // get first sensor (in z) and its reference point
+      // get first sensor (in z)
       auto firstsensor = sensors.front();
-      auto ref         = PV + ( firstsensor.z - PV.Z() ) * tvec;
-      // look if this sensor has hits and if they are in search window
+      // its associated module (to not be affected by borders between sensors)
+      auto dummy_id = ChannelID();
+      dummy_id.setSensor( firstsensor.id );
+      auto crossed_module = dummy_id.module();
+      // first sensor reference point
+      auto ref = PV + ( firstsensor.z - PV.Z() ) * tvec;
+      // look if this sensor (and associated station) has hits and if they are in search window
       for ( auto hit : hits.scalar() ) {
-        auto id = LHCb::Detector::VPChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
-        if ( id.sensor() == firstsensor.id ) {
+        auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
+        if ( id.module() == crossed_module ) {
           auto pos  = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
           auto diff = ( pos.vec3() - ref ).mag();
-          if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3()} ); }
+          if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3(), id} ); }
         }
       }
       return poshits;
     }
 
+    template <typename Composite_t, typename Vector3_t>
+    auto getMissingMomentumVars( Composite_t const& composite, Vector3_t const& direction ) {
+      using Sel::Utils::mass2;
+      using Sel::Utils::threeMomentum;
+      using std::sqrt;
+      auto const mom = threeMomentum( composite );
+      // corrected mass
+      auto const perp  = mom - direction * dot( mom, direction );
+      auto const pt    = perp.mag();
+      auto const mcorr = pt + sqrt( mass2( composite ) + pt * pt );
+      // operning angle
+      auto const cosangle = dot( mom, direction ) / mom.mag();
+      return std::make_tuple( mcorr, cosangle );
+    }
+
   } // namespace
 
   /**
@@ -134,57 +172,41 @@ namespace LHCb::Pr::Velo {
         auto const  tvec      = diff / diff.Z();
 
         // search window crosses sensors?
-        int         nInSensor = 0;
-        SensorInfos crossedSensors;
-        auto        pos = PV;
-        for ( auto sinfo : info_sensors ) {
-          if ( sinfo.z < pos.Z() ) continue;
-          if ( sinfo.z > endVertex.Z() ) break;
-          // if in range look if search window intersects with sensor
-          pos = pos + ( sinfo.z - pos.Z() ) * tvec;
-          if ( velo.sensor( sinfo.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) {
-            nInSensor++;
-            crossedSensors.push_back( {pos.Z().cast(), sinfo.id} );
-          }
-        }
+        auto crossedSensors = getCrossedSensors( velo, info_sensors, endVertex, PV, tvec );
+        int  nInSensor      = crossedSensors.size();
 
         // look for hits in relevant sensors
-        // TODO maybe check all sensors in this plane, in case the search window is at the edge?
-        auto hits_first_sensor = getHitsForFirstCrossedSensor( crossedSensors, hits, PV, tvec );
-        std::sort( hits_first_sensor.begin(), hits_first_sensor.end() );
-        bool hasHitInFirstSensor = hits_first_sensor.size() > 0;
+        auto first_hits = getHitsForFirstCrossedModule( crossedSensors, hits, PV, tvec, m_max_distance );
+        std::sort( first_hits.begin(), first_hits.end() );
+        bool hasHitInFirstModule = first_hits.size() > 0;
 
         // determine mother particle direction based on closest hit in first crossed sensor
         // or using SV-PV if not available
-        auto direction =
-            hasHitInFirstSensor ? ( hits_first_sensor.front().position - PV ).normalize() : diff.normalize();
+        auto direction = hasHitInFirstModule ? ( first_hits.front().position - PV ).normalize() : diff.normalize();
 
         // obtain corrected mass and opening angle
-        using Sel::Utils::mass2;
-        using Sel::Utils::threeMomentum;
-        using std::sqrt;
-        auto const mom      = threeMomentum( *composite );
-        auto const perp     = mom - direction * dot( mom, direction );
-        auto const pt       = perp.mag();
-        auto const mcorr    = pt + sqrt( mass2( *composite ) + pt * pt );
-        auto const cosangle = dot( mom, direction ) / mom.mag();
+        auto [mcorr, cosangle] = getMissingMomentumVars( *composite, direction );
 
         // statistics
         m_nInSensor += nInSensor;
-        if ( nInSensor > 0 ) m_nHitCandidates += hits_first_sensor.size();
+        if ( nInSensor > 0 ) m_nHitCandidates += first_hits.size();
 
         // set output
         output_table.emplace_back<SIMDWrapper::InstructionSet::Scalar>().set( composite->key(), nInSensor,
-                                                                              hasHitInFirstSensor, mcorr, cosangle );
+                                                                              hasHitInFirstModule, mcorr, cosangle );
       }
 
       return output_table;
     }
 
   private:
+    // properties
+    Gaudi::Property<float> m_max_distance{this, "MaxDistanceFromPVSV", 0.5 * Gaudi::Units::mm,
+                                          "maximum hit distance from PV-SV search window"};
+
     // statistics monitoring
-    mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "number of sensors crossed"};
-    mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "number of hit candidates per crossed sensor"};
+    mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "Nb of sensors crossed"};
+    mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "Nb hit candidates per sensor"};
   };
 
 } // namespace LHCb::Pr::Velo
-- 
GitLab


From d840280b0848d54960600ea5f8c83719d399c38a Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Tue, 25 Oct 2022 14:49:52 +0200
Subject: [PATCH 03/14] large update; tracking (info) more general, saved in
 soa

---
 .../PrKernel/PrVeloHeavyFlavourUtils.h        | 104 +++++++
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   | 294 ++++++++++--------
 2 files changed, 271 insertions(+), 127 deletions(-)
 create mode 100644 Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h

diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
new file mode 100644
index 00000000000..5be01d6168b
--- /dev/null
+++ b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
@@ -0,0 +1,104 @@
+/*****************************************************************************\
+* (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.                                       *
+\*****************************************************************************/
+
+#include "Event/Particle.h"
+#include "Event/Particle_v2.h"
+#include "Event/PrTracksTag.h"
+#include "Event/PrimaryVertices.h"
+#include "Event/RelationTable.h"
+#include "VPDet/DeVP.h"
+
+namespace LHCb::Pr::Velo {
+
+  namespace HeavyFlavourTracking {
+
+    // main heavy flavour track storage
+    namespace Tag {
+      // hit info
+      struct LHCbID : Event::lhcbid_field {};
+      struct HitPosition : Event::vec3_field {};
+      struct Hits : Event::vector_field<Event::struct_field<LHCbID, HitPosition>> {};
+
+      // crossed sensor info
+      struct nSensors : Event::int_field {};
+
+      // geometry
+      struct OriginState : Event::state_field {};
+
+      // track quality info
+      struct Chi2 : Event::float_field {};
+
+      // main storage: relating track to associated composite and PV in full SoA event model
+      template <typename Composites, typename PVs>
+      using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>;
+
+      // (hopefully temporary) fix with SOACollection with keys/indices to non SoA objects
+      struct CompositeKey : Event::int_field {};
+      struct PrimaryVertexKey : Event::int_field {};
+
+      template <typename T>
+      using Track_t = Event::SOACollection<T, CompositeKey, PrimaryVertexKey, Hits, nSensors, OriginState, Chi2>;
+    } // namespace Tag
+
+    namespace Proxy {
+
+      // base class for proxy of for a track 'object' in the track container
+      template <typename ProxyBase>
+      struct TrackProxy : ProxyBase {
+        using ProxyBase::Proxy;
+        // hit info
+        [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); }
+        [[nodiscard]] auto hitID( std::size_t i ) const {
+          return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>();
+        }
+        [[nodiscard]] auto hitPosition( std::size_t i ) const {
+          return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>();
+        }
+        // crossed sensor info
+        [[nodiscard]] auto nSensors() const { return this->template field<Tag::nSensors>(); }
+        // geometry
+        [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); }
+        // track quality info
+        [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); }
+      };
+
+    } // namespace Proxy
+
+    namespace v3 {
+      using Composites = LHCb::Event::Composites;
+      using PVs        = LHCb::Event::PV::PrimaryVertexContainer; // FIXME
+
+      template <typename Composites, typename PVs>
+      struct HeavyFlavourTracks : Tag::RelationTable2D_t<Composites, PVs> {
+        using Tag::RelationTable2D_t<Composites, PVs>::RelationTable2D;
+
+        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType>
+        using proxy_type = Proxy::TrackProxy<
+            typename Tag::RelationTable2D_t<Composites, PVs>::proxy_type<simd, behaviour, ContainerType>>;
+      };
+    } // namespace v3
+
+    namespace v1 {
+      using Composites = LHCb::Particles;
+      using PVs        = LHCb::Event::PV::PrimaryVertexContainer;
+
+      struct HeavyFlavourTracks : Tag::Track_t<HeavyFlavourTracks> {
+        using base_t = typename Tag::Track_t<HeavyFlavourTracks>;
+        using base_t::base_t;
+
+        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType>
+        using proxy_type = Proxy::TrackProxy<typename Event::Proxy<simd, behaviour, ContainerType>>;
+      };
+    } // namespace v1
+
+  } // namespace HeavyFlavourTracking
+
+} // namespace LHCb::Pr::Velo
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index 81751211bed..b66c44d3dab 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -11,63 +11,78 @@
 
 #include "DetDesc/GenericConditionAccessorHolder.h"
 #include "DetDesc/IConditionDerivationMgr.h"
-#include "Event/Particle.h"
 #include "Event/PrHits.h"
-#include "Event/PrimaryVertices.h"
-#include "Event/RawEvent.h"
-#include "Event/RelationTable.h"
 #include "GaudiAlg/Transformer.h"
 #include "Kernel/LHCbID.h"
 #include "LHCbMath/MatVec.h"
+#include "PrKernel/PrVeloHeavyFlavourUtils.h"
 #include "SelKernel/Utilities.h"
 #include "SelKernel/VertexRelation.h"
-#include "VPDet/DeVP.h"
 #include "VPDet/VPDetPaths.h"
-#include <vector>
+#include "VeloKalmanHelpers.h"
 
 namespace LHCb::Pr::Velo {
 
   namespace {
 
-    namespace VeloInfo {
-      int constexpr NSensors = 208;
-    } // namespace VeloInfo
-
-    namespace TrackInfo {
-      struct nInSensor : LHCb::Event::int_field {};
-      struct hasHitInFirstSensor : LHCb::Event::int_field {};
-      struct mCorr : LHCb::Event::float_field {};
-      struct cosAngle : LHCb::Event::float_field {};
-    } // namespace TrackInfo
-
-    // in/outputs
-    using Composites = LHCb::Particles;
-    using VeloHits   = LHCb::Pr::VP::Hits;
-    using PVs        = LHCb::Event::PV::PrimaryVertexContainer;
-    using Output     = LHCb::Event::RelationTable1D<Composites, TrackInfo::nInSensor, TrackInfo::hasHitInFirstSensor,
-                                                TrackInfo::mCorr, TrackInfo::cosAngle>;
+    // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA)
+    using namespace LHCb::Pr::Velo::HeavyFlavourTracking::v1;
+    using VeloHits = LHCb::Pr::VP::Hits;
+    using Output   = HeavyFlavourTracks;
 
     // miscellaneous
-    using HitPos    = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>;
-    using SensorID  = LHCb::Detector::VPChannelID::SensorID;
-    using ChannelID = LHCb::Detector::VPChannelID;
+    namespace TrackTag = LHCb::Pr::Velo::HeavyFlavourTracking::Tag;
+    using FType        = SIMDWrapper::scalar::float_v;
+    using Position     = LHCb::LinAlg::Vec3<FType>;
+    using SensorID     = LHCb::Detector::VPChannelID::SensorID;
+    using ChannelID    = LHCb::Detector::VPChannelID;
 
     // to store relevant sensor information
     struct SensorInfo {
-      float       z;
-      SensorID    id;
+      FType    z;
+      SensorID id;
+      // sorting
       friend bool operator<( SensorInfo const& lhs, SensorInfo const& rhs ) {
         return std::pair{lhs.z, lhs.id} < std::pair{rhs.z, rhs.id};
       }
+      // module information
+      auto module() const {
+        auto cid = ChannelID();
+        cid.setSensor( this->id );
+        return cid.module();
+      }
     };
     using SensorInfos = std::vector<SensorInfo>;
 
+    class Sensors {
+    private:
+      SensorInfos m_sensors;
+
+    public:
+      // constructors
+      Sensors() = default;
+      // conditions related constructor
+      Sensors( DeVP const& velo ) {
+        m_sensors = SensorInfos{};
+        m_sensors.reserve( velo.numberSensors() );
+        auto getSensorZs = [&]( auto const& sensor ) {
+          m_sensors.push_back( {(float)sensor.z(), sensor.sensorNumber()} );
+        };
+        velo.runOnAllSensors( getSensorZs );
+        std::sort( m_sensors.begin(), m_sensors.end() );
+      };
+      // main access
+      SensorInfos const& operator()() const { return m_sensors; }
+    };
+
     // to store relevant hit information
     struct HitInfo {
-      float       d;
-      HitPos      position;
+      FType       d;
+      Position    position;
       ChannelID   id;
-      friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.d < rhs.d; }
+      friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) {
+        return std::pair{lhs.position.z.cast(), lhs.d.cast()} < std::pair{rhs.position.z.cast(), rhs.d.cast()};
+      }
     };
     using HitInfos = std::vector<HitInfo>;
 
@@ -89,126 +104,151 @@ namespace LHCb::Pr::Velo {
     }
 
     template <typename Position3, typename Vector3>
-    HitInfos getHitsForFirstCrossedModule( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV,
-                                           Vector3 const& tvec, float max_distance ) {
+    HitInfos getHitsForCrossedModules( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start,
+                                       Vector3 const& tvec, FType max_distance ) {
       HitInfos poshits;
       if ( !sensors.size() ) return poshits;
-      // get first sensor (in z)
-      auto firstsensor = sensors.front();
-      // its associated module (to not be affected by borders between sensors)
-      auto dummy_id = ChannelID();
-      dummy_id.setSensor( firstsensor.id );
-      auto crossed_module = dummy_id.module();
-      // first sensor reference point
-      auto ref = PV + ( firstsensor.z - PV.Z() ) * tvec;
-      // look if this sensor (and associated station) has hits and if they are in search window
-      for ( auto hit : hits.scalar() ) {
-        auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
-        if ( id.module() == crossed_module ) {
-          auto pos  = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
-          auto diff = ( pos.vec3() - ref ).mag();
-          if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3(), id} ); }
+      poshits.reserve( sensors.size() );
+      // starting point search
+      auto pos = start;
+      // check sensors
+      for ( auto const sensor : sensors ) {
+        // its associated module (to not be affected by borders between sensors)
+        auto crossed_module = sensor.module();
+        // go to sensor location
+        pos = pos + ( sensor.z - pos.z.cast() ) * tvec;
+        // look if this sensor (and associated station) has hits and if they are in search window
+        for ( auto hit : hits.scalar() ) {
+          auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
+          if ( id.module() == crossed_module ) {
+            auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
+            auto diff   = ( hitpos.vec3() - pos ).mag();
+            if ( diff < max_distance ) poshits.push_back( {diff.cast(), hitpos.vec3(), id} );
+          }
         }
       }
       return poshits;
     }
 
-    template <typename Composite_t, typename Vector3_t>
-    auto getMissingMomentumVars( Composite_t const& composite, Vector3_t const& direction ) {
-      using Sel::Utils::mass2;
-      using Sel::Utils::threeMomentum;
-      using std::sqrt;
-      auto const mom = threeMomentum( composite );
-      // corrected mass
-      auto const perp  = mom - direction * dot( mom, direction );
-      auto const pt    = perp.mag();
-      auto const mcorr = pt + sqrt( mass2( composite ) + pt * pt );
-      // operning angle
-      auto const cosangle = dot( mom, direction ) / mom.mag();
-      return std::make_tuple( mcorr, cosangle );
-    }
-
   } // namespace
 
   /**
-   * @class VeloHeavyFlavourTracking
+   * @class HeavyFlavourTrackFinder
    * Search for hits of charged b/c hadrons based on primary and secondary vertices
    *
    * @author Maarten van Veghel (Nikhef, RUG)
    */
 
-  class VeloHeavyFlavourTracking
-      : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const& ),
-                                              LHCb::DetDesc::usesConditions<DeVP>> {
+  class HeavyFlavourTrackFinder
+      : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const&,
+                                                      Sensors const& ),
+                                              LHCb::DetDesc::usesConditions<DeVP, Sensors>> {
   public:
     // standard constructor
-    VeloHeavyFlavourTracking( const std::string& name, ISvcLocator* pSvcLocator )
-        : Transformer( name, pSvcLocator,
-                       {KeyValue{"Composites", ""}, KeyValue{"PVs", ""}, KeyValue{"Hits", ""},
-                        KeyValue{"DeVP", LHCb::Det::VP::det_path}},
-                       {KeyValue{"OutputLocation", ""}} ) {}
-
-    // main execution
-    Output operator()( Composites const& composites, PVs const& pvs, VeloHits const& hits,
-                       DeVP const& velo ) const override {
-      // get z positions of sensors (should be in initialize TODO)
-      SensorInfos info_sensors( VeloInfo::NSensors );
-      auto        getSensorZs = [&info_sensors]( auto const& sensor ) {
-        info_sensors.push_back( {(float)sensor.z(), sensor.sensorNumber()} );
-      };
-      velo.runOnAllSensors( getSensorZs );
-      std::sort( info_sensors.begin(), info_sensors.end() );
-
-      // output
-      Output output_table( &composites );
-      output_table.reserve( composites.size() );
-
-      for ( auto const composite : composites ) {
-        // get SV and PV information
-        using Sel::Utils::endVertexPos;
-        auto const& endVertex = endVertexPos( *composite );
-        auto const& PV        = endVertexPos( Sel::getBestPV( *composite, pvs ) );
-        auto const  diff      = endVertex - PV;
-        auto const  tvec      = diff / diff.Z();
-
-        // search window crosses sensors?
-        auto crossedSensors = getCrossedSensors( velo, info_sensors, endVertex, PV, tvec );
-        int  nInSensor      = crossedSensors.size();
-
-        // look for hits in relevant sensors
-        auto first_hits = getHitsForFirstCrossedModule( crossedSensors, hits, PV, tvec, m_max_distance );
-        std::sort( first_hits.begin(), first_hits.end() );
-        bool hasHitInFirstModule = first_hits.size() > 0;
-
-        // determine mother particle direction based on closest hit in first crossed sensor
-        // or using SV-PV if not available
-        auto direction = hasHitInFirstModule ? ( first_hits.front().position - PV ).normalize() : diff.normalize();
-
-        // obtain corrected mass and opening angle
-        auto [mcorr, cosangle] = getMissingMomentumVars( *composite, direction );
-
-        // statistics
-        m_nInSensor += nInSensor;
-        if ( nInSensor > 0 ) m_nHitCandidates += first_hits.size();
-
-        // set output
-        output_table.emplace_back<SIMDWrapper::InstructionSet::Scalar>().set( composite->key(), nInSensor,
-                                                                              hasHitInFirstModule, mcorr, cosangle );
-      }
+    HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator );
 
-      return output_table;
-    }
+    // initialize
+    StatusCode initialize() override;
+
+    // main function/operator
+    Output operator()( Composites const&, PVs const&, VeloHits const&, DeVP const&, Sensors const& ) const override;
 
   private:
     // properties
     Gaudi::Property<float> m_max_distance{this, "MaxDistanceFromPVSV", 0.5 * Gaudi::Units::mm,
-                                          "maximum hit distance from PV-SV search window"};
+                                          "maximum hit distance from PV-SV line (search window)"};
 
     // statistics monitoring
     mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "Nb of sensors crossed"};
-    mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "Nb hit candidates per sensor"};
-  };
+    mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "Nb hit candidates"};
+  }; // namespace LHCb::Pr::Velo
+
+  DECLARE_COMPONENT_WITH_ID( HeavyFlavourTrackFinder, "VeloHeavyFlavourTrackFinder" )
+
+  // ============================= IMPLEMENTATION ===================================
+
+  // constructor
+  HeavyFlavourTrackFinder::HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator )
+      : Transformer( name, pSvcLocator,
+                     {KeyValue( "Composites", "" ), KeyValue( "PVs", "" ), KeyValue( "Hits", "" ),
+                      KeyValue( "DeVP", LHCb::Det::VP::det_path ),
+#ifdef USE_DD4HEP
+                      KeyValue( "Sensors", {"/world:AlgorithmSpecific-" + name + "-sensorinfos"} )},
+#else
+                      KeyValue( "Sensors", {"AlgorithmSpecific-" + name + "-sensorinfos"} )},
+#endif
+
+                     {KeyValue{"OutputLocation", ""}} ) {
+  }
+
+  // initialization
+  StatusCode LHCb::Pr::Velo::HeavyFlavourTrackFinder::initialize() {
+    auto sc = Transformer::initialize().andThen( [&] {
+      addConditionDerivation<Sensors( DeVP const& )>( {LHCb::Det::VP::det_path}, inputLocation<Sensors>() );
+    } );
+    return sc;
+  }
+
+  // main execution
+  Output LHCb::Pr::Velo::HeavyFlavourTrackFinder::operator()( Composites const& composites, PVs const& pvs,
+                                                              VeloHits const& velohits, DeVP const& velo,
+                                                              Sensors const& sensors ) const {
+    // output
+    Output output{}; // composites.zipIdentifier()};
+    output.reserve( composites.size() );
+
+    for ( auto const composite : composites ) {
+      // set output
+      auto track = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>();
+      track.field<TrackTag::CompositeKey>().set( composite->key() );
+
+      // primary vertex
+      using Sel::Utils::endVertexPos;
+      auto const& PV    = Sel::getBestPV( *composite, pvs );
+      auto const& PVpos = endVertexPos( PV );
+      track.field<TrackTag::CompositeKey>().set( PV.key() );
+
+      // secondary vertex
+      auto const& endVertex = endVertexPos( *composite );
+      auto const  diff      = endVertex - PVpos;
+      auto const  tvec      = diff / diff.Z();
+
+      // crossed sensors in PV-SV line
+      auto crossedSensors = getCrossedSensors( velo, sensors(), endVertex, PVpos, tvec );
+      track.field<TrackTag::nSensors>().set( crossedSensors.size() );
+
+      // look for hits in relevant sensors and sort in z and distance to search window
+      auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, PVpos, tvec, m_max_distance.value() );
+      std::sort( found_hits.begin(), found_hits.end() );
+
+      // store sorted hits
+      auto saved_hits = track.field<TrackTag::Hits>();
+      auto nhits      = found_hits.size();
+      saved_hits.resize( nhits );
+      for ( long unsigned int i = 0; i < nhits; i++ ) {
+        auto found_hit = found_hits[i];
+        saved_hits[i].field<TrackTag::LHCbID>().set( LHCbID( found_hit.id ) );
+        saved_hits[i].field<TrackTag::HitPosition>().set( found_hit.position );
+      }
 
-} // namespace LHCb::Pr::Velo
+      // determine mother particle direction based on closest hit
+      // in first crossed sensor or using SV-PV if not available
+      auto direction    = found_hits.size() ? found_hits.front().position - PVpos : diff;
+      auto origin_state = track.field<TrackTag::OriginState>();
+      origin_state.setPosition( PVpos );
+      origin_state.setDirection( direction.x / direction.z, direction.y / direction.z );
+      origin_state.setQOverP( 1. / composite->p() );
+
+      // fitting
+      // TODO use stuff from VeloKalmanHelpers
+      track.field<TrackTag::Chi2>().set( -1.f );
+
+      // statistics
+      m_nInSensor += crossedSensors.size();
+      m_nHitCandidates += nhits;
+    }
 
-DECLARE_COMPONENT_WITH_ID( LHCb::Pr::Velo::VeloHeavyFlavourTracking, "VeloHeavyFlavourTracking" )
+    return output;
+  }
+
+} // namespace LHCb::Pr::Velo
-- 
GitLab


From 562404700ab34cda9329ded7af9f080f620bc7af Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Tue, 25 Oct 2022 22:37:24 +0200
Subject: [PATCH 04/14] hit duplicate removal

---
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 63 ++++++++++++---------
 1 file changed, 36 insertions(+), 27 deletions(-)

diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index b66c44d3dab..b96369a569c 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -77,26 +77,29 @@ namespace LHCb::Pr::Velo {
 
     // to store relevant hit information
     struct HitInfo {
-      FType       d;
-      Position    position;
-      ChannelID   id;
+      FType     d;
+      Position  position;
+      ChannelID id;
+      // sorting in z first, then distance to search window
       friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) {
         return std::pair{lhs.position.z.cast(), lhs.d.cast()} < std::pair{rhs.position.z.cast(), rhs.d.cast()};
       }
+      // equal just by unique channel id
+      friend bool operator==( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.id == rhs.id; }
     };
     using HitInfos = std::vector<HitInfo>;
 
     // functions
     template <typename Position3, typename Vector3>
-    SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& endVertex,
-                                   Position3 const& PV, Vector3 const& tvec ) {
+    SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& start,
+                                   Position3 const& end, Vector3 const& slopes ) {
       SensorInfos crossedSensors;
-      auto        pos = PV;
+      auto        pos = start;
       for ( auto sensor : sensors ) {
         if ( sensor.z < pos.Z() ) continue;
-        if ( sensor.z > endVertex.Z() ) break;
+        if ( sensor.z > end.Z() ) break;
         // if in range look if search window intersects with sensor
-        pos = pos + ( sensor.z - pos.Z() ) * tvec;
+        pos = pos + ( sensor.z - pos.Z() ) * slopes;
         if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) )
           crossedSensors.push_back( {pos.Z().cast(), sensor.id} );
       }
@@ -105,10 +108,10 @@ namespace LHCb::Pr::Velo {
 
     template <typename Position3, typename Vector3>
     HitInfos getHitsForCrossedModules( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start,
-                                       Vector3 const& tvec, FType max_distance ) {
+                                       Vector3 const& slopes, FType max_distance ) {
       HitInfos poshits;
       if ( !sensors.size() ) return poshits;
-      poshits.reserve( sensors.size() );
+      poshits.reserve( 2 * sensors.size() );
       // starting point search
       auto pos = start;
       // check sensors
@@ -116,17 +119,24 @@ namespace LHCb::Pr::Velo {
         // its associated module (to not be affected by borders between sensors)
         auto crossed_module = sensor.module();
         // go to sensor location
-        pos = pos + ( sensor.z - pos.z.cast() ) * tvec;
-        // look if this sensor (and associated station) has hits and if they are in search window
+        pos = pos + ( sensor.z - pos.z.cast() ) * slopes;
+        // look if this sensor (and associated module) has hits and if they are in search window
         for ( auto hit : hits.scalar() ) {
           auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
           if ( id.module() == crossed_module ) {
             auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
-            auto diff   = ( hitpos.vec3() - pos ).mag();
-            if ( diff < max_distance ) poshits.push_back( {diff.cast(), hitpos.vec3(), id} );
+            // make sure they are in the same plane up to some tolerance
+            // (not all sensors in a module are exactly)
+            if ( abs( hitpos.vec3().z.cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) {
+              auto d = ( hitpos.vec3() - pos ).mag();
+              if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos.vec3(), id} );
+            }
           }
         }
       }
+      // sort and remove duplicates (can arise with only module/plane (not sensor) checking)
+      std::sort( poshits.begin(), poshits.end() );
+      poshits.erase( std::unique( poshits.begin(), poshits.end() ), poshits.end() );
       return poshits;
     }
 
@@ -194,7 +204,7 @@ namespace LHCb::Pr::Velo {
                                                               VeloHits const& velohits, DeVP const& velo,
                                                               Sensors const& sensors ) const {
     // output
-    Output output{}; // composites.zipIdentifier()};
+    Output output{};
     output.reserve( composites.size() );
 
     for ( auto const composite : composites ) {
@@ -204,22 +214,21 @@ namespace LHCb::Pr::Velo {
 
       // primary vertex
       using Sel::Utils::endVertexPos;
-      auto const& PV    = Sel::getBestPV( *composite, pvs );
-      auto const& PVpos = endVertexPos( PV );
+      auto const& PV = Sel::getBestPV( *composite, pvs );
+      auto const& pv = endVertexPos( PV );
       track.field<TrackTag::CompositeKey>().set( PV.key() );
 
       // secondary vertex
-      auto const& endVertex = endVertexPos( *composite );
-      auto const  diff      = endVertex - PVpos;
-      auto const  tvec      = diff / diff.Z();
+      auto const& end    = endVertexPos( *composite );
+      auto const  fd     = end - pv;
+      auto const  slopes = fd / fd.Z();
 
       // crossed sensors in PV-SV line
-      auto crossedSensors = getCrossedSensors( velo, sensors(), endVertex, PVpos, tvec );
+      auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes );
       track.field<TrackTag::nSensors>().set( crossedSensors.size() );
 
-      // look for hits in relevant sensors and sort in z and distance to search window
-      auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, PVpos, tvec, m_max_distance.value() );
-      std::sort( found_hits.begin(), found_hits.end() );
+      // look for hits in relevant sensors (duplicated remoded and sorted in z and distance to search window)
+      auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() );
 
       // store sorted hits
       auto saved_hits = track.field<TrackTag::Hits>();
@@ -233,14 +242,14 @@ namespace LHCb::Pr::Velo {
 
       // determine mother particle direction based on closest hit
       // in first crossed sensor or using SV-PV if not available
-      auto direction    = found_hits.size() ? found_hits.front().position - PVpos : diff;
+      auto direction    = found_hits.size() ? found_hits.front().position - pv : fd;
       auto origin_state = track.field<TrackTag::OriginState>();
-      origin_state.setPosition( PVpos );
+      origin_state.setPosition( pv );
       origin_state.setDirection( direction.x / direction.z, direction.y / direction.z );
       origin_state.setQOverP( 1. / composite->p() );
 
       // fitting
-      // TODO use stuff from VeloKalmanHelpers
+      // TODO (using VeloKalmanHelpers)
       track.field<TrackTag::Chi2>().set( -1.f );
 
       // statistics
-- 
GitLab


From e17338cbe7adde486d55a84b57cac4572dc0ddfc Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Tue, 1 Nov 2022 18:54:41 +0100
Subject: [PATCH 05/14] added mc checker of hf tracking

---
 .../PrKernel/PrVeloHeavyFlavourUtils.h        |  83 +++---
 Pr/PrMCTools/CMakeLists.txt                   |   2 +
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 249 ++++++++++++++++++
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   |   4 +-
 4 files changed, 291 insertions(+), 47 deletions(-)
 create mode 100644 Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp

diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
index 5be01d6168b..941c308233b 100644
--- a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
+++ b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
@@ -36,66 +36,59 @@ namespace LHCb::Pr::Velo {
       // track quality info
       struct Chi2 : Event::float_field {};
 
-      // main storage: relating track to associated composite and PV in full SoA event model
-      template <typename Composites, typename PVs>
-      using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>;
-
       // (hopefully temporary) fix with SOACollection with keys/indices to non SoA objects
-      struct CompositeKey : Event::int_field {};
-      struct PrimaryVertexKey : Event::int_field {};
+      struct CompositeIndex : Event::int_field {};
+      struct PrimaryVertexIndex : Event::int_field {};
 
       template <typename T>
-      using Track_t = Event::SOACollection<T, CompositeKey, PrimaryVertexKey, Hits, nSensors, OriginState, Chi2>;
-    } // namespace Tag
-
-    namespace Proxy {
-
-      // base class for proxy of for a track 'object' in the track container
-      template <typename ProxyBase>
-      struct TrackProxy : ProxyBase {
-        using ProxyBase::Proxy;
-        // hit info
-        [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); }
-        [[nodiscard]] auto hitID( std::size_t i ) const {
-          return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>();
-        }
-        [[nodiscard]] auto hitPosition( std::size_t i ) const {
-          return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>();
-        }
-        // crossed sensor info
-        [[nodiscard]] auto nSensors() const { return this->template field<Tag::nSensors>(); }
-        // geometry
-        [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); }
-        // track quality info
-        [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); }
-      };
-
-    } // namespace Proxy
-
-    namespace v3 {
-      using Composites = LHCb::Event::Composites;
-      using PVs        = LHCb::Event::PV::PrimaryVertexContainer; // FIXME
+      using Track_t = Event::SOACollection<T, CompositeIndex, PrimaryVertexIndex, Hits, nSensors, OriginState, Chi2>;
 
+      /*
+      // main storage: relating track to associated composite and PV in full SoA event model
       template <typename Composites, typename PVs>
-      struct HeavyFlavourTracks : Tag::RelationTable2D_t<Composites, PVs> {
-        using Tag::RelationTable2D_t<Composites, PVs>::RelationTable2D;
-
-        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType>
-        using proxy_type = Proxy::TrackProxy<
-            typename Tag::RelationTable2D_t<Composites, PVs>::proxy_type<simd, behaviour, ContainerType>>;
-      };
-    } // namespace v3
+      using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>;
+      */
+    } // namespace Tag
 
     namespace v1 {
       using Composites = LHCb::Particles;
+      using Composite  = LHCb::Particle;
       using PVs        = LHCb::Event::PV::PrimaryVertexContainer;
 
       struct HeavyFlavourTracks : Tag::Track_t<HeavyFlavourTracks> {
         using base_t = typename Tag::Track_t<HeavyFlavourTracks>;
         using base_t::base_t;
 
+        template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType>
+        struct TrackProxy : Event::Proxy<simd, behaviour, ContainerType> {
+          using base_t = typename Event::Proxy<simd, behaviour, ContainerType>;
+          using base_t::Proxy;
+
+          // composite / pv info
+          [[nodiscard]] auto compositeIndex() const { return this->template get<Tag::CompositeIndex>(); }
+          [[nodiscard]] auto primaryVertexIndex() const { return this->template get<Tag::PrimaryVertexIndex>(); }
+
+          // hit info
+          [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); }
+          [[nodiscard]] auto hitID( std::size_t i ) const {
+            return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>();
+          }
+          [[nodiscard]] auto hitPosition( std::size_t i ) const {
+            return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>();
+          }
+
+          // crossed sensor info
+          [[nodiscard]] auto nSensors() const { return this->template get<Tag::nSensors>(); }
+
+          // geometry
+          [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); }
+
+          // track quality info
+          [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); }
+        };
+
         template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType>
-        using proxy_type = Proxy::TrackProxy<typename Event::Proxy<simd, behaviour, ContainerType>>;
+        using proxy_type = TrackProxy<simd, behaviour, ContainerType>;
       };
     } // namespace v1
 
diff --git a/Pr/PrMCTools/CMakeLists.txt b/Pr/PrMCTools/CMakeLists.txt
index 62b61171ae9..4cae29bd6a7 100644
--- a/Pr/PrMCTools/CMakeLists.txt
+++ b/Pr/PrMCTools/CMakeLists.txt
@@ -38,6 +38,7 @@ gaudi_add_module(PrMCTools
         src/PrVPHitsChecker.cpp
         src/PrUTHitsChecker.cpp
         src/PrFTHitsChecker.cpp
+        src/PrVeloHeavyFlavourTrackingChecker.cpp
         src/VPClusterEfficiency.cpp
         src/PrGhostDumper.cpp
         src/PrDownstreamDumper.cpp
@@ -70,6 +71,7 @@ gaudi_add_module(PrMCTools
         Rec::PrFitParamsLib
         Rec::PrKernel
         Rec::TrackInterfacesLib
+        Rec::SelKernelLib
         ROOT::Core
         ROOT::RIO
         ROOT::Tree
diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
new file mode 100644
index 00000000000..503e2afedff
--- /dev/null
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -0,0 +1,249 @@
+/*****************************************************************************\
+* (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.                                       *
+\*****************************************************************************/
+#include "DetDesc/GenericConditionAccessorHolder.h"
+#include "Event/LinksByKey.h"
+#include "Event/MCHit.h"
+#include "Event/MCParticle.h"
+#include "Event/RawEvent.h"
+#include "Event/VPDigit.h"
+#include "Event/VPFullCluster.h"
+#include "GaudiAlg/Consumer.h"
+#include "GaudiAlg/GaudiTupleAlg.h"
+#include "Linker/LinkedFrom.h"
+#include "Linker/LinkedTo.h"
+#include "PrKernel/PrVeloHeavyFlavourUtils.h"
+
+/** @class PrVeloHeavyFlavourTrackingChecker PrVeloHeavyFlavourTrackingChecker.cpp
+ *
+ * For efficiency checker of VeloHeavyFlavourTrackFinder on MC
+ *
+ */
+
+namespace {
+  // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA)
+  using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks;
+  using Composites         = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composites;
+  using Composite          = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite;
+  using MCParticles        = LHCb::MCParticles;
+  using MCHits             = LHCb::MCHits;
+  using Links              = LHCb::LinksByKey;
+
+  // helper classes
+  struct DecayMCParticles {
+    LHCb::MCParticle const*                B     = NULL;
+    LHCb::MCParticle const*                tau   = NULL;
+    std::array<LHCb::MCParticle const*, 3> pions = {NULL, NULL, NULL};
+  };
+
+  // mcparticle filter
+  std::optional<DecayMCParticles> has_b2taunu_tau23pions( LHCb::MCParticle const* mcp ) {
+    if ( mcp->particleID().hasBottom() ) {
+      auto vtxs = mcp->endVertices();
+      if ( vtxs.size() ) {
+        auto sv   = vtxs.at( vtxs.size() - 1 );
+        auto daus = sv->products();
+        if ( daus.size() >= 2 ) {
+          bool   has_taup  = false;
+          bool   has_ntau  = false;
+          size_t n_photons = 0;
+          for ( LHCb::MCParticle const* dau : daus ) {
+            auto dID = dau->particleID().abspid();
+            if ( dID == 15 ) has_taup = true;
+            if ( dID == 16 ) has_ntau = true;
+            if ( dID == 22 ) n_photons += 1;
+          }
+          if ( has_taup && has_ntau && ( n_photons == ( daus.size() - 2 ) ) ) {
+            auto parts = DecayMCParticles();
+            parts.B    = mcp;
+            for ( LHCb::MCParticle const* dau : daus ) {
+              if ( dau->particleID().abspid() == 15 ) parts.tau = dau;
+            }
+            auto tau_vtxs = parts.tau->endVertices();
+            if ( tau_vtxs.size() ) {
+              auto tv     = tau_vtxs.at( tau_vtxs.size() - 1 );
+              int  i_pion = 0;
+              for ( LHCb::MCParticle const* tdau : tv->products() ) {
+                if ( tdau->particleID().abspid() == 211 ) {
+                  parts.pions[i_pion] = tdau;
+                  i_pion += 1;
+                }
+              }
+            }
+            return parts;
+          } // getting all relevant MCParticles
+        }
+      }
+    }
+    return std::nullopt;
+  }
+
+  // match MCParticle to Composites
+  template <typename LinkMap>
+  std::optional<Composite const*> matchToComposites( DecayMCParticles const& mcps, Composites const& composites,
+                                                     LinkMap const& map, MCParticles const& mcparts ) {
+    for ( Composite const* comp : composites ) {
+      size_t nmatches = 0;
+      for ( Composite const* dau : comp->daughtersVector() ) {
+        if ( !dau->proto() ) continue;
+        for ( LHCb::MCParticle const* pion : mcps.pions ) {
+          map.applyToLinks( dau->proto()->track()->key(), [&nmatches, &mcparts,
+                                                           &pion]( unsigned int, unsigned int mcPartKey, float ) {
+            if ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) nmatches += 1;
+          } );
+        }
+      }
+      // assuming there is only one matching composite
+      if ( nmatches == mcps.pions.size() ) return comp;
+    }
+    return std::nullopt;
+  }
+
+  // match MCParticles to MCHits
+  auto getMCHits( DecayMCParticles const& mcps, MCHits const& mchits ) {
+    std::vector<LHCb::MCHit const*> hits;
+    // which ones come from B or tau?
+    for ( LHCb::MCHit const* mchit : mchits ) {
+      auto hitmcp = mchit->mcParticle();
+      if ( hitmcp == mcps.B || hitmcp == mcps.tau ) hits.push_back( mchit );
+    }
+    // order them in z
+    auto ordering = []( LHCb::MCHit const* i1, LHCb::MCHit const* i2 ) { return i1->entry().z() < i2->entry().z(); };
+    std::sort( hits.begin(), hits.end(), ordering );
+    return hits;
+  }
+
+  // get track
+  std::optional<int> getHeavyFlavourTrack( HeavyFlavourTracks const& tracks, int key ) {
+    for ( auto track : tracks.scalar() ) {
+      if ( track.compositeIndex() == key ) return track.indices().cast();
+    }
+    return std::nullopt;
+  }
+
+} // namespace
+
+class PrVeloHeavyFlavourTrackingChecker
+    : public Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&,
+                                               MCHits const&, Links const&, Links const& ),
+                                         LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>> {
+public:
+  using base_type = Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&,
+                                                      MCHits const&, Links const&, Links const& ),
+                                                LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>>;
+  using KeyValue  = base_type::KeyValue;
+
+  // standard constructor
+  PrVeloHeavyFlavourTrackingChecker( std::string const& name, ISvcLocator* pSvc )
+      : base_type{name,
+                  pSvc,
+                  {KeyValue{"MCParticles", LHCb::MCParticleLocation::Default}, KeyValue{"Composites", ""},
+                   KeyValue{"HeavyFlavourTracks", ""}, KeyValue{"MCHits", "/Event/MC/VP/Hits"}, KeyValue{"VPLinks", ""},
+                   KeyValue{"TrackLinks", ""}}} {}
+
+  // main execution
+  void operator()( MCParticles const&, Composites const&, HeavyFlavourTracks const&, MCHits const&, Links const&,
+                   Links const& ) const override;
+
+private:
+  // monitoring
+  mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""};
+};
+
+DECLARE_COMPONENT_WITH_ID( PrVeloHeavyFlavourTrackingChecker, "PrVeloHeavyFlavourTrackingChecker" )
+
+//// implementation /////
+
+// main execution
+void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, Composites const& composites,
+                                                    HeavyFlavourTracks const& tracks, MCHits const& mchits,
+                                                    Links const& vplinks, Links const& tracklinks ) const {
+  // setup ntuple
+  auto tuple = this->nTuple( "MCParticleTuple" );
+
+  // table linking a LHCb::MCHit* to std::vector<LHCb::Detector::VPChannelID>> and vice versa
+  std::map<LHCb::MCHit const*, std::vector<unsigned int>> channelIDForMCHit;
+  vplinks.applyToAllLinks( [&channelIDForMCHit, &mchits]( unsigned int channelID, unsigned int mcHitKey, float ) {
+    channelIDForMCHit[mchits[mcHitKey]].emplace_back( channelID );
+  } );
+
+  // fill for all tracks
+  for ( auto const mcp : mcparts ) {
+    if ( !mcp ) continue;
+    // check if it is the right MCParticle
+    auto check_mcps = has_b2taunu_tau23pions( mcp );
+    if ( !check_mcps.has_value() ) continue;
+    auto mcps = check_mcps.value();
+    // get MCHits for the candidate
+    auto bhits       = getMCHits( mcps, mchits );
+    auto first_mchit = bhits.size() ? bhits.front()->midPoint() : Gaudi::XYZPoint( 0., 0., 0. );
+    // first link MCParticle to daughters (pions) of composite
+    auto composite = matchToComposites( mcps, composites, tracklinks, mcparts );
+    auto has_comp  = composite.has_value();
+    // get heavy flavour track
+    auto itrack    = has_comp ? getHeavyFlavourTrack( tracks, composite.value()->key() ) : std::nullopt;
+    auto has_track = has_comp ? itrack.has_value() : false;
+    auto has_hits  = has_track ? (int)tracks.scalar()[itrack.value()].nHits() > 0 : false;
+    //// fill ntuple with info
+    std::string base = "MCP_";
+    // mc truth
+    auto sc = tuple->column( base + "TRUEID", mcp->particleID().pid() );
+    sc &= tuple->column( base + "P", mcp->momentum().P() );
+    // mc hits
+    sc &= tuple->column( base + "has_VP_MCHits", (bool)bhits.size() );
+    sc &= tuple->column( base + "VP_first_MCHit_x", first_mchit.x() );
+    sc &= tuple->column( base + "VP_first_MCHit_y", first_mchit.y() );
+    sc &= tuple->column( base + "VP_first_MCHit_z", first_mchit.z() );
+    // composite info
+    sc &= tuple->column( base + "has_Composite", has_comp );
+    std::string compbase = base + "Composite_";
+    auto        endvtx   = has_comp ? composite.value()->endVertex()->position() : Gaudi::XYZPoint( 0., 0., 0. );
+    sc &= tuple->column( compbase + "EndVertex_x", endvtx.x() );
+    sc &= tuple->column( compbase + "EndVertex_y", endvtx.y() );
+    sc &= tuple->column( compbase + "EndVertex_z", endvtx.z() );
+    auto compmom =
+        has_comp ? composite.value()->momentum() : Gaudi::LorentzVector( 0., 0., 0., -1 * Gaudi::Units::GeV );
+    sc &= tuple->column( compbase + "P", compmom.R() );
+    sc &= tuple->column( compbase + "PT", compmom.Pt() );
+    sc &= tuple->column( compbase + "eta", compmom.Eta() );
+    sc &= tuple->column( compbase + "phi", compmom.Phi() );
+    // track info
+    std::string hfbase = base + "HeavyFlavourTracking_";
+    sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks.scalar()[itrack.value()].nHits() : -1 );
+    // track mc matching info
+    bool fh_from_mcps = false;
+    int  fh_trueid    = 0;
+    if ( has_hits ) {
+      float              max_weight = 0.f;
+      LHCb::MCHit const* mchit      = nullptr;
+      unsigned int       hid        = tracks.scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID();
+      vplinks.applyToLinks( hid, [&max_weight, &mchits, &mchit]( unsigned int, unsigned int mcHitKey, float weight ) {
+        if ( weight > max_weight ) {
+          max_weight = weight;
+          mchit      = static_cast<LHCb::MCHit const*>( mchits.containedObject( mcHitKey ) );
+        }
+      } );
+      if ( mchit ) {
+        fh_trueid    = mchit->mcParticle()->particleID().pid();
+        fh_from_mcps = ( mchit->mcParticle() == mcps.B ) || ( mchit->mcParticle() == mcps.tau );
+      }
+    }
+    sc &= tuple->column( hfbase + "Track_first_hit_x",
+                         has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. );
+    sc &= tuple->column( hfbase + "Track_first_hit_y",
+                         has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. );
+    sc &= tuple->column( hfbase + "Track_first_hit_z",
+                         has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. );
+    sc &= tuple->column( hfbase + "Track_first_hit_from_MCPs", fh_from_mcps );
+    sc &= tuple->column( hfbase + "Track_first_hit_MC_TRUEID", fh_trueid );
+    // write result
+    sc.andThen( [&] { return tuple->write(); } ).ignore();
+  }
+}
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index b96369a569c..5fab6af5c4b 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -210,13 +210,13 @@ namespace LHCb::Pr::Velo {
     for ( auto const composite : composites ) {
       // set output
       auto track = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>();
-      track.field<TrackTag::CompositeKey>().set( composite->key() );
+      track.field<TrackTag::CompositeIndex>().set( composite->key() );
 
       // primary vertex
       using Sel::Utils::endVertexPos;
       auto const& PV = Sel::getBestPV( *composite, pvs );
       auto const& pv = endVertexPos( PV );
-      track.field<TrackTag::CompositeKey>().set( PV.key() );
+      track.field<TrackTag::PrimaryVertexIndex>().set( PV.key() );
 
       // secondary vertex
       auto const& end    = endVertexPos( *composite );
-- 
GitLab


From 8fee8c4381110ab43bae8fe81cc1e7adfdccfee2 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Tue, 8 Nov 2022 14:16:42 +0100
Subject: [PATCH 06/14] hf track mc checker running on all mcparticles

---
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 227 +++++++++++++++---
 1 file changed, 190 insertions(+), 37 deletions(-)

diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index 503e2afedff..a7caf773073 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -13,13 +13,17 @@
 #include "Event/MCHit.h"
 #include "Event/MCParticle.h"
 #include "Event/RawEvent.h"
+#include "Event/Track_v3.h"
 #include "Event/VPDigit.h"
 #include "Event/VPFullCluster.h"
+#include "Functors/Common.h"
 #include "GaudiAlg/Consumer.h"
 #include "GaudiAlg/GaudiTupleAlg.h"
 #include "Linker/LinkedFrom.h"
 #include "Linker/LinkedTo.h"
 #include "PrKernel/PrVeloHeavyFlavourUtils.h"
+#include "SelKernel/Utilities.h"
+#include "SelKernel/VertexRelation.h"
 
 /** @class PrVeloHeavyFlavourTrackingChecker PrVeloHeavyFlavourTrackingChecker.cpp
  *
@@ -32,15 +36,24 @@ namespace {
   using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks;
   using Composites         = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composites;
   using Composite          = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite;
+  using PVs                = LHCb::Pr::Velo::HeavyFlavourTracking::v1::PVs;
   using MCParticles        = LHCb::MCParticles;
   using MCHits             = LHCb::MCHits;
   using Links              = LHCb::LinksByKey;
 
   // helper classes
   struct DecayMCParticles {
-    LHCb::MCParticle const*                B     = NULL;
-    LHCb::MCParticle const*                tau   = NULL;
-    std::array<LHCb::MCParticle const*, 3> pions = {NULL, NULL, NULL};
+    LHCb::MCParticle const*                B        = NULL;
+    LHCb::MCParticle const*                tau      = NULL;
+    std::array<LHCb::MCParticle const*, 3> pions    = {NULL, NULL, NULL};
+    std::array<Composite const*, 3>        recpions = {NULL, NULL, NULL};
+  };
+
+  struct WeightedHit {
+    unsigned int id;
+    float        weight;
+    // reverse sorting (high to low)
+    friend bool operator<( WeightedHit& lhs, WeightedHit& rhs ) { return lhs.weight > rhs.weight; }
   };
 
   // mcparticle filter
@@ -85,18 +98,50 @@ namespace {
     return std::nullopt;
   }
 
+  // check if MCParticle has match to track
+  template <typename LinkMap>
+  bool hasTrack( LHCb::MCParticle const* mcp, LinkMap const& map, float min_weight = 0.7 ) {
+    bool has_track = false;
+    map.applyToAllLinks( [&has_track, &mcp, &min_weight]( unsigned int, unsigned int mcPartKey, float weight ) {
+      if ( ( (unsigned int)mcp->key() == mcPartKey ) && weight > min_weight ) has_track = true;
+    } );
+    return has_track;
+  }
+
+  // match MCParticle daughters (pions from tau) to Tracks (general), 'reconstructed' variable
+  template <typename LinkMap>
+  bool matchToTracks( DecayMCParticles const& mcps, LinkMap const& map, MCParticles const& mcparts,
+                      float min_weight = 0.7 ) {
+    size_t nmatches = 0;
+    for ( LHCb::MCParticle const* pion : mcps.pions ) {
+      map.applyToAllLinks(
+          [&pion, &mcps, &mcparts, &min_weight, &nmatches]( unsigned int, unsigned int mcPartKey, float weight ) {
+            if ( ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) &&
+                 ( weight > min_weight ) )
+              nmatches += 1;
+          } );
+    }
+    return nmatches == mcps.pions.size();
+  }
+
   // match MCParticle to Composites
   template <typename LinkMap>
-  std::optional<Composite const*> matchToComposites( DecayMCParticles const& mcps, Composites const& composites,
-                                                     LinkMap const& map, MCParticles const& mcparts ) {
+  std::optional<Composite const*> matchToComposites( DecayMCParticles& mcps, Composites const& composites,
+                                                     LinkMap const& map, MCParticles const& mcparts,
+                                                     float min_weight = 0.7 ) {
     for ( Composite const* comp : composites ) {
       size_t nmatches = 0;
       for ( Composite const* dau : comp->daughtersVector() ) {
         if ( !dau->proto() ) continue;
-        for ( LHCb::MCParticle const* pion : mcps.pions ) {
-          map.applyToLinks( dau->proto()->track()->key(), [&nmatches, &mcparts,
-                                                           &pion]( unsigned int, unsigned int mcPartKey, float ) {
-            if ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) nmatches += 1;
+        for ( int i = 0; i < (int)mcps.pions.size(); i++ ) {
+          LHCb::MCParticle const* pion = mcps.pions[i];
+          map.applyToLinks( dau->proto()->track()->key(), [&i, &dau, &mcps, &nmatches, &mcparts, &pion, &min_weight](
+                                                              unsigned int, unsigned int mcPartKey, float weight ) {
+            if ( ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) &&
+                 ( weight >= min_weight ) ) {
+              mcps.recpions[i] = dau;
+              nmatches += 1;
+            }
           } );
         }
       }
@@ -116,7 +161,7 @@ namespace {
     }
     // order them in z
     auto ordering = []( LHCb::MCHit const* i1, LHCb::MCHit const* i2 ) { return i1->entry().z() < i2->entry().z(); };
-    std::sort( hits.begin(), hits.end(), ordering );
+    if ( hits.size() > 1 ) std::sort( hits.begin(), hits.end(), ordering );
     return hits;
   }
 
@@ -128,31 +173,41 @@ namespace {
     return std::nullopt;
   }
 
+  // composite observables
+  float getMcorr( Composite const* composite, Gaudi::XYZVector const& direction ) {
+    auto mom   = composite->momentum().Vect();
+    auto mass2 = composite->momentum().mass2();
+    auto pt    = ( mom - direction * mom.Dot( direction ) ).r();
+    return pt + std::sqrt( mass2 + pt * pt );
+  }
+
 } // namespace
 
 class PrVeloHeavyFlavourTrackingChecker
-    : public Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&,
-                                               MCHits const&, Links const&, Links const& ),
+    : public Gaudi::Functional::Consumer<void( MCParticles const&, PVs const&, MCHits const&, Links const&,
+                                               Links const& ),
                                          LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>> {
 public:
-  using base_type = Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&,
-                                                      MCHits const&, Links const&, Links const& ),
-                                                LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>>;
-  using KeyValue  = base_type::KeyValue;
+  using base_type =
+      Gaudi::Functional::Consumer<void( MCParticles const&, PVs const&, MCHits const&, Links const&, Links const& ),
+                                  LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>>;
+  using KeyValue = base_type::KeyValue;
 
   // standard constructor
   PrVeloHeavyFlavourTrackingChecker( std::string const& name, ISvcLocator* pSvc )
       : base_type{name,
                   pSvc,
-                  {KeyValue{"MCParticles", LHCb::MCParticleLocation::Default}, KeyValue{"Composites", ""},
-                   KeyValue{"HeavyFlavourTracks", ""}, KeyValue{"MCHits", "/Event/MC/VP/Hits"}, KeyValue{"VPLinks", ""},
-                   KeyValue{"TrackLinks", ""}}} {}
+                  {KeyValue{"MCParticles", LHCb::MCParticleLocation::Default}, KeyValue{"PVs", ""},
+                   KeyValue{"MCHits", "/Event/MC/VP/Hits"}, KeyValue{"VPLinks", ""}, KeyValue{"TrackLinks", ""}}} {}
 
   // main execution
-  void operator()( MCParticles const&, Composites const&, HeavyFlavourTracks const&, MCHits const&, Links const&,
-                   Links const& ) const override;
+  void operator()( MCParticles const&, PVs const&, MCHits const&, Links const&, Links const& ) const override;
 
 private:
+  // data that might not be there
+  DataObjectReadHandle<Composites>         m_composites{this, "Composites", ""};
+  DataObjectReadHandle<HeavyFlavourTracks> m_tracks{this, "HeavyFlavourTracks", ""};
+
   // monitoring
   mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""};
 };
@@ -162,17 +217,17 @@ DECLARE_COMPONENT_WITH_ID( PrVeloHeavyFlavourTrackingChecker, "PrVeloHeavyFlavou
 //// implementation /////
 
 // main execution
-void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, Composites const& composites,
-                                                    HeavyFlavourTracks const& tracks, MCHits const& mchits,
+void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, PVs const& pvs, MCHits const& mchits,
                                                     Links const& vplinks, Links const& tracklinks ) const {
   // setup ntuple
   auto tuple = this->nTuple( "MCParticleTuple" );
 
   // table linking a LHCb::MCHit* to std::vector<LHCb::Detector::VPChannelID>> and vice versa
-  std::map<LHCb::MCHit const*, std::vector<unsigned int>> channelIDForMCHit;
-  vplinks.applyToAllLinks( [&channelIDForMCHit, &mchits]( unsigned int channelID, unsigned int mcHitKey, float ) {
-    channelIDForMCHit[mchits[mcHitKey]].emplace_back( channelID );
-  } );
+  std::map<LHCb::MCHit const*, std::vector<WeightedHit>> channelIDForMCHit;
+  vplinks.applyToAllLinks(
+      [&channelIDForMCHit, &mchits]( unsigned int channelID, unsigned int mcHitKey, float weight ) {
+        channelIDForMCHit[mchits[mcHitKey]].emplace_back( WeightedHit{channelID, weight} );
+      } );
 
   // fill for all tracks
   for ( auto const mcp : mcparts ) {
@@ -184,46 +239,117 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     // get MCHits for the candidate
     auto bhits       = getMCHits( mcps, mchits );
     auto first_mchit = bhits.size() ? bhits.front()->midPoint() : Gaudi::XYZPoint( 0., 0., 0. );
-    // first link MCParticle to daughters (pions) of composite
-    auto composite = matchToComposites( mcps, composites, tracklinks, mcparts );
-    auto has_comp  = composite.has_value();
+    // look at long tracks and see if pions match to all (as in 'reconstructed')
+    auto reconstructed = matchToTracks( mcps, tracklinks, mcparts, 0.7 );
+    // link MCParticle to daughters (pions) of composite
+    auto composite =
+        m_composites.exist() ? matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 ) : std::nullopt;
+    auto has_comp = composite.has_value();
     // get heavy flavour track
-    auto itrack    = has_comp ? getHeavyFlavourTrack( tracks, composite.value()->key() ) : std::nullopt;
+    auto const* tracks = m_tracks.getIfExists();
+    auto itrack    = ( has_comp && tracks ) ? getHeavyFlavourTrack( *tracks, composite.value()->key() ) : std::nullopt;
     auto has_track = has_comp ? itrack.has_value() : false;
-    auto has_hits  = has_track ? (int)tracks.scalar()[itrack.value()].nHits() > 0 : false;
+    auto has_hits  = has_track ? (int)tracks->scalar()[itrack.value()].nHits() > 0 : false;
+    int  pv_key    = has_track ? (int)tracks->scalar()[itrack.value()].primaryVertexIndex().cast() : -1;
     //// fill ntuple with info
     std::string base = "MCP_";
     // mc truth
     auto sc = tuple->column( base + "TRUEID", mcp->particleID().pid() );
     sc &= tuple->column( base + "P", mcp->momentum().P() );
+    sc &= tuple->column( base + "PT", mcp->momentum().Pt() );
+    sc &= tuple->column( base + "eta", mcp->momentum().Eta() );
+    sc &= tuple->column( base + "phi", mcp->momentum().Phi() );
+    sc &= tuple->column( base + "tx", mcp->momentum().Px() / mcp->momentum().Pz() );
+    sc &= tuple->column( base + "ty", mcp->momentum().Py() / mcp->momentum().Pz() );
+    auto ovtx = mcps.B->originVertex()->position();
+    sc &= tuple->column( base + "OriginVertex_x", ovtx.x() );
+    sc &= tuple->column( base + "OriginVertex_y", ovtx.y() );
+    sc &= tuple->column( base + "OriginVertex_z", ovtx.z() );
+    auto bvtx = mcps.B->endVertices().size() ? mcps.B->endVertices().at( mcps.B->endVertices().size() - 1 )->position()
+                                             : Gaudi::XYZPoint( 0., 0., -100 * Gaudi::Units::m );
+    sc &= tuple->column( base + "EndVertex_x", bvtx.x() );
+    sc &= tuple->column( base + "EndVertex_y", bvtx.y() );
+    sc &= tuple->column( base + "EndVertex_z", bvtx.z() );
+    auto tvtx = mcps.tau->endVertices().size()
+                    ? mcps.tau->endVertices().at( mcps.tau->endVertices().size() - 1 )->position()
+                    : Gaudi::XYZPoint( 0., 0., -110 * Gaudi::Units::m );
+    sc &= tuple->column( base + "Tau_EndVertex_x", tvtx.x() );
+    sc &= tuple->column( base + "Tau_EndVertex_y", tvtx.y() );
+    sc &= tuple->column( base + "Tau_EndVertex_z", tvtx.z() );
+    auto pionsmom = mcps.pions[0]->momentum() + mcps.pions[1]->momentum() + mcps.pions[2]->momentum();
+    sc &= tuple->column( base + "Pions_P", pionsmom.P() );
+    sc &= tuple->column( base + "Pions_PT", pionsmom.Pt() );
+    sc &= tuple->column( base + "Pions_Mass", pionsmom.mass() );
+    for ( int i = 0; i < (int)mcps.pions.size(); i++ ) {
+      std::string pbase = "Pion_" + std::to_string( i ) + "_";
+      auto        pion  = mcps.pions[i];
+      sc &= tuple->column( base + pbase + "TRUEID", pion->particleID().pid() );
+      sc &= tuple->column( base + pbase + "P", pion->momentum().P() );
+      sc &= tuple->column( base + pbase + "PT", pion->momentum().Pt() );
+      sc &= tuple->column( base + pbase + "has_Track", hasTrack( pion, tracklinks ) );
+    }
     // mc hits
     sc &= tuple->column( base + "has_VP_MCHits", (bool)bhits.size() );
     sc &= tuple->column( base + "VP_first_MCHit_x", first_mchit.x() );
     sc &= tuple->column( base + "VP_first_MCHit_y", first_mchit.y() );
     sc &= tuple->column( base + "VP_first_MCHit_z", first_mchit.z() );
+    // associated VP hits
+    std::vector<WeightedHit> hits_for_mchit =
+        bhits.size() ? channelIDForMCHit[bhits.front()] : std::vector<WeightedHit>();
+    if ( hits_for_mchit.size() > 1 ) std::sort( hits_for_mchit.begin(), hits_for_mchit.end() );
+    sc &= tuple->column( base + "VP_first_MCHit_has_VPHit", hits_for_mchit.size() > 0 );
+    sc &=
+        tuple->column( base + "VP_first_MCHit_VPHit_ChannelID", hits_for_mchit.size() ? hits_for_mchit.front().id : 0 );
     // composite info
+    sc &= tuple->column( base + "reconstructed", reconstructed );
     sc &= tuple->column( base + "has_Composite", has_comp );
     std::string compbase = base + "Composite_";
     auto        endvtx   = has_comp ? composite.value()->endVertex()->position() : Gaudi::XYZPoint( 0., 0., 0. );
+    float       vtxchi2  = has_comp ? composite.value()->endVertex()->chi2() : -1.f;
     sc &= tuple->column( compbase + "EndVertex_x", endvtx.x() );
     sc &= tuple->column( compbase + "EndVertex_y", endvtx.y() );
     sc &= tuple->column( compbase + "EndVertex_z", endvtx.z() );
+    sc &= tuple->column( compbase + "EndVertex_chi2", vtxchi2 );
     auto compmom =
         has_comp ? composite.value()->momentum() : Gaudi::LorentzVector( 0., 0., 0., -1 * Gaudi::Units::GeV );
     sc &= tuple->column( compbase + "P", compmom.R() );
     sc &= tuple->column( compbase + "PT", compmom.Pt() );
     sc &= tuple->column( compbase + "eta", compmom.Eta() );
     sc &= tuple->column( compbase + "phi", compmom.Phi() );
+    sc &= tuple->column( compbase + "Mass", compmom.mass() );
+    // composite daughters info
+    using Functors::Common::ImpactParameter;
+    using Functors::Common::Position;
+    using Functors::Common::ToLinAlg;
+    using Sel::getBestPV;
+    using Sel::Utils::impactParameterChi2;
+    for ( int i = 0; i < (int)mcps.pions.size(); i++ ) {
+      std::string pbase = "Pion_" + std::to_string( i ) + "_";
+      auto        pion  = mcps.recpions[i];
+      float       ip    = pion ? ImpactParameter{}(
+                            ToLinAlg{}( Position{}( has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ) ), pion )
+                            .cast()
+                      : -1.f;
+      float ipchi2 =
+          pion ? impactParameterChi2( pion, has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ).cast() : -1.f;
+      sc &= tuple->column( compbase + pbase + "match", (bool)pion );
+      sc &= tuple->column( compbase + pbase + "P", pion ? pion->momentum().R() : 0. );
+      sc &= tuple->column( compbase + pbase + "PT", pion ? pion->momentum().Pt() : 0. );
+      sc &= tuple->column( compbase + pbase + "IP", ip );
+      sc &= tuple->column( compbase + pbase + "IPChi2", ipchi2 );
+      sc &= tuple->column( compbase + pbase + "PIDK",
+                           pion ? pion->proto()->info( LHCb::ProtoParticle::additionalInfo::CombDLLk, 0. ) : 0. );
+    }
     // track info
     std::string hfbase = base + "HeavyFlavourTracking_";
-    sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks.scalar()[itrack.value()].nHits() : -1 );
+    sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks->scalar()[itrack.value()].nHits() : -1 );
     // track mc matching info
     bool fh_from_mcps = false;
     int  fh_trueid    = 0;
     if ( has_hits ) {
       float              max_weight = 0.f;
       LHCb::MCHit const* mchit      = nullptr;
-      unsigned int       hid        = tracks.scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID();
+      unsigned int       hid        = tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID();
       vplinks.applyToLinks( hid, [&max_weight, &mchits, &mchit]( unsigned int, unsigned int mcHitKey, float weight ) {
         if ( weight > max_weight ) {
           max_weight = weight;
@@ -236,13 +362,40 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
       }
     }
     sc &= tuple->column( hfbase + "Track_first_hit_x",
-                         has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. );
+                         has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. );
     sc &= tuple->column( hfbase + "Track_first_hit_y",
-                         has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. );
+                         has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. );
     sc &= tuple->column( hfbase + "Track_first_hit_z",
-                         has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. );
+                         has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. );
+    sc &= tuple->column( hfbase + "Track_first_hit_LHCbID",
+                         has_hits ? (unsigned int)tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID() : 0 );
     sc &= tuple->column( hfbase + "Track_first_hit_from_MCPs", fh_from_mcps );
     sc &= tuple->column( hfbase + "Track_first_hit_MC_TRUEID", fh_trueid );
+    // more advanced variables for analysis
+    auto  truedir    = mcps.B->momentum().Vect().Unit();
+    auto  mcorr_true = has_comp ? getMcorr( composite.value(), truedir ) : 0.f;
+    float mcorr_fh   = 0.f;
+    float mcorr_ps   = 0.f;
+    auto  slopes     = Gaudi::XYZVector( 0., 0., 1. );
+    auto  pvpos      = Gaudi::XYZPoint( 0., 0., -10 * Gaudi::Units::m );
+    if ( has_track ) {
+      auto track = tracks->scalar()[itrack.value()];
+      auto state = track.originState();
+      slopes.SetXYZ( state.tx().cast(), state.ty().cast(), 1. );
+      auto fhdir = slopes.Unit();
+      pvpos.SetXYZ( state.x().cast(), state.y().cast(), state.z().cast() );
+      auto pvsv = ( composite.value()->endVertex()->position() - pvpos ).Unit();
+      mcorr_fh  = getMcorr( composite.value(), fhdir );
+      mcorr_ps  = getMcorr( composite.value(), pvsv );
+    }
+    sc &= tuple->column( hfbase + "tx", slopes.x() );
+    sc &= tuple->column( hfbase + "ty", slopes.y() );
+    sc &= tuple->column( hfbase + "true_dir_mcorr", mcorr_true );
+    sc &= tuple->column( hfbase + "first_hit_mcorr", mcorr_fh );
+    sc &= tuple->column( hfbase + "pvsv_mcorr", mcorr_ps );
+    sc &= tuple->column( hfbase + "PV_x", pvpos.x() );
+    sc &= tuple->column( hfbase + "PV_y", pvpos.y() );
+    sc &= tuple->column( hfbase + "PV_z", pvpos.z() );
     // write result
     sc.andThen( [&] { return tuple->write(); } ).ignore();
   }
-- 
GitLab


From f6ce5ae8e031b57b16de2b3585014a4872785382 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <mveghel@cern.ch>
Date: Wed, 9 Nov 2022 18:49:39 +0100
Subject: [PATCH 07/14] added btracking functor

---
 Phys/FunctorCore/include/Functors/Composite.h | 39 +++++++++++++++++++
 .../include/Functors/JIT_includes.h           |  1 +
 Phys/FunctorCore/python/Functors/__init__.py  | 20 ++++++++++
 .../tests/src/InstantiateFunctors.cpp         |  6 +++
 .../PrKernel/PrVeloHeavyFlavourUtils.h        | 27 ++++++++++++-
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 11 ++----
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   | 16 ++++----
 7 files changed, 103 insertions(+), 17 deletions(-)

diff --git a/Phys/FunctorCore/include/Functors/Composite.h b/Phys/FunctorCore/include/Functors/Composite.h
index 284ba53cc60..cde1e593fb8 100755
--- a/Phys/FunctorCore/include/Functors/Composite.h
+++ b/Phys/FunctorCore/include/Functors/Composite.h
@@ -321,4 +321,43 @@ namespace Functors::Composite {
     }
   };
 
+  namespace BTracking {
+    /** @brief Calculate the corrected mass using heavy flavour tracking info. */
+    struct CorrectedMass : public Function {
+      static constexpr auto name() { return "BTracking::CorrectedMass"; }
+
+      template <typename HeavyFlavourTracks, typename Composite>
+      auto operator()( HeavyFlavourTracks const& tracks, Composite const& composite ) const {
+
+        using FType = SIMDWrapper::scalar::float_v;
+        using Sel::Utils::mass2;
+        using Sel::Utils::threeMomentum;
+        using std::sqrt;
+        auto const& comp = Sel::Utils::deref_if_ptr( composite );
+
+        // first get the right heavy flavour track
+        auto get_track = [&tracks]( int index ) -> std::optional<int> {
+          for ( auto track : tracks.scalar() ) {
+            if ( track.compositeIndex() == index ) return track.indices().cast();
+          }
+          return std::nullopt;
+        };
+        auto itrack = get_track( composite.key() );
+        if ( !itrack.has_value() ) return FType{0};
+
+        // direction using first hit (or end vertex if not available, as saved in heavy flavour track)
+        auto const state = tracks.scalar()[itrack.value()].originState();
+        auto const d     = LHCb::LinAlg::Vec<FType, 3>{state.tx(), state.ty(), 1.};
+
+        // Get the pT variable that we need
+        auto const mom  = threeMomentum( comp );
+        auto const perp = mom - d * ( dot( mom, d ) / d.mag2() );
+        auto const pt   = perp.mag();
+
+        // Calculate the corrected mass
+        return pt + sqrt( mass2( comp ) + pt * pt );
+      }
+    };
+  } // namespace BTracking
+
 } // namespace Functors::Composite
diff --git a/Phys/FunctorCore/include/Functors/JIT_includes.h b/Phys/FunctorCore/include/Functors/JIT_includes.h
index a83d656b988..8e3554e4d50 100644
--- a/Phys/FunctorCore/include/Functors/JIT_includes.h
+++ b/Phys/FunctorCore/include/Functors/JIT_includes.h
@@ -65,5 +65,6 @@
 #include "SelKernel/VertexRelation.h"
 // PrKernel
 #include "PrKernel/PrSelection.h"
+#include "PrKernel/PrVeloHeavyFlavourUtils.h"
 // RelTables
 #include "Event/TableView.h"
diff --git a/Phys/FunctorCore/python/Functors/__init__.py b/Phys/FunctorCore/python/Functors/__init__.py
index a542a25e872..63e58f5cebc 100644
--- a/Phys/FunctorCore/python/Functors/__init__.py
+++ b/Phys/FunctorCore/python/Functors/__init__.py
@@ -1372,6 +1372,14 @@ MVA technique to establish pion mass hypothesis which keep in consideration the
 General information on PID could be found at the following `TWiki page <https://twiki.cern.ch/twiki/bin/view/LHCbPhysics/ChargedPID>`_
 """)
 
+_BPVCORRM = Functor(
+    '_BPVCORRM', 'Composite::CorrectedMass',
+    'Compute the corrected mass of the composite using the associated [primary] vertex.'
+)
+_BTRACKING_BPVCORRM = Functor(
+    '_BTRACKING_BPVCORRM', 'Composite::BTracking::CorrectedMass',
+    'Compute the corrected mass of the composite using the associated heavy flavour track.'
+)
 
 def BPVCORRM(Vertices: DataHandle):
     """ Compute the corrected mass of the composite using the :py:func:`~BPV`.
@@ -1405,6 +1413,18 @@ def BPVCORRMERR(Vertices: DataHandle):
     return _BPVCORRMERR.bind(BPV(Vertices), FORWARDARGS)
 
 
+def BTRACKING_BPVCORRM(HeavyFlavourTracks: DataHandle):
+    """ Compute the corrected mass of the composite using the associated heavy flavour track.
+
+    Functor's call operator expects a composite like object
+
+    Args:
+        Vertices: DataHandle of the primary vertices
+
+    """
+    return _BTRACKING_BPVCORRM.bind(TES(HeavyFlavourTracks), FORWARDARGS)
+
+
 VTX_FDCHI2 = Functor(
     'VTX_FDCHI2', 'Composite::FlightDistanceChi2ToVertex',
     '''Return the flight distance chi2 w.r.t. the given vertex.''')
diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
index 7b4b936f371..7abeacb2a0e 100644
--- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
+++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
@@ -208,6 +208,12 @@ using cmerr5 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHC
 // FIXME shouldn't these work?
 // using cmerr6 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHCb::RecVertex, LHCb::Particle const&>;
 
+//
+// BTracking::CorrectedMass( heavyflavourtracks, composite )
+//
+using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass,
+                                   LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks, LHCb::Particle const&>;
+
 //
 // MassWithHypotheses( vertices, composite )
 //
diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
index 941c308233b..eb5c7eff50b 100644
--- a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
+++ b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
@@ -24,7 +24,7 @@ namespace LHCb::Pr::Velo {
     namespace Tag {
       // hit info
       struct LHCbID : Event::lhcbid_field {};
-      struct HitPosition : Event::vec3_field {};
+      struct HitPosition : Event::Vec_field<3> {};
       struct Hits : Event::vector_field<Event::struct_field<LHCbID, HitPosition>> {};
 
       // crossed sensor info
@@ -55,10 +55,22 @@ namespace LHCb::Pr::Velo {
       using Composite  = LHCb::Particle;
       using PVs        = LHCb::Event::PV::PrimaryVertexContainer;
 
-      struct HeavyFlavourTracks : Tag::Track_t<HeavyFlavourTracks> {
+      class HeavyFlavourTracks : public Tag::Track_t<HeavyFlavourTracks> {
+      public:
         using base_t = typename Tag::Track_t<HeavyFlavourTracks>;
+        using base_t::allocator_type;
         using base_t::base_t;
 
+        HeavyFlavourTracks( UniqueIDGenerator const& unique_id_gen,
+                            Zipping::ZipFamilyNumber zip_identifier = Zipping::generateZipIdentifier(),
+                            allocator_type           alloc          = {} )
+            : base_t{std::move( zip_identifier ), std::move( alloc )}, m_unique_id_gen_tag{unique_id_gen.tag()} {}
+
+        // special constructor for zipping
+        HeavyFlavourTracks( Zipping::ZipFamilyNumber zip_identifier, HeavyFlavourTracks const& other )
+            : base_t{std::move( zip_identifier ), other} {}
+
+        // proxy
         template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType>
         struct TrackProxy : Event::Proxy<simd, behaviour, ContainerType> {
           using base_t = typename Event::Proxy<simd, behaviour, ContainerType>;
@@ -89,6 +101,17 @@ namespace LHCb::Pr::Velo {
 
         template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType>
         using proxy_type = TrackProxy<simd, behaviour, ContainerType>;
+
+        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour>
+        using reference = LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, HeavyFlavourTracks>>;
+        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour>
+        using const_reference = const LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, const HeavyFlavourTracks>>;
+
+        auto const& unique_id_gen_tag() const { return m_unique_id_gen_tag; }
+
+      private:
+        // keep the identifier of the generator used to build this container
+        boost::uuids::uuid m_unique_id_gen_tag;
       };
     } // namespace v1
 
diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index a7caf773073..94d3418ee9c 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -34,8 +34,8 @@
 namespace {
   // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA)
   using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks;
-  using Composites         = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composites;
   using Composite          = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite;
+  using Composites         = Composite::Selection;
   using PVs                = LHCb::Pr::Velo::HeavyFlavourTracking::v1::PVs;
   using MCParticles        = LHCb::MCParticles;
   using MCHits             = LHCb::MCHits;
@@ -326,12 +326,9 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     for ( int i = 0; i < (int)mcps.pions.size(); i++ ) {
       std::string pbase = "Pion_" + std::to_string( i ) + "_";
       auto        pion  = mcps.recpions[i];
-      float       ip    = pion ? ImpactParameter{}(
-                            ToLinAlg{}( Position{}( has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ) ), pion )
-                            .cast()
-                      : -1.f;
-      float ipchi2 =
-          pion ? impactParameterChi2( pion, has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ).cast() : -1.f;
+      float       ip =
+          ( pion && has_track ) ? ImpactParameter{}( ToLinAlg{}( Position{}( pvs[pv_key] ) ), pion ).cast() : -1.f;
+      float ipchi2 = ( pion && has_track ) ? impactParameterChi2( pion, pvs[pv_key] ).cast() : -1.f;
       sc &= tuple->column( compbase + pbase + "match", (bool)pion );
       sc &= tuple->column( compbase + pbase + "P", pion ? pion->momentum().R() : 0. );
       sc &= tuple->column( compbase + pbase + "PT", pion ? pion->momentum().Pt() : 0. );
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index 5fab6af5c4b..24c75c8b8fe 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -33,7 +33,7 @@ namespace LHCb::Pr::Velo {
     // miscellaneous
     namespace TrackTag = LHCb::Pr::Velo::HeavyFlavourTracking::Tag;
     using FType        = SIMDWrapper::scalar::float_v;
-    using Position     = LHCb::LinAlg::Vec3<FType>;
+    using Position     = LHCb::LinAlg::Vec<FType, 3>;
     using SensorID     = LHCb::Detector::VPChannelID::SensorID;
     using ChannelID    = LHCb::Detector::VPChannelID;
 
@@ -82,7 +82,7 @@ namespace LHCb::Pr::Velo {
       ChannelID id;
       // sorting in z first, then distance to search window
       friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) {
-        return std::pair{lhs.position.z.cast(), lhs.d.cast()} < std::pair{rhs.position.z.cast(), rhs.d.cast()};
+        return std::pair{lhs.position.z().cast(), lhs.d.cast()} < std::pair{rhs.position.z().cast(), rhs.d.cast()};
       }
       // equal just by unique channel id
       friend bool operator==( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.id == rhs.id; }
@@ -119,7 +119,7 @@ namespace LHCb::Pr::Velo {
         // its associated module (to not be affected by borders between sensors)
         auto crossed_module = sensor.module();
         // go to sensor location
-        pos = pos + ( sensor.z - pos.z.cast() ) * slopes;
+        pos = pos + ( sensor.z - pos.z().cast() ) * slopes;
         // look if this sensor (and associated module) has hits and if they are in search window
         for ( auto hit : hits.scalar() ) {
           auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
@@ -127,9 +127,9 @@ namespace LHCb::Pr::Velo {
             auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
             // make sure they are in the same plane up to some tolerance
             // (not all sensors in a module are exactly)
-            if ( abs( hitpos.vec3().z.cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) {
-              auto d = ( hitpos.vec3() - pos ).mag();
-              if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos.vec3(), id} );
+            if ( abs( hitpos.z().cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) {
+              auto d = ( hitpos - pos ).mag();
+              if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos, id} );
             }
           }
         }
@@ -214,7 +214,7 @@ namespace LHCb::Pr::Velo {
 
       // primary vertex
       using Sel::Utils::endVertexPos;
-      auto const& PV = Sel::getBestPV( *composite, pvs );
+      auto const& PV = Sel::getBestPV( Sel::Utils::deref_if_ptr( composite ), pvs );
       auto const& pv = endVertexPos( PV );
       track.field<TrackTag::PrimaryVertexIndex>().set( PV.key() );
 
@@ -245,7 +245,7 @@ namespace LHCb::Pr::Velo {
       auto direction    = found_hits.size() ? found_hits.front().position - pv : fd;
       auto origin_state = track.field<TrackTag::OriginState>();
       origin_state.setPosition( pv );
-      origin_state.setDirection( direction.x / direction.z, direction.y / direction.z );
+      origin_state.setDirection( direction.x() / direction.z(), direction.y() / direction.z() );
       origin_state.setQOverP( 1. / composite->p() );
 
       // fitting
-- 
GitLab


From 77bcba1e1c4252553254ff5b192080f3f9bcb550 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Tue, 25 Jul 2023 19:03:06 +0200
Subject: [PATCH 08/14] move to v1 event model

---
 Phys/FunctorCore/include/Functors/Composite.h |  17 +--
 .../include/Functors/JIT_includes.h           |   1 +
 .../tests/src/InstantiateFunctors.cpp         |   4 +-
 .../PrKernel/PrVeloHeavyFlavourUtils.h        | 120 ------------------
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp |  94 +++++++-------
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   | 116 +++++++++--------
 6 files changed, 121 insertions(+), 231 deletions(-)
 delete mode 100644 Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h

diff --git a/Phys/FunctorCore/include/Functors/Composite.h b/Phys/FunctorCore/include/Functors/Composite.h
index cde1e593fb8..9a00ace18f7 100755
--- a/Phys/FunctorCore/include/Functors/Composite.h
+++ b/Phys/FunctorCore/include/Functors/Composite.h
@@ -326,8 +326,8 @@ namespace Functors::Composite {
     struct CorrectedMass : public Function {
       static constexpr auto name() { return "BTracking::CorrectedMass"; }
 
-      template <typename HeavyFlavourTracks, typename Composite>
-      auto operator()( HeavyFlavourTracks const& tracks, Composite const& composite ) const {
+      template <typename Composite2TrackRelations, typename Composite>
+      auto operator()( Composite2TrackRelations const& relations, Composite const& composite ) const {
 
         using FType = SIMDWrapper::scalar::float_v;
         using Sel::Utils::mass2;
@@ -336,17 +336,12 @@ namespace Functors::Composite {
         auto const& comp = Sel::Utils::deref_if_ptr( composite );
 
         // first get the right heavy flavour track
-        auto get_track = [&tracks]( int index ) -> std::optional<int> {
-          for ( auto track : tracks.scalar() ) {
-            if ( track.compositeIndex() == index ) return track.indices().cast();
-          }
-          return std::nullopt;
-        };
-        auto itrack = get_track( composite.key() );
-        if ( !itrack.has_value() ) return FType{0};
+        auto track_range = relations.relations( &comp );
+        if ( track_range.empty() ) return FType{0};
+        auto track = track_range.front().to();
 
         // direction using first hit (or end vertex if not available, as saved in heavy flavour track)
-        auto const state = tracks.scalar()[itrack.value()].originState();
+        auto const state = track->firstState();
         auto const d     = LHCb::LinAlg::Vec<FType, 3>{state.tx(), state.ty(), 1.};
 
         // Get the pT variable that we need
diff --git a/Phys/FunctorCore/include/Functors/JIT_includes.h b/Phys/FunctorCore/include/Functors/JIT_includes.h
index 8e3554e4d50..b1510005c2d 100644
--- a/Phys/FunctorCore/include/Functors/JIT_includes.h
+++ b/Phys/FunctorCore/include/Functors/JIT_includes.h
@@ -68,3 +68,4 @@
 #include "PrKernel/PrVeloHeavyFlavourUtils.h"
 // RelTables
 #include "Event/TableView.h"
+#include "Relations/Relation2D.h"
diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
index 7abeacb2a0e..6c08a94cff3 100644
--- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
+++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
@@ -209,10 +209,10 @@ using cmerr5 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHC
 // using cmerr6 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHCb::RecVertex, LHCb::Particle const&>;
 
 //
-// BTracking::CorrectedMass( heavyflavourtracks, composite )
+// BTracking::CorrectedMass( heavyflavourtrackrelations, composite )
 //
 using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass,
-                                   LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks, LHCb::Particle const&>;
+                                   LHCb::Relation2D<LHCb::Particle, LHCb::Event::v1::Track>, LHCb::Particle const&>;
 
 //
 // MassWithHypotheses( vertices, composite )
diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
deleted file mode 100644
index eb5c7eff50b..00000000000
--- a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h
+++ /dev/null
@@ -1,120 +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.                                       *
-\*****************************************************************************/
-
-#include "Event/Particle.h"
-#include "Event/Particle_v2.h"
-#include "Event/PrTracksTag.h"
-#include "Event/PrimaryVertices.h"
-#include "Event/RelationTable.h"
-#include "VPDet/DeVP.h"
-
-namespace LHCb::Pr::Velo {
-
-  namespace HeavyFlavourTracking {
-
-    // main heavy flavour track storage
-    namespace Tag {
-      // hit info
-      struct LHCbID : Event::lhcbid_field {};
-      struct HitPosition : Event::Vec_field<3> {};
-      struct Hits : Event::vector_field<Event::struct_field<LHCbID, HitPosition>> {};
-
-      // crossed sensor info
-      struct nSensors : Event::int_field {};
-
-      // geometry
-      struct OriginState : Event::state_field {};
-
-      // track quality info
-      struct Chi2 : Event::float_field {};
-
-      // (hopefully temporary) fix with SOACollection with keys/indices to non SoA objects
-      struct CompositeIndex : Event::int_field {};
-      struct PrimaryVertexIndex : Event::int_field {};
-
-      template <typename T>
-      using Track_t = Event::SOACollection<T, CompositeIndex, PrimaryVertexIndex, Hits, nSensors, OriginState, Chi2>;
-
-      /*
-      // main storage: relating track to associated composite and PV in full SoA event model
-      template <typename Composites, typename PVs>
-      using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>;
-      */
-    } // namespace Tag
-
-    namespace v1 {
-      using Composites = LHCb::Particles;
-      using Composite  = LHCb::Particle;
-      using PVs        = LHCb::Event::PV::PrimaryVertexContainer;
-
-      class HeavyFlavourTracks : public Tag::Track_t<HeavyFlavourTracks> {
-      public:
-        using base_t = typename Tag::Track_t<HeavyFlavourTracks>;
-        using base_t::allocator_type;
-        using base_t::base_t;
-
-        HeavyFlavourTracks( UniqueIDGenerator const& unique_id_gen,
-                            Zipping::ZipFamilyNumber zip_identifier = Zipping::generateZipIdentifier(),
-                            allocator_type           alloc          = {} )
-            : base_t{std::move( zip_identifier ), std::move( alloc )}, m_unique_id_gen_tag{unique_id_gen.tag()} {}
-
-        // special constructor for zipping
-        HeavyFlavourTracks( Zipping::ZipFamilyNumber zip_identifier, HeavyFlavourTracks const& other )
-            : base_t{std::move( zip_identifier ), other} {}
-
-        // proxy
-        template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType>
-        struct TrackProxy : Event::Proxy<simd, behaviour, ContainerType> {
-          using base_t = typename Event::Proxy<simd, behaviour, ContainerType>;
-          using base_t::Proxy;
-
-          // composite / pv info
-          [[nodiscard]] auto compositeIndex() const { return this->template get<Tag::CompositeIndex>(); }
-          [[nodiscard]] auto primaryVertexIndex() const { return this->template get<Tag::PrimaryVertexIndex>(); }
-
-          // hit info
-          [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); }
-          [[nodiscard]] auto hitID( std::size_t i ) const {
-            return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>();
-          }
-          [[nodiscard]] auto hitPosition( std::size_t i ) const {
-            return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>();
-          }
-
-          // crossed sensor info
-          [[nodiscard]] auto nSensors() const { return this->template get<Tag::nSensors>(); }
-
-          // geometry
-          [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); }
-
-          // track quality info
-          [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); }
-        };
-
-        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType>
-        using proxy_type = TrackProxy<simd, behaviour, ContainerType>;
-
-        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour>
-        using reference = LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, HeavyFlavourTracks>>;
-        template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour>
-        using const_reference = const LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, const HeavyFlavourTracks>>;
-
-        auto const& unique_id_gen_tag() const { return m_unique_id_gen_tag; }
-
-      private:
-        // keep the identifier of the generator used to build this container
-        boost::uuids::uuid m_unique_id_gen_tag;
-      };
-    } // namespace v1
-
-  } // namespace HeavyFlavourTracking
-
-} // namespace LHCb::Pr::Velo
diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index 94d3418ee9c..711dbf61ead 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -13,7 +13,7 @@
 #include "Event/MCHit.h"
 #include "Event/MCParticle.h"
 #include "Event/RawEvent.h"
-#include "Event/Track_v3.h"
+#include "Event/Track.h"
 #include "Event/VPDigit.h"
 #include "Event/VPFullCluster.h"
 #include "Functors/Common.h"
@@ -21,7 +21,7 @@
 #include "GaudiAlg/GaudiTupleAlg.h"
 #include "Linker/LinkedFrom.h"
 #include "Linker/LinkedTo.h"
-#include "PrKernel/PrVeloHeavyFlavourUtils.h"
+#include "Relations/Relation2D.h"
 #include "SelKernel/Utilities.h"
 #include "SelKernel/VertexRelation.h"
 
@@ -32,14 +32,16 @@
  */
 
 namespace {
-  // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA)
-  using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks;
-  using Composite          = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite;
-  using Composites         = Composite::Selection;
-  using PVs                = LHCb::Pr::Velo::HeavyFlavourTracking::v1::PVs;
-  using MCParticles        = LHCb::MCParticles;
-  using MCHits             = LHCb::MCHits;
-  using Links              = LHCb::LinksByKey;
+  // in/outputs
+  using Composite    = LHCb::Particle;
+  using Composites   = LHCb::Particle::Selection;
+  using PVs          = LHCb::Event::PV::PrimaryVertexContainer;
+  using Track        = LHCb::Event::v1::Track;
+  using Tracks       = Track::Container;
+  using P2TRelations = LHCb::Relation2D<LHCb::Particle, Track>;
+  using MCParticles  = LHCb::MCParticles;
+  using MCHits       = LHCb::MCHits;
+  using Links        = LHCb::LinksByKey;
 
   // helper classes
   struct DecayMCParticles {
@@ -165,14 +167,6 @@ namespace {
     return hits;
   }
 
-  // get track
-  std::optional<int> getHeavyFlavourTrack( HeavyFlavourTracks const& tracks, int key ) {
-    for ( auto track : tracks.scalar() ) {
-      if ( track.compositeIndex() == key ) return track.indices().cast();
-    }
-    return std::nullopt;
-  }
-
   // composite observables
   float getMcorr( Composite const* composite, Gaudi::XYZVector const& direction ) {
     auto mom   = composite->momentum().Vect();
@@ -205,8 +199,8 @@ public:
 
 private:
   // data that might not be there
-  DataObjectReadHandle<Composites>         m_composites{this, "Composites", ""};
-  DataObjectReadHandle<HeavyFlavourTracks> m_tracks{this, "HeavyFlavourTracks", ""};
+  DataObjectReadHandle<Composites>   m_composites{this, "Composites", ""};
+  DataObjectReadHandle<P2TRelations> m_relations{this, "Composite2HeavyFlavourTrackRelations", ""};
 
   // monitoring
   mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""};
@@ -242,15 +236,20 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     // look at long tracks and see if pions match to all (as in 'reconstructed')
     auto reconstructed = matchToTracks( mcps, tracklinks, mcparts, 0.7 );
     // link MCParticle to daughters (pions) of composite
-    auto composite =
-        m_composites.exist() ? matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 ) : std::nullopt;
+    std::optional<Composite const*> composite = std::nullopt;
+    if ( m_composites.exist() ) {
+      composite = matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 );
+    }
     auto has_comp = composite.has_value();
     // get heavy flavour track
-    auto const* tracks = m_tracks.getIfExists();
-    auto itrack    = ( has_comp && tracks ) ? getHeavyFlavourTrack( *tracks, composite.value()->key() ) : std::nullopt;
-    auto has_track = has_comp ? itrack.has_value() : false;
-    auto has_hits  = has_track ? (int)tracks->scalar()[itrack.value()].nHits() > 0 : false;
-    int  pv_key    = has_track ? (int)tracks->scalar()[itrack.value()].primaryVertexIndex().cast() : -1;
+    auto const*  relations = m_relations.getIfExists();
+    Track const* track     = nullptr;
+    if ( has_comp && relations ) {
+      auto trackrange = relations->relations( composite.value() );
+      if ( trackrange.size() ) track = trackrange.front();
+    }
+    auto has_track = track ? true : false;
+    auto has_hits  = has_track ? track->nHits() > 0 : false;
     //// fill ntuple with info
     std::string base = "MCP_";
     // mc truth
@@ -324,11 +323,15 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     using Sel::getBestPV;
     using Sel::Utils::impactParameterChi2;
     for ( int i = 0; i < (int)mcps.pions.size(); i++ ) {
-      std::string pbase = "Pion_" + std::to_string( i ) + "_";
-      auto        pion  = mcps.recpions[i];
-      float       ip =
-          ( pion && has_track ) ? ImpactParameter{}( ToLinAlg{}( Position{}( pvs[pv_key] ) ), pion ).cast() : -1.f;
-      float ipchi2 = ( pion && has_track ) ? impactParameterChi2( pion, pvs[pv_key] ).cast() : -1.f;
+      std::string pbase  = "Pion_" + std::to_string( i ) + "_";
+      auto        pion   = mcps.recpions[i];
+      float       ip     = -1.f;
+      float       ipchi2 = -1.f;
+      if ( pion && has_track ) {
+        auto pv = getBestPV( pion, pvs );
+        ip      = ImpactParameter{}( ToLinAlg{}( Position{}( pv ) ), pion ).cast();
+        ipchi2  = impactParameterChi2( pion, pv ).cast();
+      }
       sc &= tuple->column( compbase + pbase + "match", (bool)pion );
       sc &= tuple->column( compbase + pbase + "P", pion ? pion->momentum().R() : 0. );
       sc &= tuple->column( compbase + pbase + "PT", pion ? pion->momentum().Pt() : 0. );
@@ -339,14 +342,16 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     }
     // track info
     std::string hfbase = base + "HeavyFlavourTracking_";
-    sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks->scalar()[itrack.value()].nHits() : -1 );
+    sc &= tuple->column( hfbase + "Track_nSensors",
+                         has_track ? track->info( LHCb::Event::Enum::Track::AdditionalInfo::VeloNSensors, 0 ) : -1. );
+    sc &= tuple->column( hfbase + "Track_nHits", has_track ? track->nHits() : -1 );
     // track mc matching info
     bool fh_from_mcps = false;
     int  fh_trueid    = 0;
     if ( has_hits ) {
       float              max_weight = 0.f;
       LHCb::MCHit const* mchit      = nullptr;
-      unsigned int       hid        = tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID();
+      unsigned int       hid        = track->lhcbIDs().front().vpID();
       vplinks.applyToLinks( hid, [&max_weight, &mchits, &mchit]( unsigned int, unsigned int mcHitKey, float weight ) {
         if ( weight > max_weight ) {
           max_weight = weight;
@@ -358,14 +363,14 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
         fh_from_mcps = ( mchit->mcParticle() == mcps.B ) || ( mchit->mcParticle() == mcps.tau );
       }
     }
-    sc &= tuple->column( hfbase + "Track_first_hit_x",
-                         has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. );
-    sc &= tuple->column( hfbase + "Track_first_hit_y",
-                         has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. );
-    sc &= tuple->column( hfbase + "Track_first_hit_z",
-                         has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. );
+    decltype( track->position() ) fh_pos( 0., 0., 0. );
+    LHCb::State const* fh_state = track ? track->stateAt( LHCb::State::Location::FirstMeasurement ) : nullptr;
+    if ( fh_state ) fh_pos = fh_state->position();
+    sc &= tuple->column( hfbase + "Track_first_hit_x", fh_pos.x() );
+    sc &= tuple->column( hfbase + "Track_first_hit_y", fh_pos.y() );
+    sc &= tuple->column( hfbase + "Track_first_hit_z", fh_pos.z() );
     sc &= tuple->column( hfbase + "Track_first_hit_LHCbID",
-                         has_hits ? (unsigned int)tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID() : 0 );
+                         has_hits ? (unsigned int)track->lhcbIDs().front().vpID() : 0 );
     sc &= tuple->column( hfbase + "Track_first_hit_from_MCPs", fh_from_mcps );
     sc &= tuple->column( hfbase + "Track_first_hit_MC_TRUEID", fh_trueid );
     // more advanced variables for analysis
@@ -376,11 +381,10 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     auto  slopes     = Gaudi::XYZVector( 0., 0., 1. );
     auto  pvpos      = Gaudi::XYZPoint( 0., 0., -10 * Gaudi::Units::m );
     if ( has_track ) {
-      auto track = tracks->scalar()[itrack.value()];
-      auto state = track.originState();
-      slopes.SetXYZ( state.tx().cast(), state.ty().cast(), 1. );
+      auto state = track->firstState();
+      slopes.SetXYZ( state.tx(), state.ty(), 1. );
       auto fhdir = slopes.Unit();
-      pvpos.SetXYZ( state.x().cast(), state.y().cast(), state.z().cast() );
+      pvpos.SetXYZ( state.x(), state.y(), state.z() );
       auto pvsv = ( composite.value()->endVertex()->position() - pvpos ).Unit();
       mcorr_fh  = getMcorr( composite.value(), fhdir );
       mcorr_ps  = getMcorr( composite.value(), pvsv );
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index 24c75c8b8fe..ce3f506d4ec 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -12,30 +12,33 @@
 #include "DetDesc/GenericConditionAccessorHolder.h"
 #include "DetDesc/IConditionDerivationMgr.h"
 #include "Event/PrHits.h"
+#include "Event/Track.h"
 #include "GaudiAlg/Transformer.h"
 #include "Kernel/LHCbID.h"
 #include "LHCbMath/MatVec.h"
-#include "PrKernel/PrVeloHeavyFlavourUtils.h"
+#include "Relations/Relation2D.h"
 #include "SelKernel/Utilities.h"
 #include "SelKernel/VertexRelation.h"
+#include "VPDet/DeVP.h"
 #include "VPDet/VPDetPaths.h"
-#include "VeloKalmanHelpers.h"
 
 namespace LHCb::Pr::Velo {
 
   namespace {
 
-    // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA)
-    using namespace LHCb::Pr::Velo::HeavyFlavourTracking::v1;
-    using VeloHits = LHCb::Pr::VP::Hits;
-    using Output   = HeavyFlavourTracks;
+    // in/outputs
+    using VeloHits  = LHCb::Pr::VP::Hits;
+    using PVs       = LHCb::Event::PV::PrimaryVertexContainer;
+    using Track     = LHCb::Event::v1::Track;
+    using Tracks    = Track::Container;
+    using Relations = LHCb::Relation2D<LHCb::Particle, Track>;
+    using Output    = std::tuple<Tracks, Relations>;
 
     // miscellaneous
-    namespace TrackTag = LHCb::Pr::Velo::HeavyFlavourTracking::Tag;
-    using FType        = SIMDWrapper::scalar::float_v;
-    using Position     = LHCb::LinAlg::Vec<FType, 3>;
-    using SensorID     = LHCb::Detector::VPChannelID::SensorID;
-    using ChannelID    = LHCb::Detector::VPChannelID;
+    using FType     = SIMDWrapper::scalar::float_v;
+    using Position  = LHCb::LinAlg::Vec<FType, 3>;
+    using SensorID  = LHCb::Detector::VPChannelID::SensorID;
+    using ChannelID = LHCb::Detector::VPChannelID;
 
     // to store relevant sensor information
     struct SensorInfo {
@@ -150,9 +153,9 @@ namespace LHCb::Pr::Velo {
    */
 
   class HeavyFlavourTrackFinder
-      : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const&,
-                                                      Sensors const& ),
-                                              LHCb::DetDesc::usesConditions<DeVP, Sensors>> {
+      : public Gaudi::Functional::MultiTransformer<Output( LHCb::Particle::Range const&, PVs const&, VeloHits const&,
+                                                           DeVP const&, Sensors const& ),
+                                                   LHCb::DetDesc::usesConditions<DeVP, Sensors>> {
   public:
     // standard constructor
     HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator );
@@ -161,7 +164,8 @@ namespace LHCb::Pr::Velo {
     StatusCode initialize() override;
 
     // main function/operator
-    Output operator()( Composites const&, PVs const&, VeloHits const&, DeVP const&, Sensors const& ) const override;
+    Output operator()( LHCb::Particle::Range const&, PVs const&, VeloHits const&, DeVP const&,
+                       Sensors const& ) const override;
 
   private:
     // properties
@@ -179,44 +183,39 @@ namespace LHCb::Pr::Velo {
 
   // constructor
   HeavyFlavourTrackFinder::HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator )
-      : Transformer( name, pSvcLocator,
-                     {KeyValue( "Composites", "" ), KeyValue( "PVs", "" ), KeyValue( "Hits", "" ),
-                      KeyValue( "DeVP", LHCb::Det::VP::det_path ),
-#ifdef USE_DD4HEP
-                      KeyValue( "Sensors", {"/world:AlgorithmSpecific-" + name + "-sensorinfos"} )},
-#else
-                      KeyValue( "Sensors", {"AlgorithmSpecific-" + name + "-sensorinfos"} )},
-#endif
-
-                     {KeyValue{"OutputLocation", ""}} ) {
-  }
+      : MultiTransformer( name, pSvcLocator,
+                          {KeyValue( "Composites", "" ), KeyValue( "PVs", "" ), KeyValue( "Hits", "" ),
+                           KeyValue( "DeVP", LHCb::Det::VP::det_path ),
+                           KeyValue( "Sensors", {"AlgorithmSpecific-" + name + "-sensorinfos"} )},
+                          {KeyValue{"OutputTracks", ""}, KeyValue{"OutputRelations", ""}} ) {}
 
   // initialization
   StatusCode LHCb::Pr::Velo::HeavyFlavourTrackFinder::initialize() {
-    auto sc = Transformer::initialize().andThen( [&] {
+    auto sc = MultiTransformer::initialize().andThen( [&] {
       addConditionDerivation<Sensors( DeVP const& )>( {LHCb::Det::VP::det_path}, inputLocation<Sensors>() );
     } );
     return sc;
   }
 
   // main execution
-  Output LHCb::Pr::Velo::HeavyFlavourTrackFinder::operator()( Composites const& composites, PVs const& pvs,
+  Output LHCb::Pr::Velo::HeavyFlavourTrackFinder::operator()( LHCb::Particle::Range const& composites, PVs const& pvs,
                                                               VeloHits const& velohits, DeVP const& velo,
                                                               Sensors const& sensors ) const {
     // output
-    Output output{};
-    output.reserve( composites.size() );
+    auto result               = Output{};
+    auto& [tracks, relations] = result;
+    tracks.reserve( composites.size() );
+    relations.reserve( composites.size() );
 
-    for ( auto const composite : composites ) {
+    for ( LHCb::Particle const* composite : composites ) {
       // set output
-      auto track = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>();
-      track.field<TrackTag::CompositeIndex>().set( composite->key() );
+      using namespace LHCb::Event::Enum::Track;
+      auto track = new Track( History::PrVeloHeavyFlavour, Type::Velo, PatRecStatus::PatRecIDs );
 
       // primary vertex
       using Sel::Utils::endVertexPos;
-      auto const& PV = Sel::getBestPV( Sel::Utils::deref_if_ptr( composite ), pvs );
-      auto const& pv = endVertexPos( PV );
-      track.field<TrackTag::PrimaryVertexIndex>().set( PV.key() );
+      auto const& pv_obj = Sel::getBestPV( *composite, pvs );
+      auto const& pv     = endVertexPos( pv_obj );
 
       // secondary vertex
       auto const& end    = endVertexPos( *composite );
@@ -225,39 +224,50 @@ namespace LHCb::Pr::Velo {
 
       // crossed sensors in PV-SV line
       auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes );
-      track.field<TrackTag::nSensors>().set( crossedSensors.size() );
+      track->addInfo( AdditionalInfo::VeloNSensors, crossedSensors.size() );
 
       // look for hits in relevant sensors (duplicated remoded and sorted in z and distance to search window)
       auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() );
+      auto nhits      = found_hits.size();
 
       // store sorted hits
-      auto saved_hits = track.field<TrackTag::Hits>();
-      auto nhits      = found_hits.size();
-      saved_hits.resize( nhits );
-      for ( long unsigned int i = 0; i < nhits; i++ ) {
-        auto found_hit = found_hits[i];
-        saved_hits[i].field<TrackTag::LHCbID>().set( LHCbID( found_hit.id ) );
-        saved_hits[i].field<TrackTag::HitPosition>().set( found_hit.position );
-      }
+      std::vector<LHCbID> found_ids;
+      found_ids.reserve( nhits );
+      for ( auto const found_hit : found_hits ) { found_ids.push_back( LHCbID( found_hit.id ) ); }
+      track->setLhcbIDs( found_ids );
 
       // determine mother particle direction based on closest hit
       // in first crossed sensor or using SV-PV if not available
       auto direction    = found_hits.size() ? found_hits.front().position - pv : fd;
-      auto origin_state = track.field<TrackTag::OriginState>();
-      origin_state.setPosition( pv );
-      origin_state.setDirection( direction.x() / direction.z(), direction.y() / direction.z() );
-      origin_state.setQOverP( 1. / composite->p() );
+      auto origin_state = LHCb::State();
+      origin_state.setState( pv.x().cast(), pv.y().cast(), pv.z().cast(), ( direction.x() / direction.z() ).cast(),
+                             ( direction.y() / direction.z() ).cast(), composite->charge() / composite->p() );
+      origin_state.setLocation( LHCb::State::Location::ClosestToBeam );
+      track->addToStates( origin_state );
+
+      // in case of found hit, add appropriate state
+      if ( found_hits.size() ) {
+        auto fh = found_hits.front();
+        auto fm_state =
+            LHCb::State( origin_state.stateVector(), fh.position.z().cast(), LHCb::State::Location::FirstMeasurement );
+        fm_state.setX( fh.position.x().cast() );
+        fm_state.setY( fh.position.y().cast() );
+        track->addToStates( fm_state );
+      }
 
-      // fitting
-      // TODO (using VeloKalmanHelpers)
-      track.field<TrackTag::Chi2>().set( -1.f );
+      // save track
+      tracks.insert( track );
+      relations.i_push( composite, track );
 
       // statistics
       m_nInSensor += crossedSensors.size();
       m_nHitCandidates += nhits;
     }
 
-    return output;
+    // mandatory sorting of relation table
+    relations.i_sort();
+
+    return result;
   }
 
 } // namespace LHCb::Pr::Velo
-- 
GitLab


From f11f4dda19f1d5789ee2de5c502e257ef8889242 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Thu, 27 Jul 2023 15:31:44 +0200
Subject: [PATCH 09/14] add tolerance to sensor lookup for search cylinder, not
 line

---
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 11 ++++-------
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   | 19 +++++++++++--------
 2 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index 711dbf61ead..fb878247a6c 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -34,10 +34,9 @@
 namespace {
   // in/outputs
   using Composite    = LHCb::Particle;
-  using Composites   = LHCb::Particle::Selection;
+  using Composites   = LHCb::Particle::Range;
   using PVs          = LHCb::Event::PV::PrimaryVertexContainer;
   using Track        = LHCb::Event::v1::Track;
-  using Tracks       = Track::Container;
   using P2TRelations = LHCb::Relation2D<LHCb::Particle, Track>;
   using MCParticles  = LHCb::MCParticles;
   using MCHits       = LHCb::MCHits;
@@ -236,11 +235,9 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     // look at long tracks and see if pions match to all (as in 'reconstructed')
     auto reconstructed = matchToTracks( mcps, tracklinks, mcparts, 0.7 );
     // link MCParticle to daughters (pions) of composite
-    std::optional<Composite const*> composite = std::nullopt;
-    if ( m_composites.exist() ) {
-      composite = matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 );
-    }
-    auto has_comp = composite.has_value();
+    auto composites = m_composites.get();
+    auto composite  = matchToComposites( mcps, composites, tracklinks, mcparts, 0.7 );
+    auto has_comp   = composite.has_value();
     // get heavy flavour track
     auto const*  relations = m_relations.getIfExists();
     Track const* track     = nullptr;
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index ce3f506d4ec..dd7dd0b21e3 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -95,15 +95,18 @@ namespace LHCb::Pr::Velo {
     // functions
     template <typename Position3, typename Vector3>
     SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& start,
-                                   Position3 const& end, Vector3 const& slopes ) {
+                                   Position3 const& end, Vector3 const& slopes, float margin ) {
       SensorInfos crossedSensors;
-      auto        pos = start;
+      auto        pos               = start;
+      auto        margin_normalized = margin / sqrt( slopes.mag2() - 1 );
       for ( auto sensor : sensors ) {
         if ( sensor.z < pos.Z() ) continue;
         if ( sensor.z > end.Z() ) break;
-        // if in range look if search window intersects with sensor
-        pos = pos + ( sensor.z - pos.Z() ) * slopes;
-        if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) )
+        // if in range look if search window intersects with sensor (up to some margin deeper into VELO)
+        pos                  = pos + ( sensor.z - pos.Z() ) * slopes;
+        auto pos_with_margin = pos + margin_normalized * slopes;
+        if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) |
+             velo.sensor( sensor.id ).isInsideSensor( pos_with_margin.X().cast(), pos_with_margin.Y().cast() ) )
           crossedSensors.push_back( {pos.Z().cast(), sensor.id} );
       }
       return crossedSensors;
@@ -222,11 +225,11 @@ namespace LHCb::Pr::Velo {
       auto const  fd     = end - pv;
       auto const  slopes = fd / fd.Z();
 
-      // crossed sensors in PV-SV line
-      auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes );
+      // crossed sensors in PV-SV line (or cylinder, with non-zero margin; should be same as for hits)
+      auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, m_max_distance.value() );
       track->addInfo( AdditionalInfo::VeloNSensors, crossedSensors.size() );
 
-      // look for hits in relevant sensors (duplicated remoded and sorted in z and distance to search window)
+      // look for hits in relevant sensors (duplicates removed and sorted in z and distance to search window)
       auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() );
       auto nhits      = found_hits.size();
 
-- 
GitLab


From 08c83c83f6dbb776dedb0e2cb78d8d0c7094600e Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Thu, 27 Jul 2023 17:37:42 +0200
Subject: [PATCH 10/14] adds Particle maker for containing heavy flavour tracks

---
 .../include/Functors/JIT_includes.h           |   1 -
 Phys/ParticleMaker/CMakeLists.txt             |   1 +
 .../ParticleWithHeavyFlavourTrackMaker.cpp    | 104 ++++++++++++++++++
 3 files changed, 105 insertions(+), 1 deletion(-)
 create mode 100644 Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp

diff --git a/Phys/FunctorCore/include/Functors/JIT_includes.h b/Phys/FunctorCore/include/Functors/JIT_includes.h
index b1510005c2d..1b05e840f21 100644
--- a/Phys/FunctorCore/include/Functors/JIT_includes.h
+++ b/Phys/FunctorCore/include/Functors/JIT_includes.h
@@ -65,7 +65,6 @@
 #include "SelKernel/VertexRelation.h"
 // PrKernel
 #include "PrKernel/PrSelection.h"
-#include "PrKernel/PrVeloHeavyFlavourUtils.h"
 // RelTables
 #include "Event/TableView.h"
 #include "Relations/Relation2D.h"
diff --git a/Phys/ParticleMaker/CMakeLists.txt b/Phys/ParticleMaker/CMakeLists.txt
index 67420524060..4d39be7209d 100644
--- a/Phys/ParticleMaker/CMakeLists.txt
+++ b/Phys/ParticleMaker/CMakeLists.txt
@@ -23,6 +23,7 @@ gaudi_add_module(ParticleMaker
         src/ParticleMassDTFMonitor.cpp
         src/ParticleMassMonitor.cpp
         src/ParticleWithBremMaker.cpp
+        src/ParticleWithHeavyFlavourTrackMaker.cpp
         src/V0FromDstMaker.cpp
         src/ParticlesEmptyProducer.cpp
     LINK
diff --git a/Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp b/Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp
new file mode 100644
index 00000000000..3307e55dbc8
--- /dev/null
+++ b/Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp
@@ -0,0 +1,104 @@
+/*****************************************************************************\
+* (c) Copyright 2000-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 "Event/Particle.h"
+#include "Event/ProtoParticle.h"
+#include "Event/Track.h"
+#include "Kernel/IParticlePropertySvc.h"
+#include "Kernel/ParticleProperty.h"
+#include "LHCbAlgs/Transformer.h"
+#include "Relations/Relation2D.h"
+
+namespace {
+  using Track         = LHCb::Event::v1::Track;
+  using RelationTable = LHCb::Relation2D<LHCb::Particle, Track>;
+  using Output        = std::tuple<LHCb::Particles, LHCb::ProtoParticles>;
+} // namespace
+
+class ParticleWithHeavyFlavourTrackMaker
+    : public LHCb::Algorithm::MultiTransformer<Output( LHCb::Particle::Range const&, RelationTable const& )> {
+public:
+  ParticleWithHeavyFlavourTrackMaker( const std::string&, ISvcLocator* );
+
+  StatusCode initialize() override;
+
+  Output operator()( LHCb::Particle::Range const&, RelationTable const& ) const override;
+
+private:
+  // properties
+  Gaudi::Property<std::string> m_pid{this, "ParticleID", "B+"};
+
+  // particle property information
+  LHCb::ParticleProperty const* m_partprop     = nullptr;
+  LHCb::ParticleProperty const* m_antipartprop = nullptr;
+
+  // services
+  ServiceHandle<LHCb::IParticlePropertySvc> m_particlePropertySvc{this, "ParticleProperty",
+                                                                  "LHCb::ParticlePropertySvc"};
+};
+
+DECLARE_COMPONENT( ParticleWithHeavyFlavourTrackMaker )
+
+ParticleWithHeavyFlavourTrackMaker::ParticleWithHeavyFlavourTrackMaker( const std::string& name,
+                                                                        ISvcLocator*       pSvcLocator )
+    : MultiTransformer( name, pSvcLocator, {KeyValue{"InputComposites", ""}, KeyValue{"Composite2TrackRelations", ""}},
+                        {KeyValue{"OutputParticles", ""}, KeyValue{"OutputProtoParticles", ""}} ) {}
+
+StatusCode ParticleWithHeavyFlavourTrackMaker::initialize() {
+  return MultiTransformer::initialize().andThen( [&] {
+    // Get particle properties from ParticlePropertySvc
+    m_partprop = m_particlePropertySvc->find( m_pid );
+    if ( !m_partprop ) {
+      throw GaudiException( "Could not find ParticleProperty for " + m_pid.value(),
+                            "ParticleWithHeavyFlavourTrackMaker::initialize", StatusCode::FAILURE );
+    }
+    m_antipartprop = m_partprop->antiParticle();
+    return StatusCode::SUCCESS;
+  } );
+}
+
+Output ParticleWithHeavyFlavourTrackMaker::operator()( LHCb::Particle::Range const& in_particles,
+                                                       RelationTable const&         relation_table ) const {
+  auto result                       = Output{};
+  auto& [out_particles, out_protos] = result;
+
+  for ( auto const* in_particle : in_particles ) {
+    // heavy flavour track
+    auto         track_range = relation_table.relations( in_particle );
+    Track const* track       = track_range.empty() ? NULL : track_range.front().to();
+
+    // create proto from heavy flavour track
+    auto* out_proto = new LHCb::ProtoParticle();
+    if ( track ) out_proto->setTrack( track );
+    out_protos.add( out_proto );
+
+    // create new particle containing proto/track and daughters (in_particle)
+    auto* out_particle = new LHCb::Particle(
+        ( out_proto->charge() == m_partprop->charge() ? m_partprop : m_antipartprop )->particleID() );
+    out_particle->setProto( out_proto );
+
+    // add information from daughter composite
+    out_particle->addToDaughters( in_particle );
+
+    out_particle->setEndVertex( in_particle->endVertex() );
+    out_particle->setReferencePoint( in_particle->referencePoint() );
+
+    out_particle->setMomentum( in_particle->momentum() );
+    out_particle->setMeasuredMass( in_particle->measuredMass() );
+
+    out_particle->setPosCovMatrix( in_particle->posCovMatrix() );
+    out_particle->setMomCovMatrix( in_particle->momCovMatrix() );
+    out_particle->setPosMomCovMatrix( in_particle->posMomCovMatrix() );
+
+    out_particles.add( out_particle );
+  }
+
+  return result;
+}
-- 
GitLab


From 33161ecf2583fe1e72e555183ee9a4c4fc38f58f Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Mon, 14 Aug 2023 16:57:42 +0200
Subject: [PATCH 11/14] added lastmeasurement state; added sensor and hit count
 functors

---
 Phys/FunctorCore/include/Functors/Composite.h | 33 +++++++---
 Phys/FunctorCore/include/Functors/TrackLike.h | 10 +++
 Phys/FunctorCore/python/Functors/__init__.py  | 66 +++++++++++++++----
 .../tests/src/InstantiateFunctors.cpp         |  4 +-
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp |  3 +-
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   | 31 ++++++---
 6 files changed, 115 insertions(+), 32 deletions(-)

diff --git a/Phys/FunctorCore/include/Functors/Composite.h b/Phys/FunctorCore/include/Functors/Composite.h
index 9a00ace18f7..8bd68583260 100755
--- a/Phys/FunctorCore/include/Functors/Composite.h
+++ b/Phys/FunctorCore/include/Functors/Composite.h
@@ -322,35 +322,48 @@ namespace Functors::Composite {
   };
 
   namespace BTracking {
+
+    template <typename Composite2TrackRelations, typename Composite>
+    auto getTrack( Composite2TrackRelations const& relations, Composite const& composite ) {
+      auto const& comp        = Sel::Utils::deref_if_ptr( composite );
+      auto        track_range = relations.relations( &comp );
+      return track_range.empty() ? std::nullopt : Functors::Optional{track_range.front().to()};
+    }
+
+    /** @brief Get heavy flavour track associated to composite. */
+    struct Track : public Function {
+      static constexpr auto name() { return "BTracking::Track"; }
+
+      template <typename Composite2TrackRelations, typename Composite>
+      auto operator()( Composite2TrackRelations const& relations, Composite const& composite ) const {
+        return getTrack( relations, composite );
+      }
+    };
+
     /** @brief Calculate the corrected mass using heavy flavour tracking info. */
     struct CorrectedMass : public Function {
       static constexpr auto name() { return "BTracking::CorrectedMass"; }
 
       template <typename Composite2TrackRelations, typename Composite>
       auto operator()( Composite2TrackRelations const& relations, Composite const& composite ) const {
-
         using FType = SIMDWrapper::scalar::float_v;
         using Sel::Utils::mass2;
         using Sel::Utils::threeMomentum;
-        using std::sqrt;
-        auto const& comp = Sel::Utils::deref_if_ptr( composite );
 
-        // first get the right heavy flavour track
-        auto track_range = relations.relations( &comp );
-        if ( track_range.empty() ) return FType{0};
-        auto track = track_range.front().to();
+        auto track = getTrack( relations, composite );
+        if ( !track.has_value() ) return FType{0};
 
         // direction using first hit (or end vertex if not available, as saved in heavy flavour track)
-        auto const state = track->firstState();
+        auto const state = track.value()->firstState();
         auto const d     = LHCb::LinAlg::Vec<FType, 3>{state.tx(), state.ty(), 1.};
 
         // Get the pT variable that we need
-        auto const mom  = threeMomentum( comp );
+        auto const mom  = threeMomentum( composite );
         auto const perp = mom - d * ( dot( mom, d ) / d.mag2() );
         auto const pt   = perp.mag();
 
         // Calculate the corrected mass
-        return pt + sqrt( mass2( comp ) + pt * pt );
+        return pt + sqrt( mass2( composite ) + pt * pt );
       }
     };
   } // namespace BTracking
diff --git a/Phys/FunctorCore/include/Functors/TrackLike.h b/Phys/FunctorCore/include/Functors/TrackLike.h
index e955aa16f0b..ea3c1d59a93 100644
--- a/Phys/FunctorCore/include/Functors/TrackLike.h
+++ b/Phys/FunctorCore/include/Functors/TrackLike.h
@@ -773,6 +773,16 @@ namespace Functors::Track {
     }
   };
 
+  /** @brief Number of expected Velo clusters from VELO 3D pattern recognition
+   */
+  struct nPRVelo3DExpect : public Function {
+    static constexpr auto name() { return "nPRVelo3DExpect"; }
+    template <typename Data>
+    auto operator()( Data const& d ) const {
+      return (int)Sel::Utils::deref_if_ptr( d ).info( LHCb::Event::Enum::Track::AdditionalInfo::nPRVelo3DExpect, -1. );
+    }
+  };
+
   /**
    * @brief MC_Reconstructed, return the reconstructed category
    * for MC Particle. The input of this functor must be reconstructed
diff --git a/Phys/FunctorCore/python/Functors/__init__.py b/Phys/FunctorCore/python/Functors/__init__.py
index 63e58f5cebc..ee5e0ec89a0 100644
--- a/Phys/FunctorCore/python/Functors/__init__.py
+++ b/Phys/FunctorCore/python/Functors/__init__.py
@@ -1031,6 +1031,13 @@ NFTHITS = Functor(
 
                    Functor's call operator expects a track like object.
                    """)
+NPRVELO3DEXPECT = Functor(
+    'NPRVELO3DEXPECT', "Track::nPRVelo3DExpect",
+    """Number of expected Velo clusters from VELO 3D pattern recognition.
+
+        Functor's call operator expects a track like object.
+        """)
+
 ###
 TRACKHISTORY = CAST_TO_INT @ Functor(
     'TRACKHISTORY', "Track::History",
@@ -1372,14 +1379,6 @@ MVA technique to establish pion mass hypothesis which keep in consideration the
 General information on PID could be found at the following `TWiki page <https://twiki.cern.ch/twiki/bin/view/LHCbPhysics/ChargedPID>`_
 """)
 
-_BPVCORRM = Functor(
-    '_BPVCORRM', 'Composite::CorrectedMass',
-    'Compute the corrected mass of the composite using the associated [primary] vertex.'
-)
-_BTRACKING_BPVCORRM = Functor(
-    '_BTRACKING_BPVCORRM', 'Composite::BTracking::CorrectedMass',
-    'Compute the corrected mass of the composite using the associated heavy flavour track.'
-)
 
 def BPVCORRM(Vertices: DataHandle):
     """ Compute the corrected mass of the composite using the :py:func:`~BPV`.
@@ -1413,16 +1412,61 @@ def BPVCORRMERR(Vertices: DataHandle):
     return _BPVCORRMERR.bind(BPV(Vertices), FORWARDARGS)
 
 
-def BTRACKING_BPVCORRM(HeavyFlavourTracks: DataHandle):
+def BTRACKING_TRACK(HeavyFlavourTrackRelations: DataHandle):
+    """ Returns heavy flavour track object associated to composite.
+
+    Functor's call operator expects a composite like object
+
+    Args:
+        HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track
+
+    """
+    _BTRACKING_TRACK = Functor('BTRACKING_TRACK',
+                               'Composite::BTracking::Track',
+                               'Heavy flavour track associated to composite.')
+
+    return _BTRACKING_TRACK.bind(TES(HeavyFlavourTrackRelations), FORWARDARGS)
+
+
+def BTRACKING_NHITS(HeavyFlavourTrackRelations: DataHandle):
+    """ Number of hits on heavy flavour track.
+
+    Functor's call operator expects a composite like object
+
+    Args:
+        HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track
+
+    """
+    return NHITS @ BTRACKING_TRACK(HeavyFlavourTrackRelations)
+
+
+def BTRACKING_NPRVELO3DEXPECT(HeavyFlavourTrackRelations: DataHandle):
+    """ Number of crossed VELO sensors according to pattern recognition window.
+
+    Functor's call operator expects a composite like object
+
+    Args:
+        HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track
+
+    """
+    return NPRVELO3DEXPECT @ BTRACKING_TRACK(HeavyFlavourTrackRelations)
+
+
+def BTRACKING_BPVCORRM(HeavyFlavourTrackRelations: DataHandle):
     """ Compute the corrected mass of the composite using the associated heavy flavour track.
 
     Functor's call operator expects a composite like object
 
     Args:
-        Vertices: DataHandle of the primary vertices
+        HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track
 
     """
-    return _BTRACKING_BPVCORRM.bind(TES(HeavyFlavourTracks), FORWARDARGS)
+    _BTRACKING_BPVCORRM = Functor(
+        'BTRACKING_BPVCORRM', 'Composite::BTracking::CorrectedMass',
+        'Compute the corrected mass of the composite using the associated heavy flavour track.'
+    )
+    return _BTRACKING_BPVCORRM.bind(
+        TES(HeavyFlavourTrackRelations), FORWARDARGS)
 
 
 VTX_FDCHI2 = Functor(
diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
index 6c08a94cff3..451e0856b2e 100644
--- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
+++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp
@@ -211,7 +211,9 @@ using cmerr5 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHC
 //
 // BTracking::CorrectedMass( heavyflavourtrackrelations, composite )
 //
-using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass,
+using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::Track,
+                                   LHCb::Relation2D<LHCb::Particle, LHCb::Event::v1::Track>, LHCb::Particle const&>;
+using cmbt2 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass,
                                    LHCb::Relation2D<LHCb::Particle, LHCb::Event::v1::Track>, LHCb::Particle const&>;
 
 //
diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index fb878247a6c..73fc4a1cf23 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -340,7 +340,8 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
     // track info
     std::string hfbase = base + "HeavyFlavourTracking_";
     sc &= tuple->column( hfbase + "Track_nSensors",
-                         has_track ? track->info( LHCb::Event::Enum::Track::AdditionalInfo::VeloNSensors, 0 ) : -1. );
+                         has_track ? track->info( LHCb::Event::Enum::Track::AdditionalInfo::nPRVelo3DExpect, -1. )
+                                   : -2. );
     sc &= tuple->column( hfbase + "Track_nHits", has_track ? track->nHits() : -1 );
     // track mc matching info
     bool fh_from_mcps = false;
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index dd7dd0b21e3..3a132041878 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -211,10 +211,6 @@ namespace LHCb::Pr::Velo {
     relations.reserve( composites.size() );
 
     for ( LHCb::Particle const* composite : composites ) {
-      // set output
-      using namespace LHCb::Event::Enum::Track;
-      auto track = new Track( History::PrVeloHeavyFlavour, Type::Velo, PatRecStatus::PatRecIDs );
-
       // primary vertex
       using Sel::Utils::endVertexPos;
       auto const& pv_obj = Sel::getBestPV( *composite, pvs );
@@ -227,7 +223,11 @@ namespace LHCb::Pr::Velo {
 
       // crossed sensors in PV-SV line (or cylinder, with non-zero margin; should be same as for hits)
       auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, m_max_distance.value() );
-      track->addInfo( AdditionalInfo::VeloNSensors, crossedSensors.size() );
+
+      // build track (output; also save if no sensors crossed, as search window info is still persisted this way)
+      using namespace LHCb::Event::Enum::Track;
+      auto track = new Track( History::PrVeloHeavyFlavour, Type::Velo, PatRecStatus::PatRecIDs );
+      track->addInfo( AdditionalInfo::nPRVelo3DExpect, crossedSensors.size() );
 
       // look for hits in relevant sensors (duplicates removed and sorted in z and distance to search window)
       auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() );
@@ -241,21 +241,34 @@ namespace LHCb::Pr::Velo {
 
       // determine mother particle direction based on closest hit
       // in first crossed sensor or using SV-PV if not available
-      auto direction    = found_hits.size() ? found_hits.front().position - pv : fd;
+      auto direction    = nhits > 0 ? found_hits.front().position - pv : fd;
+      auto qoverp       = composite->charge() / composite->p();
       auto origin_state = LHCb::State();
       origin_state.setState( pv.x().cast(), pv.y().cast(), pv.z().cast(), ( direction.x() / direction.z() ).cast(),
-                             ( direction.y() / direction.z() ).cast(), composite->charge() / composite->p() );
+                             ( direction.y() / direction.z() ).cast(), qoverp );
       origin_state.setLocation( LHCb::State::Location::ClosestToBeam );
       track->addToStates( origin_state );
 
-      // in case of found hit, add appropriate state
-      if ( found_hits.size() ) {
+      // in case of found hit(s), add appropriate states
+      if ( nhits > 0 ) {
         auto fh = found_hits.front();
         auto fm_state =
             LHCb::State( origin_state.stateVector(), fh.position.z().cast(), LHCb::State::Location::FirstMeasurement );
         fm_state.setX( fh.position.x().cast() );
         fm_state.setY( fh.position.y().cast() );
         track->addToStates( fm_state );
+
+        // in case more hits, add last measurement state
+        if ( nhits > 1 ) {
+          auto lh           = found_hits.back();
+          auto lh_direction = lh.position - pv;
+          auto last_state   = LHCb::State();
+          last_state.setState( lh.position.x().cast(), lh.position.y().cast(), lh.position.z().cast(),
+                               ( lh_direction.x() / lh_direction.z() ).cast(),
+                               ( lh_direction.y() / lh_direction.z() ).cast(), qoverp );
+          last_state.setLocation( LHCb::State::Location::LastMeasurement );
+          track->addToStates( last_state );
+        }
       }
 
       // save track
-- 
GitLab


From 46a62e7c8c0c217e08e0836518a345883a299a18 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Wed, 16 Aug 2023 17:23:43 +0200
Subject: [PATCH 12/14] adding efficiency counters

---
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 85 ++++++++++++++++++-
 1 file changed, 83 insertions(+), 2 deletions(-)

diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index 73fc4a1cf23..ca36d6d97d7 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -202,7 +202,22 @@ private:
   DataObjectReadHandle<P2TRelations> m_relations{this, "Composite2HeavyFlavourTrackRelations", ""};
 
   // monitoring
-  mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_hftrack_global_eff{this,
+                                                                      "Eff: has HF track | sel'd FS & HF reco'ible"};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_hftrack_pure_eff{this,
+                                                                    "HF track: MC-matched | sel'd FS & HF reco'ible"};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_reconstructible_hf_track{this, "HF track: reco'ible"};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_reconstructed_hf_track{this, "HF track: reco'd | reco'ible"};
+
+  mutable Gaudi::Accumulators::StatCounter<int>  m_hftrack_nhits{this, "HF track: #hits | sel'd FS & HF track"};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_hftrack_purity{this, "HF track: hit purity | sel'd FS & HF track"};
+
+  mutable Gaudi::Accumulators::BinomialCounter<> m_vphit_for_mchit_eff{this, "Eff: VP hit for VP MCHit"};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_hit_on_track_for_vphit_with_mchit_eff{
+      this, "Eff: true hit on track | sel'd FS & HF track"};
+
+  mutable Gaudi::Accumulators::BinomialCounter<> m_reconstructed_eff{this, "Final state (FS): reco'd"};
+  mutable Gaudi::Accumulators::BinomialCounter<> m_final_state_eff{this, "Final state (FS): sel'd | reco'd"};
 };
 
 DECLARE_COMPONENT_WITH_ID( PrVeloHeavyFlavourTrackingChecker, "PrVeloHeavyFlavourTrackingChecker" )
@@ -246,7 +261,73 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts,
       if ( trackrange.size() ) track = trackrange.front();
     }
     auto has_track = track ? true : false;
-    auto has_hits  = has_track ? track->nHits() > 0 : false;
+    int  n_hits    = has_track ? track->nHits() : 0;
+    auto has_hits  = n_hits > 0;
+
+    //// counters
+
+    // efficiencies of final state reconstruction and selection
+    m_reconstructed_eff += reconstructed;
+    if ( reconstructed ) m_final_state_eff += has_comp;
+
+    // check MCHits associated to B+/tau+ if they are reconstructible
+    // and if so look if they appear on the reconstructed track in case
+    // the final state is reconstructed, selected and matched and has associated HF track
+    bool has_reconstructible_hf_track = false;
+    bool has_reconstructed_hf_track   = false;
+    for ( auto const bhit : bhits ) {
+      auto ids_for_mchit       = channelIDForMCHit.find( bhit );
+      auto has_vphit_for_mchit = ids_for_mchit != channelIDForMCHit.end();
+      m_vphit_for_mchit_eff += has_vphit_for_mchit;
+      if ( !has_vphit_for_mchit ) continue;
+      has_reconstructible_hf_track = true;
+      auto ids                     = ids_for_mchit->second;
+      int  n_matches               = 0;
+      if ( has_hits ) {
+        auto track_lhcbids = track->lhcbIDs();
+        n_matches          = std::count_if( ids.begin(), ids.end(), [&track_lhcbids]( auto whit ) {
+          return std::count_if( track_lhcbids.begin(), track_lhcbids.end(),
+                                [&whit]( auto lhcbid ) { return lhcbid.vpID() == whit.id; } ) > 0;
+        } );
+      }
+      if ( reconstructed && has_comp && has_hits ) {
+        has_reconstructed_hf_track |= n_matches > 0;
+        m_hit_on_track_for_vphit_with_mchit_eff += n_matches > 0;
+      }
+    }
+
+    float track_purity = -1.f;
+    if ( reconstructed && has_comp && has_hits ) {
+      int n_matched_hits = 0;
+      for ( auto const lhcbid : track->lhcbIDs() ) {
+        unsigned int channelid = lhcbid.vpID();
+        int          n_match_per_hit =
+            bhits.empty() ? 0
+                          : std::count_if( bhits.begin(), bhits.end(), [&channelIDForMCHit, &channelid]( auto mchit ) {
+                              auto map_iter = channelIDForMCHit.find( mchit );
+                              int  n_match_per_hit_per_mchit =
+                                  ( map_iter == channelIDForMCHit.end() )
+                                      ? 0
+                                      : std::count_if( map_iter->second.begin(), map_iter->second.end(),
+                                                       [&channelid]( auto whit ) { return whit.id == channelid; } );
+                              return n_match_per_hit_per_mchit > 0;
+                            } );
+        if ( n_match_per_hit > 0 ) n_matched_hits++;
+        m_hftrack_purity += n_match_per_hit > 0;
+      }
+      track_purity = n_matched_hits / (float)n_hits;
+      m_hftrack_nhits += n_hits;
+    }
+
+    m_reconstructible_hf_track += has_reconstructible_hf_track;
+    if ( has_reconstructible_hf_track ) {
+      m_reconstructed_hf_track += has_reconstructed_hf_track;
+      if ( reconstructed && has_comp ) {
+        m_hftrack_global_eff += has_hits;
+        m_hftrack_pure_eff += has_hits && track_purity >= 0.7;
+      }
+    }
+
     //// fill ntuple with info
     std::string base = "MCP_";
     // mc truth
-- 
GitLab


From f815da3cec8a0f55a55d10ba35af7720a06227be Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Thu, 24 Aug 2023 17:08:22 +0200
Subject: [PATCH 13/14] added control of charge setting

---
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index 3a132041878..6f861b214ba 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -174,6 +174,7 @@ namespace LHCb::Pr::Velo {
     // properties
     Gaudi::Property<float> m_max_distance{this, "MaxDistanceFromPVSV", 0.5 * Gaudi::Units::mm,
                                           "maximum hit distance from PV-SV line (search window)"};
+    Gaudi::Property<bool>  m_revert_charge{this, "RevertCharge", false, "Revert charge with respect to composite"};
 
     // statistics monitoring
     mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "Nb of sensors crossed"};
@@ -242,7 +243,8 @@ namespace LHCb::Pr::Velo {
       // determine mother particle direction based on closest hit
       // in first crossed sensor or using SV-PV if not available
       auto direction    = nhits > 0 ? found_hits.front().position - pv : fd;
-      auto qoverp       = composite->charge() / composite->p();
+      auto charge       = composite->charge() != 0 ? ( m_revert_charge ? -1 : 1 ) * composite->charge() : 1;
+      auto qoverp       = charge / composite->p();
       auto origin_state = LHCb::State();
       origin_state.setState( pv.x().cast(), pv.y().cast(), pv.z().cast(), ( direction.x() / direction.z() ).cast(),
                              ( direction.y() / direction.z() ).cast(), qoverp );
-- 
GitLab


From 4170bbfc06736b09865d6cfd9bf27a8cc765e331 Mon Sep 17 00:00:00 2001
From: Maarten Van Veghel <maarten.vanveghel@cern.ch>
Date: Wed, 1 Nov 2023 18:47:29 +0100
Subject: [PATCH 14/14] sensor search with distance to sensor (full cylinder
 assumption)

---
 .../tests/options/test_functors_tested.py     |  9 ++--
 .../src/PrVeloHeavyFlavourTrackingChecker.cpp |  2 +-
 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp   | 47 +++++++------------
 3 files changed, 23 insertions(+), 35 deletions(-)

diff --git a/Phys/FunctorCore/tests/options/test_functors_tested.py b/Phys/FunctorCore/tests/options/test_functors_tested.py
index 5eedfb9d66e..8348dc2dbd5 100644
--- a/Phys/FunctorCore/tests/options/test_functors_tested.py
+++ b/Phys/FunctorCore/tests/options/test_functors_tested.py
@@ -65,10 +65,11 @@ exceptions = [
     'Functors::Track::Flag', 'Functors::Track::TX', 'Functors::Track::TY',
     'Functors::Track::history', 'Functors::Track::nDoF',
     'Functors::Track::nFTHits', 'Functors::Track::nUTHits',
-    'Functors::Track::nVPHits', 'Functors::Common::ForwardArgs<>( ',
-    'Functors::Common::ForwardArgs<0>( ', 'Functors::Common::ForwardArgs<1>( ',
-    'Functors::Common::ForwardArgs<2>( ', 'Functors::Common::UnitVector',
-    'Functors::Common::TES', 'Functors::Common::Call', 'Functors::Common::Dot',
+    'Functors::Track::nVPHits', 'Functors::Track::nPRVelo3DExpect',
+    'Functors::Common::ForwardArgs<>( ', 'Functors::Common::ForwardArgs<0>( ',
+    'Functors::Common::ForwardArgs<1>( ', 'Functors::Common::ForwardArgs<2>( ',
+    'Functors::Common::UnitVector', 'Functors::Common::TES',
+    'Functors::Common::Call', 'Functors::Common::Dot',
     'Functors::Common::NormedDot', 'Functors::Common::BestPV',
     'Functors::Common::ToLinAlg', 'Functors::Common::Phi_Coordinate',
     'Functors::Examples::PlusN',
diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
index ca36d6d97d7..f270550f811 100644
--- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
+++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp
@@ -116,7 +116,7 @@ namespace {
     size_t nmatches = 0;
     for ( LHCb::MCParticle const* pion : mcps.pions ) {
       map.applyToAllLinks(
-          [&pion, &mcps, &mcparts, &min_weight, &nmatches]( unsigned int, unsigned int mcPartKey, float weight ) {
+          [&pion, &mcparts, &min_weight, &nmatches]( unsigned int, unsigned int mcPartKey, float weight ) {
             if ( ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) &&
                  ( weight > min_weight ) )
               nmatches += 1;
diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
index 6f861b214ba..986db9a3fc2 100644
--- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
+++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp
@@ -48,12 +48,6 @@ namespace LHCb::Pr::Velo {
       friend bool operator<( SensorInfo const& lhs, SensorInfo const& rhs ) {
         return std::pair{lhs.z, lhs.id} < std::pair{rhs.z, rhs.id};
       }
-      // module information
-      auto module() const {
-        auto cid = ChannelID();
-        cid.setSensor( this->id );
-        return cid.module();
-      }
     };
     using SensorInfos = std::vector<SensorInfo>;
 
@@ -95,25 +89,24 @@ namespace LHCb::Pr::Velo {
     // functions
     template <typename Position3, typename Vector3>
     SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& start,
-                                   Position3 const& end, Vector3 const& slopes, float margin ) {
+                                   Position3 const& end, Vector3 const& slopes, float max_distance_squared ) {
       SensorInfos crossedSensors;
-      auto        pos               = start;
-      auto        margin_normalized = margin / sqrt( slopes.mag2() - 1 );
+      auto        pos = start;
       for ( auto sensor : sensors ) {
         if ( sensor.z < pos.Z() ) continue;
         if ( sensor.z > end.Z() ) break;
-        // if in range look if search window intersects with sensor (up to some margin deeper into VELO)
-        pos                  = pos + ( sensor.z - pos.Z() ) * slopes;
-        auto pos_with_margin = pos + margin_normalized * slopes;
-        if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) |
-             velo.sensor( sensor.id ).isInsideSensor( pos_with_margin.X().cast(), pos_with_margin.Y().cast() ) )
-          crossedSensors.push_back( {pos.Z().cast(), sensor.id} );
+        // if in range look if search window intersects with sensor
+        pos = pos + ( sensor.z - pos.Z() ) * slopes;
+        auto dist =
+            velo.sensor( sensor.id )
+                .closestDistanceToSensor( ROOT::Math::XYZPoint( pos.X().cast(), pos.Y().cast(), pos.Z().cast() ) );
+        if ( dist.perp2() < max_distance_squared ) crossedSensors.push_back( {pos.Z().cast(), sensor.id} );
       }
       return crossedSensors;
     }
 
     template <typename Position3, typename Vector3>
-    HitInfos getHitsForCrossedModules( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start,
+    HitInfos getHitsForCrossedSensors( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start,
                                        Vector3 const& slopes, FType max_distance ) {
       HitInfos poshits;
       if ( !sensors.size() ) return poshits;
@@ -122,25 +115,19 @@ namespace LHCb::Pr::Velo {
       auto pos = start;
       // check sensors
       for ( auto const sensor : sensors ) {
-        // its associated module (to not be affected by borders between sensors)
-        auto crossed_module = sensor.module();
         // go to sensor location
         pos = pos + ( sensor.z - pos.z().cast() ) * slopes;
-        // look if this sensor (and associated module) has hits and if they are in search window
+        // look if this sensor has hits and if they are in the search window
         for ( auto hit : hits.scalar() ) {
           auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() );
-          if ( id.module() == crossed_module ) {
+          if ( id.sensor() == sensor.id ) {
             auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>();
-            // make sure they are in the same plane up to some tolerance
-            // (not all sensors in a module are exactly)
-            if ( abs( hitpos.z().cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) {
-              auto d = ( hitpos - pos ).mag();
-              if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos, id} );
-            }
+            auto dist   = ( hitpos - pos ).rho();
+            if ( dist < max_distance ) poshits.push_back( {dist.cast(), hitpos, id} );
           }
         }
       }
-      // sort and remove duplicates (can arise with only module/plane (not sensor) checking)
+      // sort and remove duplicates (shouldn't be possible in this logic, but leave out of paranoia)
       std::sort( poshits.begin(), poshits.end() );
       poshits.erase( std::unique( poshits.begin(), poshits.end() ), poshits.end() );
       return poshits;
@@ -222,8 +209,8 @@ namespace LHCb::Pr::Velo {
       auto const  fd     = end - pv;
       auto const  slopes = fd / fd.Z();
 
-      // crossed sensors in PV-SV line (or cylinder, with non-zero margin; should be same as for hits)
-      auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, m_max_distance.value() );
+      // crossed sensors in PV-SV cylinder (with m_max_distance as radius)
+      auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, std::pow( m_max_distance, 2 ) );
 
       // build track (output; also save if no sensors crossed, as search window info is still persisted this way)
       using namespace LHCb::Event::Enum::Track;
@@ -231,7 +218,7 @@ namespace LHCb::Pr::Velo {
       track->addInfo( AdditionalInfo::nPRVelo3DExpect, crossedSensors.size() );
 
       // look for hits in relevant sensors (duplicates removed and sorted in z and distance to search window)
-      auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() );
+      auto found_hits = getHitsForCrossedSensors( crossedSensors, velohits, pv, slopes, m_max_distance.value() );
       auto nhits      = found_hits.size();
 
       // store sorted hits
-- 
GitLab