From 635a32a03ddbe94df6e54f18e2f2e99486768126 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Tue, 8 Apr 2025 10:14:32 +0200
Subject: [PATCH 01/10] Load magnetic field grid into texture memory

---
 CMakeLists.txt                                |  11 ++
 .../include/Dumpers/Identifiers.h             |   9 +-
 .../common/include/MagneticField.cuh          |  91 ++++++++++++
 .../non_event_data/include/Consumers.h        |  15 +-
 .../non_event_data/src/MagneticField.cpp      | 131 +++++++++++++++++-
 integration/non_event_data/src/Updater.cpp    |   3 +-
 main/src/RegisterConsumers.cpp                |  12 +-
 stream/sequence/include/Constants.cuh         |   6 +
 8 files changed, 269 insertions(+), 9 deletions(-)
 create mode 100644 device/event_model/common/include/MagneticField.cuh

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 415ccec0e40..f0738cd2dbb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -167,6 +167,17 @@ if(CALLGRIND_PROFILE)
   add_compile_definitions(CALLGRIND_PROFILE)
 endif()
 
+# Choice of memory backend for magnetic field
+option(MAGFIELD_USE_TEXTURE "Use cuda texture memory for magnetic field" ON)
+if(NOT TARGET_DEVICE STREQUAL "CUDA")
+  # if target is not cuda, fallback to regular memory
+  set(MAGFIELD_USE_TEXTURE OFF)
+endif()
+message(STATUS "Magfield use texture: " ${MAGFIELD_USE_TEXTURE})
+if(MAGFIELD_USE_TEXTURE)
+  add_compile_definitions(MAGFIELD_USE_TEXTURE)
+endif()
+
 set(WITH_Allen_PRIVATE_DEPENDENCIES TRUE)
 include(AllenDependencies)
 
diff --git a/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h b/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h
index 152ab3d6025..352bca5456d 100644
--- a/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h
+++ b/Dumpers/BinaryDumpers/include/Dumpers/Identifiers.h
@@ -66,7 +66,14 @@ namespace Allen {
       inline static std::string const id = "CrossingAngles";
     };
 
-    /** @class UTLookupTables
+    /** @class MagneticFieldPolarity
+     *  Identifier for the magnetic field non-event data for Allen
+     */
+    struct MagneticFieldPolarity : Identifier {
+      inline static std::string const id = "MagneticFieldPolarity";
+    };
+
+    /** @class MagneticField
      *  Identifier for the magnetic field non-event data for Allen
      */
     struct MagneticField : Identifier {
diff --git a/device/event_model/common/include/MagneticField.cuh b/device/event_model/common/include/MagneticField.cuh
new file mode 100644
index 00000000000..9fa1d592775
--- /dev/null
+++ b/device/event_model/common/include/MagneticField.cuh
@@ -0,0 +1,91 @@
+/*****************************************************************************\
+* (c) Copyright 2025 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+namespace MagneticField {
+  struct Magfield {
+#ifdef MAGFIELD_USE_TEXTURE
+    // https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#texture-and-surface-memory
+
+    __device__ float3 fieldVectorLinearInterpolation([[maybe_unused]] float3 pos) const
+    {
+#ifdef __CUDA_ARCH__
+      const float x = (pos.x - minX) * invDx;
+      const float y = (pos.y - minY) * invDy;
+      const float z = (pos.z - minZ) * invDz;
+      return {tex3D<float>(tex_Bx, x, y, z), tex3D<float>(tex_Bx, x, y, z), tex3D<float>(tex_Bx, x, y, z)};
+#else
+      return {0, 0, 0};
+#endif
+    }
+
+    cudaArray_t array_Bx;
+    cudaArray_t array_By;
+    cudaArray_t array_Bz;
+
+    cudaTextureObject_t tex_Bx;
+    cudaTextureObject_t tex_By;
+    cudaTextureObject_t tex_Bz;
+#else
+    __device__ float3 fieldVectorLinearInterpolation(float3 pos) const
+    {
+      const float x = (pos.x - minX) * invDx;
+      const float y = (pos.y - minY) * invDy;
+      const float z = (pos.z - minZ) * invDz;
+
+      const int ix = (int) x;
+      const int iy = (int) y;
+      const int iz = (int) z;
+
+      if (ix < 0 || iy < 0 || iz < 0 || ix >= Nx - 1 || iy >= Ny - 1 || iz >= Nz - 1) {
+        return {0.f, 0.f, 0.f};
+      }
+
+      const int i000 = Nx * (Ny * iz + iy) + ix;
+      const int i001 = Nx * (Ny * (iz + 1) + iy) + ix;
+      const int i010 = Nx * (Ny * iz + iy + 1) + ix;
+      const int i011 = Nx * (Ny * (iz + 1) + iy + 1) + ix;
+
+      const float h1x = x - ix;
+      const float h1y = y - iy;
+      const float h1z = z - iz;
+
+      const float h0x = 1.0f - h1x;
+      const float h0y = 1.0f - h1y;
+      const float h0z = 1.0f - h1z;
+
+      const float h00 = h0x * h0y;
+      const float h01 = h0x * h1y;
+      const float h10 = h1x * h0y;
+      const float h11 = h1x * h1y;
+
+      return {(h0z * (h00 * Bx[i000] + h10 * Bx[i000 + 1] + h01 * Bx[i010] + h11 * Bx[i010 + 1]) +
+               h1z * (h00 * Bx[i001] + h10 * Bx[i001 + 1] + h01 * Bx[i011] + h11 * Bx[i011 + 1])),
+
+              (h0z * (h00 * By[i000] + h10 * By[i000 + 1] + h01 * By[i010] + h11 * By[i010 + 1]) +
+               h1z * (h00 * By[i001] + h10 * By[i001 + 1] + h01 * By[i011] + h11 * By[i011 + 1])),
+
+              (h0z * (h00 * Bz[i000] + h10 * Bz[i000 + 1] + h01 * Bz[i010] + h11 * Bz[i010 + 1]) +
+               h1z * (h00 * Bz[i001] + h10 * Bz[i001 + 1] + h01 * Bz[i011] + h11 * Bz[i011 + 1]))};
+    }
+
+    float* Bx;
+    float* By;
+    float* Bz;
+
+    int Nx, Ny, Nz;
+#endif
+
+    float invDx, invDy, invDz;
+    float minX, minY, minZ;
+  };
+
+} // namespace MagneticField
\ No newline at end of file
diff --git a/integration/non_event_data/include/Consumers.h b/integration/non_event_data/include/Consumers.h
index 044789c4703..e723b28c283 100644
--- a/integration/non_event_data/include/Consumers.h
+++ b/integration/non_event_data/include/Consumers.h
@@ -124,9 +124,10 @@ namespace Consumers {
   private:
     std::reference_wrapper<Constants> m_constants;
   };
-  struct MagneticField final : public Allen::NonEventData::Consumer {
+  
+  struct MagneticFieldPolarity final : public Allen::NonEventData::Consumer {
   public:
-    MagneticField(gsl::span<float>&, std::vector<float>&);
+    MagneticFieldPolarity(gsl::span<float>&, std::vector<float>&);
 
     void consume(std::vector<char> const& data) override;
 
@@ -135,6 +136,16 @@ namespace Consumers {
     std::reference_wrapper<std::vector<float>> m_host_magnet_polarity;
   };
 
+  struct MagneticField final : public Allen::NonEventData::Consumer {
+  public:
+    MagneticField(Constants& constants);
+
+    void consume(std::vector<char> const& data) override;
+
+  private:
+    std::reference_wrapper<Constants> m_constants;
+  };
+
   struct MuonGeometry final : public Allen::NonEventData::Consumer {
   public:
     MuonGeometry(std::vector<char>& host_geometry_raw, char*& dev_geometry_raw, Muon::MuonGeometry*& muon_geometry);
diff --git a/integration/non_event_data/src/MagneticField.cpp b/integration/non_event_data/src/MagneticField.cpp
index 14be298aa85..d82ea8597af 100644
--- a/integration/non_event_data/src/MagneticField.cpp
+++ b/integration/non_event_data/src/MagneticField.cpp
@@ -13,20 +13,21 @@
 #include <BackendCommon.h>
 #include <Common.h>
 #include <Consumers.h>
+#include "MagneticField.cuh"
 
 namespace {
   using std::string;
   using std::to_string;
 } // namespace
 
-Consumers::MagneticField::MagneticField(
+Consumers::MagneticFieldPolarity::MagneticFieldPolarity(
   gsl::span<float>& dev_magnet_polarity,
   std::vector<float>& host_magnet_polarity) :
   m_dev_magnet_polarity {dev_magnet_polarity},
   m_host_magnet_polarity {host_magnet_polarity}
 {}
 
-void Consumers::MagneticField::consume(std::vector<char> const& data)
+void Consumers::MagneticFieldPolarity::consume(std::vector<char> const& data)
 {
   auto& dev_magnet_polarity = m_dev_magnet_polarity.get();
   auto& host_magnet_polarity = m_host_magnet_polarity.get();
@@ -45,3 +46,129 @@ void Consumers::MagneticField::consume(std::vector<char> const& data)
   Allen::memcpy(dev_magnet_polarity.data(), data.data(), data.size(), Allen::memcpyHostToDevice);
   Allen::memcpy(host_magnet_polarity.data(), data.data(), data.size(), Allen::memcpyHostToHost);
 }
+
+Consumers::MagneticField::MagneticField(Constants& constants) : m_constants {constants} {}
+
+void Consumers::MagneticField::consume(std::vector<char> const& data)
+{
+  // Parse header
+  const float* invDxyz = reinterpret_cast<const float*>(data.data());
+  const int* Nxyz = reinterpret_cast<const int*>(data.data() + sizeof(float) * 4);
+  const float* min = reinterpret_cast<const float*>(data.data() + sizeof(float) * 8);
+  const float* B = reinterpret_cast<const float*>(data.data() + sizeof(float) * 12);
+
+  printf("Loading magfield grid, nbins x,y,z : (%d, %d, %d)\n", Nxyz[0], Nxyz[1], Nxyz[2]);
+  printf(
+    "dx, xmin, xmax: (%f, %f, %f)\n",
+    (double) (1.f / invDxyz[0]),
+    (double) min[0],
+    (double) (min[0] + (Nxyz[0] - 1) / invDxyz[0]));
+  printf(
+    "dy, ymin, ymax: (%f, %f, %f)\n",
+    (double) (1.f / invDxyz[1]),
+    (double) min[1],
+    (double) (min[1] + (Nxyz[1] - 1) / invDxyz[1]));
+  printf(
+    "dz, zmin, zmax: (%f, %f, %f)\n",
+    (double) (1.f / invDxyz[2]),
+    (double) min[2],
+    (double) (min[2] + (Nxyz[2] - 1) / invDxyz[2]));
+
+  // Reorder B field
+  std::size_t N = Nxyz[0] * Nxyz[1] * Nxyz[2];
+  std::vector<float> Bx, By, Bz;
+  Bx.reserve(N);
+  By.reserve(N);
+  Bz.reserve(N);
+  for (std::size_t i = 0; i < N * 4; i += 4) {
+    Bx.emplace_back(B[i]);
+    By.emplace_back(B[i + 1]);
+    Bz.emplace_back(B[i + 2]);
+  }
+
+  auto& magnetic_field = m_constants.get().magnetic_field;
+  if (magnetic_field == nullptr) {
+    magnetic_field = new ::MagneticField::Magfield();
+    magnetic_field->invDx = invDxyz[0];
+    magnetic_field->invDy = invDxyz[1];
+    magnetic_field->invDz = invDxyz[2];
+    magnetic_field->minX = min[0];
+    magnetic_field->minY = min[1];
+    magnetic_field->minZ = min[2];
+
+#ifdef MAGFIELD_USE_TEXTURE
+    // Allocate CUDA arrays in device memory
+    cudaChannelFormatDesc chDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
+    cudaMallocArray(&magnetic_field->array_Bx, &chDesc, Nxyz[0], Nxyz[1], Nxyz[2]);
+    cudaMallocArray(&magnetic_field->array_By, &chDesc, Nxyz[0], Nxyz[1], Nxyz[2]);
+    cudaMallocArray(&magnetic_field->array_Bz, &chDesc, Nxyz[0], Nxyz[1], Nxyz[2]);
+
+    // Specify texture object parameters
+    cudaTextureDesc texDesc;
+    memset(&texDesc, 0, sizeof(texDesc));
+    texDesc.addressMode[0] = cudaAddressModeClamp;
+    texDesc.addressMode[1] = cudaAddressModeClamp;
+    texDesc.addressMode[2] = cudaAddressModeClamp;
+    texDesc.filterMode = cudaFilterModeLinear;
+    texDesc.readMode = cudaReadModeElementType;
+    texDesc.normalizedCoords = 0;
+
+    // Specify textures
+    cudaResourceDesc resDesc;
+    std::memset(&resDesc, 0, sizeof(resDesc));
+    resDesc.resType = cudaResourceTypeArray;
+
+    resDesc.res.array.array = magnetic_field->array_Bx;
+    magnetic_field->tex_Bx = 0;
+    cudaCreateTextureObject(&magnetic_field->tex_Bx, &resDesc, &texDesc, NULL);
+
+    resDesc.res.array.array = magnetic_field->array_By;
+    magnetic_field->tex_By = 0;
+    cudaCreateTextureObject(&magnetic_field->tex_By, &resDesc, &texDesc, NULL);
+
+    resDesc.res.array.array = magnetic_field->array_Bz;
+    magnetic_field->tex_Bz = 0;
+    cudaCreateTextureObject(&magnetic_field->tex_Bz, &resDesc, &texDesc, NULL);
+#else
+    magnetic_field->Nx = Nxyz[0];
+    magnetic_field->Ny = Nxyz[1];
+    magnetic_field->Nz = Nxyz[2];
+
+    Allen::malloc((void**) &magnetic_field->Bx, sizeof(float) * N);
+    Allen::malloc((void**) &magnetic_field->By, sizeof(float) * N);
+    Allen::malloc((void**) &magnetic_field->Bz, sizeof(float) * N);
+#endif
+  }
+
+  // FIXME need to check the size of data is as expected
+
+#ifdef MAGFIELD_USE_TEXTURE
+
+  cudaMemcpy3DParms cpyParms;
+  memset(&cpyParms, 0, sizeof(cpyParms));
+  cpyParms.srcPtr.pitch = sizeof(float);
+  cpyParms.srcPtr.xsize = Nxyz[0];
+  cpyParms.srcPtr.ysize = Nxyz[1];
+  cpyParms.kind = cudaMemcpyHostToDevice;
+  cpyParms.extent.width = Nxyz[0];
+  cpyParms.extent.height = Nxyz[1];
+  cpyParms.extent.depth = Nxyz[2];
+
+  cpyParms.srcPtr.ptr = Bx.data();
+  cpyParms.dstArray = magnetic_field->array_Bx;
+  cudaMemcpy3D(&cpyParms);
+
+  cpyParms.srcPtr.ptr = By.data();
+  cpyParms.dstArray = magnetic_field->array_By;
+  cudaMemcpy3D(&cpyParms);
+
+  cpyParms.srcPtr.ptr = Bz.data();
+  cpyParms.dstArray = magnetic_field->array_Bz;
+  cudaMemcpy3D(&cpyParms);
+
+#else
+  Allen::memcpy(magnetic_field->Bx, Bx.data(), sizeof(float) * N, Allen::memcpyHostToDevice);
+  Allen::memcpy(magnetic_field->By, By.data(), sizeof(float) * N, Allen::memcpyHostToDevice);
+  Allen::memcpy(magnetic_field->Bz, Bz.data(), sizeof(float) * N, Allen::memcpyHostToDevice);
+#endif
+}
diff --git a/integration/non_event_data/src/Updater.cpp b/integration/non_event_data/src/Updater.cpp
index 55c9a9907ed..65be48b3e8d 100644
--- a/integration/non_event_data/src/Updater.cpp
+++ b/integration/non_event_data/src/Updater.cpp
@@ -56,7 +56,8 @@ namespace Allen {
       tuple producers {tuple {NonEventData::VeloGeometry {}, std::string("velo_geometry.bin")},
                        tuple {NonEventData::UTBoards {}, std::string("ut_boards.bin")},
                        tuple {NonEventData::Beamline {}, std::string("beamline.bin")},
-                       tuple {NonEventData::MagneticField {}, std::string("polarity.bin")},
+                       tuple {NonEventData::MagneticFieldPolarity {}, std::string("polarity.bin")},
+                       tuple {NonEventData::MagneticField {}, std::string("magfield.bin")},
                        tuple {NonEventData::UTGeometry {}, std::string("ut_geometry.bin")},
                        tuple {NonEventData::UTLookupTables {}, std::string("ut_tables.bin")},
                        tuple {NonEventData::SciFiGeometry {}, std::string("scifi_geometry.bin")},
diff --git a/main/src/RegisterConsumers.cpp b/main/src/RegisterConsumers.cpp
index 4d320bc3722..4eeee89dfd6 100644
--- a/main/src/RegisterConsumers.cpp
+++ b/main/src/RegisterConsumers.cpp
@@ -104,9 +104,15 @@ void register_consumers(
       },
       BankTypes::Rich2));
 
-  const auto unconditional_consumers =
-    std::make_tuple(std::make_tuple(Allen::NonEventData::MagneticField {}, [&constants]() {
-      return std::make_unique<Consumers::MagneticField>(constants.dev_magnet_polarity, constants.host_magnet_polarity);
+  const auto unconditional_consumers = std::make_tuple(
+    std::make_tuple(
+      Allen::NonEventData::MagneticFieldPolarity {},
+      [&constants]() {
+        return std::make_unique<Consumers::MagneticFieldPolarity>(
+          constants.dev_magnet_polarity, constants.host_magnet_polarity);
+      }),
+    std::make_tuple(Allen::NonEventData::MagneticField {}, [&constants]() {
+      return std::make_unique<Consumers::MagneticField>(constants);
     }));
 
   for_each(consumers, [updater, requested_banks](const auto& c) {
diff --git a/stream/sequence/include/Constants.cuh b/stream/sequence/include/Constants.cuh
index e1b3518a849..9abbc552dbe 100644
--- a/stream/sequence/include/Constants.cuh
+++ b/stream/sequence/include/Constants.cuh
@@ -51,6 +51,9 @@ namespace Rich::Future::DAQ::Allen {
 namespace UT::Constants {
   struct UTLayerGeometry;
 }
+namespace MagneticField {
+  struct Magfield;
+} // namespace MagneticField
 
 /**
  * @brief Struct intended as a singleton with constants defined on GPU.
@@ -100,6 +103,9 @@ struct Constants {
   gsl::span<float> dev_magnet_polarity;
   std::vector<float> host_magnet_polarity;
 
+  // Magnetic field
+  MagneticField::Magfield* magnetic_field = nullptr;
+
   // Looking forward
   LookingForward::Constants* host_looking_forward_constants;
 
-- 
GitLab


From c6a3c61e7fcf1aa7046f59a6c524acc5dbd63a42 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Tue, 8 Apr 2025 16:12:52 +0200
Subject: [PATCH 02/10] Add parabolic extrapolator example

---
 configuration/python/AllenConf/HLT1.py        |  8 +++
 .../python/AllenConf/hlt1_reconstruction.py   | 10 ++-
 .../secondary_vertex_reconstruction.py        | 11 ++-
 .../common/include/MagneticField.cuh          |  4 +-
 device/event_model/common/include/States.cuh  |  2 +
 .../ParKalman/include/ExtrapolateStates.cuh   | 38 ++++++++++
 .../include/ParabolicExtrapolator.cuh         | 57 +++++++++++++++
 .../kalman/ParKalman/src/ExtrapolateStates.cu | 69 +++++++++++++++++++
 .../non_event_data/src/MagneticField.cpp      | 24 ++++---
 9 files changed, 209 insertions(+), 14 deletions(-)
 create mode 100644 device/kalman/ParKalman/include/ExtrapolateStates.cuh
 create mode 100644 device/kalman/ParKalman/include/ParabolicExtrapolator.cuh
 create mode 100644 device/kalman/ParKalman/src/ExtrapolateStates.cu

diff --git a/configuration/python/AllenConf/HLT1.py b/configuration/python/AllenConf/HLT1.py
index 64a2736507a..da99bb613f4 100644
--- a/configuration/python/AllenConf/HLT1.py
+++ b/configuration/python/AllenConf/HLT1.py
@@ -1569,6 +1569,14 @@ def setup_hlt1_node(enablePhysics=True,
             NodeLogic.NONLAZY_AND,
             force_order=False)
 
+    if with_fullKF:
+        # Added for testing
+        hlt1_node = CompositeNode(
+            "AllenExtrapolatedStates",
+            [hlt1_node, reconstructed_objects["extrapolated_states"]],
+            NodeLogic.NONLAZY_AND,
+            force_order=True)
+
     if with_rich:
         hlt1_node = CompositeNode(
             "AllenWithRich", [
diff --git a/configuration/python/AllenConf/hlt1_reconstruction.py b/configuration/python/AllenConf/hlt1_reconstruction.py
index 306df863bb5..e79d39e8110 100644
--- a/configuration/python/AllenConf/hlt1_reconstruction.py
+++ b/configuration/python/AllenConf/hlt1_reconstruction.py
@@ -19,7 +19,7 @@ from AllenConf.primary_vertex_reconstruction import make_pvs
 from AllenConf.secondary_vertex_reconstruction import (
     make_kalman_velo_only, make_basic_particles, fit_secondary_vertices,
     make_sv_track_pairs, make_sv_pairs, make_generic_sv_pairs,
-    make_three_body_svs, make_kalman_long)
+    make_three_body_svs, make_kalman_long, make_extrapolated_states)
 from AllenConf.jet_reconstruction import make_cone_jets
 from AllenConf.validators import (
     velo_validation, veloUT_validation, seeding_validation,
@@ -159,6 +159,14 @@ def hlt1_reconstruction(algorithm_name='',
         kalman_long_tracks = make_kalman_long(long_tracks, pvs, muonID)
         KF_long_track = kalman_long_tracks
         output.update({"kalman_long_track": kalman_long_tracks})
+
+        # Added for testing magfield extrapolation:
+        output.update({
+            "extrapolated_states":
+            make_extrapolated_states(
+                KF_long_track["dev_kalman_states_view"],
+                long_tracks["host_number_of_reconstructed_scifi_tracks"])
+        })
     else:
         kalman_velo_only = make_kalman_velo_only(long_tracks, pvs, muonID)
         output.update({"kalman_velo_only": kalman_velo_only})
diff --git a/configuration/python/AllenConf/secondary_vertex_reconstruction.py b/configuration/python/AllenConf/secondary_vertex_reconstruction.py
index 677f5d2a59a..bad112d0b87 100644
--- a/configuration/python/AllenConf/secondary_vertex_reconstruction.py
+++ b/configuration/python/AllenConf/secondary_vertex_reconstruction.py
@@ -13,7 +13,7 @@ from AllenCore.algorithms import (
     make_long_track_particles_t, filter_tracks_t, fit_secondary_vertices_t,
     empty_lepton_id_t, sv_combiner_t, filter_svs_t, filter_two_svs_t,
     generic_sv_combiner_t, calc_max_combos_t, filter_sv_track_t,
-    kalman_filter_t, combine_sv_track_t, flatten_svs_t)
+    kalman_filter_t, combine_sv_track_t, flatten_svs_t, extrapolate_states_t)
 from AllenConf.utils import initialize_number_of_events, mep_layout
 from AllenConf.velo_reconstruction import run_velo_kalman_filter
 from AllenCore.generator import make_algorithm
@@ -48,6 +48,15 @@ def make_kalman_long(long_tracks, pvs, is_muon_result,
         "dev_kalman_states_view": kalman.dev_kalman_states_view_t
     }
 
+def make_extrapolated_states(states, n_states, step_dz=100.0, n_steps=100):
+    extrapolator = make_algorithm(
+        extrapolate_states_t,
+        name='extrapolate_states_{hash}',
+        host_number_of_input_states_t=n_states,
+        dev_kalman_states_view_t=states,
+        n_steps=n_steps,
+        step_dz=step_dz)
+    return extrapolator.dev_states_t
 
 def make_kalman_velo_only(long_tracks,
                           pvs,
diff --git a/device/event_model/common/include/MagneticField.cuh b/device/event_model/common/include/MagneticField.cuh
index 9fa1d592775..231d50af94f 100644
--- a/device/event_model/common/include/MagneticField.cuh
+++ b/device/event_model/common/include/MagneticField.cuh
@@ -21,9 +21,9 @@ namespace MagneticField {
       const float x = (pos.x - minX) * invDx;
       const float y = (pos.y - minY) * invDy;
       const float z = (pos.z - minZ) * invDz;
-      return {tex3D<float>(tex_Bx, x, y, z), tex3D<float>(tex_Bx, x, y, z), tex3D<float>(tex_Bx, x, y, z)};
+      return {tex3D<float>(tex_Bx, x, y, z), tex3D<float>(tex_By, x, y, z), tex3D<float>(tex_Bz, x, y, z)};
 #else
-      return {0, 0, 0};
+      return {1, 2, 3};
 #endif
     }
 
diff --git a/device/event_model/common/include/States.cuh b/device/event_model/common/include/States.cuh
index 43c7e46e29a..55c75466434 100644
--- a/device/event_model/common/include/States.cuh
+++ b/device/event_model/common/include/States.cuh
@@ -261,6 +261,8 @@ namespace Allen {
           assert(track_index < m_size);
           return KalmanState {m_base_pointer, m_offset + track_index, m_total_number_of_tracks};
         }
+
+        __host__ __device__ unsigned total_number_of_states() const { return m_total_number_of_tracks; }
       };
 
       struct SecondaryVertex {
diff --git a/device/kalman/ParKalman/include/ExtrapolateStates.cuh b/device/kalman/ParKalman/include/ExtrapolateStates.cuh
new file mode 100644
index 00000000000..6d42efec7b9
--- /dev/null
+++ b/device/kalman/ParKalman/include/ExtrapolateStates.cuh
@@ -0,0 +1,38 @@
+/*****************************************************************************\
+* (c) Copyright 2018-2020 CERN for the benefit of the LHCb Collaboration      *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+#include "AlgorithmTypes.cuh"
+#include "States.cuh"
+#include "ParabolicExtrapolator.cuh"
+
+namespace extrapolate_states {
+  struct Parameters {
+    HOST_INPUT(host_number_of_input_states_t, unsigned) host_number_of_input_states;
+    MASK_INPUT(dev_event_list_t) dev_event_list;
+    DEVICE_INPUT(dev_kalman_states_view_t, Allen::Views::Physics::KalmanStates) dev_kalman_states_view;
+    DEVICE_OUTPUT(dev_states_t, Extrapolators::State) dev_states;
+  };
+
+  struct extrapolate_states_t : public DeviceAlgorithm, Parameters {
+    void set_arguments_size(ArgumentReferences<Parameters>, const RuntimeOptions&, const Constants&) const;
+
+    void operator()(
+      const ArgumentReferences<Parameters>&,
+      const RuntimeOptions&,
+      const Constants&,
+      const Allen::Context&) const;
+
+  private:
+    Allen::Property<float> m_step_dz {this, "step_dz", 100.0, "size of a step in z (in mm)"};
+    Allen::Property<unsigned> m_n_steps {this, "n_steps", 100, "number of steps"};
+  };
+} // namespace extrapolate_states
\ No newline at end of file
diff --git a/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh b/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh
new file mode 100644
index 00000000000..6a01ce7f904
--- /dev/null
+++ b/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh
@@ -0,0 +1,57 @@
+/*****************************************************************************\
+* (c) Copyright 2025 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+#include <MagneticField.cuh>
+
+namespace Extrapolators {
+  struct State {
+    float x, y, z, tx, ty, qop;
+  };
+
+  template<typename ftype = float>
+  struct ParabolicExtrapolator {
+    // Implementation taken from:
+    // https://gitlab.cern.ch/lhcb/Rec/-/blob/master/Tr/TrackExtrapolators/src/TrackParabolicExtrapolator.cpp
+
+    __device__ static void propagate(State& state, ftype dz, const MagneticField::Magfield& field)
+    {
+      // Get the B field at midpoint
+      const ftype xMid = state.x + state.tx * dz * ftype {0.5};
+      const ftype yMid = state.y + state.ty * dz * ftype {0.5};
+      const float3 P = make_float3(xMid, yMid, state.z + dz * ftype {0.5});
+      const float3 B = field.fieldVectorLinearInterpolation(P);
+
+      const ftype tx = state.tx;
+      const ftype ty = state.ty;
+      const ftype norm = fsqrt(ftype {1.} + tx * tx + ty * ty);
+
+      const ftype ax = norm * (ty * (tx * B.x + B.z) - ((ftype {1.} + tx * tx) * B.y));
+      const ftype ay = norm * (-tx * (ty * B.y + B.z) + ((ftype {1.} + ty * ty) * B.x));
+
+      const ftype fact = eplus * c_light * dz * state.qop;
+
+      // Update the state parameters (exact extrapolation)
+      state.x += dz * (tx + ftype {0.5} * ax * fact);
+      state.y += dz * (ty + ftype {0.5} * ay * fact);
+      state.z += dz;
+      state.tx += ax * fact;
+      state.ty += ay * fact;
+    }
+
+  private:
+    static constexpr ftype c_light = 2.99792458e+8 * 1000. / 1.e+9;
+    static constexpr ftype eplus = 1.;
+
+    __device__ static float fsqrt(float x) { return sqrtf(x); }
+    __device__ static double fsqrt(double x) { return sqrt(x); }
+  };
+} // namespace Extrapolators
\ No newline at end of file
diff --git a/device/kalman/ParKalman/src/ExtrapolateStates.cu b/device/kalman/ParKalman/src/ExtrapolateStates.cu
new file mode 100644
index 00000000000..df99fae27f6
--- /dev/null
+++ b/device/kalman/ParKalman/src/ExtrapolateStates.cu
@@ -0,0 +1,69 @@
+/*****************************************************************************\
+* (c) Copyright 2025 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), 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 "ExtrapolateStates.cuh"
+
+INSTANTIATE_ALGORITHM(extrapolate_states::extrapolate_states_t)
+
+void extrapolate_states::extrapolate_states_t::set_arguments_size(
+  ArgumentReferences<Parameters> arguments,
+  const RuntimeOptions&,
+  const Constants&) const
+{
+  auto n_states = first<host_number_of_input_states_t>(arguments);
+  set_size<dev_states_t>(arguments, n_states);
+}
+
+__global__ void extrapolate_states_kernel(
+  extrapolate_states::Parameters parameters,
+  float dz,
+  unsigned n_steps,
+  const MagneticField::Magfield field)
+{
+
+  const auto* input_states = parameters.dev_kalman_states_view.data();
+  const unsigned n_states = input_states->total_number_of_states();
+
+  for (unsigned i = blockIdx.x * blockDim.x + threadIdx.x; i < n_states; i += blockDim.x * gridDim.x) {
+    const auto& input = input_states->state(i);
+
+    Extrapolators::State output;
+    output.x = input.x();
+    output.y = input.y();
+    output.z = input.z();
+    output.tx = input.tx();
+    output.ty = input.ty();
+    output.qop = input.qop();
+
+    for (unsigned s = 0; s < n_steps; s++) {
+      Extrapolators::ParabolicExtrapolator<float>::propagate(output, dz, field);
+
+      /*if (threadIdx.x == 0 && blockIdx.x == 0) {
+        float3 B = field.fieldVectorLinearInterpolation(make_float3(output.x, output.y, output.z));
+        printf("[%g, %g, %g, %g, %g, %g], // Bx=%g By=%g Bz=%g\n",
+        output.x, output.y, output.z, output.tx, output.ty, output.qop, B.x, B.y, B.z);
+      }*/
+    }
+
+    parameters.dev_states[i] = output;
+  }
+}
+
+void extrapolate_states::extrapolate_states_t::operator()(
+  const ArgumentReferences<Parameters>& arguments,
+  const RuntimeOptions&,
+  const Constants& constants,
+  const Allen::Context& context) const
+{
+  const int block_dim = 256;
+  const int grid_dim = (first<host_number_of_input_states_t>(arguments) + block_dim - 1) / block_dim;
+  global_function(extrapolate_states_kernel)(grid_dim, block_dim, context)(
+    arguments, m_step_dz, m_n_steps, *constants.magnetic_field);
+}
diff --git a/integration/non_event_data/src/MagneticField.cpp b/integration/non_event_data/src/MagneticField.cpp
index d82ea8597af..0b0873a6ef0 100644
--- a/integration/non_event_data/src/MagneticField.cpp
+++ b/integration/non_event_data/src/MagneticField.cpp
@@ -99,9 +99,13 @@ void Consumers::MagneticField::consume(std::vector<char> const& data)
 #ifdef MAGFIELD_USE_TEXTURE
     // Allocate CUDA arrays in device memory
     cudaChannelFormatDesc chDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
-    cudaMallocArray(&magnetic_field->array_Bx, &chDesc, Nxyz[0], Nxyz[1], Nxyz[2]);
-    cudaMallocArray(&magnetic_field->array_By, &chDesc, Nxyz[0], Nxyz[1], Nxyz[2]);
-    cudaMallocArray(&magnetic_field->array_Bz, &chDesc, Nxyz[0], Nxyz[1], Nxyz[2]);
+    cudaExtent extent;
+    extent.width = Nxyz[0];
+    extent.height = Nxyz[1];
+    extent.depth = Nxyz[2];
+    cudaCheck(cudaMalloc3DArray(&magnetic_field->array_Bx, &chDesc, extent));
+    cudaCheck(cudaMalloc3DArray(&magnetic_field->array_By, &chDesc, extent));
+    cudaCheck(cudaMalloc3DArray(&magnetic_field->array_Bz, &chDesc, extent));
 
     // Specify texture object parameters
     cudaTextureDesc texDesc;
@@ -120,15 +124,15 @@ void Consumers::MagneticField::consume(std::vector<char> const& data)
 
     resDesc.res.array.array = magnetic_field->array_Bx;
     magnetic_field->tex_Bx = 0;
-    cudaCreateTextureObject(&magnetic_field->tex_Bx, &resDesc, &texDesc, NULL);
+    cudaCheck(cudaCreateTextureObject(&magnetic_field->tex_Bx, &resDesc, &texDesc, NULL));
 
     resDesc.res.array.array = magnetic_field->array_By;
     magnetic_field->tex_By = 0;
-    cudaCreateTextureObject(&magnetic_field->tex_By, &resDesc, &texDesc, NULL);
+    cudaCheck(cudaCreateTextureObject(&magnetic_field->tex_By, &resDesc, &texDesc, NULL));
 
     resDesc.res.array.array = magnetic_field->array_Bz;
     magnetic_field->tex_Bz = 0;
-    cudaCreateTextureObject(&magnetic_field->tex_Bz, &resDesc, &texDesc, NULL);
+    cudaCheck(cudaCreateTextureObject(&magnetic_field->tex_Bz, &resDesc, &texDesc, NULL));
 #else
     magnetic_field->Nx = Nxyz[0];
     magnetic_field->Ny = Nxyz[1];
@@ -146,7 +150,7 @@ void Consumers::MagneticField::consume(std::vector<char> const& data)
 
   cudaMemcpy3DParms cpyParms;
   memset(&cpyParms, 0, sizeof(cpyParms));
-  cpyParms.srcPtr.pitch = sizeof(float);
+  cpyParms.srcPtr.pitch = Nxyz[0] * sizeof(float);
   cpyParms.srcPtr.xsize = Nxyz[0];
   cpyParms.srcPtr.ysize = Nxyz[1];
   cpyParms.kind = cudaMemcpyHostToDevice;
@@ -156,15 +160,15 @@ void Consumers::MagneticField::consume(std::vector<char> const& data)
 
   cpyParms.srcPtr.ptr = Bx.data();
   cpyParms.dstArray = magnetic_field->array_Bx;
-  cudaMemcpy3D(&cpyParms);
+  cudaCheck(cudaMemcpy3D(&cpyParms));
 
   cpyParms.srcPtr.ptr = By.data();
   cpyParms.dstArray = magnetic_field->array_By;
-  cudaMemcpy3D(&cpyParms);
+  cudaCheck(cudaMemcpy3D(&cpyParms));
 
   cpyParms.srcPtr.ptr = Bz.data();
   cpyParms.dstArray = magnetic_field->array_Bz;
-  cudaMemcpy3D(&cpyParms);
+  cudaCheck(cudaMemcpy3D(&cpyParms));
 
 #else
   Allen::memcpy(magnetic_field->Bx, Bx.data(), sizeof(float) * N, Allen::memcpyHostToDevice);
-- 
GitLab


From 35f3e7830123fcd5a3e0e53a2ce1d32cd97bcff1 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Wed, 9 Apr 2025 17:14:34 +0200
Subject: [PATCH 03/10] Implement Runge-Kutta extrapolator familly

---
 .../ParKalman/include/ButcherTableau.cuh      | 91 +++++++++++++++++++
 .../ParKalman/include/ExtrapolateStates.cuh   |  1 +
 .../ParKalman/include/ExtrapolatorCommon.cuh  | 46 ++++++++++
 .../include/ParabolicExtrapolator.cuh         | 33 ++-----
 .../include/RungeKuttaExtrapolator.cuh        | 43 +++++++++
 .../kalman/ParKalman/src/ExtrapolateStates.cu | 18 ++--
 6 files changed, 199 insertions(+), 33 deletions(-)
 create mode 100644 device/kalman/ParKalman/include/ButcherTableau.cuh
 create mode 100644 device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
 create mode 100644 device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh

diff --git a/device/kalman/ParKalman/include/ButcherTableau.cuh b/device/kalman/ParKalman/include/ButcherTableau.cuh
new file mode 100644
index 00000000000..ce22b96e9f2
--- /dev/null
+++ b/device/kalman/ParKalman/include/ButcherTableau.cuh
@@ -0,0 +1,91 @@
+/*****************************************************************************\
+* (c) Copyright 2025 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+// This is a workaround to cuda not allowing static constexpr data member on device
+// Hidden in macros to make the definitions more readable
+// It checks the number of values match N_stages
+#define a_table(...)                                                                                                 \
+  __device__ static constexpr ftype a(int i, int j)                                                                  \
+  {                                                                                                                  \
+    ftype a[] {__VA_ARGS__};                                                                                         \
+    static_assert(N_stages * (N_stages - 1) / 2 == 0 || N_stages * (N_stages - 1) / 2 == sizeof(a) / sizeof(ftype)); \
+    return a[i * (i - 1) / 2 + j];                                                                                   \
+  }
+#define b_table(...)                         \
+  __device__ static constexpr ftype b(int i) \
+  {                                          \
+    ftype b[N_stages] {__VA_ARGS__};         \
+    return b[i];                             \
+  }
+#define b_star_table(...)                         \
+  __device__ static constexpr ftype b_star(int i) \
+  {                                               \
+    ftype b_star[N_stages] {__VA_ARGS__};         \
+    return b_star[i];                             \
+  }
+
+namespace ButcherTableau {
+  // c1 |
+  // c2 | a1
+  // c3 | a2 a3
+  // .. | a4 a5 a6
+  // cs | ...
+  // ------------------
+  //    | b1 b2 B3 ... bs
+
+  template<typename ftype = float>
+  struct Euler {
+    static constexpr unsigned N_stages = 1;
+    a_table(0.) b_table(1.)
+  };
+
+  template<typename ftype = float>
+  struct HeunEuler {
+    static constexpr unsigned N_stages = 2;
+    a_table(1.) b_table(1. / 2., 1. / 2.)
+  };
+
+  template<typename ftype = float>
+  struct RK4 {
+    static constexpr unsigned N_stages = 4;
+    a_table(
+      1. / 2., //
+      0.,
+      1. / 2., //
+      0.,
+      0.,
+      1. //
+      ) b_table(1. / 6., 1. / 3., 1. / 3., 1. / 6.)
+  };
+
+  template<typename ftype = float>
+  struct CashKarp {
+    static constexpr unsigned N_stages = 6;
+    a_table(
+      1. / 5., //
+      3. / 40.,
+      9. / 40., //
+      3. / 10.,
+      -9. / 10.,
+      6. / 5., //
+      -11. / 54.,
+      5. / 2.,
+      -70. / 27.,
+      35. / 27., //
+      1631. / 55296.,
+      175. / 512.,
+      575. / 13824.,
+      44275. / 110592.,
+      253. / 4096.) b_table(37. / 378., 0., 250. / 621., 125. / 594., 0., 512. / 1771.)
+      b_star_table(2825. / 27648., 0., 18575. / 48384., 13525. / 55296., 277. / 14336., 1. / 4.)
+  };
+} // namespace ButcherTableau
\ No newline at end of file
diff --git a/device/kalman/ParKalman/include/ExtrapolateStates.cuh b/device/kalman/ParKalman/include/ExtrapolateStates.cuh
index 6d42efec7b9..acba44ef92e 100644
--- a/device/kalman/ParKalman/include/ExtrapolateStates.cuh
+++ b/device/kalman/ParKalman/include/ExtrapolateStates.cuh
@@ -13,6 +13,7 @@
 #include "AlgorithmTypes.cuh"
 #include "States.cuh"
 #include "ParabolicExtrapolator.cuh"
+#include "RungeKuttaExtrapolator.cuh"
 
 namespace extrapolate_states {
   struct Parameters {
diff --git a/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh b/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
new file mode 100644
index 00000000000..13245b74bec
--- /dev/null
+++ b/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
@@ -0,0 +1,46 @@
+/*****************************************************************************\
+* (c) Copyright 2025 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+#include <MagneticField.cuh>
+
+namespace Extrapolators {
+  constexpr float c_light = 2.99792458e+8 * 1000. / 1.e+9; // 299.792458 (float=299.792449951171875)
+  constexpr float eplus = 1.;
+
+  struct State {
+    float x {0}, y {0}, z {0}, tx {0}, ty {0}, qop {0};
+
+    struct Derivative {
+      float dx, dy, dz, dtx, dty;
+      __device__ friend Derivative operator*(const Derivative& a, float b)
+      {
+        return {a.dx * b, a.dy * b, a.dz * b, a.dtx * b, a.dty * b};
+      }
+    };
+
+    // dState / Dz
+    __device__ friend Derivative derivative(const State& v, const float3& B)
+    {
+      const auto tx2 = v.tx * v.tx;
+      const auto ty2 = v.ty * v.ty;
+      const auto norm = sqrtf(1.f + tx2 + ty2);
+      const auto ax = norm * (v.ty * (v.tx * B.x + B.z) - (1.f + tx2) * B.y);
+      const auto ay = norm * (-v.tx * (v.ty * B.y + B.z) + (1.f + ty2) * B.x);
+      return {v.tx, v.ty, 1.f, v.qop * ax, v.qop * ay};
+    }
+
+    __device__ friend State operator+(const State& a, const Derivative& b)
+    {
+      return {a.x + b.dx, a.y + b.dy, a.z + b.dz, a.tx + b.dtx, a.ty + b.dty};
+    }
+  };
+} // namespace Extrapolators
\ No newline at end of file
diff --git a/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh b/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh
index 6a01ce7f904..ae70f7465df 100644
--- a/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh
+++ b/device/kalman/ParKalman/include/ParabolicExtrapolator.cuh
@@ -10,13 +10,10 @@
 \*****************************************************************************/
 #pragma once
 
+#include <ExtrapolatorCommon.cuh>
 #include <MagneticField.cuh>
 
 namespace Extrapolators {
-  struct State {
-    float x, y, z, tx, ty, qop;
-  };
-
   template<typename ftype = float>
   struct ParabolicExtrapolator {
     // Implementation taken from:
@@ -27,31 +24,17 @@ namespace Extrapolators {
       // Get the B field at midpoint
       const ftype xMid = state.x + state.tx * dz * ftype {0.5};
       const ftype yMid = state.y + state.ty * dz * ftype {0.5};
-      const float3 P = make_float3(xMid, yMid, state.z + dz * ftype {0.5});
+      const ftype zMid = state.z + dz * ftype {0.5};
+      const float3 P = make_float3(xMid, yMid, zMid);
       const float3 B = field.fieldVectorLinearInterpolation(P);
 
-      const ftype tx = state.tx;
-      const ftype ty = state.ty;
-      const ftype norm = fsqrt(ftype {1.} + tx * tx + ty * ty);
-
-      const ftype ax = norm * (ty * (tx * B.x + B.z) - ((ftype {1.} + tx * tx) * B.y));
-      const ftype ay = norm * (-tx * (ty * B.y + B.z) + ((ftype {1.} + ty * ty) * B.x));
-
-      const ftype fact = eplus * c_light * dz * state.qop;
-
       // Update the state parameters (exact extrapolation)
-      state.x += dz * (tx + ftype {0.5} * ax * fact);
-      state.y += dz * (ty + ftype {0.5} * ay * fact);
+      const auto d = derivative(state, B);
+      state.x += dz * (d.dx + ftype {0.5} * dz * d.dtx);
+      state.y += dz * (d.dy + ftype {0.5} * dz * d.dty);
       state.z += dz;
-      state.tx += ax * fact;
-      state.ty += ay * fact;
+      state.tx += dz * d.dtx;
+      state.ty += dz * d.dty;
     }
-
-  private:
-    static constexpr ftype c_light = 2.99792458e+8 * 1000. / 1.e+9;
-    static constexpr ftype eplus = 1.;
-
-    __device__ static float fsqrt(float x) { return sqrtf(x); }
-    __device__ static double fsqrt(double x) { return sqrt(x); }
   };
 } // namespace Extrapolators
\ No newline at end of file
diff --git a/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh b/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
new file mode 100644
index 00000000000..2d6be8bc1db
--- /dev/null
+++ b/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
@@ -0,0 +1,43 @@
+/*****************************************************************************\
+* (c) Copyright 2025 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+#include <ButcherTableau.cuh>
+#include <ExtrapolatorCommon.cuh>
+#include <MagneticField.cuh>
+
+namespace Extrapolators {
+  template<typename ftype = float, typename Table = ButcherTableau::CashKarp<ftype>>
+  struct RungeKuttaExtrapolator {
+    // Implementation taken from:
+    // https://gitlab.cern.ch/lhcb/Rec/-/blob/master/Tr/TrackExtrapolators/src/TrackRungeKuttaExtrapolator.cpp
+
+    __device__ static void propagate(State& state, ftype dz, const MagneticField::Magfield& field)
+    {
+      State::Derivative k[Table::N_stages];
+#pragma unroll
+      for (int stage = 0; stage < Table::N_stages; stage++) {
+        State s = state;
+#pragma unroll
+        for (int i = 0; i < stage - 1; i++) {
+          s = s + k[i] * Table::a(stage, i);
+        }
+        float3 B = field.fieldVectorLinearInterpolation(make_float3(s.x, s.y, s.z));
+        k[stage] = derivative(state, B) * dz;
+      }
+
+#pragma unroll
+      for (int i = 0; i < Table::N_stages; i++) {
+        state = state + k[i] * Table::b(i);
+      }
+    }
+  };
+} // namespace Extrapolators
\ No newline at end of file
diff --git a/device/kalman/ParKalman/src/ExtrapolateStates.cu b/device/kalman/ParKalman/src/ExtrapolateStates.cu
index df99fae27f6..9f1bfe6b207 100644
--- a/device/kalman/ParKalman/src/ExtrapolateStates.cu
+++ b/device/kalman/ParKalman/src/ExtrapolateStates.cu
@@ -34,16 +34,18 @@ __global__ void extrapolate_states_kernel(
   for (unsigned i = blockIdx.x * blockDim.x + threadIdx.x; i < n_states; i += blockDim.x * gridDim.x) {
     const auto& input = input_states->state(i);
 
-    Extrapolators::State output;
-    output.x = input.x();
-    output.y = input.y();
-    output.z = input.z();
-    output.tx = input.tx();
-    output.ty = input.ty();
-    output.qop = input.qop();
+    Extrapolators::State output {
+      input.x(),
+      input.y(),
+      input.z(),
+      input.tx(),
+      input.ty(),
+      input.qop() * Extrapolators::c_light * Extrapolators::eplus // Don't forget to do this conversion !
+    };
 
     for (unsigned s = 0; s < n_steps; s++) {
-      Extrapolators::ParabolicExtrapolator<float>::propagate(output, dz, field);
+      // Extrapolators::ParabolicExtrapolator<float>::propagate(output, dz, field);
+      Extrapolators::RungeKuttaExtrapolator<float, ButcherTableau::CashKarp<float>>::propagate(output, dz, field);
 
       /*if (threadIdx.x == 0 && blockIdx.x == 0) {
         float3 B = field.fieldVectorLinearInterpolation(make_float3(output.x, output.y, output.z));
-- 
GitLab


From 4d23c90945b5117c791739fbad030e74d83e487b Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Tue, 15 Apr 2025 12:09:46 +0200
Subject: [PATCH 04/10] Provide all vector types for CPU backend

---
 backend/include/CPUBackend.h | 101 ++++++++---------------------------
 1 file changed, 21 insertions(+), 80 deletions(-)

diff --git a/backend/include/CPUBackend.h b/backend/include/CPUBackend.h
index 8a5a1833ec1..2a3fd6369a6 100644
--- a/backend/include/CPUBackend.h
+++ b/backend/include/CPUBackend.h
@@ -72,86 +72,27 @@ unsigned inline __ballot_sync(unsigned mask, int predicate) { return predicate &
   auto _dynamic_shared_memory_buffer = std::vector<_type>(_config.dynamic_shared_memory_size() / sizeof(_type)); \
   auto _instance = _dynamic_shared_memory_buffer.data();
 
-struct alignas(unsigned long long int) uint2 {
-  unsigned int x;
-  unsigned int y;
-};
-
-inline uint2 make_uint2(unsigned int x, unsigned int y)
-{
-  uint2 out;
-  out.x = x;
-  out.y = y;
-  return out;
-}
-
-struct int2 {
-  int x;
-  int y;
-};
-
-struct char2 {
-  char x;
-  char y;
-};
-
-struct uint4 {
-  unsigned int x;
-  unsigned int y;
-  unsigned int z;
-  unsigned int w;
-};
-
-struct alignas(unsigned int) ushort2 {
-  unsigned short x;
-  unsigned short y;
-};
-
-inline ushort2 make_ushort2(ushort x, ushort y)
-{
-  ushort2 out;
-  out.x = x;
-  out.y = y;
-  return out;
-}
-
-struct ushort4 {
-  unsigned short x;
-  unsigned short y;
-  unsigned short z;
-  unsigned short w;
-};
-
-struct short2 {
-  short x;
-  short y;
-};
-
-inline short2 make_short2(short x, short y)
-{
-  short2 out;
-  out.x = x;
-  out.y = y;
-  return out;
-}
-struct float3 {
-  float x;
-  float y;
-  float z;
-};
-
-struct float2 {
-  float x;
-  float y;
-};
-struct float4 {
-  float x;
-  float y;
-  float z;
-  float w;
-};
-
-struct dim3 {
+#define struct1(type)                                   \
+  struct alignas(1 * sizeof(type)) type##1 { type x; }; \
+  inline type##1 make_##type##1(type x) { return {x}; }
+#define struct2(type)                                      \
+  struct alignas(2 * sizeof(type)) type##2 { type x, y; }; \
+  inline type##2 make_##type##2(type x, type y) { return {x, y}; }
+#define struct3(type)                                         \
+  struct alignas(1 * sizeof(type)) type##3 { type x, y, z; }; \
+  inline type##3 make_##type##3(type x, type y, type z) { return {x, y, z}; }
+#define struct4(type)                                            \
+  struct alignas(4 * sizeof(type)) type##4 { type x, y, z, w; }; \
+  inline type##4 make_##type##4(type x, type y, type z, type w) { return {x, y, z, w}; }
+#define vectype(type) struct1(type) struct2(type) struct3(type) struct4(type)
+
+vectype(char) using uchar = unsigned char;
+vectype(uchar) vectype(short) vectype(ushort) vectype(int) vectype(uint) vectype(long) using ulong = unsigned long;
+vectype(ulong) using longlong = long long;
+vectype(longlong) using ulonglong = unsigned long long;
+vectype(ulonglong) vectype(float) vectype(double)
+
+  struct dim3 {
   unsigned int x = 1;
   unsigned int y = 1;
   unsigned int z = 1;
-- 
GitLab


From 36f01fda73055109c87eb305dab8b6920f343756 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Tue, 15 Apr 2025 12:10:17 +0200
Subject: [PATCH 05/10] Fix cpu warnings

---
 device/kalman/ParKalman/include/ButcherTableau.cuh     | 10 +++++-----
 .../ParKalman/include/RungeKuttaExtrapolator.cuh       |  2 +-
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/device/kalman/ParKalman/include/ButcherTableau.cuh b/device/kalman/ParKalman/include/ButcherTableau.cuh
index ce22b96e9f2..40d54d99720 100644
--- a/device/kalman/ParKalman/include/ButcherTableau.cuh
+++ b/device/kalman/ParKalman/include/ButcherTableau.cuh
@@ -44,19 +44,19 @@ namespace ButcherTableau {
 
   template<typename ftype = float>
   struct Euler {
-    static constexpr unsigned N_stages = 1;
+    static constexpr int N_stages = 1;
     a_table(0.) b_table(1.)
   };
 
   template<typename ftype = float>
   struct HeunEuler {
-    static constexpr unsigned N_stages = 2;
+    static constexpr int N_stages = 2;
     a_table(1.) b_table(1. / 2., 1. / 2.)
   };
 
   template<typename ftype = float>
   struct RK4 {
-    static constexpr unsigned N_stages = 4;
+    static constexpr int N_stages = 4;
     a_table(
       1. / 2., //
       0.,
@@ -69,7 +69,7 @@ namespace ButcherTableau {
 
   template<typename ftype = float>
   struct CashKarp {
-    static constexpr unsigned N_stages = 6;
+    static constexpr int N_stages = 6;
     a_table(
       1. / 5., //
       3. / 40.,
@@ -88,4 +88,4 @@ namespace ButcherTableau {
       253. / 4096.) b_table(37. / 378., 0., 250. / 621., 125. / 594., 0., 512. / 1771.)
       b_star_table(2825. / 27648., 0., 18575. / 48384., 13525. / 55296., 277. / 14336., 1. / 4.)
   };
-} // namespace ButcherTableau
\ No newline at end of file
+} // namespace ButcherTableau
diff --git a/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh b/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
index 2d6be8bc1db..5ed6d5f44bf 100644
--- a/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
+++ b/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
@@ -40,4 +40,4 @@ namespace Extrapolators {
       }
     }
   };
-} // namespace Extrapolators
\ No newline at end of file
+} // namespace Extrapolators
-- 
GitLab


From 8eec9e735f4406fa35f744b5e26a29b61357556a Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Thu, 10 Jul 2025 14:16:27 +0200
Subject: [PATCH 06/10] Fix CPU pragma unroll warnings, implement error in RK

---
 backend/include/BackendCommon.h               | 11 +++++++++-
 .../include/hybrid_seeding_helpers.cuh        |  7 ++----
 .../looking_forward/src/LFTripletSeeding.cu   |  7 ++----
 device/downstream/common/DownstreamHelper.cuh |  4 +---
 .../downstream/common/DownstreamV2Helper.cuh  |  8 ++-----
 .../common/include/MagneticField.cuh          | 10 ++++-----
 .../ParKalman/include/ButcherTableau.cuh      |  2 +-
 .../ParKalman/include/ExtrapolatorCommon.cuh  | 17 +++++++++++++-
 .../include/RungeKuttaExtrapolator.cuh        | 11 ++++++----
 .../kalman/ParKalman/src/ExtrapolateStates.cu |  3 ++-
 .../track_matching/match/src/TrackMatching.cu |  3 ---
 .../mva_models/include/SingleLayerFCNN.cuh    | 22 ++++++-------------
 .../non_event_data/include/Consumers.h        |  2 +-
 13 files changed, 56 insertions(+), 51 deletions(-)

diff --git a/backend/include/BackendCommon.h b/backend/include/BackendCommon.h
index 6d200e59392..4616bee3d32 100644
--- a/backend/include/BackendCommon.h
+++ b/backend/include/BackendCommon.h
@@ -39,6 +39,15 @@
 #include "CUDABackend.h"
 #endif
 
+#define DO_PRAGMA_(x) _Pragma(#x)
+#define DO_PRAGMA(x) DO_PRAGMA_(x)
+
+#if defined(__clang__) or defined(__NVCC__)
+#define UNROLL(n) DO_PRAGMA(unroll)
+#elif defined(__GNUC__)
+#define UNROLL(n) DO_PRAGMA(GCC unroll n)
+#endif
+
 #if defined(DEVICE_COMPILER)
 namespace Allen {
   namespace device {
@@ -231,4 +240,4 @@ namespace Allen {
     __builtin_unreachable();
 #endif
   }
-} // namespace Allen
\ No newline at end of file
+} // namespace Allen
diff --git a/device/SciFi/hybridseeding/include/hybrid_seeding_helpers.cuh b/device/SciFi/hybridseeding/include/hybrid_seeding_helpers.cuh
index c7ecc8feaf7..dd6854dad99 100644
--- a/device/SciFi/hybridseeding/include/hybrid_seeding_helpers.cuh
+++ b/device/SciFi/hybridseeding/include/hybrid_seeding_helpers.cuh
@@ -11,6 +11,7 @@
 #pragma once
 
 #include <tuple>
+#include <BackendCommon.h>
 #include "States.cuh"
 #include "SciFiEventModel.cuh"
 
@@ -26,11 +27,7 @@ namespace hybrid_seeding {
     unsigned int size = array_size;
 
     // Unroll 10 time to cover arrays of size max 1024
-#if defined(__clang__) or defined(__NVCC__)
-#pragma unroll
-#elif defined(__GNUC__)
-#pragma GCC unroll 10
-#endif
+    UNROLL(10)
     for (unsigned int step = 0; step < 10; step++) {
       unsigned int half = size / 2;
       low += (array[low + half] < needle) * (size - half);
diff --git a/device/SciFi/looking_forward/src/LFTripletSeeding.cu b/device/SciFi/looking_forward/src/LFTripletSeeding.cu
index b3fa5e03486..5b348582e51 100644
--- a/device/SciFi/looking_forward/src/LFTripletSeeding.cu
+++ b/device/SciFi/looking_forward/src/LFTripletSeeding.cu
@@ -8,6 +8,7 @@
 * granted to it by virtue of its status as an Intergovernmental Organization  *
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
+#include "BackendCommon.h"
 #include "LFTripletSeeding.cuh"
 #include "LookingForwardTools.cuh"
 #include "BinarySearch.cuh"
@@ -65,11 +66,7 @@ __device__ unsigned int binary_search_leftmost_unrolled(const T* array, const un
   unsigned int size = array_size;
 
   // Unroll 9 time to cover arrays of size max 512
-#if defined(__clang__) or defined(__NVCC__)
-#pragma unroll
-#elif defined(__GNUC__)
-#pragma GCC unroll 9
-#endif
+  UNROLL(9)
   for (unsigned int step = 0; step < 9; step++) {
     unsigned int half = size / 2;
     low += (array[low + half] < needle) * (size - half);
diff --git a/device/downstream/common/DownstreamHelper.cuh b/device/downstream/common/DownstreamHelper.cuh
index 284b7e87205..cdd45b053b0 100644
--- a/device/downstream/common/DownstreamHelper.cuh
+++ b/device/downstream/common/DownstreamHelper.cuh
@@ -27,9 +27,7 @@ namespace Downstream {
     __device__ inline float find_root(f_t const& f, fprime_t const& fprime, const float x0)
     {
       float x = x0;
-#if (defined(TARGET_DEVICE_CUDA) && defined(__CUDACC__))
-#pragma unroll
-#endif
+      UNROLL(5)
       for (unsigned i = 0; i < NumIteration; i++) {
         const auto h = f(x) / fprime(x);
         x -= h;
diff --git a/device/downstream/common/DownstreamV2Helper.cuh b/device/downstream/common/DownstreamV2Helper.cuh
index 6fd2653c54f..c9be0adc826 100644
--- a/device/downstream/common/DownstreamV2Helper.cuh
+++ b/device/downstream/common/DownstreamV2Helper.cuh
@@ -149,14 +149,10 @@ namespace Downstream::Helpers {
   {
     // Fetch sector idx
     unsigned sector_idx = 0;
-#if defined(__clang__) or defined(__NVCC__)
-#pragma unroll
-#elif defined(__GNUC__)
-#pragma GCC unroll 16
-#endif
+    UNROLL(16)
     for (unsigned i = 0; i < 16; i++) {
       sector_idx += (hit_idx >= selected_offsets[i]);
     }
     return sector_idx;
   }
-} // namespace Downstream::Helpers
\ No newline at end of file
+} // namespace Downstream::Helpers
diff --git a/device/event_model/common/include/MagneticField.cuh b/device/event_model/common/include/MagneticField.cuh
index 231d50af94f..1d7618299d6 100644
--- a/device/event_model/common/include/MagneticField.cuh
+++ b/device/event_model/common/include/MagneticField.cuh
@@ -15,17 +15,17 @@ namespace MagneticField {
 #ifdef MAGFIELD_USE_TEXTURE
     // https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#texture-and-surface-memory
 
-    __device__ float3 fieldVectorLinearInterpolation([[maybe_unused]] float3 pos) const
-    {
 #ifdef __CUDA_ARCH__
+    __device__ float3 fieldVectorLinearInterpolation(float3 pos) const
+    {
       const float x = (pos.x - minX) * invDx;
       const float y = (pos.y - minY) * invDy;
       const float z = (pos.z - minZ) * invDz;
       return {tex3D<float>(tex_Bx, x, y, z), tex3D<float>(tex_By, x, y, z), tex3D<float>(tex_Bz, x, y, z)};
+    }
 #else
-      return {1, 2, 3};
+    __device__ float3 fieldVectorLinearInterpolation(float3 pos) const;
 #endif
-    }
 
     cudaArray_t array_Bx;
     cudaArray_t array_By;
@@ -88,4 +88,4 @@ namespace MagneticField {
     float minX, minY, minZ;
   };
 
-} // namespace MagneticField
\ No newline at end of file
+} // namespace MagneticField
diff --git a/device/kalman/ParKalman/include/ButcherTableau.cuh b/device/kalman/ParKalman/include/ButcherTableau.cuh
index 40d54d99720..f3dd5606b87 100644
--- a/device/kalman/ParKalman/include/ButcherTableau.cuh
+++ b/device/kalman/ParKalman/include/ButcherTableau.cuh
@@ -51,7 +51,7 @@ namespace ButcherTableau {
   template<typename ftype = float>
   struct HeunEuler {
     static constexpr int N_stages = 2;
-    a_table(1.) b_table(1. / 2., 1. / 2.)
+    a_table(1.) b_table(1., 0.) b_star_table(1. / 2., 1. / 2.)
   };
 
   template<typename ftype = float>
diff --git a/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh b/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
index 13245b74bec..6ab468007cc 100644
--- a/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
+++ b/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
@@ -27,6 +27,21 @@ namespace Extrapolators {
       }
     };
 
+    struct Error {
+      float x {0}, y {0}, tx {0}, ty {0};
+      __device__ void clear()
+      {
+        x = 0;
+        y = 0;
+        tx = 0;
+        ty = 0;
+      }
+      __device__ friend Error operator+(const Error& a, const Derivative& b)
+      {
+        return {a.x + b.dx, a.y + b.dy, a.tx + b.dtx, a.ty + b.dty};
+      }
+    };
+
     // dState / Dz
     __device__ friend Derivative derivative(const State& v, const float3& B)
     {
@@ -43,4 +58,4 @@ namespace Extrapolators {
       return {a.x + b.dx, a.y + b.dy, a.z + b.dz, a.tx + b.dtx, a.ty + b.dty};
     }
   };
-} // namespace Extrapolators
\ No newline at end of file
+} // namespace Extrapolators
diff --git a/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh b/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
index 5ed6d5f44bf..157e9d7ac4f 100644
--- a/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
+++ b/device/kalman/ParKalman/include/RungeKuttaExtrapolator.cuh
@@ -10,6 +10,7 @@
 \*****************************************************************************/
 #pragma once
 
+#include <BackendCommon.h>
 #include <ButcherTableau.cuh>
 #include <ExtrapolatorCommon.cuh>
 #include <MagneticField.cuh>
@@ -20,13 +21,13 @@ namespace Extrapolators {
     // Implementation taken from:
     // https://gitlab.cern.ch/lhcb/Rec/-/blob/master/Tr/TrackExtrapolators/src/TrackRungeKuttaExtrapolator.cpp
 
-    __device__ static void propagate(State& state, ftype dz, const MagneticField::Magfield& field)
+    __device__ static void propagate(State& state, State::Error& err, ftype dz, const MagneticField::Magfield& field)
     {
       State::Derivative k[Table::N_stages];
-#pragma unroll
+      UNROLL(10)
       for (int stage = 0; stage < Table::N_stages; stage++) {
         State s = state;
-#pragma unroll
+        UNROLL(10)
         for (int i = 0; i < stage - 1; i++) {
           s = s + k[i] * Table::a(stage, i);
         }
@@ -34,8 +35,10 @@ namespace Extrapolators {
         k[stage] = derivative(state, B) * dz;
       }
 
-#pragma unroll
+      err.clear();
+      UNROLL(10)
       for (int i = 0; i < Table::N_stages; i++) {
+        err = err + k[i] * (Table::b(i) - Table::b_star(i));
         state = state + k[i] * Table::b(i);
       }
     }
diff --git a/device/kalman/ParKalman/src/ExtrapolateStates.cu b/device/kalman/ParKalman/src/ExtrapolateStates.cu
index 9f1bfe6b207..53327ca727b 100644
--- a/device/kalman/ParKalman/src/ExtrapolateStates.cu
+++ b/device/kalman/ParKalman/src/ExtrapolateStates.cu
@@ -45,7 +45,8 @@ __global__ void extrapolate_states_kernel(
 
     for (unsigned s = 0; s < n_steps; s++) {
       // Extrapolators::ParabolicExtrapolator<float>::propagate(output, dz, field);
-      Extrapolators::RungeKuttaExtrapolator<float, ButcherTableau::CashKarp<float>>::propagate(output, dz, field);
+      Extrapolators::State::Error err;
+      Extrapolators::RungeKuttaExtrapolator<float, ButcherTableau::CashKarp<float>>::propagate(output, err, dz, field);
 
       /*if (threadIdx.x == 0 && blockIdx.x == 0) {
         float3 B = field.fieldVectorLinearInterpolation(make_float3(output.x, output.y, output.z));
diff --git a/device/track_matching/match/src/TrackMatching.cu b/device/track_matching/match/src/TrackMatching.cu
index ba992417b52..8b74a36ba14 100644
--- a/device/track_matching/match/src/TrackMatching.cu
+++ b/device/track_matching/match/src/TrackMatching.cu
@@ -407,9 +407,6 @@ __global__ void track_matching::track_matching_add_ut_hits(
   // Try up to MaxNumIteration times.
   // If a candidate overflows, try using a tighter threshold.
   // If it still exceeds after MaxNumIteration times, simply reset the event to 0.
-#if (defined(TARGET_DEVICE_CUDA) && defined(__CUDACC__))
-#pragma unroll
-#endif
   for (unsigned iteration = 0; iteration < MaxNumIteration; iteration++) {
     bool overflow = false;
     for (unsigned layer = 0; layer < UT::Constants::n_layers; layer++) {
diff --git a/device/utils/mva_models/include/SingleLayerFCNN.cuh b/device/utils/mva_models/include/SingleLayerFCNN.cuh
index f80e2033afb..827d09a6039 100644
--- a/device/utils/mva_models/include/SingleLayerFCNN.cuh
+++ b/device/utils/mva_models/include/SingleLayerFCNN.cuh
@@ -134,24 +134,18 @@ template<unsigned num_input, unsigned num_node>
 __device__ inline float Allen::MVAModels::DeviceSingleLayerFCNN<num_input, num_node>::evaluate(float* input) const
 {
   using ModelType = Allen::MVAModels::DeviceSingleLayerFCNN<num_input, num_node>;
-// Data preprocessing
-#if (defined(TARGET_DEVICE_CUDA) && defined(__CUDACC__))
-#pragma unroll
-#endif
+  // Data preprocessing
+  UNROLL(100)
   for (unsigned i = 0; i < ModelType::nInput; i++) {
     // input[i] = __fdividef(input[i] - input_mean[i], input_std[i]);
     input[i] = (input[i] - mean[i]) / std[i];
   }
   float h1[ModelType::nNode] = {0.f};
 
-// First layer
-#if (defined(TARGET_DEVICE_CUDA) && defined(__CUDACC__))
-#pragma unroll
-#endif
+  // First layer
+  UNROLL(100)
   for (unsigned i = 0; i < ModelType::nNode; i++) {
-#if (defined(TARGET_DEVICE_CUDA) && defined(__CUDACC__))
-#pragma unroll
-#endif
+    UNROLL(100)
     for (unsigned j = 0; j < ModelType::nInput; j++) {
       h1[i] += input[j] * weights1[i][j];
     }
@@ -160,9 +154,7 @@ __device__ inline float Allen::MVAModels::DeviceSingleLayerFCNN<num_input, num_n
 
   // Output layer
   float output = 0.f;
-#if (defined(TARGET_DEVICE_CUDA) && defined(__CUDACC__))
-#pragma unroll
-#endif
+  UNROLL(100)
   for (unsigned i = 0; i < ModelType::nNode; i++) {
     output += h1[i] * weights2[i];
   }
@@ -170,4 +162,4 @@ __device__ inline float Allen::MVAModels::DeviceSingleLayerFCNN<num_input, num_n
   output = ActivateFunction::sigmoid(output + bias2);
 
   return output;
-}
\ No newline at end of file
+}
diff --git a/integration/non_event_data/include/Consumers.h b/integration/non_event_data/include/Consumers.h
index e723b28c283..4b058ee7f3a 100644
--- a/integration/non_event_data/include/Consumers.h
+++ b/integration/non_event_data/include/Consumers.h
@@ -124,7 +124,7 @@ namespace Consumers {
   private:
     std::reference_wrapper<Constants> m_constants;
   };
-  
+
   struct MagneticFieldPolarity final : public Allen::NonEventData::Consumer {
   public:
     MagneticFieldPolarity(gsl::span<float>&, std::vector<float>&);
-- 
GitLab


From f93cd637b89739af322dc43b0ae71e733ef361e8 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Thu, 10 Jul 2025 14:26:34 +0200
Subject: [PATCH 07/10] Add symbolic links to magfield.bin on cvmfs

---
 input/allen_geometries/SMOG2_exp2024/magfield.bin                | 1 +
 .../magfield.bin                                                 | 1 +
 .../geometry_dddb-20230313_sim-20230626-vc-md100/magfield.bin    | 1 +
 .../geometry_dddb-20231010_sim-20231116-vc-md100/magfield.bin    | 1 +
 .../geometry_dddb-20231017_sim-20231017-vc-md100/magfield.bin    | 1 +
 .../magfield.bin                                                 | 1 +
 .../magfield.bin                                                 | 1 +
 .../allen_geometries/geometry_run3_2024.Q1.2-v00.00/magfield.bin | 1 +
 .../allen_geometries/geometry_sim-20231017-vc-mu100/magfield.bin | 1 +
 .../magfield.bin                                                 | 1 +
 .../upgrade_DC19_01_MinBiasMD_new_UT_geometry/magfield.bin       | 1 +
 11 files changed, 11 insertions(+)
 create mode 120000 input/allen_geometries/SMOG2_exp2024/magfield.bin
 create mode 120000 input/allen_geometries/geometry_dddb-20180815_sim-20180530-vc-md100_new_UT_geometry/magfield.bin
 create mode 120000 input/allen_geometries/geometry_dddb-20230313_sim-20230626-vc-md100/magfield.bin
 create mode 120000 input/allen_geometries/geometry_dddb-20231010_sim-20231116-vc-md100/magfield.bin
 create mode 120000 input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100/magfield.bin
 create mode 120000 input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100_new_SciFi_geometry/magfield.bin
 create mode 120000 input/allen_geometries/geometry_dddb-20240427_sim10-2024.Q1.2-v1.1-md100/magfield.bin
 create mode 120000 input/allen_geometries/geometry_run3_2024.Q1.2-v00.00/magfield.bin
 create mode 120000 input/allen_geometries/geometry_sim-20231017-vc-mu100/magfield.bin
 create mode 120000 input/allen_geometries/upgrade-magdown-sim10-up08-11102202-digi_new_UT_geometry/magfield.bin
 create mode 120000 input/allen_geometries/upgrade_DC19_01_MinBiasMD_new_UT_geometry/magfield.bin

diff --git a/input/allen_geometries/SMOG2_exp2024/magfield.bin b/input/allen_geometries/SMOG2_exp2024/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/SMOG2_exp2024/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_dddb-20180815_sim-20180530-vc-md100_new_UT_geometry/magfield.bin b/input/allen_geometries/geometry_dddb-20180815_sim-20180530-vc-md100_new_UT_geometry/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_dddb-20180815_sim-20180530-vc-md100_new_UT_geometry/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_dddb-20230313_sim-20230626-vc-md100/magfield.bin b/input/allen_geometries/geometry_dddb-20230313_sim-20230626-vc-md100/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_dddb-20230313_sim-20230626-vc-md100/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_dddb-20231010_sim-20231116-vc-md100/magfield.bin b/input/allen_geometries/geometry_dddb-20231010_sim-20231116-vc-md100/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_dddb-20231010_sim-20231116-vc-md100/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100/magfield.bin b/input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100_new_SciFi_geometry/magfield.bin b/input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100_new_SciFi_geometry/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_dddb-20231017_sim-20231017-vc-md100_new_SciFi_geometry/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_dddb-20240427_sim10-2024.Q1.2-v1.1-md100/magfield.bin b/input/allen_geometries/geometry_dddb-20240427_sim10-2024.Q1.2-v1.1-md100/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_dddb-20240427_sim10-2024.Q1.2-v1.1-md100/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_run3_2024.Q1.2-v00.00/magfield.bin b/input/allen_geometries/geometry_run3_2024.Q1.2-v00.00/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/geometry_run3_2024.Q1.2-v00.00/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/geometry_sim-20231017-vc-mu100/magfield.bin b/input/allen_geometries/geometry_sim-20231017-vc-mu100/magfield.bin
new file mode 120000
index 00000000000..56df97de512
--- /dev/null
+++ b/input/allen_geometries/geometry_sim-20231017-vc-mu100/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.up.bin
\ No newline at end of file
diff --git a/input/allen_geometries/upgrade-magdown-sim10-up08-11102202-digi_new_UT_geometry/magfield.bin b/input/allen_geometries/upgrade-magdown-sim10-up08-11102202-digi_new_UT_geometry/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/upgrade-magdown-sim10-up08-11102202-digi_new_UT_geometry/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
diff --git a/input/allen_geometries/upgrade_DC19_01_MinBiasMD_new_UT_geometry/magfield.bin b/input/allen_geometries/upgrade_DC19_01_MinBiasMD_new_UT_geometry/magfield.bin
new file mode 120000
index 00000000000..6cc91b95fcf
--- /dev/null
+++ b/input/allen_geometries/upgrade_DC19_01_MinBiasMD_new_UT_geometry/magfield.bin
@@ -0,0 +1 @@
+/cvmfs/lhcb.cern.ch/lib/lhcb/DBASE/FieldMap/v8r1/cdf/field.v8r1.down.bin
\ No newline at end of file
-- 
GitLab


From ab6fa9aa8b90cd4e88e058404eb30c85bc5b93e6 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Thu, 21 Aug 2025 16:01:26 +0200
Subject: [PATCH 08/10] Comment extrapolated states test

---
 configuration/python/AllenConf/HLT1.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/configuration/python/AllenConf/HLT1.py b/configuration/python/AllenConf/HLT1.py
index da99bb613f4..beee26aac44 100644
--- a/configuration/python/AllenConf/HLT1.py
+++ b/configuration/python/AllenConf/HLT1.py
@@ -1568,14 +1568,15 @@ def setup_hlt1_node(enablePhysics=True,
             "AllenWithLumi", [hlt1_node, lumi_with_prefilter],
             NodeLogic.NONLAZY_AND,
             force_order=False)
-
+    """
     if with_fullKF:
-        # Added for testing
+        # Added for testing magnetic field
         hlt1_node = CompositeNode(
             "AllenExtrapolatedStates",
             [hlt1_node, reconstructed_objects["extrapolated_states"]],
             NodeLogic.NONLAZY_AND,
             force_order=True)
+    """
 
     if with_rich:
         hlt1_node = CompositeNode(
-- 
GitLab


From f4a60ef85e5c1f23d85ea924f0b40a038fafbc58 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Fri, 22 Aug 2025 14:53:06 +0200
Subject: [PATCH 09/10] Add magnetic field dumper

---
 Dumpers/BinaryDumpers/CMakeLists.txt          |  1 +
 .../BinaryDumpers/src/DumpMagneticField.cpp   | 36 +++-----
 .../src/DumpMagneticFieldPolarity.cpp         | 86 +++++++++++++++++++
 Rec/Allen/python/Allen/config.py              | 19 ++--
 4 files changed, 111 insertions(+), 31 deletions(-)
 create mode 100644 Dumpers/BinaryDumpers/src/DumpMagneticFieldPolarity.cpp

diff --git a/Dumpers/BinaryDumpers/CMakeLists.txt b/Dumpers/BinaryDumpers/CMakeLists.txt
index 7675512a65e..f9cfc428b34 100644
--- a/Dumpers/BinaryDumpers/CMakeLists.txt
+++ b/Dumpers/BinaryDumpers/CMakeLists.txt
@@ -47,6 +47,7 @@ gaudi_add_module(BinaryDumpersModule
                    src/DumpFTGeometry.cpp
                    src/DumpFTHits.cpp
                    src/DumpMagneticField.cpp
+                   src/DumpMagneticFieldPolarity.cpp
                    src/DumpMuonCommonHits.cpp
                    src/DumpMuonCoords.cpp
                    src/DumpMuonGeometry.cpp
diff --git a/Dumpers/BinaryDumpers/src/DumpMagneticField.cpp b/Dumpers/BinaryDumpers/src/DumpMagneticField.cpp
index 1b60bc08306..833ac1efab8 100644
--- a/Dumpers/BinaryDumpers/src/DumpMagneticField.cpp
+++ b/Dumpers/BinaryDumpers/src/DumpMagneticField.cpp
@@ -1,5 +1,5 @@
 /*****************************************************************************\
-* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      *
+* (c) Copyright 2000-2025 CERN for the benefit of the LHCb Collaboration      *
 *                                                                             *
 * This software is distributed under the terms of the Apache License          *
 * version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
@@ -9,9 +9,6 @@
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
 
-#ifndef DUMPMAGNETICFIELD_H
-#define DUMPMAGNETICFIELD_H 1
-
 #include <tuple>
 #include <vector>
 
@@ -29,23 +26,16 @@
 
 #include "Dumper.h"
 
-/** @class DumpMagneticField
- *  Dump Magnetic Field Polarity.
- *
- *  @author Nabil Garroum
- *  @date   2022-04-21
- *  This Class dumps Data for Magnetic Field polarity using DD4HEP and Gaudi Algorithm
- *  This Class uses a detector description
- *  This Class is basically an instation of a Gaudi algorithm with specific inputs and outputs:
- *  The role of this class is to get data from TES to Allen for the Magnetic Field
- */
 namespace Dumpers {
   struct MagneticField {
-    MagneticField(std::vector<char>& data, const DeMagnet& magField)
+    MagneticField(std::vector<char>& data, const DeMagnet& deMagnet)
     {
+      const auto* grid = deMagnet.fieldGrid();
       DumpUtils::Writer output {};
-      float polarity = magField.isDown() ? -1.f : 1.f;
-      output.write(polarity);
+      output.write(grid->invDXYZ());
+      output.write(grid->sizeXYZ());
+      output.write(grid->minXYZ());
+      output.write(grid->field());
       data = output.buffer();
     }
   };
@@ -70,22 +60,20 @@ DECLARE_COMPONENT(DumpMagneticField)
 // Add the multitransformer call , which keyvalues for Magnetic Field ?
 
 DumpMagneticField::DumpMagneticField(const std::string& name, ISvcLocator* svcLoc) :
-  Dumper(name, svcLoc, {KeyValue {"MagneticFieldLocation", location(name, "Polarity")}})
+  Dumper(name, svcLoc, {KeyValue {"MagneticFieldLocation", location(name, "Magfield")}})
 {}
 
 StatusCode DumpMagneticField::initialize()
 {
   return Dumper::initialize().andThen([&] {
-    register_producer(Allen::NonEventData::MagneticField::id, "polarity", m_data);
+    register_producer(Allen::NonEventData::MagneticField::id, "magfield", m_data);
     addConditionDerivation(
-      {LHCb::Det::Magnet::det_path}, inputLocation<Dumpers::MagneticField>(), [&](DeMagnet const& magField) {
-        auto polarity = Dumpers::MagneticField {m_data, magField};
+      {LHCb::Det::Magnet::det_path}, inputLocation<Dumpers::MagneticField>(), [&](DeMagnet const& deMagnet) {
+        auto magfield = Dumpers::MagneticField {m_data, deMagnet};
         dump();
-        return polarity;
+        return magfield;
       });
   });
 }
 
 void DumpMagneticField::operator()(const Dumpers::MagneticField&) const {}
-
-#endif // DUMPMAGNETICFIELD_H
diff --git a/Dumpers/BinaryDumpers/src/DumpMagneticFieldPolarity.cpp b/Dumpers/BinaryDumpers/src/DumpMagneticFieldPolarity.cpp
new file mode 100644
index 00000000000..e44df4e02da
--- /dev/null
+++ b/Dumpers/BinaryDumpers/src/DumpMagneticFieldPolarity.cpp
@@ -0,0 +1,86 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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 <tuple>
+#include <vector>
+
+// Gaudi
+#include <GaudiAlg/Transformer.h>
+#include <GaudiAlg/FunctionalUtilities.h>
+
+// LHCb
+#include <Magnet/DeMagnet.h>
+#include <DetDesc/GenericConditionAccessorHolder.h>
+
+// Include files
+#include <Dumpers/Utils.h>
+#include <Dumpers/Identifiers.h>
+
+#include "Dumper.h"
+
+/** @class DumpMagneticFieldPolarity
+ *  Dump Magnetic Field Polarity.
+ *
+ *  @author Nabil Garroum
+ *  @date   2022-04-21
+ *  This Class dumps Data for Magnetic Field polarity using DD4HEP and Gaudi Algorithm
+ *  This Class uses a detector description
+ *  This Class is basically an instation of a Gaudi algorithm with specific inputs and outputs:
+ *  The role of this class is to get data from TES to Allen for the Magnetic Field
+ */
+namespace Dumpers {
+  struct MagneticFieldPolarity {
+    MagneticFieldPolarity(std::vector<char>& data, const DeMagnet& magField)
+    {
+      DumpUtils::Writer output {};
+      float polarity = magField.isDown() ? -1.f : 1.f;
+      output.write(polarity);
+      data = output.buffer();
+    }
+  };
+} // namespace Dumpers
+
+class DumpMagneticFieldPolarity final : public Allen::Dumpers::Dumper<
+                                          void(Dumpers::MagneticFieldPolarity const&),
+                                          LHCb::Algorithm::Traits::usesConditions<Dumpers::MagneticFieldPolarity>> {
+public:
+  DumpMagneticFieldPolarity(const std::string& name, ISvcLocator* svcLoc);
+
+  void operator()(const Dumpers::MagneticFieldPolarity&) const override;
+
+  StatusCode initialize() override;
+
+private:
+  std::vector<char> m_data;
+};
+
+DECLARE_COMPONENT(DumpMagneticFieldPolarity)
+
+// Add the multitransformer call , which keyvalues for Magnetic Field ?
+
+DumpMagneticFieldPolarity::DumpMagneticFieldPolarity(const std::string& name, ISvcLocator* svcLoc) :
+  Dumper(name, svcLoc, {KeyValue {"MagneticFieldPolarityLocation", location(name, "Polarity")}})
+{}
+
+StatusCode DumpMagneticFieldPolarity::initialize()
+{
+  return Dumper::initialize().andThen([&] {
+    register_producer(Allen::NonEventData::MagneticFieldPolarity::id, "polarity", m_data);
+    addConditionDerivation(
+      {LHCb::Det::Magnet::det_path}, inputLocation<Dumpers::MagneticFieldPolarity>(), [&](DeMagnet const& magField) {
+        auto polarity = Dumpers::MagneticFieldPolarity {m_data, magField};
+        dump();
+        return polarity;
+      });
+  });
+}
+
+void DumpMagneticFieldPolarity::operator()(const Dumpers::MagneticFieldPolarity&) const {}
diff --git a/Rec/Allen/python/Allen/config.py b/Rec/Allen/python/Allen/config.py
index ed58528b995..572b7ea16ea 100755
--- a/Rec/Allen/python/Allen/config.py
+++ b/Rec/Allen/python/Allen/config.py
@@ -17,9 +17,10 @@ from PyConf import configurable
 from PyConf.control_flow import CompositeNode, NodeLogic
 from PyConf.application import ApplicationOptions, configure_input, configure
 from PyConf.Algorithms import (
-    DumpBeamline, DumpCaloGeometry, DumpMagneticField, DumpVPGeometry,
-    DumpCrossingAngles, DumpFTGeometry, DumpUTGeometry, DumpUTLookupTables,
-    DumpMuonGeometry, DumpMuonTable, AllenODINProducer, DumpRichPDMDBMapping,
+    DumpBeamline, DumpCaloGeometry, DumpMagneticField,
+    DumpMagneticFieldPolarity, DumpVPGeometry, DumpCrossingAngles,
+    DumpFTGeometry, DumpUTGeometry, DumpUTLookupTables, DumpMuonGeometry,
+    DumpMuonTable, AllenODINProducer, DumpRichPDMDBMapping,
     DumpRichCableMapping)
 from DDDB.CheckDD4Hep import UseDD4Hep
 from PyConf.reading import get_generator_BeamParameters
@@ -109,8 +110,10 @@ def setup_allen_non_event_data_service(allen_event_loop=False,
                     'ut_tables')],
             'ECal': [(DumpCaloGeometry, 'DeviceCaloGeometry', {},
                       'ecal_geometry')],
-            'Magnet': [(DumpMagneticField, 'DeviceMagneticField', {},
-                        'polarity')],
+            'Magnet':
+            [(DumpMagneticField, 'DeviceMagneticField', {}, 'magfield'),
+             (DumpMagneticFieldPolarity, 'DeviceMagneticFieldPolarity', {},
+              'polarity')],
             'FTCluster': [(DumpFTGeometry, 'DeviceFTGeometry', {},
                            'scifi_geometry')],
             'Muon': [(DumpMuonGeometry, 'DeviceMuonGeometry', {},
@@ -132,8 +135,10 @@ def setup_allen_non_event_data_service(allen_event_loop=False,
                     'ut_tables')],
             'ECal': [(DumpCaloGeometry, 'DeviceCaloGeometry', {},
                       'ecal_geometry')],
-            'Magnet': [(DumpMagneticField, 'DeviceMagneticField', {},
-                        'polarity')],
+            'Magnet':
+            [(DumpMagneticField, 'DeviceMagneticField', {}, 'magfield'),
+             (DumpMagneticFieldPolarity, 'DeviceMagneticFieldPolarity', {},
+              'polarity')],
             'FTCluster': [(DumpFTGeometry, 'DeviceFTGeometry', {},
                            'scifi_geometry')],
             'Muon': [(DumpMuonGeometry, 'DeviceMuonGeometry', {},
-- 
GitLab


From d8b10ab07d3402f603ea5ab04a96d2d47cca62e8 Mon Sep 17 00:00:00 2001
From: Arthur Hennequin <arthur.hennequin@cern.ch>
Date: Tue, 14 Oct 2025 11:17:58 +0200
Subject: [PATCH 10/10] Fix qop initialization in State operator+

---
 .../python/AllenConf/secondary_vertex_reconstruction.py         | 2 ++
 device/kalman/ParKalman/include/ExtrapolatorCommon.cuh          | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/configuration/python/AllenConf/secondary_vertex_reconstruction.py b/configuration/python/AllenConf/secondary_vertex_reconstruction.py
index bad112d0b87..f98352fefc2 100644
--- a/configuration/python/AllenConf/secondary_vertex_reconstruction.py
+++ b/configuration/python/AllenConf/secondary_vertex_reconstruction.py
@@ -48,6 +48,7 @@ def make_kalman_long(long_tracks, pvs, is_muon_result,
         "dev_kalman_states_view": kalman.dev_kalman_states_view_t
     }
 
+
 def make_extrapolated_states(states, n_states, step_dz=100.0, n_steps=100):
     extrapolator = make_algorithm(
         extrapolate_states_t,
@@ -58,6 +59,7 @@ def make_extrapolated_states(states, n_states, step_dz=100.0, n_steps=100):
         step_dz=step_dz)
     return extrapolator.dev_states_t
 
+
 def make_kalman_velo_only(long_tracks,
                           pvs,
                           is_muon_result,
diff --git a/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh b/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
index 6ab468007cc..251c701e515 100644
--- a/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
+++ b/device/kalman/ParKalman/include/ExtrapolatorCommon.cuh
@@ -55,7 +55,7 @@ namespace Extrapolators {
 
     __device__ friend State operator+(const State& a, const Derivative& b)
     {
-      return {a.x + b.dx, a.y + b.dy, a.z + b.dz, a.tx + b.dtx, a.ty + b.dty};
+      return {a.x + b.dx, a.y + b.dy, a.z + b.dz, a.tx + b.dtx, a.ty + b.dty, a.qop};
     }
   };
 } // namespace Extrapolators
-- 
GitLab