diff --git a/CMakeLists.txt b/CMakeLists.txt
index ffeb915be4f7434edbf30e29b3b0d88ebc7720d8..06ac84cf1956cfd632ab44f45b9740f3ffa21c8f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -215,7 +215,7 @@ include_directories(cuda/SciFi/looking_forward_sbt/extend_tracks_uv/include)
 include_directories(cuda/SciFi/looking_forward_sbt/quality_filter/include)
 include_directories(cuda/SciFi/looking_forward_sbt/quality_filter_x/include)
 include_directories(cuda/SciFi/looking_forward_sbt/search_uv_windows/include)
-include_directories(cuda/SciFi/PrForward/include)
+include_directories(cuda/SciFi/classifiers/include)
 include_directories(cuda/SciFi/consolidate/include)
 include_directories(cuda/muon/common/include)
 include_directories(cuda/utils/prefix_sum/include)
@@ -227,7 +227,6 @@ include_directories(checker/tracking/include)
 include_directories(checker/pv/include)
 include_directories(checker/selections/include)
 include_directories(stream/sequence/include)
-include_directories(x86/SciFi/PrForward/include)
 include_directories(x86/SciFi/LookingForward/include)
 include_directories(x86/SciFi/MomentumForward/include)
 include_directories(cuda/UT/UTDecoding/include)
diff --git a/checker/tracking/CMakeLists.txt b/checker/tracking/CMakeLists.txt
index fa931478c6ae8b4ce968a69ef99f6a51acfcfb10..7ff1dcecfc3d0fa08e13846443e0482cc73dfbf7 100644
--- a/checker/tracking/CMakeLists.txt
+++ b/checker/tracking/CMakeLists.txt
@@ -9,7 +9,6 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/UT/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/SciFi/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/kalman/ParKalman/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/beamlinePV/include)
diff --git a/cuda/SciFi/CMakeLists.txt b/cuda/SciFi/CMakeLists.txt
index a387c9b81560bcc73489e1bd4b46af874e0e75ae..fcffc61c06ce11283d2e535c1f9265ca3a19cbd6 100644
--- a/cuda/SciFi/CMakeLists.txt
+++ b/cuda/SciFi/CMakeLists.txt
@@ -19,7 +19,8 @@ file(GLOB scifi_lf_quality_filter "looking_forward_sbt/quality_filter/src/*cu")
 file(GLOB scifi_lf_quality_filter_x "looking_forward_sbt/quality_filter_x/src/*cu")
 file(GLOB scifi_lf_search_uv_windows "looking_forward_sbt/search_uv_windows/src/*cu")
 file(GLOB scifi_lf_fit "looking_forward_sbt/fit/src/*cu")
-file(GLOB scifi_prforward "PrForward/src/*cu")
+file(GLOB scifi_utils "utils/src/*cu")
+file(GLOB scifi_classifiers "classifiers/src/*cu")
 file(GLOB scifi_consolidate "consolidate/src/*cu")
 
 include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
@@ -61,12 +62,14 @@ include_directories(looking_forward_sbt/quality_filter/include)
 include_directories(looking_forward_sbt/quality_filter_x/include)
 include_directories(looking_forward_sbt/search_uv_windows/include)
 include_directories(looking_forward_sbt/fit/include)
-include_directories(PrForward/include)
+include_directories(utils/include)
+include_directories(classifiers/include)
 include_directories(consolidate/include)
 
 add_library(SciFi STATIC
   ${scifi_common}
   ${scifi_preprocessing}
+  ${scifi_utils}
   ${scifi_lf_common}
   ${scifi_lf_calculate_first_layer_window}
   ${scifi_lf_calculate_second_layer_window}
@@ -86,7 +89,7 @@ add_library(SciFi STATIC
   ${scifi_lf_quality_filter_x}
   ${scifi_lf_search_uv_windows}
   ${scifi_lf_fit}
-  ${scifi_prforward}
+  ${scifi_classifiers}
   ${scifi_consolidate}
 )
 
diff --git a/cuda/SciFi/PrForward/.gitignore b/cuda/SciFi/PrForward/.gitignore
deleted file mode 100644
index 378eac25d311703f3f2cd456d8036da525cd0366..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/.gitignore
+++ /dev/null
@@ -1 +0,0 @@
-build
diff --git a/cuda/SciFi/PrForward/README.md b/cuda/SciFi/PrForward/README.md
deleted file mode 100644
index cc8a2c28798afad75d8007a0ba8fe7eb7450dd9f..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/README.md
+++ /dev/null
@@ -1 +0,0 @@
-# Implementation of the forward based on Rec v23r3
diff --git a/cuda/SciFi/PrForward/include/FindStereoHits.cuh b/cuda/SciFi/PrForward/include/FindStereoHits.cuh
deleted file mode 100644
index c14347c56570a8a7c943353e1ca17425b7cf92b3..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/FindStereoHits.cuh
+++ /dev/null
@@ -1,49 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <array>
-#include <vector>
-#include <algorithm>
-#include <fstream>
-#include "SciFiDefinitions.cuh"
-#include "SciFiEventModel.cuh"
-#include "TrackUtils.cuh"
-#include "HitUtils.cuh"
-
-/**
-   Functions related to selecting hits on the uv planes,
-   which match to the VeloUT input track
- */
-
-__host__ __device__ void collectStereoHits(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track& track,
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars,
-  const SciFi::Tracking::Arrays* constArrays,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits);
-
-__host__ __device__ bool selectStereoHits(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track& track,
-  const SciFi::Tracking::Arrays* constArrays,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits,
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars_cur);
-
-__host__ __device__ bool addHitsOnEmptyStereoLayers(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track& track,
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits,
-  const SciFi::Tracking::Arrays* constArrays,
-  PlaneCounter& planeCounter,
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars_cur);
diff --git a/cuda/SciFi/PrForward/include/FindXHits.cuh b/cuda/SciFi/PrForward/include/FindXHits.cuh
deleted file mode 100644
index 136f83fd08c8a1c92d23120f81523a7dd64fee0b..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/FindXHits.cuh
+++ /dev/null
@@ -1,120 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <array>
-#include <vector>
-#include <algorithm>
-#include <fstream>
-#include "SciFiDefinitions.cuh"
-#include "PrForwardConstants.cuh"
-#include "UTDefinitions.cuh"
-#include "TrackUtils.cuh"
-#include "HitUtils.cuh"
-#include "LinearFitting.cuh"
-#include "ReferencePlaneProjection.cuh"
-#include "SciFiEventModel.cuh"
-
-#include "LookingForwardUtils.h"
-
-/**
-   Functions related to selecting hits on the x planes,
-   which match to the VeloUT input track
- */
-__host__ void collectAllXHits_proto_p(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state,
-  const MiniState& UT_state,
-  const float qOverP,
-  int side,
-  std::array<int, 2 * 6>& windows_x,
-  std::array<int, 2 * 6>& windows_uv,
-  std::array<float, 4 * 6>& parameters_uv,
-  const SciFiWindowsParams& window_params,
-  const std::array<int, 12> true_scifi_indices_per_layer);
-
-__host__ void x_limits_from_dxRef(
-  const SciFi::Tracking::Arrays* constArrays,
-  const MiniState& velo_state,
-  const float InvPz,
-  const float p,
-  const float tx2,
-  const float ty2,
-  const bool wSignTreatment,
-  float& xBoundOnRef,
-  float& xBoundOnRefWS);
-
-__host__ void collectAllXHits_proto(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& UT_state,
-  const float qOverP,
-  int side,
-  std::array<int, 2 * 6>& windows_x,
-  std::array<int, 2 * 6>& windows_uv,
-  std::array<float, 4 * 6>& parameters_uv);
-
-__host__ __device__ void collectAllXHits(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  int allXHits[SciFi::Tracking::max_x_hits],
-  int& n_x_hits,
-  float coordX[SciFi::Tracking::max_x_hits],
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state,
-  const float qop,
-  int side);
-
-__host__ __device__ void improveXCluster(
-  int& it2,
-  const int it1,
-  const int itEnd,
-  const int n_x_hits,
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const float xWindow,
-  const SciFi::Tracking::HitSearchCuts& pars,
-  PlaneCounter& planeCounter,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits);
-
-__host__ __device__ void selectXCandidates(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  int allXHits[SciFi::Tracking::max_x_hits],
-  int& n_x_hits,
-  bool usedHits[SciFi::Tracking::max_x_hits],
-  float coordX[SciFi::Tracking::max_x_hits],
-  SciFi::Tracking::Track candidate_tracks[SciFi::Constants::max_tracks],
-  int& n_candidate_tracks,
-  const float zRef_track,
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars,
-  const SciFi::Tracking::Arrays* constArrays,
-  int side,
-  const bool secondLoop);
-
-__host__ __device__ bool addHitsOnEmptyXLayers(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  bool fullFit,
-  int coordToFit[SciFi::Tracking::max_coordToFit],
-  int& n_coordToFit,
-  const SciFi::Tracking::Arrays* constArrays,
-  PlaneCounter& planeCounter,
-  SciFi::Tracking::HitSearchCuts& pars_cur,
-  int side);
diff --git a/cuda/SciFi/PrForward/include/HitUtils.cuh b/cuda/SciFi/PrForward/include/HitUtils.cuh
deleted file mode 100644
index dad26033d92ea5b485bdd4b50d16696aa333784e..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/HitUtils.cuh
+++ /dev/null
@@ -1,220 +0,0 @@
-#pragma once
-
-#include "SciFiDefinitions.cuh"
-#include "PrForwardConstants.cuh"
-#include "SciFiEventModel.cuh"
-
-/**
-   Helper functions related to properties of hits on planes
- */
-
-// Helper used to keep track of how many x / stereo hits per lane have
-// been added to a candidate track
-struct PlaneCounter {
-  int planeList[SciFi::Constants::n_layers] = {0};
-  unsigned int nbDifferent = 0;
-
-  __host__ __device__ inline void addHit(int plane)
-  {
-    assert(plane < SciFi::Constants::n_layers);
-    nbDifferent += (int) ((planeList[plane] += 1) == 1);
-  }
-
-  __host__ __device__ inline void removeHit(int plane)
-  {
-    assert(plane < SciFi::Constants::n_layers);
-    nbDifferent -= ((int) ((planeList[plane] -= 1) == 0));
-  }
-
-  __host__ __device__ inline int nbInPlane(int plane) const
-  {
-    assert(plane < SciFi::Constants::n_layers);
-    return planeList[plane];
-  }
-
-  __host__ __device__ inline int nbSingle() const
-  {
-    int single = 0;
-    for (int i = 0; i < SciFi::Constants::n_layers; ++i) {
-      single += planeList[i] == 1 ? 1 : 0;
-    }
-    return single;
-  }
-
-  __host__ __device__ inline void clear()
-  {
-    nbDifferent = 0;
-    for (int i = 0; i < SciFi::Constants::n_layers; ++i) {
-      planeList[i] = 0;
-    }
-  }
-};
-
-__host__ __device__ void countPlanesOfXHits(
-  PlaneCounter& planeCounter,
-  const int it1,
-  const int it2,
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits);
-
-__host__ __device__ void countUnusedXHitsOnPlanes(
-  PlaneCounter& lplaneCounter,
-  const int itWindowStart,
-  const int itWindowEnd,
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits);
-
-__host__ __device__ void addXHitsForCandidateWithTooFewPlanes(
-  int& itWindowStart,
-  int& itWindowEnd,
-  const int it2,
-  const int itEnd,
-  float& minInterval,
-  PlaneCounter& lplaneCounter,
-  const int nPlanes,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  int& best,
-  int& bestEnd,
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits);
-
-__host__ __device__ void collectXHitsToFit(
-  const int it1,
-  const int it2,
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  bool usedHits[SciFi::Tracking::max_x_hits],
-  int coordToFit[SciFi::Tracking::max_x_hits],
-  int& n_coordToFit,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  float& xAtRef);
-
-__host__ __device__ int findBestXHitOnEmptyLayer(
-  const int itEnd,
-  const int itH,
-  const SciFi::Hits& scifi_hits,
-  const float maxX,
-  const float xPred);
-
-template<int N>
-__host__ __device__ void sortHitsByKey(float* keys, int n, int* hits)
-{
-  // find permutations
-  uint permutations[N];
-  assert(n <= N);
-  for (int i = 0; i < n; ++i) {
-    uint position = 0;
-    for (int j = 0; j < n; ++j) {
-      // sort keys in ascending order
-      int sort_result = -1;
-      if (keys[i] > keys[j]) sort_result = 1;
-      if (keys[i] == keys[j]) sort_result = 0;
-      position += sort_result > 0 || (sort_result == 0 && i > j);
-    }
-    permutations[position] = i;
-  }
-
-  // apply permutations, store hits in temporary container
-  int hits_tmp[N];
-  float keys_tmp[N];
-  for (int i = 0; i < n; ++i) {
-    const int index = permutations[i];
-    hits_tmp[i] = hits[index];
-    keys_tmp[i] = keys[index];
-  }
-
-  // copy hits back to original container
-  for (int i = 0; i < n; ++i) {
-    hits[i] = hits_tmp[i];
-    keys[i] = keys_tmp[i];
-  }
-}
-
-// check that val is within [min, max]
-__host__ __device__ inline bool isInside(float val, const float min, const float max)
-{
-  return (val > min) && (val < max);
-}
-
-// get lowest index where range[index] > value, within [start,end] of range
-__host__ __device__ inline int getLowerBound(float range[], float value, int start, int end)
-{
-  int i = start;
-  for (; i < end; i++) {
-    if (range[i] > value) break;
-  }
-  return i;
-}
-
-// match stereo hits to x hits
-__host__ __device__ bool matchStereoHit(
-  const int itUV1,
-  const int uv_zone_offset_end,
-  const SciFi::Hits& scifi_hits,
-  const float xMinUV,
-  const float xMaxUV);
-
-__host__ __device__ bool matchStereoHitWithTriangle(
-  const int itUV2,
-  const int triangle_zone_offset_end,
-  const float yInZone,
-  const SciFi::Hits& scifi_hits,
-  const float xMinUV,
-  const float xMaxUV,
-  const int side);
-
-__host__ __device__ void removeOutlier(
-  const SciFi::Hits& scifi_hits,
-  PlaneCounter& planeCounter,
-  int* coordToFit,
-  int& n_coordToFit,
-  const int worst);
-
-__host__ __device__ void findStereoHitsWithinXTol(
-  const int itBegin,
-  const int itEnd,
-  const SciFi::Hits& scifi_hits,
-  const float yZone,
-  const float xPred,
-  const float dxTol,
-  const bool triangleSearch,
-  const float dxDySign,
-  int& n_stereoHits,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits]);
-
-__host__ __device__ void findStereoHitClusterByDx(
-  PlaneCounter& planeCounter,
-  int& endRange,
-  const SciFi::Tracking::HitSearchCuts& pars,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  const int n_stereoHits,
-  const SciFi::Hits& scifi_hits,
-  float& sumCoord,
-  int& first_hit);
-
-__host__ __device__ void cleanStereoHitCluster(
-  int& beginRange,
-  int& endRange,
-  const int n_stereoHits,
-  const int stereoHits[SciFi::Tracking::max_stereo_hits],
-  const float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  float& sumCoord,
-  PlaneCounter& planeCounter,
-  const SciFi::Hits& scifi_hits);
-
-__host__ __device__ int findBestStereoHitOnEmptyLayer(
-  const int itBegin,
-  const int itEnd,
-  const SciFi::Hits& scifi_hits,
-  const float yZone,
-  const float xPred,
-  const float dxTol,
-  const bool triangleSearch);
diff --git a/cuda/SciFi/PrForward/include/LinearFitting.cuh b/cuda/SciFi/PrForward/include/LinearFitting.cuh
deleted file mode 100644
index 486f7b509059d9810c72920d2206137b337ff268..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/LinearFitting.cuh
+++ /dev/null
@@ -1,77 +0,0 @@
-#pragma once
-
-#include "SciFiDefinitions.cuh"
-#include "TrackUtils.cuh"
-#include "HitUtils.cuh"
-#include "SciFiEventModel.cuh"
-
-#include <cmath>
-
-namespace SciFi {
-  namespace Tracking {
-
-    struct LineFitterPars {
-      float m_z0 = 0.;
-      float m_c0 = 0.;
-      float m_tc = 0.;
-
-      float m_s0 = 0.;
-      float m_sz = 0.;
-      float m_sz2 = 0.;
-      float m_sc = 0.;
-      float m_scz = 0.;
-    };
-  } // namespace Tracking
-} // namespace SciFi
-
-__host__ __device__ void incrementLineFitParameters(
-  SciFi::Tracking::LineFitterPars& parameters,
-  const SciFi::Hits& scifi_hits,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int it);
-
-__host__ __device__ void fitHitsFromSingleHitPlanes(
-  const int it1,
-  const int it2,
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int n_x_hits,
-  const PlaneCounter planeCounter,
-  SciFi::Tracking::LineFitterPars& lineFitParameters,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  int otherHits[SciFi::Constants::n_layers][SciFi::Tracking::max_other_hits],
-  int nOtherHits[SciFi::Constants::n_layers]);
-
-__host__ __device__ void addAndFitHitsFromMultipleHitPlanes(
-  const int nOtherHits[SciFi::Constants::n_layers],
-  SciFi::Tracking::LineFitterPars& lineFitParameters,
-  const SciFi::Hits& scifi_hits,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int otherHits[SciFi::Constants::n_layers][SciFi::Tracking::max_other_hits]);
-
-__host__ __device__ float getLineFitDistance(
-  SciFi::Tracking::LineFitterPars& parameters,
-  const SciFi::Hits& scifi_hits,
-  float coordX[SciFi::Tracking::max_x_hits],
-  int allXHits[SciFi::Tracking::max_x_hits],
-  int it);
-
-__host__ __device__ float getLineFitChi2(
-  SciFi::Tracking::LineFitterPars& parameters,
-  const SciFi::Hits& scifi_hits,
-  float coordX[SciFi::Tracking::max_x_hits],
-  int allXHits[SciFi::Tracking::max_x_hits],
-  int it);
-
-__host__ __device__ void solveLineFit(SciFi::Tracking::LineFitterPars& parameters);
-
-__host__ __device__ void fastLinearFit(
-  const SciFi::Hits& scifi_hits,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  int coordToFit[SciFi::Tracking::max_coordToFit],
-  int& n_coordToFit,
-  PlaneCounter planeCounter,
-  SciFi::Tracking::HitSearchCuts& pars_cur);
diff --git a/cuda/SciFi/PrForward/include/ParabolaFitting.cuh b/cuda/SciFi/PrForward/include/ParabolaFitting.cuh
deleted file mode 100644
index 0cfb8d66da7fa8ea281d66977eff4327ec9ce10d..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/ParabolaFitting.cuh
+++ /dev/null
@@ -1,18 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <array>
-#include <vector>
-#include <algorithm>
-#include <fstream>
-
-#include "SciFiDefinitions.cuh"
-#include "TrackUtils.cuh"
-#include "SciFiEventModel.cuh"
-
-__host__ __device__ int fitParabola(
-  int* coordToFit,
-  const int n_coordToFit,
-  const SciFi::Hits& scifi_hits,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  const bool xFit);
diff --git a/cuda/SciFi/PrForward/include/PrForward.cuh b/cuda/SciFi/PrForward/include/PrForward.cuh
deleted file mode 100644
index 376a16fe8fed90d80f15fa59a5efebca11228635..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/PrForward.cuh
+++ /dev/null
@@ -1,59 +0,0 @@
-#pragma once
-
-#include "PrForwardTools.cuh"
-#include "Handler.cuh"
-// #include "ArgumentsCommon.cuh"
-#include "ArgumentsVelo.cuh"
-#include "ArgumentsUT.cuh"
-#include "ArgumentsSciFi.cuh"
-
-/** @class PrForward PrForward.h
- *
- *  - InputTracksName: Input location for VeloUT tracks
- *  - OutputTracksName: Output location for Forward tracks
- *  Based on code written by
- *  2012-03-20 : Olivier Callot
- *  2013-03-15 : Thomas Nikodem
- *  2015-02-13 : Sevda Esen [additional search in the triangles by Marian Stahl]
- *  2016-03-09 : Thomas Nikodem [complete restructuring]
- *  2018-08    : Vava Gligorov [extract code from Rec, make compile within GPU framework
- *  2018-09    : Dorothea vom Bruch [convert to CUDA, runs on GPU]
- */
-
-__global__ void scifi_pr_forward(
-  uint32_t* dev_scifi_hits,
-  const uint32_t* dev_scifi_hit_count,
-  const int* dev_atomics_velo,
-  const uint* dev_velo_track_hit_number,
-  const char* dev_velo_states,
-  const int* dev_atomics_ut,
-  const char* dev_ut_track_hits,
-  const uint* dev_ut_track_hit_number,
-  const float* dev_ut_qop,
-  const uint* dev_ut_track_velo_indices,
-  SciFi::TrackHits* dev_scifi_tracks,
-  int* dev_atomics_scifi,
-  const SciFi::Tracking::TMVA* dev_tmva1,
-  const SciFi::Tracking::TMVA* dev_tmva2,
-  const SciFi::Tracking::Arrays* dev_constArrays,
-  const float* dev_magnet_polarity,
-  const char* dev_scifi_geometry,
-  const float* dev_inv_clus_res);
-
-ALGORITHM(
-  scifi_pr_forward,
-  scifi_pr_forward_t,
-  ARGUMENTS(
-    dev_scifi_hits,
-    dev_scifi_hit_count,
-    dev_atomics_velo,
-    dev_velo_track_hit_number,
-    dev_velo_states,
-    dev_atomics_ut,
-    dev_ut_track_hits,
-    dev_ut_track_hit_number,
-    dev_ut_qop,
-    dev_ut_track_velo_indices,
-    dev_scifi_tracks,
-    dev_scifi_selected_track_indices,
-    dev_atomics_scifi))
diff --git a/cuda/SciFi/PrForward/include/PrForwardConstants.cuh b/cuda/SciFi/PrForward/include/PrForwardConstants.cuh
deleted file mode 100644
index 4a20d34ba4f4459317685350faec1813c16d1318..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/PrForwardConstants.cuh
+++ /dev/null
@@ -1,236 +0,0 @@
-#pragma once
-
-/**
-   Contains constants needed for the forward tracking
-   - cut values
-   - geometry descriptions
-   - parameterizations
-
-   12/09/2018: cut values are those defined in:
-   https://gitlab.cern.ch/lhcb/Rec/blob/master/Tf/TrackSys/python/TrackSys/Configuration.py
-   https://gitlab.cern.ch/lhcb/Rec/blob/master/Tf/TrackSys/python/TrackSys/RecoUpgradeTracking.py
-
-   for the RecoFastTrackingStage, using the default values of ConfigHLT1 (master branch of Rec)
-
- */
-
-#include "VeloEventModel.cuh"
-#include "SystemOfUnits.h"
-#include "SciFiDefinitions.cuh"
-#include <cassert>
-
-namespace SciFi {
-
-  namespace Tracking {
-
-    constexpr int max_candidate_tracks = 5;   // max # of candidate tracks from x hits only
-    constexpr int max_tracks_second_loop = 5; // same as above, but for second loop
-    constexpr int max_selected_tracks = max_candidate_tracks + max_tracks_second_loop;
-    constexpr int max_x_hits = 500;     // max # of hits in all x layers
-    constexpr int max_other_hits = 5;   // max # of hits from x planes with more than 1 hit
-    constexpr int max_stereo_hits = 25; // max # of hits in all stereo layers
-    constexpr int max_coordToFit = 15;  // only for x layers
-    constexpr int max_scifi_hits = 20;  // for x and u/v layers
-
-    constexpr int nTrackParams = 9;
-
-    constexpr int TMVA_Nvars = 7;
-    constexpr int TMVA_Nlayers = 5;
-
-    // Formerly PrParameters
-    struct HitSearchCuts {
-      __host__ __device__ HitSearchCuts(
-        unsigned int minXHits_,
-        float maxXWindow_,
-        float maxXWindowSlope_,
-        float maxXGap_,
-        unsigned int minStereoHits_) :
-        minXHits {minXHits_},
-        maxXWindow {maxXWindow_}, maxXWindowSlope {maxXWindowSlope_}, maxXGap {maxXGap_}, minStereoHits {minStereoHits_}
-      {}
-      const unsigned int minXHits;
-      const float maxXWindow;
-      const float maxXWindowSlope;
-      const float maxXGap;
-      unsigned int minStereoHits;
-    };
-
-    // dump a bunch of options here
-    constexpr float deltaQuality = 0.1; // Difference in quality btw two tracks which share hits when clone killing
-    constexpr float cloneFraction =
-      0.4; // The fraction of shared SciFi hits btw two tracks to trigger the clone killing
-
-    constexpr float yTolUVSearch = 11. * Gaudi::Units::mm;
-    constexpr float tolY = 5. * Gaudi::Units::mm;
-    constexpr float tolYSlope = 0.002 * Gaudi::Units::mm;
-    constexpr float maxChi2LinearFit = 100.;
-    constexpr float maxChi2XProjection = 15.;
-    constexpr float maxChi2PerDoF = 7.;
-
-    constexpr float tolYMag = 10. * Gaudi::Units::mm;
-    constexpr float tolYMagSlope = 0.015;
-    constexpr float minYGap = 0.4 * Gaudi::Units::mm;
-
-    constexpr unsigned int minTotalHits = 10;
-    constexpr float maxChi2StereoLinear = 60.;
-    constexpr float maxChi2Stereo = 4.5;
-
-    // first loop Hough Cluster search
-    constexpr unsigned int minXHits = 5;
-    constexpr float maxXWindow = 1. * Gaudi::Units::mm; // 1.2 * Gaudi::Units::mm  ;
-    constexpr float maxXWindowSlope = 0.002 * Gaudi::Units::mm;
-    constexpr float maxXGap = 1. * Gaudi::Units::mm; // 1.2 * Gaudi::Units::mm  ;
-    constexpr unsigned int minSingleHits = 2;
-
-    // second loop Hough Cluster search
-    constexpr bool secondLoop = true;
-    constexpr unsigned int minXHits_2nd = 4;
-    constexpr float maxXWindow_2nd = 1.5 * Gaudi::Units::mm;
-    constexpr float maxXWindowSlope_2nd = 0.002 * Gaudi::Units::mm;
-    constexpr float maxXGap_2nd = 0.5 * Gaudi::Units::mm;
-
-    // collectX search
-    constexpr float minPt = 400 * Gaudi::Units::MeV; // 500 * Gaudi::Units::MeV ;
-    // stereo hit matching
-    constexpr float tolYCollectX = 3.5 * Gaudi::Units::mm;        // 4.1* Gaudi::Units::mm ;
-    constexpr float tolYSlopeCollectX = 0.001 * Gaudi::Units::mm; // 0.0018 * Gaudi::Units::mm ;
-    constexpr float tolYTriangleSearch = 20.f;
-    // veloUT momentum estimate
-    constexpr bool useMomentumEstimate = true;
-    constexpr bool useWrongSignWindow = true;
-    constexpr float wrongSignPT = 2000. * Gaudi::Units::MeV;
-    // Track Quality NN
-    constexpr float maxQuality = 0.9;
-    constexpr float deltaQuality_NN = 0.1;
-
-    // parameterizations
-    constexpr float byParams = -0.667996;
-    constexpr float cyParams = -3.68424e-05;
-
-    // z Reference plane
-    constexpr float zReference = 8520. * Gaudi::Units::mm; // in T2
-    constexpr float zRefInv = 1.f / zReference;
-
-    // TODO: CHECK THESE VALUES USING FRAMEWORK
-    constexpr float xLim_Max = 3300.;
-    constexpr float yLim_Max = 2500.;
-    constexpr float xLim_Min = -3300.;
-    constexpr float yLim_Min = -25.;
-
-    // TO BE READ FROM XML EVENTUALLY
-    // constexpr float magscalefactor = -1;
-    constexpr int zoneoffsetpar = 6;
-
-    struct Arrays {
-      // Returns whether the current layer is an X plane
-      const bool is_x_plane[12] {1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1};
-
-      // the Magnet Parametrization
-      // parameterized in offset [0], (slope difference due to kick)^2 [1],
-      // tx^2 [2], ty^2 [3]
-      const float zMagnetParams[4] = {5212.38, 406.609, -1102.35, -498.039};
-
-      // more Parametrizations
-      const float xParams[2] = {18.6195, -5.55793};
-
-      // momentum Parametrization
-      const float momentumParams[6] = {1.21014, 0.637339, -0.200292, 0.632298, 3.23793, -27.0259};
-
-      // covariance values
-      const float covarianceValues[5] = {4.0, 400.0, 4.e-6, 1.e-4, 0.1};
-
-      // definition of zones
-      // access upper with offset of 6
-      const int zoneoffsetpar = 6;
-      const int xZones[12] = {0, 6, 8, 14, 16, 22, 1, 7, 9, 15, 17, 23};
-      const int uvZones[12] = {2, 4, 10, 12, 18, 20, 3, 5, 11, 13, 19, 21};
-
-      // ASSORTED GEOMETRY VALUES, eventually read this from some xml
-      const float xZone_zPos[6] = {7826., 8036., 8508., 8718., 9193., 9403.};
-      const float uvZone_zPos[12] =
-        {7896., 7966., 8578., 8648., 9263., 9333., 7896., 7966., 8578., 8648., 9263., 9333.};
-      const float uvZone_dxdy[12] = {0.0874892,
-                                     -0.0874892,
-                                     0.0874892,
-                                     -0.0874892,
-                                     0.0874892,
-                                     -0.0874892,
-                                     0.0874892,
-                                     -0.0874892,
-                                     0.0874892,
-                                     -0.0874892,
-                                     0.0874892,
-                                     -0.0874892};
-      const float Zone_dzdy[24] = {0.0036010};
-
-      // this is used by looking_forward_sbt maybe this is not the right place to put it
-      const float uv_dx[6] = {1.6739478541449213,
-                              1.6738495069872612,
-                              1.935683825160498,
-                              1.9529279746403518,
-                              2.246931985749485,
-                              2.2797556995480273};
-    };
-
-    // parameters for extrapolating from EndVelo to ZReference
-    constexpr float xExtParams[6] = {4.08934e+06f, 6.31187e+08f, 131.999f, -1433.64f, -325.055f, 3173.52f};
-    // Params for momentum dependent search window estimate
-    // upper window to include 98% of hits(can't be too greedy
-    // here or the window size would explode)
-    constexpr float pUp[3] = {1.46244e+02f, 5.15348e+02f, -4.17237e-05f};
-    // lower window, the same to include 98% of hits
-    constexpr float pLo[3] = {5.00000e+01f, 9.61409e+02f, -1.31317e-04f};
-
-    // Track object used for finding tracks, not the final container for storing the tracks
-    struct Track {
-      int hit_indices[max_scifi_hits];
-      float qop;
-      int hitsNum = 0;
-      float quality;
-      float chi2;
-      // [0]: xRef
-      // [1]: (xRef-xMag)/(zRef-zMag)
-      // [2]: xParams[0] * dSlope
-      // [3]: xParams[1] * dSlope
-      // [4]: y param
-      // [5]: y param
-      // [6]: y param
-      // [7]: chi2
-      // [8]: nDoF
-      float trackParams[SciFi::Tracking::nTrackParams];
-
-      __host__ __device__ void addHit(int hit)
-      {
-        assert(hitsNum < max_scifi_hits - 1);
-        hit_indices[hitsNum++] = hit;
-      }
-
-      __host__ __device__ void set_qop(float _qop) { qop = _qop; }
-
-      __host__ __device__ float x(const float z) const
-      {
-        float dz = z - zReference;
-        return trackParams[0] + dz * (trackParams[1] + dz * (trackParams[2] + dz * trackParams[3]));
-      }
-
-      __host__ __device__ float xSlope(const float z) const
-      {
-        float dz = z - zReference;
-        return trackParams[1] + dz * (2.f * trackParams[2] + 3.f * dz * trackParams[3]);
-      }
-
-      __host__ __device__ float y(const float z) const
-      {
-        float dz = z - zReference;
-        return trackParams[4] + dz * (trackParams[5] + dz * trackParams[6]);
-      }
-
-      __host__ __device__ float ySlope(const float z) const
-      {
-        float dz = z - zReference;
-        return trackParams[5] + dz * 2.f * trackParams[6];
-      }
-    };
-
-  } // namespace Tracking
-} // namespace SciFi
diff --git a/cuda/SciFi/PrForward/include/PrForwardTools.cuh b/cuda/SciFi/PrForward/include/PrForwardTools.cuh
deleted file mode 100644
index a4c924a45bfab2f62feef14eb5f1cbc10c1c5e86..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/PrForwardTools.cuh
+++ /dev/null
@@ -1,88 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <array>
-#include <vector>
-#include <algorithm>
-#include <fstream>
-
-#include <cassert>
-#include "Logger.h"
-#include "SystemOfUnits.h"
-#include "TMVA_Forward_1.cuh"
-#include "TMVA_Forward_2.cuh"
-#include "SciFiDefinitions.cuh"
-#include "VeloDefinitions.cuh"
-#include "UTDefinitions.cuh"
-#include "PrForwardConstants.cuh"
-#include "TrackUtils.cuh"
-#include "LinearFitting.cuh"
-#include "HitUtils.cuh"
-#include "FindXHits.cuh"
-#include "FindStereoHits.cuh"
-#include "VeloEventModel.cuh"
-#include "VeloConsolidated.cuh"
-#include "UTConsolidated.cuh"
-#include "SciFiEventModel.cuh"
-#include "States.cuh"
-
-__host__ __device__ void find_forward_tracks(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const float qop_ut,
-  const int i_veloUT_track,
-  SciFi::TrackHits* outputTracks,
-  uint* n_forward_tracks,
-  const int ut_event_number_of_tracks,
-  const SciFi::Tracking::TMVA* tmva1,
-  const SciFi::Tracking::TMVA* tmva2,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state);
-
-__host__ __device__ void selectFullCandidates(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track* candidate_tracks,
-  int& n_candidate_tracks,
-  SciFi::Tracking::Track* selected_tracks,
-  int& n_selected_tracks,
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  MiniState velo_state,
-  const float VeloUT_qOverP,
-  SciFi::Tracking::HitSearchCuts& pars_cur,
-  const SciFi::Tracking::TMVA* tmva1,
-  const SciFi::Tracking::TMVA* tmva2,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const bool secondLoop);
-
-__host__ __device__ SciFi::TrackHits makeTrack(SciFi::Tracking::Track track);
-
-template<class T>
-__host__ __device__ void sort_tracks(SciFi::Tracking::Track* tracks, const int n, const T& sort_function)
-{
-  // find permutations based on sort_function
-  uint permutations[SciFi::Tracking::max_selected_tracks];
-  for (int i = 0; i < n; ++i) {
-    uint position = 0;
-    for (int j = 0; j < n; ++j) {
-      const int sort_result = sort_function(tracks[i], tracks[j]);
-      position += sort_result > 0 || (sort_result == 0 && i > j);
-    }
-    permutations[position] = i;
-  }
-
-  // apply permutations, store tracks in temporary container
-  SciFi::Tracking::Track tracks_tmp[SciFi::Tracking::max_selected_tracks];
-  for (int i = 0; i < n; ++i) {
-    const int index = permutations[i];
-    tracks_tmp[i] = tracks[index];
-  }
-
-  // copy tracks back to original container
-  for (int i = 0; i < n; ++i) {
-    tracks[i] = tracks_tmp[i];
-  }
-}
diff --git a/cuda/SciFi/PrForward/include/ReferencePlaneProjection.cuh b/cuda/SciFi/PrForward/include/ReferencePlaneProjection.cuh
deleted file mode 100644
index b1e549f5948fd03bc8dcb890e81fd05b5d17d210..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/ReferencePlaneProjection.cuh
+++ /dev/null
@@ -1,27 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <array>
-#include <vector>
-#include <algorithm>
-#include <fstream>
-
-#include "SciFiDefinitions.cuh"
-#include "TrackUtils.cuh"
-#include "SciFiEventModel.cuh"
-
-/**
-   Project x hits onto reference plane
-*/
-
-__host__ __device__ void xAtRef_SamePlaneHits(
-  const SciFi::Hits& scifi_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int n_x_hits,
-  float coordX[SciFi::Tracking::max_x_hits],
-  const float xParams_seed[4],
-  const SciFi::Tracking::Arrays* constArrays,
-  MiniState velo_state,
-  const float zMag,
-  int itH,
-  int itEnd);
diff --git a/cuda/SciFi/PrForward/include/TrackUtils.cuh b/cuda/SciFi/PrForward/include/TrackUtils.cuh
deleted file mode 100644
index 2e4d33889c37837b37b4f88867620bb0aa625ef8..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/include/TrackUtils.cuh
+++ /dev/null
@@ -1,78 +0,0 @@
-#pragma once
-
-#include <cmath>
-
-#include "SciFiDefinitions.cuh"
-#include "PrForwardConstants.cuh"
-#include "HitUtils.cuh"
-#include "ParabolaFitting.cuh"
-#include "SciFiEventModel.cuh"
-
-/**
-   Helper functions related to track properties
- */
-
-// extrapolate x position from given state to z
-__host__ __device__ float xFromVelo(const float z, const MiniState& velo_state);
-
-// extrapolate y position from given state to z
-__host__ __device__ float yFromVelo(const float z, const MiniState& velo_state);
-
-__host__ __device__ float evalCubicParameterization(const float params[4], float z);
-
-__host__ __device__ float evalParameterizationX(const float* params, float z);
-
-__host__ __device__ float evalParameterizationY(const float* params, float z);
-
-__host__ __device__ bool lowerByQuality(SciFi::Tracking::Track t1, SciFi::Tracking::Track t2);
-
-__host__ __device__ inline float straightLinePropagationFromReferencePlane(const float params[4], float z)
-{
-  float dz = z - SciFi::Tracking::zReference;
-  return params[0] + params[1] * dz;
-}
-
-__host__ __device__ inline float straightLinePropagationFromReferencePlane(const float x0, const float tx, float z)
-{
-  float dz = z - SciFi::Tracking::zReference;
-  return x0 + tx * dz;
-}
-
-__host__ __device__ void getTrackParameters(
-  float xAtRef,
-  const MiniState& velo_state,
-  const SciFi::Tracking::Arrays* constArrays,
-  float trackParams[SciFi::Tracking::nTrackParams]);
-
-__host__ __device__ float calcqOverP(
-  float bx,
-  const SciFi::Tracking::Arrays* constArrays,
-  const MiniState& velo_state,
-  const float magnet_polarity);
-
-__host__ __device__ float zMagnet(const MiniState& velo_state, const SciFi::Tracking::Arrays* constArrays);
-
-__host__ __device__ float calcDxRef(float pt, const MiniState& velo_state);
-
-__host__ __device__ float
-trackToHitDistance(const float trackParameters[SciFi::Tracking::nTrackParams], const SciFi::Hits& scifi_hits, int hit);
-
-__host__ __device__ float chi2XHit(const float parsX[4], const SciFi::Hits& scifi_hits, const int hit);
-
-__host__ __device__ bool quadraticFitX(
-  const SciFi::Hits& scifi_hits,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  int coordToFit[SciFi::Tracking::max_coordToFit],
-  int& n_coordToFit,
-  PlaneCounter& planeCounter,
-  SciFi::Tracking::HitSearchCuts& pars_cur);
-
-__host__ __device__ bool fitYProjection(
-  const SciFi::Hits& scifi_hits,
-  SciFi::Tracking::Track& track,
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits,
-  PlaneCounter& planeCounter,
-  const MiniState& velo_state,
-  const SciFi::Tracking::Arrays* constArrays,
-  SciFi::Tracking::HitSearchCuts& pars_cur);
diff --git a/cuda/SciFi/PrForward/src/FindStereoHits.cu b/cuda/SciFi/PrForward/src/FindStereoHits.cu
deleted file mode 100644
index ce5a89f300ea63008884d190ec05176c5120e3ab..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/FindStereoHits.cu
+++ /dev/null
@@ -1,247 +0,0 @@
-#include "FindStereoHits.cuh"
-
-//=========================================================================
-//  Collect all hits in the stereo planes compatible with the track
-//=========================================================================
-__host__ __device__ void collectStereoHits(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track& track,
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars,
-  const SciFi::Tracking::Arrays* constArrays,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits)
-{
-  for (int zone = 0; zone < SciFi::Constants::n_layers; ++zone) {
-    // get yZone and xPred: x and y values at z of layer based on candidate track parameters
-    const float parsX[4] = {track.trackParams[0], track.trackParams[1], track.trackParams[2], track.trackParams[3]};
-    const float parsY[4] = {track.trackParams[4], track.trackParams[5], track.trackParams[6], 0.f};
-    float zZone = constArrays->uvZone_zPos[zone];
-    const float yZone = evalCubicParameterization(parsY, zZone);
-    assert(constArrays->uvZones[zone] < SciFi::Constants::n_zones);
-    zZone += constArrays->Zone_dzdy[constArrays->uvZones[zone]] * yZone; // Correct for dzDy
-    const float xPred = evalCubicParameterization(parsX, zZone);
-
-    const bool triangleSearch = fabsf(yZone) < SciFi::Tracking::tolYTriangleSearch;
-    // even zone number: if ( yZone > 0 ) continue;
-    // odd zone number: if ( -yZone > 0 ) continue;
-    // -> check for upper / lower half
-    // -> only continue if yZone is in the correct half
-    if (!triangleSearch && (2.f * float(((constArrays->uvZones[zone]) % 2) == 0) - 1.f) * yZone > 0.f) continue;
-
-    const float dxDySign = constArrays->uvZone_dxdy[zone] < 0.f ? -1.f : 1.f;
-    const float seed_x_at_zZone = xFromVelo(zZone, velo_state);
-    const float dxTol =
-      SciFi::Tracking::tolY + SciFi::Tracking::tolYSlope * (fabsf(xPred - seed_x_at_zZone) + fabsf(yZone));
-
-    // find stereo hits whose x coordinate is within xTol of the prediction
-    // from the candidate track
-    // This takes the y value (max, min) into account
-    const float lower_bound_at = -dxTol - yZone * constArrays->uvZone_dxdy[zone] + xPred;
-    int uv_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->uvZones[zone]);
-    int uv_zone_offset_end = uv_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->uvZones[zone]);
-
-    int itBegin = getLowerBound(scifi_hits.x0, lower_bound_at, uv_zone_offset_begin, uv_zone_offset_end);
-    int itEnd = uv_zone_offset_end;
-
-    findStereoHitsWithinXTol(
-      itBegin,
-      itEnd,
-      scifi_hits,
-      yZone,
-      xPred,
-      dxTol,
-      triangleSearch,
-      dxDySign,
-      n_stereoHits,
-      stereoCoords,
-      stereoHits);
-
-    if (n_stereoHits >= SciFi::Tracking::max_stereo_hits) break;
-  }
-
-  // Sort hits by coord
-  // not using thrust::sort due to temporary_buffer::allocate:: get_temporary_buffer failed" error
-  // thrust::sort_by_key(thrust::seq, stereoCoords, stereoCoords + n_stereoHits, stereoHits);
-  sortHitsByKey<SciFi::Tracking::max_stereo_hits>(stereoCoords, n_stereoHits, stereoHits);
-}
-
-//=========================================================================
-//  Fit the stereo hits
-//=========================================================================
-__host__ __device__ bool selectStereoHits(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track& track,
-  const SciFi::Tracking::Arrays* constArrays,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits,
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars)
-{
-  int bestStereoHits[SciFi::Tracking::max_stereo_hits];
-  int n_bestStereoHits = 0;
-  float originalYParams[3] = {track.trackParams[4], track.trackParams[5], track.trackParams[6]};
-  float bestYParams[3];
-  float bestMeanDy = 1e9f;
-
-  if (pars.minStereoHits > n_stereoHits) return false;
-  int endLoop = n_stereoHits - pars.minStereoHits;
-
-  PlaneCounter planeCounter;
-  for (int beginRange = 0; beginRange < endLoop; ++beginRange) {
-    planeCounter.clear();
-    int endRange = beginRange;
-
-    float sumCoord = 0.;
-    // bad hack to reproduce itereator behavior from before: *(-1) = 0
-    int first_hit;
-    if (endRange == 0)
-      first_hit = 0;
-    else
-      first_hit = endRange - 1;
-
-    findStereoHitClusterByDx(
-      planeCounter, endRange, pars, stereoCoords, stereoHits, n_stereoHits, scifi_hits, sumCoord, first_hit);
-
-    cleanStereoHitCluster(
-      beginRange, endRange, n_stereoHits, stereoHits, stereoCoords, sumCoord, planeCounter, scifi_hits);
-
-    // Now we have a candidate, lets fit him
-    // track = original; //only yparams are changed
-    track.trackParams[4] = originalYParams[0];
-    track.trackParams[5] = originalYParams[1];
-    track.trackParams[6] = originalYParams[2];
-
-    int trackStereoHits[SciFi::Tracking::max_stereo_hits];
-    int n_trackStereoHits = 0;
-    assert(endRange < n_stereoHits);
-    for (int range = beginRange; range < endRange; ++range) {
-      trackStereoHits[n_trackStereoHits++] = stereoHits[range];
-    }
-
-    // fit Y Projection of track using stereo hits
-    if (!fitYProjection(
-          scifi_hits, track, trackStereoHits, n_trackStereoHits, planeCounter, velo_state, constArrays, pars))
-      continue;
-
-    if (!addHitsOnEmptyStereoLayers(
-          scifi_hits,
-          scifi_hit_count,
-          track,
-          trackStereoHits,
-          n_trackStereoHits,
-          constArrays,
-          planeCounter,
-          velo_state,
-          pars))
-      continue;
-
-    if (n_trackStereoHits < n_bestStereoHits) continue; // number of hits most important selection criteria!
-
-    //== Calculate  dy chi2 /ndf
-    float meanDy = 0.;
-    assert(n_trackStereoHits < n_stereoHits);
-    for (int i_hit = 0; i_hit < n_trackStereoHits; ++i_hit) {
-      const int hit = trackStereoHits[i_hit];
-      const float d = trackToHitDistance(track.trackParams, scifi_hits, hit) / scifi_hits.dxdy(hit);
-      meanDy += d * d;
-    }
-    meanDy /= float(n_trackStereoHits - 1);
-
-    if (n_trackStereoHits > n_bestStereoHits || meanDy < bestMeanDy) {
-      // if same number of hits take smaller chi2
-      bestYParams[0] = track.trackParams[4];
-      bestYParams[1] = track.trackParams[5];
-      bestYParams[2] = track.trackParams[6];
-      bestMeanDy = meanDy;
-
-      n_bestStereoHits = 0;
-      for (int i_hit = 0; i_hit < n_trackStereoHits; ++i_hit) {
-        assert(n_bestStereoHits < SciFi::Tracking::max_stereo_hits);
-        bestStereoHits[n_bestStereoHits++] = trackStereoHits[i_hit];
-      }
-    }
-
-  } // beginRange loop (<endLoop)
-
-  if (n_bestStereoHits > 0) {
-    track.trackParams[4] = bestYParams[0];
-    track.trackParams[5] = bestYParams[1];
-    track.trackParams[6] = bestYParams[2];
-    assert(n_bestStereoHits < n_stereoHits);
-
-    for (int i_hit = 0; i_hit < n_bestStereoHits; ++i_hit) {
-      int hit = bestStereoHits[i_hit];
-      if (track.hitsNum >= SciFi::Tracking::max_scifi_hits) break;
-      assert(track.hitsNum < SciFi::Tracking::max_scifi_hits);
-      track.addHit(hit);
-    }
-    return true;
-  }
-  return false;
-}
-
-//=========================================================================
-//  Add hits on empty stereo layers, and refit if something was added
-//=========================================================================
-__host__ __device__ bool addHitsOnEmptyStereoLayers(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track& track,
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits,
-  const SciFi::Tracking::Arrays* constArrays,
-  PlaneCounter& planeCounter,
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars)
-{
-  // at this point pc is counting only stereo HITS!
-  if (planeCounter.nbDifferent > 5) return true;
-
-  bool added = false;
-  for (unsigned int zone = 0; zone < SciFi::Constants::n_layers; zone += 1) {
-    assert(constArrays->uvZones[zone] < SciFi::Constants::n_zones);
-    if (planeCounter.nbInPlane(constArrays->uvZones[zone] / 2) != 0) continue; // there is already one hit
-
-    float zZone = constArrays->uvZone_zPos[zone];
-
-    const float parsX[4] = {track.trackParams[0], track.trackParams[1], track.trackParams[2], track.trackParams[3]};
-    const float parsY[4] = {track.trackParams[4], track.trackParams[5], track.trackParams[6], 0.};
-
-    float yZone = evalCubicParameterization(parsY, zZone);
-    zZone = constArrays->Zone_dzdy[constArrays->uvZones[zone]] * yZone; // Correct for dzDy
-    yZone = evalCubicParameterization(parsY, zZone);
-    const float xPred = evalCubicParameterization(parsX, zZone);
-
-    const bool triangleSearch = fabsf(yZone) < SciFi::Tracking::tolYTriangleSearch;
-    // change sign of yZone depending on whether we are in the upper or lower half
-    if (!triangleSearch && (2.f * float((((constArrays->uvZones[zone]) % 2) == 0)) - 1.f) * yZone > 0.f) continue;
-
-    // only version without triangle search!
-    const float dxTol =
-      SciFi::Tracking::tolY + SciFi::Tracking::tolYSlope *
-                                (fabsf(xPred - velo_state.x + (zZone - velo_state.z) * velo_state.tx) + fabsf(yZone));
-    // -- Use a binary search to find the lower bound of the range of x values
-    // -- This takes the y value into account
-    const float lower_bound_at = -dxTol - yZone * constArrays->uvZone_dxdy[zone] + xPred;
-    int uv_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->uvZones[zone]);
-    int uv_zone_offset_end = uv_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->uvZones[zone]);
-    int itBegin = getLowerBound(scifi_hits.x0, lower_bound_at, uv_zone_offset_begin, uv_zone_offset_end);
-    int itEnd = uv_zone_offset_end;
-
-    int best = findBestStereoHitOnEmptyLayer(itBegin, itEnd, scifi_hits, yZone, xPred, dxTol, triangleSearch);
-
-    if (-1 != best) {
-      assert(n_stereoHits < SciFi::Tracking::max_stereo_hits);
-      stereoHits[n_stereoHits++] = best;
-      planeCounter.addHit(scifi_hits.planeCode(best) / 2);
-      added = true;
-    }
-  }
-  if (!added) return true;
-  return fitYProjection(scifi_hits, track, stereoHits, n_stereoHits, planeCounter, velo_state, constArrays, pars);
-}
diff --git a/cuda/SciFi/PrForward/src/FindXHits.cu b/cuda/SciFi/PrForward/src/FindXHits.cu
deleted file mode 100644
index c952be9bb087a55521597ebd797794c3595dad08..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/FindXHits.cu
+++ /dev/null
@@ -1,829 +0,0 @@
-#include "FindXHits.cuh"
-#include "BinarySearch.cuh"
-
-__host__ void collectAllXHits_proto_p(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state,
-  const MiniState& UT_state,
-  const float qOverP,
-  int side,
-  std::array<int, 2 * 6>& windows_x,
-  std::array<int, 2 * 6>& windows_uv,
-  std::array<float, 4 * 6>& parameters_uv,
-  const SciFiWindowsParams& window_params,
-  const std::array<int, 12> true_scifi_indices_per_layer)
-{
-  const float tx2 = velo_state.tx * velo_state.tx;
-  const float ty2 = velo_state.ty * velo_state.ty;
-  const float slope2 = tx2 + ty2;
-  const float pt = sqrtf(slope2 / (1.f + slope2)) / fabsf(qOverP);
-  const float p = 1.f / std::abs(qOverP);
-
-  /* OPTIMIZE: possibly use these variables for wrong sign treatment */
-  // const float q = qOverP > 0.f ? 1.f : -1.f;
-  // const bool wSignTreatment = pt > SciFi::Tracking::wrongSignPT;
-  // float zMag = zMagnet(velo_state, constArrays);
-  // const float qop_WS = sqrtf(slope2 / (1.f + slope2) ) / pt;
-  // float dxRefWS = 0.f;
-  // if ( wSignTreatment ) {
-  //   dxRefWS = 0.9f * calcDxRef(SciFi::Tracking::wrongSignPT, velo_state);
-  // }
-  // const float dir = q * SciFi::Tracking::magscalefactor * (-1.f); // needed for wrong sign treatment
-  // const float xTolWS = dx_calc(velo_state, qop_WS, window_params);
-
-  // const float q = qOverP > 0.f ? 1.f : -1.f;
-  // const float dir = q * magnet_polarity * (-1.f);
-
-  // const bool wSignTreatment = pt > SciFi::Tracking::wrongSignPT;
-  float zMag = zMagnet(velo_state, constArrays);
-  const float qop_WS = sqrtf(slope2 / (1.f + slope2)) / pt;
-  // float dxRefWS = 0.f;
-  // if ( wSignTreatment ) {
-  //   dxRefWS = 0.9f * calcDxRef(SciFi::Tracking::wrongSignPT, velo_state);
-  // }
-  const float xAtRef = xFromVelo(SciFi::Tracking::zReference, UT_state);
-  float xParams_seed[4] {xAtRef, UT_state.tx, 0, 0};
-
-  // use parametrization to propagate from UT to SciFi
-  const auto state_zRef = propagate_state_from_velo(UT_state, qOverP, 5); // zRef is between layers 4 and 5
-  const float xTol = dx_calc(velo_state, qOverP, window_params);
-  int iZoneStartingPoint = side > 0 ? constArrays->zoneoffsetpar : 0;
-
-  for (unsigned int iZone = iZoneStartingPoint; iZone < iZoneStartingPoint + constArrays->zoneoffsetpar; iZone++) {
-    assert(iZone - iZoneStartingPoint < SciFi::Constants::n_zones);
-    assert(iZone - iZoneStartingPoint < 12);
-
-    const auto izone_rel = iZone - iZoneStartingPoint;
-    const float zZone = constArrays->xZone_zPos[izone_rel];
-
-    // const int layer = constArrays->xZones[iZone] / 2;
-    const float dz_x = (zZone - SciFi::Tracking::zReference);
-    const float xInZone = scifi_propagation(state_zRef.x, UT_state.tx, qOverP, dz_x);
-    const float yInZone = yFromVelo(zZone, velo_state);
-
-    const float xInZoneStraight = evalCubicParameterization(xParams_seed, zZone);
-
-    if (side > 0) {
-      if (
-        !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
-        !isInside(yInZone, SciFi::Tracking::yLim_Min, SciFi::Tracking::yLim_Max))
-        continue;
-    }
-    else {
-      if (
-        !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
-        !isInside(yInZone, side * SciFi::Tracking::yLim_Max, side * SciFi::Tracking::yLim_Min))
-        continue;
-    }
-
-    float xMin = xInZone - xTol;
-    float xMax = xInZone + xTol;
-
-    /* OPTIMIZE: how do we take care of wrong sign tracks with momentum windows? */
-    // float xTolWS = 0.0;
-    // if (wSignTreatment) {
-    //   // xTolWS = (zZone < SciFi::Tracking::zReference) ?
-    //   //   dxRefWS * zZone / SciFi::Tracking::zReference :
-    //   //   dxRefWS * (zZone - zMag) / (SciFi::Tracking::zReference - zMag);
-    //   // if (dir > 0) {
-    //   //   xMin = xInZone - xTolWS;
-    //   // }
-    //   // else {
-    //   //   xMax = xInZone + xTolWS;
-    //   // }
-
-    //   debug_cout << "\t before WS treatment: xMin = " << xMin << ", xMax = " << xMax << ", WS = " <<
-    //   int(wSignTreatment) << ", pt = " << pt << std::endl; if (dir > 0) {
-    //     xMin = -1.f * xInZone - xTolWS;
-    //   }
-    //   else {
-    //     xMax = xInZone + xTolWS;
-    //   }
-    //   debug_cout << "\t after WS treatment: xMin = " << xMin << ", xMax = " << xMax << std::endl;
-    // }
-
-    // Get the hits within the bounds
-    assert(iZone < SciFi::Constants::n_layers);
-    assert(constArrays->xZones[iZone] < SciFi::Constants::n_zones);
-    int x_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->xZones[iZone]);
-    int x_zone_offset_end = x_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->xZones[iZone]);
-    const int itH = getLowerBound(scifi_hits.x0, xMin, x_zone_offset_begin, x_zone_offset_end);
-    const int itEnd = getLowerBound(scifi_hits.x0, xMax, x_zone_offset_begin, x_zone_offset_end);
-    assert(itH >= x_zone_offset_begin && itH <= x_zone_offset_end);
-    assert(itEnd >= x_zone_offset_begin && itEnd <= x_zone_offset_end);
-
-    windows_x[2 * izone_rel] = itH;
-    windows_x[2 * izone_rel + 1] = itEnd - itH;
-
-    // Now match the stereo hits
-    const float this_uv_z = constArrays->uvZone_zPos[izone_rel];
-    const float dz_uv = this_uv_z - SciFi::Tracking::zReference;
-    const float yInUVZone = yFromVelo(this_uv_z, velo_state);
-    const float dx = yInUVZone * constArrays->uvZone_dxdy[izone_rel];
-    const float xPredUv = scifi_propagation(state_zRef.x, UT_state.tx, qOverP, dz_uv) - dx;
-    const int uv_layer = constArrays->uvZones[iZone] / 2;
-    // To Do: study this window
-    const float xBound = 70.f * window_params.extrapolation_stddev[uv_layer];
-    const float maxDx = xBound; // * ratio;
-
-    const float xMinUV = xPredUv - maxDx;
-    const float xMaxUV = xPredUv + maxDx;
-
-    // Get bounds in UV layers
-    // do one search on the same side as the x module
-    // if we are close to y = 0, also look within a region on the other side module ("triangle search")
-    assert(constArrays->uvZones[iZone] < SciFi::Constants::n_zones);
-    const int uv_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->uvZones[iZone]);
-    const int uv_zone_offset_end =
-      uv_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone]);
-    /* OPTIMIZE: check how large the effect on the efficiency is to include the triangle hits */
-    //    const int triangleOffset = side > 0 ? -1 : 1;
-    // assert(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset] < SciFi::Constants::n_zones);
-    // const int triangle_zone_offset_begin =
-    //   scifi_hit_count.zone_offset(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset]);
-    // assert(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset] < SciFi::Constants::n_zones);
-    // const int triangle_zone_offset_end =
-    //   triangle_zone_offset_begin +
-    //   scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset]);
-    int itUVStart = getLowerBound(scifi_hits.x0, xMinUV, uv_zone_offset_begin, uv_zone_offset_end);
-    int itUVEnd = getLowerBound(scifi_hits.x0, xMaxUV, uv_zone_offset_begin, uv_zone_offset_end);
-    //    int itUV2 = getLowerBound(scifi_hits.x0, xMinUV, triangle_zone_offset_begin, triangle_zone_offset_end);
-
-    windows_uv[2 * izone_rel] = itUVStart;
-    windows_uv[2 * izone_rel + 1] = itUVEnd;
-  }
-}
-
-__host__ void collectAllXHits_proto(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state,
-  const float qOverP,
-  int side,
-  std::array<int, 2 * 6>& windows_x,
-  std::array<int, 2 * 6>& windows_uv,
-  std::array<float, 4 * 6>& parameters_uv)
-{
-  // Find size of search window on reference plane, using Velo slopes and min pT as input
-  float dxRef = 0.9f * calcDxRef(SciFi::Tracking::minPt, velo_state);
-  // find position within magnet where bending happens
-  float zMag = zMagnet(velo_state, constArrays);
-
-  const float q = qOverP > 0.f ? 1.f : -1.f;
-  const float dir = q * magnet_polarity * (-1.f);
-
-  float slope2 = velo_state.tx * velo_state.tx + velo_state.ty * velo_state.ty;
-  const float pt = std::sqrt(slope2 / (1.f + slope2)) / std::abs(qOverP);
-  const bool wSignTreatment = SciFi::Tracking::useWrongSignWindow && pt > SciFi::Tracking::wrongSignPT;
-
-  float dxRefWS = 0.f;
-  if (wSignTreatment) {
-    dxRefWS = 0.9f * calcDxRef(
-                       SciFi::Tracking::wrongSignPT,
-                       velo_state); // make windows a bit too small - FIXME check effect of this, seems wrong
-  }
-
-  int iZoneStartingPoint = side > 0 ? constArrays->zoneoffsetpar : 0;
-
-  for (unsigned int iZone = iZoneStartingPoint; iZone < iZoneStartingPoint + constArrays->zoneoffsetpar; iZone++) {
-    assert(iZone - iZoneStartingPoint < SciFi::Constants::n_zones);
-    assert(iZone - iZoneStartingPoint < 12);
-
-    const auto izone_rel = iZone - iZoneStartingPoint;
-    const float zZone = constArrays->xZone_zPos[izone_rel];
-    const float xInZone = evalCubicParameterization(xParams_seed, zZone);
-    const float yInZone = evalCubicParameterization(yParams_seed, zZone);
-
-    // Now the code checks if the x and y are in the zone limits. I am really not sure
-    // why this is done here, surely could just check if within limits for the last zone
-    // in T3 and go from there? Need to think more about this.
-    //
-    // Here for now I assume the same min/max x and y for all stations, this again needs to
-    // be read from some file blablabla although actually I suspect having some general tolerances
-    // here is anyway good enough since we are doing a straight line extrapolation in the first place
-    // check (roughly) whether the extrapolated velo track is within the current zone
-    // if (side > 0) {
-    //   if (
-    //     !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
-    //     !isInside(yInZone, SciFi::Tracking::yLim_Min, SciFi::Tracking::yLim_Max))
-    //     continue;
-    // }
-    // else {
-    //   if (
-    //     !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
-    //     !isInside(yInZone, side * SciFi::Tracking::yLim_Max, side * SciFi::Tracking::yLim_Min))
-    //     continue;
-    // }
-
-    // extrapolate dxRef (x window on reference plane) to plane of current zone
-    const float xTol = (zZone < SciFi::Tracking::zReference) ?
-                         dxRef * zZone / SciFi::Tracking::zReference :
-                         dxRef * (zZone - zMag) / (SciFi::Tracking::zReference - zMag);
-    float xMin = xInZone - xTol;
-    float xMax = xInZone + xTol;
-
-    if (SciFi::Tracking::useMomentumEstimate) { // For VeloUT tracks, suppress check if track actually has qOverP set,
-                                                // get the option right!
-      float xTolWS = 0.0;
-      if (wSignTreatment) {
-        xTolWS = (zZone < SciFi::Tracking::zReference) ?
-                   dxRefWS * zZone / SciFi::Tracking::zReference :
-                   dxRefWS * (zZone - zMag) / (SciFi::Tracking::zReference - zMag);
-      }
-      if (dir > 0) {
-        xMin = xInZone - xTolWS;
-      }
-      else {
-        xMax = xInZone + xTolWS;
-      }
-    }
-
-    // Get the hits within the bounds
-    assert(iZone < SciFi::Constants::n_layers);
-    assert(constArrays->xZones[iZone] < SciFi::Constants::n_zones);
-    int x_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->xZones[iZone]);
-    int x_zone_size = scifi_hit_count.zone_number_of_hits(constArrays->xZones[iZone]);
-    int itH = binary_search_leftmost(scifi_hits.x0 + x_zone_offset_begin, x_zone_size, xMin);
-    const int itSize = binary_search_leftmost(scifi_hits.x0 + x_zone_offset_begin + itH, x_zone_size - itH, xMax);
-    itH += x_zone_offset_begin;
-
-    windows_x[2 * izone_rel] = itH;
-    windows_x[2 * izone_rel + 1] = itSize;
-
-    // Skip making range but continue if the size is zero
-    if (itSize == 0) continue;
-
-    // Now match the stereo hits
-    const float this_uv_z = constArrays->uvZone_zPos[iZone - iZoneStartingPoint];
-    const float xInUv = evalCubicParameterization(xParams_seed, this_uv_z);
-    const float zRatio = (this_uv_z - zMag) / (zZone - zMag);
-    const float dx = yInZone * constArrays->uvZone_dxdy[iZone - iZoneStartingPoint];
-    const float xCentral = xInZone + dx;
-    const float xPredUv = xInUv + (scifi_hits.x0[itH] - xInZone) * zRatio - dx;
-    const float maxDx = SciFi::Tracking::tolYCollectX +
-                        (fabsf(scifi_hits.x0[itH] - xCentral) + fabsf(yInZone)) * SciFi::Tracking::tolYSlopeCollectX;
-    const float xMinUV = xPredUv - maxDx;
-
-    // Get bounds in UV layers
-    // do one search on the same side as the x module
-    // if we are close to y = 0, also look within a region on the other side module ("triangle search")
-    assert(constArrays->uvZones[iZone] < SciFi::Constants::n_zones);
-    const int uv_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->uvZones[iZone]);
-    const int uv_zone_size = scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone]);
-    const int itUV1 = binary_search_leftmost(scifi_hits.x0 + uv_zone_offset_begin, uv_zone_size, xMinUV);
-
-    const float xPredUVProto = xInUv - xInZone * zRatio - dx;
-    const float maxDxProto = SciFi::Tracking::tolYCollectX + fabsf(yInZone) * SciFi::Tracking::tolYSlopeCollectX;
-
-    windows_uv[2 * izone_rel] = itUV1 + uv_zone_offset_begin;
-    windows_uv[2 * izone_rel + 1] = uv_zone_size - itUV1;
-
-    parameters_uv[4 * izone_rel] = xPredUVProto;
-    parameters_uv[4 * izone_rel + 1] = zRatio;
-    parameters_uv[4 * izone_rel + 2] = maxDxProto;
-    parameters_uv[4 * izone_rel + 3] = xCentral;
-  }
-}
-
-//=========================================================================
-// From LHCb Forward tracking description
-//
-// Collect all X hits, within a window defined by the minimum Pt.
-// Better restrictions possible, if we use the momentum of the input track.
-// Ask for the presence of a stereo hit in the same biLayer compatible.
-// This reduces the efficiency. X-alone hits to be re-added later in the processing
-//
-// side = 1  -> upper y half
-// side = -1 -> lower y half
-//=========================================================================
-//
-__host__ __device__ void collectAllXHits(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  int allXHits[SciFi::Tracking::max_x_hits],
-  int& n_x_hits,
-  float coordX[SciFi::Tracking::max_x_hits],
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state,
-  const float qOverP,
-  int side)
-{
-  // Find size of search window on reference plane, using Velo slopes and min pT as input
-  float dxRef = 0.9f * calcDxRef(SciFi::Tracking::minPt, velo_state);
-  // find position within magnet where bending happens
-  float zMag = zMagnet(velo_state, constArrays);
-
-  const float q = qOverP > 0.f ? 1.f : -1.f;
-  const float dir = q * magnet_polarity * (-1.f);
-
-  float slope2 = velo_state.tx * velo_state.tx + velo_state.ty * velo_state.ty;
-  const float pt = std::sqrt(slope2 / (1.f + slope2)) / std::abs(qOverP);
-  const bool wSignTreatment = SciFi::Tracking::useWrongSignWindow && pt > SciFi::Tracking::wrongSignPT;
-
-  float dxRefWS = 0.f;
-  if (wSignTreatment) {
-    // DvB: what happens if we use the actual momentum from VeloUT here instead of a constant?
-    dxRefWS = 0.9f * calcDxRef(
-                       SciFi::Tracking::wrongSignPT,
-                       velo_state); // make windows a bit too small - FIXME check effect of this, seems wrong
-  }
-
-  int iZoneEnd[7]; // 6 x planes
-  iZoneEnd[0] = 0;
-  int cptZone = 1;
-
-  int iZoneStartingPoint = side > 0 ? constArrays->zoneoffsetpar : 0;
-
-  for (unsigned int iZone = iZoneStartingPoint; iZone < iZoneStartingPoint + constArrays->zoneoffsetpar; iZone++) {
-    assert(iZone - iZoneStartingPoint < SciFi::Constants::n_zones);
-    assert(iZone - iZoneStartingPoint < 12);
-    const float zZone = constArrays->xZone_zPos[iZone - iZoneStartingPoint];
-    const float xInZone = evalCubicParameterization(xParams_seed, zZone);
-    const float yInZone = evalCubicParameterization(yParams_seed, zZone);
-
-    // Now the code checks if the x and y are in the zone limits. I am really not sure
-    // why this is done here, surely could just check if within limits for the last zone
-    // in T3 and go from there? Need to think more about this.
-    //
-    // Here for now I assume the same min/max x and y for all stations, this again needs to
-    // be read from some file blablabla although actually I suspect having some general tolerances
-    // here is anyway good enough since we are doing a straight line extrapolation in the first place
-    // check (roughly) whether the extrapolated velo track is within the current zone
-    if (side > 0) {
-      if (
-        !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
-        !isInside(yInZone, SciFi::Tracking::yLim_Min, SciFi::Tracking::yLim_Max))
-        continue;
-    }
-    else {
-      if (
-        !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
-        !isInside(yInZone, side * SciFi::Tracking::yLim_Max, side * SciFi::Tracking::yLim_Min))
-        continue;
-    }
-
-    // extrapolate dxRef (x window on reference plane) to plane of current zone
-    const float xTol = (zZone < SciFi::Tracking::zReference) ?
-                         dxRef * zZone / SciFi::Tracking::zReference :
-                         dxRef * (zZone - zMag) / (SciFi::Tracking::zReference - zMag);
-    float xMin = xInZone - xTol;
-    float xMax = xInZone + xTol;
-
-    if (SciFi::Tracking::useMomentumEstimate) { // For VeloUT tracks, suppress check if track actually has qOverP set,
-                                                // get the option right!
-      float xTolWS = 0.0;
-      if (wSignTreatment) {
-        xTolWS = (zZone < SciFi::Tracking::zReference) ?
-                   dxRefWS * zZone / SciFi::Tracking::zReference :
-                   dxRefWS * (zZone - zMag) / (SciFi::Tracking::zReference - zMag);
-      }
-      if (dir > 0) {
-        xMin = xInZone - xTolWS;
-      }
-      else {
-        xMax = xInZone + xTolWS;
-      }
-    }
-
-    // Get the hits within the bounds
-    assert(iZone < SciFi::Constants::n_layers);
-    assert(constArrays->xZones[iZone] < SciFi::Constants::n_zones);
-    int x_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->xZones[iZone]);
-    int x_zone_offset_end = x_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->xZones[iZone]);
-    const int itH = getLowerBound(scifi_hits.x0, xMin, x_zone_offset_begin, x_zone_offset_end);
-    const int itEnd = getLowerBound(scifi_hits.x0, xMax, x_zone_offset_begin, x_zone_offset_end);
-    assert(itH >= x_zone_offset_begin && itH <= x_zone_offset_end);
-    assert(itEnd >= x_zone_offset_begin && itEnd <= x_zone_offset_end);
-
-    // Skip making range but continue if the end is before or equal to the start
-    if (!(itEnd > itH)) continue;
-
-    // Now match the stereo hits
-    const float this_uv_z = constArrays->uvZone_zPos[iZone - iZoneStartingPoint];
-    const float xInUv = evalCubicParameterization(xParams_seed, this_uv_z);
-    const float zRatio = (this_uv_z - zMag) / (zZone - zMag);
-    const float dx = yInZone * constArrays->uvZone_dxdy[iZone - iZoneStartingPoint];
-    const float xCentral = xInZone + dx;
-    const float xPredUv = xInUv + (scifi_hits.x0[itH] - xInZone) * zRatio - dx;
-    const float maxDx = SciFi::Tracking::tolYCollectX +
-                        (fabsf(scifi_hits.x0[itH] - xCentral) + fabsf(yInZone)) * SciFi::Tracking::tolYSlopeCollectX;
-    const float xMinUV = xPredUv - maxDx;
-
-    // Get bounds in UV layers
-    // do one search on the same side as the x module
-    // if we are close to y = 0, also look within a region on the other side module ("triangle search")
-    assert(constArrays->uvZones[iZone] < SciFi::Constants::n_zones);
-    const int uv_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->uvZones[iZone]);
-    const int uv_zone_offset_end =
-      uv_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone]);
-    const int triangleOffset = side > 0 ? -1 : 1;
-    assert(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset] < SciFi::Constants::n_zones);
-    const int triangle_zone_offset_begin =
-      scifi_hit_count.zone_offset(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset]);
-    assert(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset] < SciFi::Constants::n_zones);
-    const int triangle_zone_offset_end =
-      triangle_zone_offset_begin +
-      scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset]);
-    int itUV1 = getLowerBound(scifi_hits.x0, xMinUV, uv_zone_offset_begin, uv_zone_offset_end);
-    int itUV2 = getLowerBound(scifi_hits.x0, xMinUV, triangle_zone_offset_begin, triangle_zone_offset_end);
-
-    const float xPredUVProto = xInUv - xInZone * zRatio - dx;
-    const float maxDxProto = SciFi::Tracking::tolYCollectX + fabsf(yInZone) * SciFi::Tracking::tolYSlopeCollectX;
-
-    const bool withTriangleSearch = fabsf(yInZone) < SciFi::Tracking::tolYTriangleSearch;
-    for (int xHit = itH; xHit < itEnd; ++xHit) { // loop over all xHits in a layer between xMin and xMax
-      const float xPredUv = xPredUVProto + scifi_hits.x0[xHit] * zRatio;
-      const float maxDx = maxDxProto + fabsf(scifi_hits.x0[xHit] - xCentral) * SciFi::Tracking::tolYSlopeCollectX;
-      const float xMinUV = xPredUv - maxDx;
-      const float xMaxUV = xPredUv + maxDx;
-
-      if (
-        matchStereoHit(itUV1, uv_zone_offset_end, scifi_hits, xMinUV, xMaxUV) //) {
-        || (withTriangleSearch &&
-            matchStereoHitWithTriangle(itUV2, triangle_zone_offset_end, yInZone, scifi_hits, xMinUV, xMaxUV, side))) {
-        if (n_x_hits >= SciFi::Tracking::max_x_hits) break;
-        assert(n_x_hits < SciFi::Tracking::max_x_hits);
-        allXHits[n_x_hits++] = xHit;
-      }
-    }
-
-    const int iStart = iZoneEnd[cptZone - 1];
-    const int iEnd = n_x_hits;
-
-    assert(cptZone < 7);
-    iZoneEnd[cptZone++] = iEnd;
-
-    // project x of all hits to reference plane
-    // save it in coordX
-    // -> first step of 1D Hough transform,
-    // select clusters of x hits on reference plane in selectXCandidates
-    if (iStart < iEnd) {
-      xAtRef_SamePlaneHits(
-        scifi_hits, allXHits, n_x_hits, coordX, xParams_seed, constArrays, velo_state, zMag, iStart, iEnd);
-    }
-    if (n_x_hits >= SciFi::Tracking::max_x_hits) break;
-  }
-
-  // Sort hits by x on reference plane
-  // not using thrust::sort due to "temporary_buffer::allocate: get_temporary_buffer failed" error
-  // every time thrust::sort is called, cudaMalloc is called, apparently there can be trouble
-  // doing this many times
-  // thrust::sort_by_key(thrust::seq, coordX, coordX + n_x_hits, allXHits);
-  sortHitsByKey<SciFi::Tracking::max_x_hits>(coordX, n_x_hits, allXHits);
-}
-
-__host__ __device__ void improveXCluster(
-  int& it2,
-  const int it1,
-  const int itEnd,
-  const int n_x_hits,
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const float xWindow,
-  const SciFi::Tracking::HitSearchCuts& pars,
-  PlaneCounter& planeCounter,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits)
-{
-
-  int itLast = it2 - 1;
-  while (it2 < itEnd) {
-    assert(it2 < n_x_hits);
-    if (usedHits[it2]) {
-      ++it2;
-      continue;
-    }
-    // now  the first and last+1 hit exist and are not used!
-
-    // Add next hit,
-    // if there is only a small gap between the hits
-    //    or inside window and plane is still empty
-    assert(it2 < itEnd);
-    if (
-      (coordX[it2] < coordX[itLast] + pars.maxXGap) ||
-      ((coordX[it2] - coordX[it1] < xWindow) &&
-       (planeCounter.nbInPlane(scifi_hits.planeCode(allXHits[it2]) / 2) == 0))) {
-      planeCounter.addHit(scifi_hits.planeCode(allXHits[it2]) / 2);
-      itLast = it2;
-      ++it2;
-      continue;
-    }
-    // Found nothing to improve
-    else {
-      break;
-    }
-  }
-}
-
-//=========================================================================
-//  Select the zones in the allXHits array where we can have a track
-//=========================================================================
-__host__ __device__ void selectXCandidates(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  int allXHits[SciFi::Tracking::max_x_hits],
-  int& n_x_hits,
-  bool usedHits[SciFi::Tracking::max_x_hits],
-  float coordX[SciFi::Tracking::max_x_hits],
-  SciFi::Tracking::Track candidate_tracks[SciFi::Tracking::max_candidate_tracks],
-  int& n_candidate_tracks,
-  const float zRef_track,
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  const MiniState& velo_state,
-  SciFi::Tracking::HitSearchCuts& pars,
-  const SciFi::Tracking::Arrays* constArrays,
-  int side,
-  const bool secondLoop)
-{
-  if (n_x_hits < pars.minXHits) return;
-  if (secondLoop)
-    if (n_candidate_tracks >= SciFi::Tracking::max_tracks_second_loop) return;
-  if (!secondLoop)
-    if (n_candidate_tracks >= SciFi::Tracking::max_candidate_tracks) return;
-
-  int itEnd = n_x_hits;
-  const float xTrack = evalCubicParameterization(xParams_seed, SciFi::Tracking::zReference);
-  int it1 = 0;
-  int it2 = 0;
-  pars.minStereoHits = 0;
-
-  PlaneCounter planeCounter;
-
-  while (it1 < itEnd) {
-    // find next unused Hits
-    assert(it1 < n_x_hits);
-    while (it1 + pars.minXHits - 1 < itEnd && usedHits[it1]) {
-      ++it1;
-    }
-
-    it2 = it1 + pars.minXHits;
-    if (it2 > itEnd) break;
-    assert(it2 - 1 < n_x_hits);
-    while (it2 <= itEnd && usedHits[it2 - 1])
-      ++it2;
-    if (it2 > itEnd) break;
-
-    // Second part of 1D Hough transform:
-    // find a cluster of x positions on the reference plane that are close to each other
-    // TODO better xWindow calculation?? how to tune this???
-    const float xWindow = pars.maxXWindow + (fabsf(coordX[it1]) + fabsf(coordX[it1] - xTrack)) * pars.maxXWindowSlope;
-
-    if ((coordX[it2 - 1] - coordX[it1]) > xWindow) {
-      ++it1;
-      continue;
-    }
-
-    // find out which planes are present in the cluster
-    countPlanesOfXHits(planeCounter, it1, it2, n_x_hits, allXHits, usedHits, scifi_hits);
-
-    // Improve cluster (at the moment only add hits to the right)
-    // to recover inefficiencies from requiring a stereo hit to
-    // be matched to an x-hit in the collectXHits step
-    improveXCluster(it2, it1, itEnd, n_x_hits, usedHits, coordX, xWindow, pars, planeCounter, allXHits, scifi_hits);
-
-    // if not enough different planes, start again from the very beginning with next right hit
-    if (planeCounter.nbDifferent < pars.minXHits) {
-      ++it1;
-      continue;
-    }
-
-    //  Now we have a (rather) clean candidate, do best hit selection
-    SciFi::Tracking::LineFitterPars lineFitParameters;
-    lineFitParameters.m_z0 = SciFi::Tracking::zReference;
-    float xAtRef = 0.;
-    const unsigned int nbSingle = planeCounter.nbSingle();
-    int coordToFit[SciFi::Tracking::max_coordToFit];
-    int n_coordToFit = 0;
-    // 1) If there are enough planes with only one hit, do a straight
-    //    line fit through these and select good matching others
-    if (nbSingle >= SciFi::Tracking::minSingleHits && nbSingle != planeCounter.nbDifferent) {
-      // 1) we have enough single planes (thus two) to make a straight line fit
-      int otherHits[SciFi::Constants::n_layers][SciFi::Tracking::max_other_hits] = {0};
-      int nOtherHits[SciFi::Constants::n_layers] = {0};
-
-      // fit hits on planes with only one hit
-      // save the other hits in separate array
-      fitHitsFromSingleHitPlanes(
-        it1,
-        it2,
-        usedHits,
-        scifi_hits,
-        allXHits,
-        n_x_hits,
-        planeCounter,
-        lineFitParameters,
-        coordX,
-        otherHits,
-        nOtherHits);
-
-      // select best other hits (only best other hit is enough!)
-      // include them in fit
-      addAndFitHitsFromMultipleHitPlanes(nOtherHits, lineFitParameters, scifi_hits, coordX, allXHits, otherHits);
-
-      xAtRef = lineFitParameters.m_c0;
-    }
-    //  2) Try to find a cluster on the reference plane with a maximal
-    //     spread and from a minimum # of different planes
-    else {
-      // 2) Try to find a small distance containing at least 5(4) different planes
-      //    Most of the time do nothing
-      const unsigned int nPlanes = fminf(planeCounter.nbDifferent, uint {5});
-      int itWindowStart = it1;
-      int itWindowEnd = it1 + nPlanes; // pointing at last+1
-      // Hit is used, go to next unused one
-      while (itWindowEnd <= it2 && usedHits[itWindowEnd - 1] && (itWindowEnd - 1) < n_x_hits)
-        ++itWindowEnd;
-      int best = itWindowStart;
-      int bestEnd = itWindowEnd;
-
-      PlaneCounter lplaneCounter;
-      countUnusedXHitsOnPlanes(lplaneCounter, itWindowStart, itWindowEnd, n_x_hits, allXHits, usedHits, scifi_hits);
-
-      // modify itWindowStart and itWindowEnd while adding more x hits
-      // set best and bestEnd to include cluster of x hits on reference plane within minInterval
-      float minInterval = 1.f;
-      addXHitsForCandidateWithTooFewPlanes(
-        itWindowStart,
-        itWindowEnd,
-        it2,
-        itEnd,
-        minInterval,
-        lplaneCounter,
-        nPlanes,
-        coordX,
-        best,
-        bestEnd,
-        usedHits,
-        n_x_hits,
-        allXHits,
-        scifi_hits);
-
-      // TODO tune minInterval cut value
-      if (minInterval < 1.f) {
-        it1 = best;
-        it2 = bestEnd;
-      }
-
-      // Fill coordToFit and compute xAtRef (average x position on reference plane)
-      collectXHitsToFit(it1, it2, n_x_hits, allXHits, usedHits, coordToFit, n_coordToFit, coordX, xAtRef);
-
-    } // end of magical second part
-    //=== We have a candidate :)
-
-    planeCounter.clear();
-    for (int j = 0; j < n_coordToFit; ++j) {
-      planeCounter.addHit(scifi_hits.planeCode(coordToFit[j]) / 2);
-    }
-
-    // Only unused(!) hits in coordToFit now
-    bool ok = planeCounter.nbDifferent > 3;
-    float trackParameters[SciFi::Tracking::nTrackParams];
-    if (ok) {
-      getTrackParameters(xAtRef, velo_state, constArrays, trackParameters);
-      // Track described by cubic function in (z-zRef), but only first two terms are adjusted
-      // during fitting procedure -> linear fit
-      // nb: usedHits are currently not un-marked when removing ourliers
-      fastLinearFit(scifi_hits, trackParameters, coordToFit, n_coordToFit, planeCounter, pars);
-      // to do: do we have to mark these as used as well?
-      addHitsOnEmptyXLayers(
-        scifi_hits,
-        scifi_hit_count,
-        trackParameters,
-        xParams_seed,
-        yParams_seed,
-        false,
-        coordToFit,
-        n_coordToFit,
-        constArrays,
-        planeCounter,
-        pars,
-        side);
-
-      ok = planeCounter.nbDifferent > 3;
-    }
-    // == Fit and remove hits...
-    // track described by cubic function, but fitting only first three parameters here
-    // -> quadratic fit
-    // to do: remove outlier from usedHits
-    if (ok) ok = quadraticFitX(scifi_hits, trackParameters, coordToFit, n_coordToFit, planeCounter, pars);
-    if (ok) ok = trackParameters[7] / trackParameters[8] < SciFi::Tracking::maxChi2PerDoF;
-    if (ok)
-      ok = addHitsOnEmptyXLayers(
-        scifi_hits,
-        scifi_hit_count,
-        trackParameters,
-        xParams_seed,
-        yParams_seed,
-        true,
-        coordToFit,
-        n_coordToFit,
-        constArrays,
-        planeCounter,
-        pars,
-        side);
-    if (ok) {
-      // save track properties in track object
-      SciFi::Tracking::Track track;
-      track.chi2 = trackParameters[7];
-      for (int k = 0; k < 7; ++k) {
-        track.trackParams[k] = trackParameters[k];
-      }
-
-      for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-        int hit = coordToFit[i_hit];
-        assert(track.hitsNum < SciFi::Tracking::max_scifi_hits);
-        track.addHit(hit);
-      }
-      if (!secondLoop) {
-        assert(n_candidate_tracks < SciFi::Tracking::max_candidate_tracks);
-        candidate_tracks[n_candidate_tracks++] = track;
-      }
-      else if (secondLoop) {
-        assert(n_candidate_tracks < SciFi::Tracking::max_tracks_second_loop);
-        candidate_tracks[n_candidate_tracks++] = track;
-      }
-    }
-    if (secondLoop) {
-      if (n_candidate_tracks >= SciFi::Tracking::max_tracks_second_loop) break;
-      assert(n_candidate_tracks < SciFi::Tracking::max_tracks_second_loop);
-    }
-    else if (!secondLoop) {
-      if (n_candidate_tracks >= SciFi::Tracking::max_candidate_tracks) break;
-      assert(n_candidate_tracks < SciFi::Tracking::max_candidate_tracks);
-    }
-
-    ++it1;
-  }
-}
-
-__host__ __device__ bool addHitsOnEmptyXLayers(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  bool fullFit,
-  int coordToFit[SciFi::Tracking::max_coordToFit],
-  int& n_coordToFit,
-  const SciFi::Tracking::Arrays* constArrays,
-  PlaneCounter& planeCounter,
-  SciFi::Tracking::HitSearchCuts& pars,
-  int side)
-{
-  // is there an empty plane? otherwise skip here!
-  if (planeCounter.nbDifferent > 11) return true;
-  bool added = false;
-  const float x1 = trackParameters[0]; // mean xRef of this candidate
-  const float xAtRefFromSeed = evalCubicParameterization(xParams_seed, SciFi::Tracking::zReference);
-  const float xWindow = pars.maxXWindow + (fabsf(x1) + fabsf(x1 - xAtRefFromSeed)) * pars.maxXWindowSlope;
-
-  int iZoneStartingPoint = side > 0 ? constArrays->zoneoffsetpar : 0;
-
-  for (unsigned int iZone = iZoneStartingPoint; iZone < iZoneStartingPoint + constArrays->zoneoffsetpar; iZone++) {
-    assert(constArrays->xZones[iZone] / 2 < SciFi::Constants::n_layers);
-    if (planeCounter.nbInPlane(constArrays->xZones[iZone] / 2) != 0) continue;
-
-    const float parsX[4] = {trackParameters[0], trackParameters[1], trackParameters[2], trackParameters[3]};
-
-    assert(iZone - iZoneStartingPoint < SciFi::Constants::n_zones);
-    const float zZone = constArrays->xZone_zPos[iZone - iZoneStartingPoint];
-    // predicted x position on this plane based on current candidate
-    const float xPred = evalCubicParameterization(parsX, zZone);
-    const float minX = xPred - xWindow;
-    const float maxX = xPred + xWindow;
-
-    // -- Use a search to find the lower bound of the range of x values
-    assert(constArrays->xZones[iZone] < SciFi::Constants::n_zones);
-    int x_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->xZones[iZone]);
-    int x_zone_offset_end = x_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->xZones[iZone]);
-    int itH = getLowerBound(scifi_hits.x0, minX, x_zone_offset_begin, x_zone_offset_end);
-    int itEnd = x_zone_offset_end;
-
-    int best = findBestXHitOnEmptyLayer(itEnd, itH, scifi_hits, maxX, xPred);
-
-    if (best > -1) {
-      if (n_coordToFit >= SciFi::Tracking::max_coordToFit) break;
-      assert(n_coordToFit < SciFi::Tracking::max_coordToFit);
-      coordToFit[n_coordToFit++] = best; // add the best hit here
-      planeCounter.addHit(scifi_hits.planeCode(best) / 2);
-      added = true;
-    }
-  }
-  if (!added) return true;
-  if (fullFit) {
-    return quadraticFitX(scifi_hits, trackParameters, coordToFit, n_coordToFit, planeCounter, pars);
-  }
-  fastLinearFit(scifi_hits, trackParameters, coordToFit, n_coordToFit, planeCounter, pars);
-  return true;
-}
diff --git a/cuda/SciFi/PrForward/src/HitUtils.cu b/cuda/SciFi/PrForward/src/HitUtils.cu
deleted file mode 100644
index d0de6e2d17de705428994ef1308cdec16a6f5014..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/HitUtils.cu
+++ /dev/null
@@ -1,324 +0,0 @@
-#include "HitUtils.cuh"
-
-// match stereo hits to x hits
-__host__ __device__ bool matchStereoHit(
-  const int itUV1,
-  const int uv_zone_offset_end,
-  const SciFi::Hits& scifi_hits,
-  const float xMinUV,
-  const float xMaxUV)
-{
-  for (int stereoHit = itUV1; stereoHit != uv_zone_offset_end; ++stereoHit) {
-    if (scifi_hits.x0[stereoHit] > xMinUV) {
-      return (scifi_hits.x0[stereoHit] < xMaxUV);
-    }
-  }
-  return false;
-}
-
-// match stereo hits to x hits using triangle method
-__host__ __device__ bool matchStereoHitWithTriangle(
-  const int itUV2,
-  const int triangle_zone_offset_end,
-  const float yInZone,
-  const SciFi::Hits& scifi_hits,
-  const float xMinUV,
-  const float xMaxUV,
-  const int side)
-{
-
-  for (int stereoHit = itUV2; stereoHit != triangle_zone_offset_end; ++stereoHit) {
-    if (scifi_hits.x0[stereoHit] > xMinUV) {
-      // Triangle search condition depends on side
-      if (side > 0) { // upper
-        if (scifi_hits.yMax(stereoHit) > yInZone - SciFi::Tracking::yTolUVSearch) {
-          return true;
-        }
-      }
-      else { // lower
-        if (scifi_hits.yMin(stereoHit) < yInZone + SciFi::Tracking::yTolUVSearch) {
-          return true;
-        }
-      }
-    }
-  }
-  return false;
-}
-
-__host__ __device__ void removeOutlier(
-  const SciFi::Hits& scifi_hits,
-  PlaneCounter& planeCounter,
-  int* coordToFit,
-  int& n_coordToFit,
-  const int worst)
-{
-  planeCounter.removeHit(scifi_hits.planeCode(worst) / 2);
-  int coordToFit_temp[SciFi::Tracking::max_stereo_hits];
-  int i_hit_temp = 0;
-  for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-    int hit = coordToFit[i_hit];
-    if (hit != worst) coordToFit_temp[i_hit_temp++] = hit;
-  }
-  n_coordToFit = i_hit_temp;
-  for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-    coordToFit[i_hit] = coordToFit_temp[i_hit];
-  }
-}
-
-__host__ __device__ void countPlanesOfXHits(
-  PlaneCounter& planeCounter,
-  const int it1,
-  const int it2,
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits)
-{
-
-  planeCounter.clear();
-  for (int itH = it1; itH != it2; ++itH) {
-    assert(itH < n_x_hits);
-    if (!usedHits[itH]) {
-      const int plane = scifi_hits.planeCode(allXHits[itH]) / 2;
-      planeCounter.addHit(plane);
-    }
-  }
-}
-
-__host__ __device__ void countUnusedXHitsOnPlanes(
-  PlaneCounter& lplaneCounter,
-  const int itWindowStart,
-  const int itWindowEnd,
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits)
-{
-  for (int itH = itWindowStart; itH != itWindowEnd; ++itH) {
-    assert(itH < n_x_hits);
-    if (!usedHits[itH]) {
-      lplaneCounter.addHit(scifi_hits.planeCode(allXHits[itH]) / 2);
-    }
-  }
-}
-
-__host__ __device__ void addXHitsForCandidateWithTooFewPlanes(
-  int& itWindowStart,
-  int& itWindowEnd,
-  const int it2,
-  const int itEnd,
-  float& minInterval,
-  PlaneCounter& lplaneCounter,
-  const int nPlanes,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  int& best,
-  int& bestEnd,
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits)
-{
-
-  while (itWindowEnd <= it2) {
-    if (lplaneCounter.nbDifferent >= nPlanes) {
-      // have nPlanes, check x distance
-      assert(itWindowEnd - 1 < SciFi::Tracking::max_x_hits);
-      assert(itWindowStart < SciFi::Tracking::max_x_hits);
-      const float dist = coordX[itWindowEnd - 1] - coordX[itWindowStart];
-      if (dist < minInterval) {
-        minInterval = dist;
-        best = itWindowStart;
-        bestEnd = itWindowEnd;
-      }
-    }
-    else {
-      // too few planes, add one hit
-      ++itWindowEnd;
-      if (itWindowEnd > it2) break;
-      assert(itWindowEnd <= n_x_hits);
-      while (itWindowEnd <= it2 && usedHits[itWindowEnd - 1] && itWindowEnd <= n_x_hits)
-        ++itWindowEnd;
-      lplaneCounter.addHit(scifi_hits.planeCode(allXHits[itWindowEnd - 1]) / 2);
-      continue;
-    }
-    // move on to the right
-
-    lplaneCounter.removeHit(scifi_hits.planeCode(allXHits[itWindowStart]) / 2);
-    ++itWindowStart;
-    assert(itWindowStart < itEnd);
-    while (itWindowStart < itWindowEnd && usedHits[itWindowStart] && itWindowStart < n_x_hits)
-      ++itWindowStart;
-    // last hit guaranteed to be not used. Therefore there is always at least one hit to go to. No additional if
-    // required.
-  }
-}
-
-__host__ __device__ void collectXHitsToFit(
-  const int it1,
-  const int it2,
-  const int n_x_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  bool usedHits[SciFi::Tracking::max_x_hits],
-  int coordToFit[SciFi::Tracking::max_x_hits],
-  int& n_coordToFit,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  float& xAtRef)
-{
-
-  for (int itH = it1; it2 != itH; ++itH) {
-    assert(itH < n_x_hits);
-    if (!usedHits[itH]) {
-      if (n_coordToFit >= SciFi::Tracking::max_coordToFit) break;
-      coordToFit[n_coordToFit++] = allXHits[itH];
-      usedHits[itH] = true;
-      xAtRef += coordX[itH];
-    }
-  }
-  xAtRef /= ((float) n_coordToFit);
-}
-
-__host__ __device__ int findBestXHitOnEmptyLayer(
-  const int itEnd,
-  const int itBegin,
-  const SciFi::Hits& scifi_hits,
-  const float maxX,
-  const float xPred)
-{
-
-  float bestChi2 = 1.e9f;
-  int best = -1;
-  for (int itH = itBegin; itEnd != itH; ++itH) {
-    if (scifi_hits.x0[itH] > maxX) break;
-    const float d = scifi_hits.x0[itH] - xPred; // fast distance good enough at this point (?!)
-    const float chi2 = d * d * scifi_hits.w(itH);
-    if (chi2 < bestChi2) {
-      bestChi2 = chi2;
-      best = itH;
-    }
-  }
-  return best;
-}
-
-__host__ __device__ void findStereoHitsWithinXTol(
-  const int itBegin,
-  const int itEnd,
-  const SciFi::Hits& scifi_hits,
-  const float yZone,
-  const float xPred,
-  const float dxTol,
-  const bool triangleSearch,
-  const float dxDySign,
-  int& n_stereoHits,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits])
-{
-
-  for (int itH = itBegin; itEnd != itH; ++itH) {
-    const float dx = scifi_hits.x0[itH] + yZone * scifi_hits.dxdy(itH) - xPred;
-    if (dx > dxTol) break;
-    if (triangleSearch) {
-      if (yZone > scifi_hits.yMax(itH) + SciFi::Tracking::yTolUVSearch) continue;
-      if (yZone < scifi_hits.yMin(itH) - SciFi::Tracking::yTolUVSearch) continue;
-    }
-    if (n_stereoHits >= SciFi::Tracking::max_stereo_hits) break;
-    assert(n_stereoHits < SciFi::Tracking::max_stereo_hits);
-    stereoHits[n_stereoHits] = itH;
-    stereoCoords[n_stereoHits++] = dx * dxDySign;
-  }
-}
-
-// find cluster of stereo hits with certain number of hits from different planes
-// and similar dx value (stored in stereoCoords)
-__host__ __device__ void findStereoHitClusterByDx(
-  PlaneCounter& planeCounter,
-  int& endRange,
-  const SciFi::Tracking::HitSearchCuts& pars,
-  float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  const int n_stereoHits,
-  const SciFi::Hits& scifi_hits,
-  float& sumCoord,
-  int& first_hit)
-{
-
-  while (planeCounter.nbDifferent < pars.minStereoHits ||
-         stereoCoords[endRange] < stereoCoords[first_hit] + SciFi::Tracking::minYGap && endRange < n_stereoHits - 1) {
-    planeCounter.addHit(scifi_hits.planeCode(stereoHits[endRange]) / 2);
-    sumCoord += stereoCoords[endRange];
-    ++endRange;
-    if (endRange == n_stereoHits - 1) break;
-    first_hit = endRange - 1;
-  }
-}
-
-// - remove hits on planes with more than one hit if
-//   the hit is farthest away from the mean
-// - add hit if no hit on that plane is present in the cluster yet
-//   and if the cluster spread is reduced by adding the hit
-__host__ __device__ void cleanStereoHitCluster(
-  int& beginRange,
-  int& endRange,
-  const int n_stereoHits,
-  const int stereoHits[SciFi::Tracking::max_stereo_hits],
-  const float stereoCoords[SciFi::Tracking::max_stereo_hits],
-  float& sumCoord,
-  PlaneCounter& planeCounter,
-  const SciFi::Hits& scifi_hits)
-{
-
-  while (endRange < n_stereoHits - 1) {
-    const float averageCoord = sumCoord / float(endRange - beginRange);
-
-    // remove first if not single and farthest from mean
-    if (
-      planeCounter.nbInPlane(scifi_hits.planeCode(stereoHits[beginRange]) / 2) > 1 &&
-      ((averageCoord - stereoCoords[beginRange]) > 1.0f * (stereoCoords[endRange - 1] - averageCoord))) {
-
-      planeCounter.removeHit(scifi_hits.planeCode(stereoHits[beginRange]) / 2);
-      sumCoord -= stereoCoords[beginRange];
-      beginRange++;
-      continue;
-    }
-
-    if (endRange == n_stereoHits - 1) break; // already at end, cluster cannot be expanded anymore
-    // add next, if it decreases the range size and is empty
-    if ((planeCounter.nbInPlane(scifi_hits.planeCode(stereoHits[beginRange]) / 2) == 0)) {
-      if ((averageCoord - stereoCoords[beginRange] > stereoCoords[endRange] - averageCoord)) {
-        planeCounter.addHit(scifi_hits.planeCode(stereoHits[endRange]) / 2);
-        sumCoord += stereoCoords[endRange];
-        endRange++;
-        continue;
-      }
-    }
-
-    break;
-  }
-}
-
-__host__ __device__ int findBestStereoHitOnEmptyLayer(
-  const int itBegin,
-  const int itEnd,
-  const SciFi::Hits& scifi_hits,
-  const float yZone,
-  const float xPred,
-  const float dxTol,
-  const bool triangleSearch)
-{
-
-  int best = -1;
-  float bestChi2 = SciFi::Tracking::maxChi2Stereo;
-  for (int itH = itBegin; itEnd != itH; ++itH) {
-    const float dx = scifi_hits.x0[itH] + yZone * scifi_hits.dxdy(itH) - xPred;
-    if (dx > dxTol) break;
-    if (triangleSearch) {
-      if (yZone > scifi_hits.yMax(itH) + SciFi::Tracking::yTolUVSearch) continue;
-      if (yZone < scifi_hits.yMin(itH) - SciFi::Tracking::yTolUVSearch) continue;
-    }
-    const float chi2 = dx * dx * scifi_hits.w(itH);
-    if (chi2 < bestChi2) {
-      bestChi2 = chi2;
-      best = itH;
-    }
-  }
-  return best;
-}
diff --git a/cuda/SciFi/PrForward/src/LinearFitting.cu b/cuda/SciFi/PrForward/src/LinearFitting.cu
deleted file mode 100644
index f8d16309a97530b3777061cb5b5d805225a649e1..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/LinearFitting.cu
+++ /dev/null
@@ -1,179 +0,0 @@
-#include "LinearFitting.cuh"
-
-/**
-   Functions related to fitting a straight line
- */
-
-__host__ __device__ float getLineFitDistance(
-  const SciFi::Tracking::LineFitterPars& parameters,
-  const SciFi::Hits& scifi_hits,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int it)
-{
-  return coordX[it] - (parameters.m_c0 + (scifi_hits.z0[allXHits[it]] - parameters.m_z0) * parameters.m_tc);
-}
-
-__host__ __device__ float getLineFitChi2(
-  const SciFi::Tracking::LineFitterPars& parameters,
-  const SciFi::Hits& scifi_hits,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int it)
-{
-  float d = getLineFitDistance(parameters, scifi_hits, coordX, allXHits, it);
-  return d * d * coordX[it];
-}
-
-__host__ __device__ void solveLineFit(SciFi::Tracking::LineFitterPars& parameters)
-{
-  float den = (parameters.m_sz * parameters.m_sz - parameters.m_s0 * parameters.m_sz2);
-  parameters.m_c0 = (parameters.m_scz * parameters.m_sz - parameters.m_sc * parameters.m_sz2) / den;
-  parameters.m_tc = (parameters.m_sc * parameters.m_sz - parameters.m_s0 * parameters.m_scz) / den;
-}
-
-__host__ __device__ void incrementLineFitParameters(
-  SciFi::Tracking::LineFitterPars& parameters,
-  const SciFi::Hits& scifi_hits,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int it)
-{
-  float c = coordX[it];
-  const int hit = allXHits[it];
-  float w = scifi_hits.w(hit);
-  float z = scifi_hits.z0[hit] - parameters.m_z0;
-  parameters.m_s0 += w;
-  parameters.m_sz += w * z;
-  parameters.m_sz2 += w * z * z;
-  parameters.m_sc += w * c;
-  parameters.m_scz += w * c * z;
-}
-
-__host__ __device__ void fitHitsFromSingleHitPlanes(
-  const int it1,
-  const int it2,
-  const bool usedHits[SciFi::Tracking::max_x_hits],
-  const SciFi::Hits& scifi_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int n_x_hits,
-  const PlaneCounter planeCounter,
-  SciFi::Tracking::LineFitterPars& lineFitParameters,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  int otherHits[SciFi::Constants::n_layers][SciFi::Tracking::max_other_hits],
-  int nOtherHits[SciFi::Constants::n_layers])
-{
-
-  for (auto itH = it1; it2 > itH; ++itH) {
-    assert(itH < n_x_hits);
-    if (usedHits[itH]) continue;
-    int planeCode = scifi_hits.planeCode(allXHits[itH]) / 2;
-    if (planeCounter.nbInPlane(planeCode) == 1) {
-      incrementLineFitParameters(lineFitParameters, scifi_hits, coordX, allXHits, itH);
-    }
-    else {
-      if (nOtherHits[planeCode] < SciFi::Tracking::max_other_hits) {
-        assert(nOtherHits[planeCode] < SciFi::Tracking::max_other_hits);
-        otherHits[planeCode][nOtherHits[planeCode]++] = itH;
-      }
-    }
-  }
-  solveLineFit(lineFitParameters);
-}
-
-__host__ __device__ void addAndFitHitsFromMultipleHitPlanes(
-  const int nOtherHits[SciFi::Constants::n_layers],
-  SciFi::Tracking::LineFitterPars& lineFitParameters,
-  const SciFi::Hits& scifi_hits,
-  const float coordX[SciFi::Tracking::max_x_hits],
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int otherHits[SciFi::Constants::n_layers][SciFi::Tracking::max_other_hits])
-{
-
-  for (int i = 0; i < SciFi::Constants::n_layers; i++) {
-    if (nOtherHits[i] == 0) continue;
-
-    float bestChi2 = 1e9f;
-    int best = -1;
-    for (int hit = 0; hit < nOtherHits[i]; ++hit) {
-      const float chi2 = getLineFitChi2(lineFitParameters, scifi_hits, coordX, allXHits, otherHits[i][hit]);
-      if (chi2 < bestChi2) {
-        bestChi2 = chi2;
-        best = hit;
-      }
-    }
-    if (best > -1) {
-      incrementLineFitParameters(lineFitParameters, scifi_hits, coordX, allXHits, otherHits[i][best]);
-    }
-  }
-  solveLineFit(lineFitParameters);
-}
-
-// the track parameterization is cubic in (z-zRef),
-// however only the first two parametres are varied in this fit
-// -> this is a linear fit
-__host__ __device__ void fastLinearFit(
-  const SciFi::Hits& scifi_hits,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  int coordToFit[SciFi::Tracking::max_coordToFit],
-  int& n_coordToFit,
-  PlaneCounter planeCounter,
-  SciFi::Tracking::HitSearchCuts& pars)
-{
-  // re-do fit every time an outlier was removed
-  bool fit = true;
-  while (fit) {
-    //== Fit a line
-    float s0 = 0.f;
-    float sz = 0.f;
-    float sz2 = 0.f;
-    float sd = 0.f;
-    float sdz = 0.f;
-
-    for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-      int hit = coordToFit[i_hit];
-      const float parsX[4] = {trackParameters[0], trackParameters[1], trackParameters[2], trackParameters[3]};
-      const float zHit = scifi_hits.z0[hit];
-      float track_x_at_zHit = evalCubicParameterization(parsX, zHit);
-      const float d = scifi_hits.x0[hit] - track_x_at_zHit;
-      const float w = scifi_hits.w(hit);
-      const float z = zHit - SciFi::Tracking::zReference;
-      s0 += w;
-      sz += w * z;
-      sz2 += w * z * z;
-      sd += w * d;
-      sdz += w * d * z;
-    }
-    float den = (sz * sz - s0 * sz2);
-    if (!(fabsf(den) > 1e-5f)) return;
-    const float da = (sdz * sz - sd * sz2) / den;
-    const float db = (sd * sz - s0 * sdz) / den;
-    trackParameters[0] += da;
-    trackParameters[1] += db;
-    fit = false;
-
-    if (n_coordToFit < pars.minXHits) return;
-
-    int worst = n_coordToFit;
-    float maxChi2 = 0.f;
-    const bool notMultiple = planeCounter.nbDifferent == n_coordToFit;
-    // TODO how many multiple hits do we normaly have?
-    // how often do we do the right thing here?
-    // delete two hits at same time?
-    for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-      int hit = coordToFit[i_hit];
-      // This could certainly be wrapped in some helper function with a lot
-      // of passing around or copying etc...
-      const float parsX[4] = {trackParameters[0], trackParameters[1], trackParameters[2], trackParameters[3]};
-      const float chi2 = chi2XHit(parsX, scifi_hits, hit);
-      if (chi2 > maxChi2 && (notMultiple || planeCounter.nbInPlane(scifi_hits.planeCode(hit) / 2) > 1)) {
-        maxChi2 = chi2;
-        worst = i_hit;
-      }
-    }
-    if (maxChi2 > SciFi::Tracking::maxChi2LinearFit || (!notMultiple && maxChi2 > 4.f)) {
-      removeOutlier(scifi_hits, planeCounter, coordToFit, n_coordToFit, coordToFit[worst]);
-      fit = true;
-    }
-  }
-}
diff --git a/cuda/SciFi/PrForward/src/ParabolaFitting.cu b/cuda/SciFi/PrForward/src/ParabolaFitting.cu
deleted file mode 100644
index e36c0d7afad95291f804dd879f1e84b78fae2a88..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/ParabolaFitting.cu
+++ /dev/null
@@ -1,59 +0,0 @@
-#include "ParabolaFitting.cuh"
-
-__host__ __device__ int fitParabola(
-  int* coordToFit,
-  const int n_coordToFit,
-  const SciFi::Hits& scifi_hits,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  const bool xFit)
-{
-
-  //== Fit a cubic
-  float s0 = 0.f;
-  float sz = 0.f;
-  float sz2 = 0.f;
-  float sz3 = 0.f;
-  float sz4 = 0.f;
-  float sd = 0.f;
-  float sdz = 0.f;
-  float sdz2 = 0.f;
-
-  for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-    int hit = coordToFit[i_hit];
-    float d = trackToHitDistance(trackParameters, scifi_hits, hit);
-    if (!xFit) d *= -1.f / scifi_hits.dxdy(hit); // TODO multiplication much faster than division!
-    float w = scifi_hits.w(hit);
-    float z = .001f * (scifi_hits.z0[hit] - SciFi::Tracking::zReference);
-    s0 += w;
-    sz += w * z;
-    sz2 += w * z * z;
-    sz3 += w * z * z * z;
-    sz4 += w * z * z * z * z;
-    sd += w * d;
-    sdz += w * d * z;
-    sdz2 += w * d * z * z;
-  }
-  const float b1 = sz * sz - s0 * sz2;
-  const float c1 = sz2 * sz - s0 * sz3;
-  const float d1 = sd * sz - s0 * sdz;
-  const float b2 = sz2 * sz2 - sz * sz3;
-  const float c2 = sz3 * sz2 - sz * sz4;
-  const float d2 = sdz * sz2 - sz * sdz2;
-  const float den = (b1 * c2 - b2 * c1);
-  if (!(fabsf(den) > 1e-5f)) return false;
-  const float db = (d1 * c2 - d2 * c1) / den;
-  const float dc = (d2 * b1 - d1 * b2) / den;
-  const float da = (sd - db * sz - dc * sz2) / s0;
-  if (xFit) {
-    trackParameters[0] += da;
-    trackParameters[1] += db * 1.e-3f;
-    trackParameters[2] += dc * 1.e-6f;
-  }
-  else {
-    trackParameters[4] += da;
-    trackParameters[5] += db * 1.e-3f;
-    trackParameters[6] += dc * 1.e-6f;
-  }
-
-  return true;
-}
diff --git a/cuda/SciFi/PrForward/src/PrForwardTools.cu b/cuda/SciFi/PrForward/src/PrForwardTools.cu
deleted file mode 100644
index e51cb3fe5299323191ea55237881b7629a43b199..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/PrForwardTools.cu
+++ /dev/null
@@ -1,399 +0,0 @@
-#include "PrForwardTools.cuh"
-
-/* Look first in x layers, then in stereo layers for hits
-   do 1D Hough transform for x- and stereo hits
-   do global 1D Hough transform
-   use TMVAs to obtain track quality */
-__host__ __device__ void find_forward_tracks(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const float qop_ut,
-  const int i_veloUT_track,
-  SciFi::TrackHits* outputTracks,
-  uint* n_forward_tracks,
-  const int ut_event_number_of_tracks,
-  const SciFi::Tracking::TMVA* tmva1,
-  const SciFi::Tracking::TMVA* tmva2,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const MiniState& velo_state)
-{
-
-  // The LHCb framework code had a PT preselection for the VeloUT tracks
-  // here, which I am removing because this should be done explicitly through
-  // track selectors if we do it at all, not hacked inside the tracking code
-
-  const float zRef_track = SciFi::Tracking::zReference;
-  const float xAtRef = xFromVelo(zRef_track, velo_state);
-  const float xParams_seed[4] = {xAtRef, velo_state.tx, 0.f, 0.f};
-  const float yAtRef = yFromVelo(zRef_track, velo_state);
-  const float yParams_seed[4] = {yAtRef, velo_state.ty, 0.f, 0.f};
-
-  // First loop Hough cluster search, set initial search windows
-  SciFi::Tracking::HitSearchCuts pars_first {SciFi::Tracking::minXHits,
-                                             SciFi::Tracking::maxXWindow,
-                                             SciFi::Tracking::maxXWindowSlope,
-                                             SciFi::Tracking::maxXGap,
-                                             4u};
-  SciFi::Tracking::HitSearchCuts pars_second {SciFi::Tracking::minXHits_2nd,
-                                              SciFi::Tracking::maxXWindow_2nd,
-                                              SciFi::Tracking::maxXWindowSlope_2nd,
-                                              SciFi::Tracking::maxXGap_2nd,
-                                              4u};
-
-  int allXHits[2][SciFi::Tracking::max_x_hits];
-  int n_x_hits[2] = {0};
-  float coordX[2][SciFi::Tracking::max_x_hits];
-
-  if (yAtRef > -5.f)
-    collectAllXHits(
-      scifi_hits,
-      scifi_hit_count,
-      allXHits[1],
-      n_x_hits[1],
-      coordX[1],
-      xParams_seed,
-      yParams_seed,
-      constArrays,
-      magnet_polarity,
-      velo_state,
-      qop_ut,
-      1);
-  if (yAtRef < 5.f)
-    collectAllXHits(
-      scifi_hits,
-      scifi_hit_count,
-      allXHits[0],
-      n_x_hits[0],
-      coordX[0],
-      xParams_seed,
-      yParams_seed,
-      constArrays,
-      magnet_polarity,
-      velo_state,
-      qop_ut,
-      -1);
-
-  SciFi::Tracking::Track candidate_tracks[SciFi::Tracking::max_candidate_tracks];
-  int n_candidate_tracks = 0;
-  bool usedHits[2][SciFi::Tracking::max_x_hits] = {false};
-
-  if (yAtRef > -5.f)
-    selectXCandidates(
-      scifi_hits,
-      scifi_hit_count,
-      allXHits[1],
-      n_x_hits[1],
-      usedHits[1],
-      coordX[1],
-      candidate_tracks,
-      n_candidate_tracks,
-      zRef_track,
-      xParams_seed,
-      yParams_seed,
-      velo_state,
-      pars_first,
-      constArrays,
-      1,
-      false);
-  if (yAtRef < 5.f)
-    selectXCandidates(
-      scifi_hits,
-      scifi_hit_count,
-      allXHits[0],
-      n_x_hits[0],
-      usedHits[0],
-      coordX[0],
-      candidate_tracks,
-      n_candidate_tracks,
-      zRef_track,
-      xParams_seed,
-      yParams_seed,
-      velo_state,
-      pars_first,
-      constArrays,
-      -1,
-      false);
-
-  SciFi::Tracking::Track selected_tracks[SciFi::Tracking::max_selected_tracks];
-  int n_selected_tracks = 0;
-
-  selectFullCandidates(
-    scifi_hits,
-    scifi_hit_count,
-    candidate_tracks,
-    n_candidate_tracks,
-    selected_tracks,
-    n_selected_tracks,
-    xParams_seed,
-    yParams_seed,
-    velo_state,
-    qop_ut,
-    pars_first,
-    tmva1,
-    tmva2,
-    constArrays,
-    magnet_polarity,
-    false);
-
-  bool ok = false;
-  for (int i_track = 0; i_track < n_selected_tracks; ++i_track) {
-    if (selected_tracks[i_track].hitsNum > 10) ok = true;
-  }
-  assert(n_selected_tracks < SciFi::Tracking::max_selected_tracks);
-
-  SciFi::Tracking::Track candidate_tracks2[SciFi::Tracking::max_tracks_second_loop];
-  int n_candidate_tracks2 = 0;
-
-  if (!ok && SciFi::Tracking::secondLoop) { // If you found nothing begin the 2nd loop
-    if (yAtRef > -5.f)
-      selectXCandidates(
-        scifi_hits,
-        scifi_hit_count,
-        allXHits[1],
-        n_x_hits[1],
-        usedHits[1],
-        coordX[1],
-        candidate_tracks2,
-        n_candidate_tracks2,
-        zRef_track,
-        xParams_seed,
-        yParams_seed,
-        velo_state,
-        pars_second,
-        constArrays,
-        1,
-        true);
-    if (yAtRef < 5.f)
-      selectXCandidates(
-        scifi_hits,
-        scifi_hit_count,
-        allXHits[0],
-        n_x_hits[0],
-        usedHits[0],
-        coordX[0],
-        candidate_tracks2,
-        n_candidate_tracks2,
-        zRef_track,
-        xParams_seed,
-        yParams_seed,
-        velo_state,
-        pars_second,
-        constArrays,
-        -1,
-        true);
-
-    SciFi::Tracking::Track selected_tracks2[SciFi::Tracking::max_tracks_second_loop];
-    int n_selected_tracks2 = 0;
-
-    selectFullCandidates(
-      scifi_hits,
-      scifi_hit_count,
-      candidate_tracks2,
-      n_candidate_tracks2,
-      selected_tracks2,
-      n_selected_tracks2,
-      xParams_seed,
-      yParams_seed,
-      velo_state,
-      qop_ut,
-      pars_second,
-      tmva1,
-      tmva2,
-      constArrays,
-      magnet_polarity,
-      true);
-
-    for (int i_track = 0; i_track < n_selected_tracks2; ++i_track) {
-      assert(n_selected_tracks < SciFi::Tracking::max_selected_tracks);
-      selected_tracks[n_selected_tracks++] = selected_tracks2[i_track];
-    }
-
-    ok = (n_selected_tracks > 0);
-  }
-
-  if (ok || !SciFi::Tracking::secondLoop) {
-
-    if (n_selected_tracks > 1) {
-      // not using thrust::sort due to temporary_buffer::allocate:: get_temporary_buffer failed" error
-      // thrust::sort( thrust::seq, selected_tracks, selected_tracks + n_selected_tracks, lowerByQuality);
-      sort_tracks(selected_tracks, n_selected_tracks, [](SciFi::Tracking::Track t1, SciFi::Tracking::Track t2) {
-        if (t1.quality < t2.quality) return -1;
-        if (t1.quality == t2.quality) return 0;
-        return 1;
-      });
-    }
-
-    const uint event_hit_offset = scifi_hit_count.event_offset();
-    float minQuality = SciFi::Tracking::maxQuality;
-    for (int i_track = 0; i_track < n_selected_tracks; ++i_track) {
-      SciFi::Tracking::Track& track = selected_tracks[i_track];
-
-      // TODO: We would have to work out a way for the Prforward to have max 12 hits
-      if (track.hitsNum > 12) {
-        track.hitsNum = 12;
-      }
-
-      if (track.quality + SciFi::Tracking::deltaQuality < minQuality)
-        minQuality = track.quality + SciFi::Tracking::deltaQuality;
-      if (!(track.quality > minQuality)) {
-
-        SciFi::TrackHits tr = makeTrack(track);
-        tr.ut_track_index = i_veloUT_track;
-
-        // state at end of SciFi
-        const float z = SciFi::Constants::ZEndT;
-        MiniState scifi_state(track.x(z), track.y(z), z, track.xSlope(z), track.ySlope(z));
-
-        // add LHCbIDs from SciFi part of the track
-        for (int i_hit = 0; i_hit < track.hitsNum; ++i_hit) {
-          // save local hit index within event to be able to use short
-          const int local_hit_index = track.hit_indices[i_hit] - event_hit_offset;
-          tr.add_hit(local_hit_index);
-        }
-
-        if (*n_forward_tracks >= ut_event_number_of_tracks * SciFi::Constants::max_SciFi_tracks_per_UT_track)
-          printf("n_forward_tracks = %u \n", *n_forward_tracks);
-        assert(*n_forward_tracks < ut_event_number_of_tracks * SciFi::Constants::max_SciFi_tracks_per_UT_track);
-#ifndef __CUDA_ARCH__
-        outputTracks[(*n_forward_tracks)++] = tr;
-#else
-        uint n_tracks = atomicAdd(n_forward_tracks, 1);
-        assert(n_tracks < ut_event_number_of_tracks * SciFi::Constants::max_SciFi_tracks_per_UT_track);
-        outputTracks[n_tracks] = tr;
-#endif
-      }
-    }
-  }
-}
-
-// Turn SciFi::Tracking::Track into a SciFi::Track
-__host__ __device__ SciFi::TrackHits makeTrack(SciFi::Tracking::Track track)
-{
-  SciFi::TrackHits tr;
-  tr.qop = track.qop;
-  tr.quality = track.quality;
-
-  return tr;
-}
-
-//=========================================================================
-//  Create Full candidates out of xCandidates
-//  Searching for stereo hits
-//  Fit of all hits
-//  save everything in track candidate folder
-//=========================================================================
-__host__ __device__ void selectFullCandidates(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  SciFi::Tracking::Track* candidate_tracks,
-  int& n_candidate_tracks,
-  SciFi::Tracking::Track* selected_tracks,
-  int& n_selected_tracks,
-  const float xParams_seed[4],
-  const float yParams_seed[4],
-  MiniState velo_state,
-  const float VeloUT_qOverP,
-  SciFi::Tracking::HitSearchCuts& pars,
-  const SciFi::Tracking::TMVA* tmva1,
-  const SciFi::Tracking::TMVA* tmva2,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  const bool secondLoop)
-{
-
-  PlaneCounter planeCounter;
-  planeCounter.clear();
-  float mlpInput[7] = {0};
-
-  for (int i_track = 0; i_track < n_candidate_tracks; ++i_track) {
-    SciFi::Tracking::Track* cand = candidate_tracks + i_track;
-
-    pars.minStereoHits = 4;
-
-    if (cand->hitsNum + pars.minStereoHits < SciFi::Tracking::minTotalHits) {
-      pars.minStereoHits = SciFi::Tracking::minTotalHits - cand->hitsNum;
-    }
-
-    int stereoHits[SciFi::Tracking::max_stereo_hits];
-    int n_stereoHits = 0;
-    float stereoCoords[SciFi::Tracking::max_stereo_hits];
-    collectStereoHits(
-      scifi_hits, scifi_hit_count, *cand, velo_state, pars, constArrays, stereoCoords, stereoHits, n_stereoHits);
-
-    if (n_stereoHits < pars.minStereoHits) continue;
-
-    if (!selectStereoHits(
-          scifi_hits, scifi_hit_count, *cand, constArrays, stereoCoords, stereoHits, n_stereoHits, velo_state, pars))
-      continue;
-
-    planeCounter.clear();
-    for (int i_hit = 0; i_hit < cand->hitsNum; ++i_hit) {
-      int hit = cand->hit_indices[i_hit];
-      planeCounter.addHit(scifi_hits.planeCode(hit) / 2);
-    }
-
-    // make a fit of ALL hits using their x coordinate
-    if (!quadraticFitX(scifi_hits, cand->trackParams, cand->hit_indices, cand->hitsNum, planeCounter, pars)) continue;
-
-    // track has enough hits, calcualte quality and save if good enough
-    if (planeCounter.nbDifferent >= SciFi::Tracking::minTotalHits) {
-
-      const float qOverP = calcqOverP(cand->trackParams[1], constArrays, velo_state, magnet_polarity);
-      // orig params before fitting , TODO faster if only calc once?? mem usage?
-      const float xAtRef = cand->trackParams[0];
-      float dSlope = (velo_state.x + (SciFi::Tracking::zReference - velo_state.z) * velo_state.tx - xAtRef) /
-                     (SciFi::Tracking::zReference - constArrays->zMagnetParams[0]);
-      const float zMagSlope =
-        constArrays->zMagnetParams[2] * pow(velo_state.tx, 2) + constArrays->zMagnetParams[3] * pow(velo_state.ty, 2);
-      const float zMag = constArrays->zMagnetParams[0] + constArrays->zMagnetParams[1] * dSlope * dSlope + zMagSlope;
-      const float xMag = velo_state.x + (zMag - velo_state.z) * velo_state.tx;
-      const float slopeT = (xAtRef - xMag) / (SciFi::Tracking::zReference - zMag);
-      dSlope = slopeT - velo_state.tx;
-      const float dyCoef = dSlope * dSlope * velo_state.ty;
-
-      float bx = slopeT;
-      float ay = velo_state.y + (SciFi::Tracking::zReference - velo_state.z) * velo_state.ty;
-      float by = velo_state.ty + dyCoef * SciFi::Tracking::byParams;
-
-      // ay,by,bx params
-      const float ay1 = cand->trackParams[4];
-      const float by1 = cand->trackParams[5];
-      const float bx1 = cand->trackParams[1];
-
-      mlpInput[0] = planeCounter.nbDifferent;
-      mlpInput[1] = qOverP;
-      mlpInput[2] = VeloUT_qOverP - qOverP;                // veloUT - scifi
-      if (fabsf(VeloUT_qOverP) < 1e-9f) mlpInput[2] = 0.f; // no momentum estiamte
-      mlpInput[3] = pow(velo_state.tx, 2) + pow(velo_state.ty, 2);
-      mlpInput[4] = by - by1;
-      mlpInput[5] = bx - bx1;
-      mlpInput[6] = ay - ay1;
-
-      float quality = 0.f;
-      /// WARNING: if the NN classes straight out of TMVA are used, put a mutex here!
-      if (pars.minXHits > 4)
-        quality = GetMvaValue(mlpInput, tmva1); // 1st loop NN
-      else
-        quality = GetMvaValue(mlpInput, tmva2); // 2nd loop NN
-
-      quality = 1.f - quality; // backward compability
-
-      if (quality < SciFi::Tracking::maxQuality) {
-        cand->quality = quality;
-        cand->qop = qOverP;
-        if (!secondLoop)
-          assert(n_selected_tracks < SciFi::Tracking::max_selected_tracks);
-        else if (secondLoop)
-          assert(n_selected_tracks < SciFi::Tracking::max_tracks_second_loop);
-        selected_tracks[n_selected_tracks++] = *cand;
-        if (!secondLoop) {
-          if (n_selected_tracks >= SciFi::Tracking::max_selected_tracks) break;
-        }
-        else if (secondLoop) {
-          if (n_selected_tracks >= SciFi::Tracking::max_tracks_second_loop) break;
-        }
-      }
-    }
-  }
-}
diff --git a/cuda/SciFi/PrForward/src/ReferencePlaneProjection.cu b/cuda/SciFi/PrForward/src/ReferencePlaneProjection.cu
deleted file mode 100644
index 98328b72f59b6b37e30a73447e47d0cef9fcd7a1..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/ReferencePlaneProjection.cu
+++ /dev/null
@@ -1,40 +0,0 @@
-#include "ReferencePlaneProjection.cuh"
-
-// project x position of hits from one plane onto the reference plane
-// save it in coordX
-// in the c++ this is vectorized, undoing because no point before CUDA (but vectorization is obvious)
-__host__ __device__ void xAtRef_SamePlaneHits(
-  const SciFi::Hits& scifi_hits,
-  const int allXHits[SciFi::Tracking::max_x_hits],
-  const int n_x_hits,
-  float coordX[SciFi::Tracking::max_x_hits],
-  const float xParams_seed[4],
-  const SciFi::Tracking::Arrays* constArrays,
-  MiniState velo_state,
-  const float zMagSlope,
-  int itH,
-  int itEnd)
-{
-  // this is quite computationally expensive mind, should take care when porting
-  assert(itH < SciFi::Tracking::max_x_hits);
-  const float zHit = scifi_hits.z0[allXHits[itH]]; // all hits in same layer
-  const float xFromVelo_Hit = evalCubicParameterization(xParams_seed, zHit);
-  const float dSlopeDivPart = 1.f / (zHit - constArrays->zMagnetParams[0]);
-  const float dz = 1.e-3f * (zHit - SciFi::Tracking::zReference);
-
-  while (itEnd > itH) {
-    float xHit = scifi_hits.x0[allXHits[itH]];
-    // difference in slope before and after the kick
-    float dSlope = (xFromVelo_Hit - xHit) * dSlopeDivPart;
-    // update zMag now that dSlope is known
-    float zMag = zMagSlope + constArrays->zMagnetParams[1] * dSlope * dSlope;
-    float xMag = xFromVelo_Hit + velo_state.tx * (zMag - zHit);
-    // calculate x position on reference plane (save in coodX)
-    // dxCoef: account for additional bending of track due to fringe field in first station
-    // expressed by quadratic and cubic term in z
-    float dxCoef = dz * dz * (constArrays->xParams[0] + dz * constArrays->xParams[1]) * dSlope;
-    float ratio = (SciFi::Tracking::zReference - zMag) / (zHit - zMag);
-    coordX[itH] = xMag + ratio * (xHit + dxCoef - xMag);
-    itH++;
-  }
-}
diff --git a/cuda/SciFi/PrForward/src/TrackUtils.cu b/cuda/SciFi/PrForward/src/TrackUtils.cu
deleted file mode 100644
index 6b81f7d8e11f269ca6ff862b2648f6453adc41cc..0000000000000000000000000000000000000000
--- a/cuda/SciFi/PrForward/src/TrackUtils.cu
+++ /dev/null
@@ -1,269 +0,0 @@
-#include "TrackUtils.cuh"
-
-// extrapolate x position from given state to z
-__host__ __device__ float xFromVelo(const float z, const MiniState& velo_state)
-{
-  return velo_state.x + (z - velo_state.z) * velo_state.tx;
-}
-
-// extrapolate y position from given state to z
-__host__ __device__ float yFromVelo(const float z, const MiniState& velo_state)
-{
-  return velo_state.y + (z - velo_state.z) * velo_state.ty;
-}
-
-__host__ __device__ float evalCubicParameterization(const float params[4], float z)
-{
-  float dz = z - SciFi::Tracking::zReference;
-  return params[0] + (params[1] + (params[2] + params[3] * dz) * dz) * dz;
-}
-
-__host__ __device__ float evalParameterizationX(const float* params, float z)
-{
-  const float dz = z - SciFi::Tracking::zReference;
-  return params[0] + (params[1] + (params[2] + params[3] * dz) * dz) * dz;
-}
-
-__host__ __device__ float evalParameterizationY(const float* params, float z)
-{
-  const float dz = z - SciFi::Tracking::zReference;
-  return params[0] + (params[1] + params[2] * dz) * dz;
-}
-
-__host__ __device__ bool lowerByQuality(SciFi::Tracking::Track t1, SciFi::Tracking::Track t2)
-{
-  return t1.quality < t2.quality;
-}
-
-__host__ __device__ void getTrackParameters(
-  float xAtRef,
-  const MiniState& velo_state,
-  const SciFi::Tracking::Arrays* constArrays,
-  float trackParams[SciFi::Tracking::nTrackParams])
-{
-  float dSlope = (xFromVelo(SciFi::Tracking::zReference, velo_state) - xAtRef) /
-                 (SciFi::Tracking::zReference - constArrays->zMagnetParams[0]);
-  const float zMagSlope = constArrays->zMagnetParams[2] * velo_state.tx * velo_state.tx +
-                          constArrays->zMagnetParams[3] * velo_state.ty * velo_state.ty;
-  const float zMag = constArrays->zMagnetParams[0] + constArrays->zMagnetParams[1] * dSlope * dSlope + zMagSlope;
-  const float xMag = xFromVelo(zMag, velo_state);
-  const float slopeT = (xAtRef - xMag) / (SciFi::Tracking::zReference - zMag);
-  dSlope = slopeT - velo_state.tx;
-  const float dyCoef = dSlope * dSlope * velo_state.ty;
-
-  trackParams[0] = xAtRef;
-  trackParams[1] = slopeT;
-  trackParams[2] = 1.e-6f * constArrays->xParams[0] * dSlope;
-  trackParams[3] = 1.e-9f * constArrays->xParams[1] * dSlope;
-  trackParams[4] = yFromVelo(SciFi::Tracking::zReference, velo_state);
-  trackParams[5] = velo_state.ty + dyCoef * SciFi::Tracking::byParams;
-  trackParams[6] = dyCoef * SciFi::Tracking::cyParams;
-  trackParams[7] = 0.0f;
-  trackParams[8] = 0.0f; // last elements are chi2 and ndof, as float
-}
-
-__host__ __device__ float calcqOverP(
-  float bx,
-  const SciFi::Tracking::Arrays* constArrays,
-  const MiniState& velo_state,
-  const float magnet_polarity)
-{
-
-  float qop = 1.0f / Gaudi::Units::GeV;
-  const float bx2 = bx * bx;
-  const float ty2 = velo_state.ty * velo_state.ty;
-  const float coef =
-    (constArrays->momentumParams[0] + constArrays->momentumParams[1] * bx2 +
-     constArrays->momentumParams[2] * bx2 * bx2 + constArrays->momentumParams[3] * bx * velo_state.tx +
-     constArrays->momentumParams[4] * ty2 + constArrays->momentumParams[5] * ty2 * ty2);
-  const float tx2 = velo_state.tx * velo_state.tx;
-  float m_slope2 = tx2 + ty2;
-  float proj = sqrtf((1.f + m_slope2) / (1.f + tx2));
-  qop = (velo_state.tx - bx) / (coef * Gaudi::Units::GeV * proj * magnet_polarity);
-  return qop;
-}
-
-// Find z zMag position within the magnet at which the bending ("kick") occurs
-// this is parameterized based on MC
-// the second parameter([1]) is multiplied by the difference in slope before and
-// after the kick, this slope is calculated from zMag and the x position of the track
-// at the reference plane -> it is calculated iteratively later
-__host__ __device__ float zMagnet(const MiniState& velo_state, const SciFi::Tracking::Arrays* constArrays)
-{
-
-  return (
-    constArrays->zMagnetParams[0] + constArrays->zMagnetParams[2] * velo_state.tx * velo_state.tx +
-    constArrays->zMagnetParams[3] * velo_state.ty * velo_state.ty);
-}
-
-// calculate difference between straight line extrapolation and
-// where a track with wrongSignPT (2 GeV) would be on the reference plane (?)
-__host__ __device__ float calcDxRef(float pt, const MiniState& velo_state)
-{
-  const float tx2 = velo_state.tx * velo_state.tx;
-  const float ty2 = velo_state.ty * velo_state.ty;
-  float m_slope2 = tx2 + ty2;
-  return 3973000.f * sqrtf(m_slope2) / pt - 2200.f * ty2 - 1000.f * tx2; // tune this window
-}
-
-__host__ __device__ float
-trackToHitDistance(const float trackParameters[SciFi::Tracking::nTrackParams], const SciFi::Hits& scifi_hits, int hit)
-{
-  const float z_Hit =
-    scifi_hits.z0[hit] + scifi_hits.dzdy(hit) * evalParameterizationY(trackParameters + 4, scifi_hits.z0[hit]);
-  const float x_track = evalParameterizationX(trackParameters, z_Hit);
-  const float y_track = evalParameterizationY(trackParameters + 4, z_Hit);
-  return scifi_hits.x0[hit] + y_track * scifi_hits.dxdy(hit) - x_track;
-}
-
-__host__ __device__ float chi2XHit(const float parsX[4], const SciFi::Hits& scifi_hits, const int hit)
-{
-  float track_x_at_zHit = evalCubicParameterization(parsX, scifi_hits.z0[hit]);
-  float hitdist = scifi_hits.x0[hit] - track_x_at_zHit;
-  return hitdist * hitdist * scifi_hits.w(hit);
-}
-
-// the track parameterization is cubic in (z-zRef),
-// however only the first three parametres are varied in this fit only
-// -> this is a quadratic fit
-__host__ __device__ bool quadraticFitX(
-  const SciFi::Hits& scifi_hits,
-  float trackParameters[SciFi::Tracking::nTrackParams],
-  int coordToFit[SciFi::Tracking::max_coordToFit],
-  int& n_coordToFit,
-  PlaneCounter& planeCounter,
-  SciFi::Tracking::HitSearchCuts& pars)
-{
-
-  if (planeCounter.nbDifferent < pars.minXHits) return false;
-  bool doFit = true;
-  while (doFit) {
-
-    fitParabola(coordToFit, n_coordToFit, scifi_hits, trackParameters, true);
-
-    float maxChi2 = 0.f;
-    float totChi2 = 0.f;
-    int nDoF = -3; // fitted 3 parameters
-    const bool notMultiple = planeCounter.nbDifferent == n_coordToFit;
-
-    int worst = n_coordToFit;
-    for (int i_hit = 0; i_hit < n_coordToFit; ++i_hit) {
-      int hit = coordToFit[i_hit];
-      float d = trackToHitDistance(trackParameters, scifi_hits, hit);
-      float chi2 = d * d * scifi_hits.w(hit);
-      totChi2 += chi2;
-      ++nDoF;
-      if (chi2 > maxChi2 && (notMultiple || planeCounter.nbInPlane(scifi_hits.planeCode(hit) / 2) > 1)) {
-        maxChi2 = chi2;
-        worst = i_hit;
-      }
-    }
-    if (nDoF < 1) return false;
-    trackParameters[7] = totChi2;
-    trackParameters[8] = (float) nDoF;
-
-    if (worst == n_coordToFit) {
-      return true;
-    }
-    doFit = false;
-    if (totChi2 / nDoF > SciFi::Tracking::maxChi2PerDoF || maxChi2 > SciFi::Tracking::maxChi2XProjection) {
-      removeOutlier(scifi_hits, planeCounter, coordToFit, n_coordToFit, coordToFit[worst]);
-      if (planeCounter.nbDifferent < pars.minXHits + pars.minStereoHits) return false;
-      doFit = true;
-    }
-  }
-  return true;
-}
-
-__host__ __device__ bool fitYProjection(
-  const SciFi::Hits& scifi_hits,
-  SciFi::Tracking::Track& track,
-  int stereoHits[SciFi::Tracking::max_stereo_hits],
-  int& n_stereoHits,
-  PlaneCounter& planeCounter,
-  const MiniState& velo_state,
-  const SciFi::Tracking::Arrays* constArrays,
-  SciFi::Tracking::HitSearchCuts& pars)
-{
-
-  float maxChi2 = 1.e9f;
-  bool parabola = false; // first linear than parabola
-  //== Fit a line
-  const float txs = track.trackParams[0]; // simplify overgeneral c++ calculation
-  const float tsxz = velo_state.x + (SciFi::Tracking::zReference - velo_state.z) * velo_state.tx;
-  const float tolYMag = SciFi::Tracking::tolYMag + SciFi::Tracking::tolYMagSlope * fabsf(txs - tsxz);
-  const float wMag = 1.f / (tolYMag * tolYMag);
-
-  bool doFit = true;
-  while (doFit) {
-    // Use position in magnet as constrain in fit
-    // although because wMag is quite small only little influence...
-    float zMag = zMagnet(velo_state, constArrays);
-    const float tys = track.trackParams[4] + (zMag - SciFi::Tracking::zReference) * track.trackParams[5];
-    const float tsyz = velo_state.y + (zMag - velo_state.z) * velo_state.ty;
-    const float dyMag = tys - tsyz;
-    zMag -= SciFi::Tracking::zReference;
-    float s0 = wMag;
-    float sz = wMag * zMag;
-    float sz2 = wMag * zMag * zMag;
-    float sd = wMag * dyMag;
-    float sdz = wMag * dyMag * zMag;
-
-    if (parabola) {
-
-      // position in magnet not used for parabola fit, hardly any influence on efficiency
-      fitParabola(stereoHits, n_stereoHits, scifi_hits, track.trackParams, false);
-    }
-    else { // straight line fit
-
-      for (int i_hit = 0; i_hit < n_stereoHits; ++i_hit) {
-        int hit = stereoHits[i_hit];
-        const float d = -trackToHitDistance(track.trackParams, scifi_hits, hit) /
-                        scifi_hits.dxdy(hit); // TODO multiplication much faster than division!
-        const float w = scifi_hits.w(hit);
-        const float z = scifi_hits.z0[hit] - SciFi::Tracking::zReference;
-        s0 += w;
-        sz += w * z;
-        sz2 += w * z * z;
-        sd += w * d;
-        sdz += w * d * z;
-      }
-      const float den = (s0 * sz2 - sz * sz);
-      if (!(fabsf(den) > 1e-5)) {
-        return false;
-      }
-      const float da = (sd * sz2 - sdz * sz) / den;
-      const float db = (sdz * s0 - sd * sz) / den;
-      track.trackParams[4] += da;
-      track.trackParams[5] += db;
-    } // fit end, now doing outlier removal
-
-    int worst = n_stereoHits;
-    maxChi2 = 0.;
-    for (int i_hit = 0; i_hit < n_stereoHits; ++i_hit) {
-      int hit = stereoHits[i_hit];
-      float d = trackToHitDistance(track.trackParams, scifi_hits, hit);
-      float chi2 = d * d * scifi_hits.w(hit);
-      if (chi2 > maxChi2) {
-        maxChi2 = chi2;
-        worst = i_hit;
-      }
-    }
-
-    if (maxChi2 < SciFi::Tracking::maxChi2StereoLinear && !parabola) {
-      parabola = true;
-      maxChi2 = 1.e9f;
-      continue;
-    }
-
-    if (maxChi2 > SciFi::Tracking::maxChi2Stereo) {
-      removeOutlier(scifi_hits, planeCounter, stereoHits, n_stereoHits, stereoHits[worst]);
-      if (planeCounter.nbDifferent < pars.minStereoHits) {
-        return false;
-      }
-      continue;
-    }
-    break;
-  }
-  return true;
-}
diff --git a/cuda/SciFi/PrForward/src/PrForward.cu b/cuda/SciFi/README.md
similarity index 58%
rename from cuda/SciFi/PrForward/src/PrForward.cu
rename to cuda/SciFi/README.md
index 6761623622e22424c992c9f29a6f072aeab99914..b8cb83696ecd304b2e2743afb47d70c330e53824 100644
--- a/cuda/SciFi/PrForward/src/PrForward.cu
+++ b/cuda/SciFi/README.md
@@ -94,10 +94,8 @@
 //
 // *****************************************************************
 
-#include "PrForward.cuh"
-
 //-----------------------------------------------------------------------------
-// Implementation file for class : PrForward
+// This description has been copied from the x86 PrForward code whose history is
 //
 // Based on code written by :
 // 2012-03-20 : Olivier Callot
@@ -107,90 +105,3 @@
 // 2018-08    : Vava Gligorov [extract code from Rec, make compile within GPU framework
 // 2018-09    : Dorothea vom Bruch [convert to CUDA, runs on GPU]
 //-----------------------------------------------------------------------------
-
-//=============================================================================
-
-// Kernel to call Forward tracking on GPU
-// Loop over veloUT input tracks using threadIdx.x
-__global__ void scifi_pr_forward(
-  uint32_t* dev_scifi_hits,
-  const uint32_t* dev_scifi_hit_count,
-  const int* dev_atomics_velo,
-  const uint* dev_velo_track_hit_number,
-  const char* dev_velo_states,
-  const int* dev_atomics_ut,
-  const char* dev_ut_track_hits,
-  const uint* dev_ut_track_hit_number,
-  const float* dev_ut_qop,
-  const uint* dev_ut_track_velo_indices,
-  SciFi::TrackHits* dev_scifi_tracks,
-  int* dev_atomics_scifi,
-  const SciFi::Tracking::TMVA* dev_tmva1,
-  const SciFi::Tracking::TMVA* dev_tmva2,
-  const SciFi::Tracking::Arrays* dev_constArrays,
-  const float* dev_magnet_polarity,
-  const char* dev_scifi_geometry,
-  const float* dev_inv_clus_res)
-{
-  const uint number_of_events = gridDim.x;
-  const uint event_number = blockIdx.x;
-
-  // Velo consolidated types
-  const Velo::Consolidated::Tracks velo_tracks {
-    (uint*) dev_atomics_velo, (uint*) dev_velo_track_hit_number, event_number, number_of_events};
-  const Velo::Consolidated::States velo_states {(char*) dev_velo_states, velo_tracks.total_number_of_tracks};
-  const uint velo_tracks_offset_event = velo_tracks.tracks_offset(event_number);
-
-  // UT consolidated tracks
-  UT::Consolidated::Tracks ut_tracks {(uint*) dev_atomics_ut,
-                                      (uint*) dev_ut_track_hit_number,
-                                      (float*) dev_ut_qop,
-                                      (uint*) dev_ut_track_velo_indices,
-                                      event_number,
-                                      number_of_events};
-  const int n_veloUT_tracks_event = ut_tracks.number_of_tracks(event_number);
-  const int ut_event_tracks_offset = ut_tracks.tracks_offset(event_number);
-
-  // SciFi un-consolidated track types
-  SciFi::TrackHits* scifi_tracks_event =
-    dev_scifi_tracks + ut_event_tracks_offset * SciFi::Constants::max_SciFi_tracks_per_UT_track;
-  int* atomics_scifi_event = dev_atomics_scifi + event_number;
-
-  // SciFi hits
-  const uint total_number_of_hits = dev_scifi_hit_count[number_of_events * SciFi::Constants::n_mat_groups_and_mats];
-  const SciFi::HitCount scifi_hit_count {(uint32_t*) dev_scifi_hit_count, event_number};
-  const SciFi::SciFiGeometry scifi_geometry {dev_scifi_geometry};
-  SciFi::Hits scifi_hits(dev_scifi_hits, total_number_of_hits, &scifi_geometry, dev_inv_clus_res);
-
-  // initialize atomic SciFi tracks counter
-  if (threadIdx.x == 0) {
-    *atomics_scifi_event = 0;
-  }
-  __syncthreads();
-
-  // Loop over the veloUT input tracks
-  for (int i = 0; i < (n_veloUT_tracks_event + blockDim.x - 1) / blockDim.x; ++i) {
-    const int i_veloUT_track = i * blockDim.x + threadIdx.x;
-    if (i_veloUT_track < n_veloUT_tracks_event) {
-      const float qop_ut = ut_tracks.qop[i_veloUT_track];
-
-      const int i_velo_track = ut_tracks.velo_track[i_veloUT_track];
-      const uint velo_states_index = velo_tracks_offset_event + i_velo_track;
-      const MiniState velo_state = velo_states.getMiniState(velo_states_index);
-
-      find_forward_tracks(
-        scifi_hits,
-        scifi_hit_count,
-        qop_ut,
-        i_veloUT_track,
-        scifi_tracks_event,
-        (uint*) atomics_scifi_event,
-        n_veloUT_tracks_event,
-        dev_tmva1,
-        dev_tmva2,
-        dev_constArrays,
-        dev_magnet_polarity[0],
-        velo_state);
-    }
-  }
-}
diff --git a/cuda/SciFi/classifiers/README.md b/cuda/SciFi/classifiers/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..74ef5519a4e393a1aba08d16048888073bdf220b
--- /dev/null
+++ b/cuda/SciFi/classifiers/README.md
@@ -0,0 +1 @@
+# Contains all the classifiers used to discriminate between real and fake tracks after the SciFi reconstruction
diff --git a/cuda/SciFi/PrForward/include/TMVA_Forward.cuh b/cuda/SciFi/classifiers/include/TMVA_Forward.cuh
similarity index 96%
rename from cuda/SciFi/PrForward/include/TMVA_Forward.cuh
rename to cuda/SciFi/classifiers/include/TMVA_Forward.cuh
index 925181aff1cb8dad75781587d5ba87b0ad8c0649..1bf3d40cdcd1033c882335cb876565193793b938 100644
--- a/cuda/SciFi/PrForward/include/TMVA_Forward.cuh
+++ b/cuda/SciFi/classifiers/include/TMVA_Forward.cuh
@@ -5,8 +5,6 @@
 #include <string>
 #include <iostream>
 
-#include "PrForwardConstants.cuh"
-
 namespace SciFi {
 
   namespace Tracking {
diff --git a/cuda/SciFi/PrForward/include/TMVA_Forward_1.cuh b/cuda/SciFi/classifiers/include/TMVA_Forward_1.cuh
similarity index 100%
rename from cuda/SciFi/PrForward/include/TMVA_Forward_1.cuh
rename to cuda/SciFi/classifiers/include/TMVA_Forward_1.cuh
diff --git a/cuda/SciFi/PrForward/include/TMVA_Forward_2.cuh b/cuda/SciFi/classifiers/include/TMVA_Forward_2.cuh
similarity index 100%
rename from cuda/SciFi/PrForward/include/TMVA_Forward_2.cuh
rename to cuda/SciFi/classifiers/include/TMVA_Forward_2.cuh
diff --git a/cuda/SciFi/PrForward/src/TMVA_Forward.cu b/cuda/SciFi/classifiers/src/TMVA_Forward.cu
similarity index 100%
rename from cuda/SciFi/PrForward/src/TMVA_Forward.cu
rename to cuda/SciFi/classifiers/src/TMVA_Forward.cu
diff --git a/cuda/SciFi/common/include/SciFiDefinitions.cuh b/cuda/SciFi/common/include/SciFiDefinitions.cuh
index 93f8da8f1ac6402c61abf49ac59c5716785a897b..5236e9d4f02e7536a8f05ec039e0fbe71a1d9972 100644
--- a/cuda/SciFi/common/include/SciFiDefinitions.cuh
+++ b/cuda/SciFi/common/include/SciFiDefinitions.cuh
@@ -8,7 +8,6 @@
 #include "device_launch_parameters.h"
 #include "Common.h"
 #include "Logger.h"
-#include "PrForwardConstants.cuh"
 #include "States.cuh"
 #include "SciFiRaw.cuh"
 
@@ -19,6 +18,97 @@ namespace SciFi {
   // need 3 arrays (size: number_of_events) for copy_and_prefix_sum_scifi_t
   static constexpr int num_atomics = 3;
 
+  namespace Tracking {
+
+    // The base PT threshold which is common to all algorithms
+    constexpr float minPt = 500 * Gaudi::Units::MeV; 
+
+    constexpr int max_scifi_hits = 20;  // for x and u/v layers
+    constexpr int nTrackParams = 9;
+
+    constexpr float tolYMag = 10. * Gaudi::Units::mm;
+    constexpr float tolYMagSlope = 0.015;
+
+    // parameterizations
+    constexpr float byParams = -0.667996;
+    constexpr float cyParams = -3.68424e-05;
+
+    // stereo hit matching
+    constexpr float tolYCollectX = 3.5 * Gaudi::Units::mm;        // 4.1* Gaudi::Units::mm ;
+    constexpr float tolYSlopeCollectX = 0.001 * Gaudi::Units::mm; // 0.0018 * Gaudi::Units::mm ;
+
+    // veloUT momentum estimate
+    constexpr bool useMomentumEstimate = true;
+    constexpr bool useWrongSignWindow = true;
+    constexpr float wrongSignPT = 2000. * Gaudi::Units::MeV;
+
+    // z Reference plane
+    constexpr float zReference = 8520. * Gaudi::Units::mm; // in T2
+    constexpr float zRefInv = 1.f / zReference; 
+
+    // TODO: CHECK THESE VALUES USING FRAMEWORK
+    constexpr float xLim_Max = 3300.;
+    constexpr float yLim_Max = 2500.;
+    constexpr float xLim_Min = -3300.;
+    constexpr float yLim_Min = -25.;
+
+    // TO BE READ FROM XML EVENTUALLY
+    //constexpr float magscalefactor = -1;
+    constexpr int zoneoffsetpar = 6;
+
+    struct Arrays {
+      // Returns whether the current layer is an X plane
+      const bool is_x_plane [12] {1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1};
+
+      // the Magnet Parametrization
+      // parameterized in offset [0], (slope difference due to kick)^2 [1],
+      // tx^2 [2], ty^2 [3]
+      const float zMagnetParams[4] = {5212.38, 406.609, -1102.35, -498.039};
+
+      // more Parametrizations
+      const float xParams[2] = {18.6195, -5.55793};
+
+      // momentum Parametrization
+      const float momentumParams[6] = {1.21014, 0.637339, -0.200292, 0.632298, 3.23793, -27.0259};
+
+      // covariance values
+      const float covarianceValues[5] = {4.0, 400.0, 4.e-6, 1.e-4, 0.1};
+
+      // definition of zones
+      // access upper with offset of 6
+      const int zoneoffsetpar = 6;
+      const int xZones[12] = {0, 6, 8, 14, 16, 22, 1, 7, 9, 15, 17, 23};
+      const int uvZones[12] = {2, 4, 10, 12, 18, 20, 3, 5, 11, 13, 19, 21};
+
+      // ASSORTED GEOMETRY VALUES, eventually read this from some xml
+      const float xZone_zPos[6] = {7826., 8036., 8508., 8718., 9193., 9403.};
+      const float uvZone_zPos[12] =
+        {7896., 7966., 8578., 8648., 9263., 9333., 7896., 7966., 8578., 8648., 9263., 9333.};
+      const float uvZone_dxdy[12] = {0.0874892,
+                                     -0.0874892,
+                                     0.0874892,
+                                     -0.0874892,
+                                     0.0874892,
+                                     -0.0874892,
+                                     0.0874892,
+                                     -0.0874892,
+                                     0.0874892,
+                                     -0.0874892,
+                                     0.0874892,
+                                     -0.0874892};
+      const float Zone_dzdy[24] = {0.0036010};
+
+      // this is used by looking_forward_sbt maybe this is not the right place to put it
+      const float uv_dx[6] = {1.6739478541449213,
+                              1.6738495069872612,
+                              1.935683825160498,
+                              1.9529279746403518,
+                              2.246931985749485,
+                              2.2797556995480273};
+    };
+ 
+  } // namespace Tracking
+
   namespace Constants {
     // Detector description
     // There are three stations with four layers each
diff --git a/cuda/SciFi/consolidate/include/ConsolidateSciFi.cuh b/cuda/SciFi/consolidate/include/ConsolidateSciFi.cuh
index 5b0e020bf048664ad84cab2e7ffef93fd6192614..a1c667931fe6a30b1e0a93e6f41bc64c819b5df7 100644
--- a/cuda/SciFi/consolidate/include/ConsolidateSciFi.cuh
+++ b/cuda/SciFi/consolidate/include/ConsolidateSciFi.cuh
@@ -8,7 +8,6 @@
 #include "ArgumentsSciFi.cuh"
 #include "ArgumentsUT.cuh"
 #include "LFFitTools.cuh"
-#include "PrForwardConstants.cuh"
 #include "LookingForwardConstants.cuh"
 
 __global__ void consolidate_scifi_tracks(
diff --git a/cuda/SciFi/looking_forward_sbt/collect_candidates/include/LFCollectCandidates.cuh b/cuda/SciFi/looking_forward_sbt/collect_candidates/include/LFCollectCandidates.cuh
index 714fbdefea21d056279f19ff149fe92b802ce63c..ea380fe634c13c3979b4c65ace5fd8de63d3d4cc 100644
--- a/cuda/SciFi/looking_forward_sbt/collect_candidates/include/LFCollectCandidates.cuh
+++ b/cuda/SciFi/looking_forward_sbt/collect_candidates/include/LFCollectCandidates.cuh
@@ -1,6 +1,5 @@
 #pragma once
 
-#include "PrForwardConstants.cuh"
 #include "VeloConsolidated.cuh"
 #include "UTConsolidated.cuh"
 #include "SciFiEventModel.cuh"
diff --git a/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindows.cuh b/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindows.cuh
index 7d9cfc74da2af2e22b9553ed7a905d17055bb31d..bea4cd10a88196e04e1b93bdb42aa4848ea61d48 100644
--- a/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindows.cuh
+++ b/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindows.cuh
@@ -1,6 +1,5 @@
 #pragma once
 
-#include "PrForwardConstants.cuh"
 #include "VeloConsolidated.cuh"
 #include "UTConsolidated.cuh"
 #include "SciFiEventModel.cuh"
@@ -11,6 +10,7 @@
 #include "ArgumentsSciFi.cuh"
 #include "LookingForwardConstants.cuh"
 #include "LookingForwardTools.cuh"
+#include "TrackUtils.cuh"
 
 __global__ void lf_search_initial_windows(
   uint32_t* dev_scifi_hits,
diff --git a/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindowsImpl.cuh b/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindowsImpl.cuh
index 3585e9534a8bd658a38d6afaa21a73f33835935b..2850cd0296cb73b43a1548f741f37ef4c2d14044 100644
--- a/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindowsImpl.cuh
+++ b/cuda/SciFi/looking_forward_sbt/search_initial_windows/include/LFSearchInitialWindowsImpl.cuh
@@ -6,14 +6,10 @@
 #include <algorithm>
 #include <fstream>
 #include "SciFiDefinitions.cuh"
-#include "PrForwardConstants.cuh"
 #include "UTDefinitions.cuh"
-#include "TrackUtils.cuh"
-#include "HitUtils.cuh"
-#include "LinearFitting.cuh"
-#include "ReferencePlaneProjection.cuh"
 #include "SciFiEventModel.cuh"
 #include "LookingForwardUtils.h"
+#include "TrackUtils.cuh"
 
 __device__ inline float evalCubicParameterization(const float value_at_ref, const float t, const float z);
 
diff --git a/cuda/SciFi/looking_forward_sbt/triplet_keep_best/include/LFTripletKeepBest.cuh b/cuda/SciFi/looking_forward_sbt/triplet_keep_best/include/LFTripletKeepBest.cuh
index db399acff6c38853b835f40476f36666cbc90727..7ed714930fa21956ea74d01d23abdb305c3aff75 100644
--- a/cuda/SciFi/looking_forward_sbt/triplet_keep_best/include/LFTripletKeepBest.cuh
+++ b/cuda/SciFi/looking_forward_sbt/triplet_keep_best/include/LFTripletKeepBest.cuh
@@ -1,6 +1,5 @@
 #pragma once
 
-#include "PrForwardConstants.cuh"
 #include "VeloConsolidated.cuh"
 #include "UTConsolidated.cuh"
 #include "SciFiEventModel.cuh"
diff --git a/cuda/SciFi/looking_forward_sbt/triplet_seeding/include/LFTripletSeeding.cuh b/cuda/SciFi/looking_forward_sbt/triplet_seeding/include/LFTripletSeeding.cuh
index 3d7473cb80e9fa65b70fa6c57f323bcfdb091725..f5ff86367b693e197a16edc747c0ddb07dc99b20 100644
--- a/cuda/SciFi/looking_forward_sbt/triplet_seeding/include/LFTripletSeeding.cuh
+++ b/cuda/SciFi/looking_forward_sbt/triplet_seeding/include/LFTripletSeeding.cuh
@@ -1,6 +1,5 @@
 #pragma once
 
-#include "PrForwardConstants.cuh"
 #include "VeloConsolidated.cuh"
 #include "UTConsolidated.cuh"
 #include "SciFiEventModel.cuh"
diff --git a/cuda/SciFi/looking_forward_sbt/triplet_seeding/src/LFTripletSeeding.cu b/cuda/SciFi/looking_forward_sbt/triplet_seeding/src/LFTripletSeeding.cu
index 362b3c4eda47d2b165d1c23c6a60b81bdbcefc56..29316bb86555118007d348cc6dc10dfbfdfc529b 100644
--- a/cuda/SciFi/looking_forward_sbt/triplet_seeding/src/LFTripletSeeding.cu
+++ b/cuda/SciFi/looking_forward_sbt/triplet_seeding/src/LFTripletSeeding.cu
@@ -1,6 +1,5 @@
 #include "LFTripletSeeding.cuh"
 #include "LFTripletSeedingImpl.cuh"
-#include "TrackUtils.cuh"
 #include "LookingForwardTools.cuh"
 
 __global__ void lf_triplet_seeding(
diff --git a/cuda/SciFi/utils/README.md b/cuda/SciFi/utils/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..18a329b1a23c68e657c4688d575ecbabdb029917
--- /dev/null
+++ b/cuda/SciFi/utils/README.md
@@ -0,0 +1 @@
+# Contains utility functions useful in the SciFi reconstruction
diff --git a/cuda/SciFi/utils/include/HitUtils.cuh b/cuda/SciFi/utils/include/HitUtils.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..f6c4f004917cc2170acd526b3596ac213ed2a3b8
--- /dev/null
+++ b/cuda/SciFi/utils/include/HitUtils.cuh
@@ -0,0 +1,28 @@
+#pragma once
+
+#include "SciFiDefinitions.cuh"
+#include "SciFiEventModel.cuh"
+
+// match stereo hits to x hits
+__host__ __device__ bool matchStereoHit(
+  const int itUV1,
+  const int uv_zone_offset_end,
+  const SciFi::Hits& scifi_hits,
+  const float xMinUV,
+  const float xMaxUV);
+
+// check that val is within [min, max]
+__host__ __device__ inline bool isInside(float val, const float min, const float max)
+{
+  return (val > min) && (val < max);
+}
+
+// get lowest index where range[index] > value, within [start,end] of range
+__host__ __device__ inline int getLowerBound(float range[], float value, int start, int end)
+{
+  int i = start;
+  for (; i < end; i++) {
+    if (range[i] > value) break;
+  }
+  return i;
+}
diff --git a/cuda/SciFi/utils/include/TrackUtils.cuh b/cuda/SciFi/utils/include/TrackUtils.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..d50dbeaf2627e042351f8627c7007655c0e57ef4
--- /dev/null
+++ b/cuda/SciFi/utils/include/TrackUtils.cuh
@@ -0,0 +1,35 @@
+#pragma once
+
+#include <cmath>
+#include "SciFiEventModel.cuh"
+#
+/**
+   Helper functions related to track properties
+ */
+
+// extrapolate x position from given state to z
+__host__ __device__ float xFromVelo(const float z, const MiniState& velo_state);
+
+// extrapolate y position from given state to z
+__host__ __device__ float yFromVelo(const float z, const MiniState& velo_state);
+
+__host__ __device__ float evalCubicParameterization(const float params[4], float z);
+
+__host__ __device__ float zMagnet(const MiniState& velo_state, const SciFi::Tracking::Arrays* constArrays);
+
+__host__ __device__ float calcDxRef(float pt, const MiniState& velo_state);
+
+__host__ __device__ float calcqOverP(float bx, const SciFi::Tracking::Arrays* constArrays, const MiniState& velo_state, const float magnet_polarity);
+
+__host__ __device__ float evalParameterizationX(const float* params, float z);
+
+__host__ __device__ float evalParameterizationY(const float* params, float z);
+
+__host__ __device__ void getTrackParameters(
+  float xAtRef,
+  const MiniState& velo_state,
+  const SciFi::Tracking::Arrays* constArrays,
+  float trackParams[SciFi::Tracking::nTrackParams]);
+
+__host__ __device__ float
+trackToHitDistance(const float trackParameters[SciFi::Tracking::nTrackParams], const SciFi::Hits& scifi_hits, int hit);
diff --git a/cuda/SciFi/utils/src/HitUtils.cu b/cuda/SciFi/utils/src/HitUtils.cu
new file mode 100644
index 0000000000000000000000000000000000000000..722efa6bad272288969bb575de384b60ef3d4bd5
--- /dev/null
+++ b/cuda/SciFi/utils/src/HitUtils.cu
@@ -0,0 +1,17 @@
+#include "HitUtils.cuh"
+
+// match stereo hits to x hits
+__host__ __device__ bool matchStereoHit(
+  const int itUV1,
+  const int uv_zone_offset_end,
+  const SciFi::Hits& scifi_hits,
+  const float xMinUV,
+  const float xMaxUV)
+{
+  for (int stereoHit = itUV1; stereoHit != uv_zone_offset_end; ++stereoHit) {
+    if (scifi_hits.x0[stereoHit] > xMinUV) {
+      return (scifi_hits.x0[stereoHit] < xMaxUV);
+    }
+  }
+  return false;
+}
diff --git a/cuda/SciFi/utils/src/TrackUtils.cu b/cuda/SciFi/utils/src/TrackUtils.cu
new file mode 100644
index 0000000000000000000000000000000000000000..69c0066c00bb5d048e0dc3225f091eac0f8282ef
--- /dev/null
+++ b/cuda/SciFi/utils/src/TrackUtils.cu
@@ -0,0 +1,108 @@
+#include "TrackUtils.cuh"
+
+// extrapolate x position from given state to z
+__host__ __device__ float xFromVelo(const float z, const MiniState& velo_state)
+{
+  return velo_state.x + (z - velo_state.z) * velo_state.tx;
+}
+
+// extrapolate y position from given state to z
+__host__ __device__ float yFromVelo(const float z, const MiniState& velo_state)
+{
+  return velo_state.y + (z - velo_state.z) * velo_state.ty;
+}
+
+__host__ __device__ float evalCubicParameterization(const float params[4], float z)
+{
+  float dz = z - SciFi::Tracking::zReference;
+  return params[0] + (params[1] + (params[2] + params[3] * dz) * dz) * dz;
+}
+
+// Find z zMag position within the magnet at which the bending ("kick") occurs
+// this is parameterized based on MC
+// the second parameter([1]) is multiplied by the difference in slope before and
+// after the kick, this slope is calculated from zMag and the x position of the track
+// at the reference plane -> it is calculated iteratively later
+__host__ __device__ float zMagnet(const MiniState& velo_state, const SciFi::Tracking::Arrays* constArrays)
+{
+
+  return (
+    constArrays->zMagnetParams[0] + constArrays->zMagnetParams[2] * velo_state.tx * velo_state.tx +
+    constArrays->zMagnetParams[3] * velo_state.ty * velo_state.ty);
+}
+
+// calculate difference between straight line extrapolation and
+// where a track with wrongSignPT (2 GeV) would be on the reference plane (?)
+__host__ __device__ float calcDxRef(float pt, const MiniState& velo_state)
+{
+  const float tx2 = velo_state.tx * velo_state.tx;
+  const float ty2 = velo_state.ty * velo_state.ty;
+  float m_slope2 = tx2 + ty2;
+  return 3973000.f * sqrtf(m_slope2) / pt - 2200.f * ty2 - 1000.f * tx2; // tune this window
+}
+
+__host__ __device__ float calcqOverP(float bx, const SciFi::Tracking::Arrays* constArrays, const MiniState& velo_state, const float magnet_polarity)
+{
+
+  float qop = 1.0f / Gaudi::Units::GeV;
+  const float bx2 = bx * bx;
+  const float ty2 = velo_state.ty * velo_state.ty;
+  const float coef =
+    (constArrays->momentumParams[0] + constArrays->momentumParams[1] * bx2 +
+     constArrays->momentumParams[2] * bx2 * bx2 + constArrays->momentumParams[3] * bx * velo_state.tx +
+     constArrays->momentumParams[4] * ty2 + constArrays->momentumParams[5] * ty2 * ty2);
+  const float tx2 = velo_state.tx * velo_state.tx;
+  float m_slope2 = tx2 + ty2;
+  float proj = sqrtf((1.f + m_slope2) / (1.f + tx2));
+  qop = (velo_state.tx - bx) / (coef * Gaudi::Units::GeV * proj * magnet_polarity);
+  return qop;
+}
+
+__host__ __device__ float evalParameterizationX(const float* params, float z)
+{
+  const float dz = z - SciFi::Tracking::zReference;
+  return params[0] + (params[1] + (params[2] + params[3] * dz) * dz) * dz;
+}
+
+__host__ __device__ float evalParameterizationY(const float* params, float z)
+{
+  const float dz = z - SciFi::Tracking::zReference;
+  return params[0] + (params[1] + params[2] * dz) * dz;
+}
+
+__host__ __device__ void getTrackParameters(
+  float xAtRef,
+  const MiniState& velo_state,
+  const SciFi::Tracking::Arrays* constArrays,
+  float trackParams[SciFi::Tracking::nTrackParams])
+{
+  float dSlope = (xFromVelo(SciFi::Tracking::zReference, velo_state) - xAtRef) /
+                 (SciFi::Tracking::zReference - constArrays->zMagnetParams[0]);
+  const float zMagSlope = constArrays->zMagnetParams[2] * velo_state.tx * velo_state.tx +
+                          constArrays->zMagnetParams[3] * velo_state.ty * velo_state.ty;
+  const float zMag = constArrays->zMagnetParams[0] + constArrays->zMagnetParams[1] * dSlope * dSlope + zMagSlope;
+  const float xMag = xFromVelo(zMag, velo_state);
+  const float slopeT = (xAtRef - xMag) / (SciFi::Tracking::zReference - zMag);
+  dSlope = slopeT - velo_state.tx;
+  const float dyCoef = dSlope * dSlope * velo_state.ty;
+
+  trackParams[0] = xAtRef;
+  trackParams[1] = slopeT;
+  trackParams[2] = 1.e-6f * constArrays->xParams[0] * dSlope;
+  trackParams[3] = 1.e-9f * constArrays->xParams[1] * dSlope;
+  trackParams[4] = yFromVelo(SciFi::Tracking::zReference, velo_state);
+  trackParams[5] = velo_state.ty + dyCoef * SciFi::Tracking::byParams;
+  trackParams[6] = dyCoef * SciFi::Tracking::cyParams;
+  trackParams[7] = 0.0f;
+  trackParams[8] = 0.0f; // last elements are chi2 and ndof, as float
+}
+
+__host__ __device__ float
+trackToHitDistance(const float trackParameters[SciFi::Tracking::nTrackParams], const SciFi::Hits& scifi_hits, int hit)
+{
+  const float z_Hit = scifi_hits.z0[hit] + scifi_hits.dzdy(hit) * evalParameterizationY(trackParameters + 4, scifi_hits.z0[hit]);
+  const float x_track = evalParameterizationX(trackParameters, z_Hit);
+  const float y_track = evalParameterizationY(trackParameters + 4, z_Hit);
+  return scifi_hits.x0[hit] + y_track * scifi_hits.dxdy(hit) - x_track;
+}
+
diff --git a/cuda/UT/compassUT/src/UTFastFitter.cu b/cuda/UT/compassUT/src/UTFastFitter.cu
index 71a1b37e9632e4789d2f5d7731f0c1620a906e39..d5ae353bf5309c10c78711eb781d92d68b0d30e3 100644
--- a/cuda/UT/compassUT/src/UTFastFitter.cu
+++ b/cuda/UT/compassUT/src/UTFastFitter.cu
@@ -51,34 +51,34 @@ __host__ __device__ float fastfitter(
   const float zDiff = 0.001f * (zKink - UT::Constants::zMidUT);
 
   // -- This is to avoid division by zero...
-  const float pHelper = std::max(float(std::abs(best_params.qp * qpxz2p)), float(1e-9));
-  const float invP = pHelper * sqrtf(1.0f + ty * ty);
-
+  const float pHelper   = std::max( float(std::abs(best_params.qp * qpxz2p)), float(1e-9));
+  const float invP      = pHelper*sqrtf(1.0f+ty*ty);
+  
   // these resolution are semi-empirical, could be tuned and might not be correct for low momentum.
   // this is the resolution due to multiple scattering between Velo and UT
-  const float error1 = 0.14f + 10000.0f * invP;
+  const float error1    = 0.14f + 10000.0f*invP; 
   // this is the resolution due to the finite Velo resolution
-  const float error2 = 0.12f + 3000.0f * invP;
-  const float error = error1 * error1 + error2 * error2;
-  const float weight = 1.0f / error;
-
-  float mat[6] = {weight, weight * zDiff, weight * zDiff * zDiff, 0.0f, 0.0f, 0.0f};
-  float rhs[3] = {weight * xMidField, weight * xMidField * zDiff, 0.0f};
+  const float error2    = 0.12f + 3000.0f*invP;  
+  const float error     = error1*error1 + error2*error2;
+  const float weight    = 1.0f/error;
+  
+  float mat[6]          = {weight, weight * zDiff, weight * zDiff * zDiff, 0.0f, 0.0f, 0.0f};
+  float rhs[3]          = {weight * xMidField, weight * xMidField * zDiff, 0.0f };
 
   const float yyProto = velo_state.y - velo_state.ty * velo_state.z;
-
+  
   for (int i = 0; i < UT::Constants::n_layers; ++i) {
     if (best_hits[i] != -1) {
       const auto hit = best_hits[i];
-
+      
       const int plane_code = i;
       const float dxDy = ut_dxDy[plane_code];
       const float yy = yyProto + (velo_state.ty * ut_hits.zAtYEq0[hit]);
       const float ui = ut_hits.xAt(hit, yy, dxDy);
       const float dz = 0.001f * (ut_hits.zAtYEq0[hit] - UT::Constants::zMidUT);
-      const float w = ut_hits.weight[hit];
-      const float t = ut_hits.sinT(hit, dxDy);
-
+      const float w = ut_hits.weight[hit];  
+      const float t = ut_hits.sinT(hit, dxDy); 
+      
       mat[0] += w;
       mat[1] += w * dz;
       mat[2] += w * dz * dz;
@@ -91,7 +91,7 @@ __host__ __device__ float fastfitter(
       rhs[2] += w * ui * t;
     }
   }
-
+  
   const float a11 = mat[2] * mat[5] - mat[4] * mat[4];
   const float a12 = mat[4] * mat[3] - mat[1] * mat[5];
   const float a13 = mat[1] * mat[4] - mat[2] * mat[3];
diff --git a/cuda/associate/CMakeLists.txt b/cuda/associate/CMakeLists.txt
index 76f5d2e27f580c913ddb540d21bfa68a845b0806..b1575b153dd7e416bcfc0bfd5f13c4f9ed6e0b29 100644
--- a/cuda/associate/CMakeLists.txt
+++ b/cuda/associate/CMakeLists.txt
@@ -3,7 +3,6 @@ include_directories(include)
 include_directories(${CMAKE_SOURCE_DIR}/main/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/associate/include)
diff --git a/cuda/global_event_cut/CMakeLists.txt b/cuda/global_event_cut/CMakeLists.txt
index 83074cb4f61296ed3ceab18007ea9122e6a73e6a..3bc76e4076cba8c6d5439af289d3e473fa79897d 100644
--- a/cuda/global_event_cut/CMakeLists.txt
+++ b/cuda/global_event_cut/CMakeLists.txt
@@ -7,7 +7,6 @@ include_directories(include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 
 add_library(GlobalEventCut STATIC
   ${global_event_cut}
diff --git a/cuda/kalman/CMakeLists.txt b/cuda/kalman/CMakeLists.txt
index 0bc451b61a8094a6d9cd0294edc3c237a9fd85a3..74984fa38defb90c55958ea3a039e72f38bfb9e8 100644
--- a/cuda/kalman/CMakeLists.txt
+++ b/cuda/kalman/CMakeLists.txt
@@ -3,7 +3,6 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/PrVeloUT/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/UT/include)
diff --git a/cuda/muon/CMakeLists.txt b/cuda/muon/CMakeLists.txt
index cbc307e087362c3ad6822f83b0220150ffaf5bdd..883b67a6a10dd26a1f9460bf8609f8ca228fe121 100644
--- a/cuda/muon/CMakeLists.txt
+++ b/cuda/muon/CMakeLists.txt
@@ -13,7 +13,6 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/SciFi/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/main/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/utils/sorting/include)
diff --git a/cuda/selections/CMakeLists.txt b/cuda/selections/CMakeLists.txt
index b8d89779947592652af40c2f1dffb1b3073b2928..2f217ccbb0e525b106bd063ed8c341640f53f9ce 100644
--- a/cuda/selections/CMakeLists.txt
+++ b/cuda/selections/CMakeLists.txt
@@ -6,7 +6,6 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/beamlinePV/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/vertex_fit/common/include)
diff --git a/cuda/utils/CMakeLists.txt b/cuda/utils/CMakeLists.txt
index 9abefb046fd906692711b95b3cf9a9ecc3e8035f..172e35a1b1966410bfbae2627f05dc89aad40775 100644
--- a/cuda/utils/CMakeLists.txt
+++ b/cuda/utils/CMakeLists.txt
@@ -14,7 +14,6 @@ include_directories(${CMAKE_SOURCE_DIR}/stream/gear/include)
 include_directories(${CMAKE_SOURCE_DIR}/stream/setup/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/muon/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/muon/decoding/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/kalman/ParKalman/include)
diff --git a/cuda/vertex_fit/CMakeLists.txt b/cuda/vertex_fit/CMakeLists.txt
index eb692eea66f60a75dd8d585cae39836f9fe839fd..46dd03486d46f1196d71f4b5f7b75b6f3866ada5 100644
--- a/cuda/vertex_fit/CMakeLists.txt
+++ b/cuda/vertex_fit/CMakeLists.txt
@@ -7,7 +7,6 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/beamlinePV/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/main/include)
diff --git a/integration/non_event_data/CMakeLists.txt b/integration/non_event_data/CMakeLists.txt
index f71b3d377139926044a43adf0b5d10216845e660..fc4d404fce2c390e31e55e212576cc7de31de26c 100644
--- a/integration/non_event_data/CMakeLists.txt
+++ b/integration/non_event_data/CMakeLists.txt
@@ -11,7 +11,7 @@ target_include_directories(NonEventData PUBLIC
   ${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include
   ${CMAKE_SOURCE_DIR}/cuda/muon/common/include
   ${CMAKE_SOURCE_DIR}/cuda/muon/decoding/include
-  ${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include
+  ${CMAKE_SOURCE_DIR}/cuda/SciFi/classifiers/include
   ${CMAKE_SOURCE_DIR}/cuda/UT/common/include
   ${CMAKE_SOURCE_DIR}/stream/sequence/include
   ${CMAKE_SOURCE_DIR}/checker/tracking/include
diff --git a/mdf/CMakeLists.txt b/mdf/CMakeLists.txt
index 3530c43f6908314f1e7f9100f0f301a7c1fbf2c0..fa2a3beff3950df686197a43cd5c8ba003d9cb1a 100644
--- a/mdf/CMakeLists.txt
+++ b/mdf/CMakeLists.txt
@@ -58,7 +58,6 @@ target_include_directories(${name} PUBLIC
   ${CMAKE_SOURCE_DIR}/cuda/velo/common/include
   ${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include
   ${CMAKE_SOURCE_DIR}/cuda/muon/common/include
-  ${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include
   ${CMAKE_SOURCE_DIR}/cuda/UT/common/include
   ${CMAKE_SOURCE_DIR}/cuda/UT/PrVeloUT/include
   ${CMAKE_SOURCE_DIR}/checker/tracking/include)
diff --git a/stream/CMakeLists.txt b/stream/CMakeLists.txt
index 6477e2a843f1893c5dbb4df68efbe243fc753225..6f5a2a053bdc6e3af6f56141244b59549ebadb79 100644
--- a/stream/CMakeLists.txt
+++ b/stream/CMakeLists.txt
@@ -29,7 +29,8 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/search_by_triplet/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/simplified_kalman_filter/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/preprocessing/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
+include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/classifiers/include)
+include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/utils/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/looking_forward/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/looking_forward/calculate_first_layer_window/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/looking_forward/calculate_second_layer_window/include)
@@ -66,7 +67,6 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/PV/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/x86/muon/decoding/include)
 include_directories(${CMAKE_SOURCE_DIR}/x86/velo/clustering/include)
 include_directories(${CMAKE_SOURCE_DIR}/x86/UT/PrVeloUT/include)
-include_directories(${CMAKE_SOURCE_DIR}/x86/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/x86/SciFi/LookingForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/x86/SciFi/MomentumForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/x86/PV/beamlinePV/include)
@@ -145,7 +145,6 @@ target_link_libraries(Stream PRIVATE
   Hlt1
   SciFi
   x86MuonDecoding
-  x86Forward
   x86LookingForward
   x86MomentumForward
   x86VeloUT
diff --git a/stream/checkers/include/SciFiSequenceCheckers_impl.cuh b/stream/checkers/include/SciFiSequenceCheckers_impl.cuh
index a7c2dc910935857b798deb3c166250550197185a..0e45cf34c667a489144b388ca229e2d86a4b595b 100644
--- a/stream/checkers/include/SciFiSequenceCheckers_impl.cuh
+++ b/stream/checkers/include/SciFiSequenceCheckers_impl.cuh
@@ -1,5 +1,4 @@
-#include "PrForward.cuh"
-#include "RunForwardCPU.h"
+#include "PrepareTracks.h"
 
 /**
  * @brief Specialization when invoking scifi_pr_forward_t as last step.
diff --git a/stream/checkers/include/UTSequenceCheckers_impl.cuh b/stream/checkers/include/UTSequenceCheckers_impl.cuh
index d98940676c98004dd1ef7c2a8e52070eb93ae46d..c552a4244072daee8c64003f304dff40d3bb72d0 100644
--- a/stream/checkers/include/UTSequenceCheckers_impl.cuh
+++ b/stream/checkers/include/UTSequenceCheckers_impl.cuh
@@ -1,6 +1,7 @@
 #include "MomentumForwardStudies.h"
 #include "LookingForwardStudies.h"
 #include "TrackChecker.h"
+#include "PrepareTracks.h"
 
 /**
  * @brief Specialization for any Velo reconstruction algorithm invoking
diff --git a/stream/checkers/include/VeloSequenceCheckers_impl.cuh b/stream/checkers/include/VeloSequenceCheckers_impl.cuh
index 0e6473b07ed044798d274126951cd204a262469f..e008f10635ac7d5885acceb2145a5e0a1449e8d0 100644
--- a/stream/checkers/include/VeloSequenceCheckers_impl.cuh
+++ b/stream/checkers/include/VeloSequenceCheckers_impl.cuh
@@ -1,5 +1,6 @@
 #include "ConsolidateVelo.cuh"
 #include "TrackChecker.h"
+#include "PrepareTracks.h"
 
 /**
  * @brief Specialization for any Velo reconstruction algorithm invoking
diff --git a/stream/sequence/include/Constants.cuh b/stream/sequence/include/Constants.cuh
index 07a4db110e60e0b0ff375ed6bdb18c7dd7c2076b..0ca59366b97b819197d3381ddde3b50c17f783f0 100644
--- a/stream/sequence/include/Constants.cuh
+++ b/stream/sequence/include/Constants.cuh
@@ -8,7 +8,6 @@
 #include "VeloDefinitions.cuh"
 #include "ClusteringDefinitions.cuh"
 #include "ClusteringCommon.h"
-#include "PrForwardConstants.cuh"
 #include "TMVA_Forward_1.cuh"
 #include "TMVA_Forward_2.cuh"
 #include "UTDefinitions.cuh"
diff --git a/stream/setup/include/ConfiguredSequence.cuh b/stream/setup/include/ConfiguredSequence.cuh
index 343e9d16e757029b2bb950a9e947af51e5ca3d7a..e90f0c6f500d53050bb2f2d2f1fc18e3b8cfe42e 100644
--- a/stream/setup/include/ConfiguredSequence.cuh
+++ b/stream/setup/include/ConfiguredSequence.cuh
@@ -28,7 +28,6 @@
 #include "ConsolidateSciFi.cuh"
 #include "SearchWindows.cuh"
 #include "CompassUT.cuh"
-#include "PrForward.cuh"
 #include "VeloKalmanFilter.cuh"
 #include "GetSeeds.cuh"
 #include "FitSeeds.cuh"
@@ -37,7 +36,6 @@
 #include "pv_beamline_peak.cuh"
 #include "pv_beamline_multi_fitter.cuh"
 #include "pv_beamline_cleanup.cuh"
-#include "RunForwardCPU.h"
 #include "VeloPVIP.cuh"
 #include "RunMomentumForwardCPU.h"
 #include "RunBeamlinePVOnCPU.h"
diff --git a/stream/visitors/SciFi/src/CpuSciFiMomentumForwardVisitor.cu b/stream/visitors/SciFi/src/CpuSciFiMomentumForwardVisitor.cu
index ce03306b26717763465bad74aa676d6767b1699c..67a1f571d2c9c5d60e377b3ab899119d490b2c7d 100644
--- a/stream/visitors/SciFi/src/CpuSciFiMomentumForwardVisitor.cu
+++ b/stream/visitors/SciFi/src/CpuSciFiMomentumForwardVisitor.cu
@@ -1,5 +1,4 @@
 #include "SequenceVisitor.cuh"
-#include "RunForwardCPU.h"
 #include "RunMomentumForwardCPU.h"
 #include "Tools.h"
 
diff --git a/stream/visitors/SciFi/src/CpuSciFiPrForwardVisitor.cu b/stream/visitors/SciFi/src/CpuSciFiPrForwardVisitor.cu
deleted file mode 100644
index dc1e13c13e688a710113b68f580c23b922406033..0000000000000000000000000000000000000000
--- a/stream/visitors/SciFi/src/CpuSciFiPrForwardVisitor.cu
+++ /dev/null
@@ -1,91 +0,0 @@
-#include "SequenceVisitor.cuh"
-#include "RunForwardCPU.h"
-#include "Tools.h"
-
-template<>
-void SequenceVisitor::set_arguments_size<cpu_scifi_pr_forward_t>(
-  cpu_scifi_pr_forward_t::arguments_t arguments,
-  const RuntimeOptions& runtime_options,
-  const Constants& constants,
-  const HostBuffers& host_buffers)
-{
-  arguments.set_size<dev_scifi_tracks>(
-    host_buffers.host_number_of_reconstructed_ut_tracks[0] * SciFi::Constants::max_SciFi_tracks_per_UT_track);
-  arguments.set_size<dev_atomics_scifi>(host_buffers.host_number_of_selected_events[0] * SciFi::num_atomics);
-}
-
-template<>
-void SequenceVisitor::visit<cpu_scifi_pr_forward_t>(
-  cpu_scifi_pr_forward_t& state,
-  const cpu_scifi_pr_forward_t::arguments_t& arguments,
-  const RuntimeOptions& runtime_options,
-  const Constants& constants,
-  HostBuffers& host_buffers,
-  cudaStream_t& cuda_stream,
-  cudaEvent_t& cuda_generic_event)
-{
-  // Synchronize previous CUDA transmissions
-  cudaEventRecord(cuda_generic_event, cuda_stream);
-  cudaEventSynchronize(cuda_generic_event);
-
-  // Run Forward on x86 architecture
-  // ATTENTION: when using SciFi raw bank version 5,
-  // need: 2*host_buffers.host_number_of_selected_events[0]*...
-  host_buffers.host_velo_states.resize(arguments.size<dev_velo_states>());
-  host_buffers.host_scifi_hits.resize(arguments.size<dev_scifi_hits>());
-  host_buffers.host_scifi_hit_count.resize(arguments.size<dev_scifi_hit_count>());
-
-  cudaCheck(cudaMemcpy(
-    host_buffers.host_scifi_hits.data(),
-    arguments.offset<dev_scifi_hits>(),
-    arguments.size<dev_scifi_hits>(),
-    cudaMemcpyDeviceToHost));
-
-  cudaCheck(cudaMemcpy(
-    host_buffers.host_scifi_hit_count.data(),
-    arguments.offset<dev_scifi_hit_count>(),
-    arguments.size<dev_scifi_hit_count>(),
-    cudaMemcpyDeviceToHost));
-
-  cudaCheck(cudaMemcpy(
-    host_buffers.host_velo_states.data(),
-    arguments.offset<dev_velo_states>(),
-    arguments.size<dev_velo_states>(),
-    cudaMemcpyDeviceToHost));
-
-  // TODO: Maybe use this rv somewhere?
-  int rv = state.invoke(
-    host_buffers.scifi_tracks_events.data(),
-    host_buffers.host_atomics_scifi,
-    host_buffers.host_scifi_hits.data(),
-    host_buffers.host_scifi_hit_count.data(),
-    constants.host_scifi_geometry.data(),
-    constants.host_inv_clus_res,
-    host_buffers.host_atomics_velo,
-    host_buffers.host_velo_track_hit_number,
-    host_buffers.host_velo_states.data(),
-    host_buffers.host_atomics_ut,
-    host_buffers.host_ut_track_hit_number,
-    host_buffers.host_ut_qop,
-    host_buffers.host_ut_track_velo_indices,
-    host_buffers.host_number_of_selected_events[0]);
-
-  if (logger::ll.verbosityLevel >= logger::debug) {
-    for (int i = 0; i < host_buffers.host_number_of_selected_events[0]; ++i) {
-      debug_cout << "Visitor: found " << host_buffers.host_atomics_scifi[i] << " tracks in event " << i << std::endl;
-    }
-  }
-
-  // copy SciFi track to device for consolidation
-  cudaCheck(cudaMemcpy(
-    arguments.offset<dev_atomics_scifi>(),
-    host_buffers.host_atomics_scifi,
-    arguments.size<dev_atomics_scifi>(),
-    cudaMemcpyHostToDevice));
-
-  cudaCheck(cudaMemcpy(
-    arguments.offset<dev_scifi_tracks>(),
-    host_buffers.scifi_tracks_events.data(),
-    arguments.size<dev_scifi_tracks>(),
-    cudaMemcpyHostToDevice));
-}
diff --git a/stream/visitors/SciFi/src/SciFiPrForwardVisitor.cu b/stream/visitors/SciFi/src/SciFiPrForwardVisitor.cu
deleted file mode 100644
index 8c5e9935eb5650dc60a8095a670e0b27aabdfcea..0000000000000000000000000000000000000000
--- a/stream/visitors/SciFi/src/SciFiPrForwardVisitor.cu
+++ /dev/null
@@ -1,48 +0,0 @@
-#include "PrForward.cuh"
-#include "SequenceVisitor.cuh"
-
-template<>
-void SequenceVisitor::set_arguments_size<scifi_pr_forward_t>(
-  scifi_pr_forward_t::arguments_t arguments,
-  const RuntimeOptions& runtime_options,
-  const Constants& constants,
-  const HostBuffers& host_buffers)
-{
-  arguments.set_size<dev_scifi_tracks>(
-    host_buffers.host_number_of_reconstructed_ut_tracks[0] * SciFi::Constants::max_SciFi_tracks_per_UT_track);
-  arguments.set_size<dev_atomics_scifi>(host_buffers.host_number_of_selected_events[0] * SciFi::num_atomics);
-}
-
-template<>
-void SequenceVisitor::visit<scifi_pr_forward_t>(
-  scifi_pr_forward_t& state,
-  const scifi_pr_forward_t::arguments_t& arguments,
-  const RuntimeOptions& runtime_options,
-  const Constants& constants,
-  HostBuffers& host_buffers,
-  cudaStream_t& cuda_stream,
-  cudaEvent_t& cuda_generic_event)
-{
-  state.set_opts(dim3(host_buffers.host_number_of_selected_events[0]), dim3(32), cuda_stream);
-  state.set_arguments(
-    arguments.offset<dev_scifi_hits>(),
-    arguments.offset<dev_scifi_hit_count>(),
-    arguments.offset<dev_atomics_velo>(),
-    arguments.offset<dev_velo_track_hit_number>(),
-    arguments.offset<dev_velo_states>(),
-    arguments.offset<dev_atomics_ut>(),
-    arguments.offset<dev_ut_track_hits>(),
-    arguments.offset<dev_ut_track_hit_number>(),
-    arguments.offset<dev_ut_qop>(),
-    arguments.offset<dev_ut_track_velo_indices>(),
-    arguments.offset<dev_scifi_tracks>(),
-    arguments.offset<dev_atomics_scifi>(),
-    constants.dev_scifi_tmva1,
-    constants.dev_scifi_tmva2,
-    constants.dev_scifi_constArrays,
-    constants.dev_magnet_polarity,
-    constants.dev_scifi_geometry,
-    constants.dev_inv_clus_res);
-
-  state.invoke();
-}
diff --git a/x86/SciFi/CMakeLists.txt b/x86/SciFi/CMakeLists.txt
index d0bfdc8dab90daa9c4c0fcd9e27468a6fbf8bc5d..22cebb5eb65eaa9e74a0e9a53334ebc3e39afde6 100644
--- a/x86/SciFi/CMakeLists.txt
+++ b/x86/SciFi/CMakeLists.txt
@@ -1,3 +1,2 @@
-add_subdirectory(PrForward)
 add_subdirectory(LookingForward)
 add_subdirectory(MomentumForward)
diff --git a/x86/SciFi/LookingForward/CMakeLists.txt b/x86/SciFi/LookingForward/CMakeLists.txt
index 7537e4ec1591716dc3cd5a19a232abcaa0f6816b..340f5b4f9b0f323998e8d1a83d9b2c364111ea60 100644
--- a/x86/SciFi/LookingForward/CMakeLists.txt
+++ b/x86/SciFi/LookingForward/CMakeLists.txt
@@ -12,7 +12,8 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/UT/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/SciFi/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
+include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/classifiers/include)
+include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/utils/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/utils/binary_search/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/looking_forward/common/include)
 include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
diff --git a/x86/SciFi/LookingForward/include/LookingForwardConstants.h b/x86/SciFi/LookingForward/include/LookingForwardConstants.h
index 230c6dfe258139ba5d1431de36f1ef3a319af784..a28950c5dce8479521e16b92aa94c4fd47c32a16 100644
--- a/x86/SciFi/LookingForward/include/LookingForwardConstants.h
+++ b/x86/SciFi/LookingForward/include/LookingForwardConstants.h
@@ -1,7 +1,7 @@
 #pragma once
 
 #include "SystemOfUnits.h"
-#include "PrForwardConstants.cuh"
+#include "SciFiDefinitions.cuh"
 
 namespace SciFi {
   namespace LookingForward {
diff --git a/x86/SciFi/LookingForward/include/LookingForwardSbt.h b/x86/SciFi/LookingForward/include/LookingForwardSbt.h
index bfaca59032096f72686d8adb793b05551bbadb02..c358457f18dbd362095effcb269ef0bd00d7dd59 100644
--- a/x86/SciFi/LookingForward/include/LookingForwardSbt.h
+++ b/x86/SciFi/LookingForward/include/LookingForwardSbt.h
@@ -1,7 +1,7 @@
 #pragma once
 
 #include "LookingForwardUtils.h"
-#include "FindXHits.cuh"
+#include "HitUtils.cuh"
 
 // Quick class for tracklets
 struct Tracklet {
diff --git a/x86/SciFi/LookingForward/include/LookingForwardStudies.h b/x86/SciFi/LookingForward/include/LookingForwardStudies.h
index d5c321362ff32655e03cbebcd37e4f1dc3e7a8e0..5cd489bdf6323a8015f7d73eb8e762ca0cbbae37 100644
--- a/x86/SciFi/LookingForward/include/LookingForwardStudies.h
+++ b/x86/SciFi/LookingForward/include/LookingForwardStudies.h
@@ -19,6 +19,7 @@
 #include "UTConsolidated.cuh"
 #include "SciFiDefinitions.cuh"
 #include "SciFiParametrization.h"
+#include "TrackUtils.cuh"
 
 #include "LookingForwardUtils.h"
 #include "LookingForwardConstants.h"
diff --git a/x86/SciFi/LookingForward/include/LookingForwardUtils.h b/x86/SciFi/LookingForward/include/LookingForwardUtils.h
index 73a0380587bf8179014ab2079a3c383de7e46a5e..8257b48e94b93d502ef89263e0f621c08ebf26e8 100644
--- a/x86/SciFi/LookingForward/include/LookingForwardUtils.h
+++ b/x86/SciFi/LookingForward/include/LookingForwardUtils.h
@@ -9,13 +9,13 @@
 
 #include <cassert>
 
+#include "HitUtils.cuh"
 #include "SciFiDefinitions.cuh"
 #include "SciFiEventModel.cuh"
 #include "LookingForwardConstants.h"
 #include "MomentumForwardUtils.h"
 #include "BinarySearch.cuh"
 #include "SciFiParametrization.h"
-#include "TrackUtils.cuh"
 #include "LookingForwardFitting.h"
 #include "TMVA_Forward.cuh"
 #include "TMVA_Forward_1.cuh"
@@ -225,3 +225,18 @@ void single_track_quality_update(
   const SciFi::Tracking::TMVA* tmva2,
   const SciFi::Hits& scifi_hits,
   const int event_offset);
+
+__host__ void collectAllXHits_proto_p(
+  const SciFi::Hits& scifi_hits,
+  const SciFi::HitCount& scifi_hit_count,
+  const SciFi::Tracking::Arrays* constArrays,
+  const float magnet_polarity,
+  const MiniState& velo_state,
+  const MiniState& UT_state,
+  const float qOverP,
+  int side, 
+  std::array<int, 2 * 6>& windows_x,
+  std::array<int, 2 * 6>& windows_uv,
+  std::array<float, 4 * 6>& parameters_uv,
+  const SciFiWindowsParams& window_params,
+  const std::array<int, 12> true_scifi_indices_per_layer);
diff --git a/x86/SciFi/LookingForward/src/LookingForwardStudies.cpp b/x86/SciFi/LookingForward/src/LookingForwardStudies.cpp
index 8743483b6f47e2c00095931ac597f228ea053c48..c9beadccb102633a2e222a893db392c0d8de6288 100644
--- a/x86/SciFi/LookingForward/src/LookingForwardStudies.cpp
+++ b/x86/SciFi/LookingForward/src/LookingForwardStudies.cpp
@@ -1,6 +1,4 @@
 #include "LookingForwardStudies.h"
-#include "TrackUtils.cuh"
-#include "FindXHits.cuh"
 #include "LookingForwardSbt.h"
 #include "LookingForwardConstants.cuh"
 #include <numeric>
diff --git a/x86/SciFi/LookingForward/src/LookingForwardUtils.cpp b/x86/SciFi/LookingForward/src/LookingForwardUtils.cpp
index fe1b294374c30de47c9999416f4e51ab5508a0ab..a25c77b0a4fd7628d26b3f3b335713d6b1802504 100644
--- a/x86/SciFi/LookingForward/src/LookingForwardUtils.cpp
+++ b/x86/SciFi/LookingForward/src/LookingForwardUtils.cpp
@@ -924,3 +924,161 @@ void filter_tracks_with_TMVA(
     selected_tracks.push_back(tracks[best_track]);
   }
 }
+
+
+__host__ void collectAllXHits_proto_p(
+  const SciFi::Hits& scifi_hits,
+  const SciFi::HitCount& scifi_hit_count,
+  const SciFi::Tracking::Arrays* constArrays,
+  const float magnet_polarity,
+  const MiniState& velo_state,
+  const MiniState& UT_state,
+  const float qOverP,
+  int side,
+  std::array<int, 2 * 6>& windows_x,
+  std::array<int, 2 * 6>& windows_uv,
+  std::array<float, 4 * 6>& parameters_uv,
+  const SciFiWindowsParams& window_params,
+  const std::array<int, 12> true_scifi_indices_per_layer)
+{
+  const float tx2 = velo_state.tx*velo_state.tx;
+  const float ty2 = velo_state.ty*velo_state.ty;
+  const float slope2 = tx2 + ty2;
+  const float pt = sqrtf(slope2 / (1.f + slope2) ) / fabsf(qOverP);
+  const float p = 1.f / std::abs(qOverP);
+
+  /* OPTIMIZE: possibly use these variables for wrong sign treatment */
+  // const float q = qOverP > 0.f ? 1.f : -1.f;
+  // const bool wSignTreatment = pt > SciFi::Tracking::wrongSignPT;
+  // float zMag = zMagnet(velo_state, constArrays);
+  // const float qop_WS = sqrtf(slope2 / (1.f + slope2) ) / pt;
+  // float dxRefWS = 0.f;
+  // if ( wSignTreatment ) {
+  //   dxRefWS = 0.9f * calcDxRef(SciFi::Tracking::wrongSignPT, velo_state);
+  // }
+  // const float dir = q * SciFi::Tracking::magscalefactor * (-1.f); // needed for wrong sign treatment
+  // const float xTolWS = dx_calc(velo_state, qop_WS, window_params);
+
+  //const float q = qOverP > 0.f ? 1.f : -1.f;
+  //const float dir = q * magnet_polarity * (-1.f);
+
+  //const bool wSignTreatment = pt > SciFi::Tracking::wrongSignPT;
+  float zMag = zMagnet(velo_state, constArrays);
+  const float qop_WS = sqrtf(slope2 / (1.f + slope2) ) / pt;
+  // float dxRefWS = 0.f;
+  // if ( wSignTreatment ) {
+  //   dxRefWS = 0.9f * calcDxRef(SciFi::Tracking::wrongSignPT, velo_state);
+  // }
+  const float xAtRef = xFromVelo(SciFi::Tracking::zReference, UT_state);
+  float xParams_seed[4] {xAtRef, UT_state.tx, 0, 0};
+
+  // use parametrization to propagate from UT to SciFi
+  const auto state_zRef = propagate_state_from_velo(UT_state, qOverP, 5);  // zRef is between layers 4 and 5
+  const float xTol = dx_calc(velo_state, qOverP, window_params);
+  int iZoneStartingPoint = side > 0 ? constArrays->zoneoffsetpar : 0;
+
+  for (unsigned int iZone = iZoneStartingPoint; iZone < iZoneStartingPoint + constArrays->zoneoffsetpar; iZone++) {
+    assert(iZone - iZoneStartingPoint < SciFi::Constants::n_zones);
+    assert(iZone - iZoneStartingPoint < 12);
+
+    const auto izone_rel = iZone - iZoneStartingPoint;
+    const float zZone = constArrays->xZone_zPos[izone_rel];
+
+    //const int layer = constArrays->xZones[iZone] / 2;
+    const float dz_x = (zZone - SciFi::Tracking::zReference);
+    const float xInZone = scifi_propagation(state_zRef.x, UT_state.tx, qOverP, dz_x);
+    const float yInZone = yFromVelo(zZone, velo_state);
+
+    const float xInZoneStraight = evalCubicParameterization(xParams_seed, zZone);
+
+    if (side > 0) {
+      if (
+        !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
+        !isInside(yInZone, SciFi::Tracking::yLim_Min, SciFi::Tracking::yLim_Max))
+        continue;
+    }
+    else {
+      if (
+        !isInside(xInZone, SciFi::Tracking::xLim_Min, SciFi::Tracking::xLim_Max) ||
+        !isInside(yInZone, side * SciFi::Tracking::yLim_Max, side * SciFi::Tracking::yLim_Min))
+        continue;
+    }
+
+    float xMin = xInZone - xTol;
+    float xMax = xInZone + xTol;
+
+    /* OPTIMIZE: how do we take care of wrong sign tracks with momentum windows? */
+    //float xTolWS = 0.0;
+    // if (wSignTreatment) {
+    //   // xTolWS = (zZone < SciFi::Tracking::zReference) ?
+    //   //   dxRefWS * zZone / SciFi::Tracking::zReference :
+    //   //   dxRefWS * (zZone - zMag) / (SciFi::Tracking::zReference - zMag);
+    //   // if (dir > 0) {
+    //   //   xMin = xInZone - xTolWS;
+    //   // }
+    //   // else {
+    //   //   xMax = xInZone + xTolWS;
+    //   // }
+
+    //   debug_cout << "\t before WS treatment: xMin = " << xMin << ", xMax = " << xMax << ", WS = " << int(wSignTreatment) << ", pt = " << pt << std::endl;
+    //   if (dir > 0) {
+    //     xMin = -1.f * xInZone - xTolWS;
+    //   }
+    //   else {
+    //     xMax = xInZone + xTolWS;
+    //   }
+    //   debug_cout << "\t after WS treatment: xMin = " << xMin << ", xMax = " << xMax << std::endl;
+    // }
+
+    // Get the hits within the bounds
+    assert(iZone < SciFi::Constants::n_layers);
+    assert(constArrays->xZones[iZone] < SciFi::Constants::n_zones);
+    int x_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->xZones[iZone]);
+    int x_zone_offset_end = x_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->xZones[iZone]);
+    const int itH = getLowerBound(scifi_hits.x0, xMin, x_zone_offset_begin, x_zone_offset_end);
+    const int itEnd = getLowerBound(scifi_hits.x0, xMax, x_zone_offset_begin, x_zone_offset_end);
+    assert(itH >= x_zone_offset_begin && itH <= x_zone_offset_end);
+    assert(itEnd >= x_zone_offset_begin && itEnd <= x_zone_offset_end);
+
+    windows_x[2*izone_rel] = itH;
+    windows_x[2*izone_rel+1] = itEnd - itH;
+
+    // Now match the stereo hits
+    const float this_uv_z = constArrays->uvZone_zPos[izone_rel];
+    const float dz_uv = this_uv_z - SciFi::Tracking::zReference;
+    const float yInUVZone = yFromVelo(this_uv_z, velo_state);
+    const float dx = yInUVZone * constArrays->uvZone_dxdy[izone_rel];
+    const float xPredUv = scifi_propagation(state_zRef.x, UT_state.tx, qOverP, dz_uv) - dx;
+    const int uv_layer = constArrays->uvZones[iZone] / 2;
+    // To Do: study this window
+    const float xBound = 70.f * window_params.extrapolation_stddev[uv_layer];
+    const float maxDx = xBound; // * ratio;
+
+    const float xMinUV = xPredUv - maxDx;
+    const float xMaxUV = xPredUv + maxDx;
+
+    // Get bounds in UV layers
+    // do one search on the same side as the x module
+    // if we are close to y = 0, also look within a region on the other side module ("triangle search")
+    assert(constArrays->uvZones[iZone] < SciFi::Constants::n_zones);
+    const int uv_zone_offset_begin = scifi_hit_count.zone_offset(constArrays->uvZones[iZone]);
+    const int uv_zone_offset_end =
+      uv_zone_offset_begin + scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone]);
+    /* OPTIMIZE: check how large the effect on the efficiency is to include the triangle hits */
+    //    const int triangleOffset = side > 0 ? -1 : 1;
+    //assert(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset] < SciFi::Constants::n_zones);
+    // const int triangle_zone_offset_begin =
+    //   scifi_hit_count.zone_offset(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset]);
+    // assert(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset] < SciFi::Constants::n_zones);
+    // const int triangle_zone_offset_end =
+    //   triangle_zone_offset_begin +
+    //   scifi_hit_count.zone_number_of_hits(constArrays->uvZones[iZone + constArrays->zoneoffsetpar * triangleOffset]);
+    int itUVStart = getLowerBound(scifi_hits.x0, xMinUV, uv_zone_offset_begin, uv_zone_offset_end);
+    int itUVEnd = getLowerBound(scifi_hits.x0, xMaxUV, uv_zone_offset_begin, uv_zone_offset_end);
+    //    int itUV2 = getLowerBound(scifi_hits.x0, xMinUV, triangle_zone_offset_begin, triangle_zone_offset_end);
+
+    windows_uv[2*izone_rel] = itUVStart;
+    windows_uv[2*izone_rel+1] = itUVEnd;
+
+  }
+}
diff --git a/x86/SciFi/MomentumForward/CMakeLists.txt b/x86/SciFi/MomentumForward/CMakeLists.txt
index 2c33bcd81955f33be0cf843370f4b1cafba6b4a2..ed98cb7844efbbd96b138cc5e80a67f1bd9c0b21 100644
--- a/x86/SciFi/MomentumForward/CMakeLists.txt
+++ b/x86/SciFi/MomentumForward/CMakeLists.txt
@@ -11,7 +11,7 @@ include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/UT/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/SciFi/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
+include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/classifiers/include)
 include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
 
 file(GLOB x86MomentumForward "src/*cpp" )
@@ -20,8 +20,6 @@ add_library(x86MomentumForward STATIC
   ${x86MomentumForward}
 )
 
-set_property(TARGET x86Forward PROPERTY
-             CUDA_SEPARABLE_COMPILATION ON)
 set_property(TARGET x86MomentumForward PROPERTY
              CUDA_RESOLVE_DEVICE_SYMBOLS ON)
 
diff --git a/x86/SciFi/MomentumForward/include/MomentumForwardConstants.h b/x86/SciFi/MomentumForward/include/MomentumForwardConstants.h
index 73bd1c667d810c4366d235e6888bb75b60fa15eb..a1cf1569f94d3182ba6d22af208ee9d5d4e379b5 100644
--- a/x86/SciFi/MomentumForward/include/MomentumForwardConstants.h
+++ b/x86/SciFi/MomentumForward/include/MomentumForwardConstants.h
@@ -1,7 +1,7 @@
 #pragma once
 
 #include "SystemOfUnits.h"
-#include "PrForwardConstants.cuh"
+#include "SciFiDefinitions.cuh"
 
 namespace SciFi {
   namespace MomentumForward {
diff --git a/x86/SciFi/PrForward/CMakeLists.txt b/x86/SciFi/PrForward/CMakeLists.txt
deleted file mode 100644
index 45de4aada34ff14cf68a4810aed675befdeb2404..0000000000000000000000000000000000000000
--- a/x86/SciFi/PrForward/CMakeLists.txt
+++ /dev/null
@@ -1,40 +0,0 @@
-include_directories(include)
-include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/looking_forward/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/checker/tracking/include)
-include_directories(${CMAKE_SOURCE_DIR}/main/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/PrVeloUT/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/UTDecoding/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/velo/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/UT/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/event_model/SciFi/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/utils/sorting)
-include_directories(${CMAKE_SOURCE_DIR}/stream/gear/include)
-include_directories(${CMAKE_SOURCE_DIR}/stream/setup/include)
-include_directories(${CMAKE_SOURCE_DIR}/x86/SciFi/LookingForward/include)
-include_directories(${CMAKE_SOURCE_DIR}/x86/SciFi/MomentumForward/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/utils/binary_search/include)
-
-
-file(GLOB x86Forward_src "src/*cpp")
-
-add_library(x86Forward STATIC
-  ${x86Forward_src}
-)
-
-set_property(TARGET x86Forward PROPERTY
-             CUDA_SEPARABLE_COMPILATION ON)
-set_property(TARGET x86Forward PROPERTY
-             CUDA_RESOLVE_DEVICE_SYMBOLS ON)
-
-if ( ROOT_FOUND )
-  target_compile_definitions(x86Forward PUBLIC WITH_ROOT)
-  target_include_directories(x86Forward BEFORE PRIVATE
-    ${ROOT_INCLUDE_DIRS}
-  )
-endif()
diff --git a/x86/SciFi/PrForward/include/PrForwardWrapper.h b/x86/SciFi/PrForward/include/PrForwardWrapper.h
deleted file mode 100644
index 058ffd002023cf0fdd09cd52d9c23f1cc847d669..0000000000000000000000000000000000000000
--- a/x86/SciFi/PrForward/include/PrForwardWrapper.h
+++ /dev/null
@@ -1,17 +0,0 @@
-#pragma once
-
-#include "PrForwardTools.cuh"
-
-void PrForwardWrapper(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const Velo::Consolidated::States& velo_states,
-  const uint event_tracks_offset,
-  const UT::Consolidated::Tracks& ut_tracks,
-  const int n_veloUT_tracks,
-  const SciFi::Tracking::TMVA* tmva1,
-  const SciFi::Tracking::TMVA* tmva2,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  SciFi::TrackHits outputTracks[SciFi::Constants::max_tracks],
-  uint* n_forward_tracks);
diff --git a/x86/SciFi/PrForward/include/RunForwardCPU.h b/x86/SciFi/PrForward/include/RunForwardCPU.h
deleted file mode 100644
index 8252af53601ef4ff25cb519433646e12ed33f40b..0000000000000000000000000000000000000000
--- a/x86/SciFi/PrForward/include/RunForwardCPU.h
+++ /dev/null
@@ -1,47 +0,0 @@
-#pragma once
-
-#include "Common.h"
-#include "VeloDefinitions.cuh"
-#include "TrackChecker.h"
-#include "PrForwardWrapper.h"
-#include "VeloEventModel.cuh"
-#include "CpuHandler.cuh"
-#include "PrepareTracks.h"
-#include "ArgumentsCommon.cuh"
-#include "ArgumentsVelo.cuh"
-#include "ArgumentsUT.cuh"
-#include "ArgumentsSciFi.cuh"
-
-int run_forward_on_CPU(
-  SciFi::TrackHits* host_scifi_tracks_events,
-  int* host_scifi_n_tracks,
-  const uint* host_scifi_hits,
-  const uint* host_scifi_hit_count,
-  const char* host_scifi_geometry,
-  const std::array<float, 9>& host_inv_clus_res,
-  const uint* host_velo_tracks_atomics,
-  const uint* host_velo_track_hit_number,
-  const char* host_velo_states,
-  const int* host_atomics_ut,
-  const uint* host_ut_track_hit_number,
-  const float* host_ut_qop,
-  const uint* host_ut_track_velo_indices,
-  const uint number_of_events);
-
-CPU_ALGORITHM(
-  run_forward_on_CPU,
-  cpu_scifi_pr_forward_t,
-  ARGUMENTS(
-    dev_scifi_tracks,
-    dev_atomics_scifi,
-    dev_scifi_selected_track_indices,
-    dev_scifi_hits,
-    dev_scifi_hit_count,
-    dev_atomics_velo,
-    dev_velo_track_hit_number,
-    dev_velo_states,
-    dev_atomics_ut,
-    dev_ut_track_hits,
-    dev_ut_track_hit_number,
-    dev_ut_qop,
-    dev_ut_track_velo_indices))
diff --git a/x86/SciFi/PrForward/src/PrForwardWrapper.cpp b/x86/SciFi/PrForward/src/PrForwardWrapper.cpp
deleted file mode 100644
index aeebab54bce5324fdc3a3e3e232774151bd131f2..0000000000000000000000000000000000000000
--- a/x86/SciFi/PrForward/src/PrForwardWrapper.cpp
+++ /dev/null
@@ -1,39 +0,0 @@
-#include "PrForwardWrapper.h"
-
-void PrForwardWrapper(
-  const SciFi::Hits& scifi_hits,
-  const SciFi::HitCount& scifi_hit_count,
-  const Velo::Consolidated::States& velo_states,
-  const uint event_tracks_offset,
-  const UT::Consolidated::Tracks& ut_tracks,
-  const int n_veloUT_tracks,
-  const SciFi::Tracking::TMVA* tmva1,
-  const SciFi::Tracking::TMVA* tmva2,
-  const SciFi::Tracking::Arrays* constArrays,
-  const float magnet_polarity,
-  SciFi::TrackHits outputTracks[SciFi::Constants::max_tracks],
-  uint* n_forward_tracks)
-{
-  // Loop over the veloUT input tracks
-  for (int i_veloUT_track = 0; i_veloUT_track < n_veloUT_tracks; ++i_veloUT_track) {
-    const float qop_ut = ut_tracks.qop[i_veloUT_track];
-
-    const int i_velo_track = ut_tracks.velo_track[i_veloUT_track];
-    const uint velo_states_index = event_tracks_offset + i_velo_track;
-    const MiniState velo_state = velo_states.get(velo_states_index);
-
-    find_forward_tracks(
-      scifi_hits,
-      scifi_hit_count,
-      qop_ut,
-      i_veloUT_track,
-      outputTracks,
-      n_forward_tracks,
-      n_veloUT_tracks,
-      tmva1,
-      tmva2,
-      constArrays,
-      magnet_polarity,
-      velo_state);
-  }
-}
diff --git a/x86/SciFi/PrForward/src/RunForwardCPU.cpp b/x86/SciFi/PrForward/src/RunForwardCPU.cpp
deleted file mode 100644
index 5e2c4383bb836005d0193d8c167e04a0af0c136c..0000000000000000000000000000000000000000
--- a/x86/SciFi/PrForward/src/RunForwardCPU.cpp
+++ /dev/null
@@ -1,158 +0,0 @@
-#include "RunForwardCPU.h"
-
-#ifdef WITH_ROOT
-#include "TH1D.h"
-#include "TFile.h"
-#include "TTree.h"
-#endif
-
-int run_forward_on_CPU(
-  SciFi::TrackHits* host_scifi_tracks,
-  int* host_scifi_n_tracks,
-  const uint* host_scifi_hits,
-  const uint* host_scifi_hit_count,
-  const char* host_scifi_geometry,
-  const std::array<float, 9>& host_inv_clus_res,
-  const uint* host_velo_tracks_atomics,
-  const uint* host_velo_track_hit_number,
-  const char* host_velo_states,
-  const int* host_atomics_ut,
-  const uint* host_ut_track_hit_number,
-  const float* host_ut_qop,
-  const uint* host_ut_track_velo_indices,
-  const uint number_of_events)
-{
-  // #ifdef WITH_ROOT
-  //   // Histograms only for checking and debugging
-  //   TFile* f = new TFile("../output/scifi.root", "RECREATE");
-  //   TTree* t_Forward_tracks = new TTree("Forward_tracks", "Forward_tracks");
-  //   TTree* t_statistics = new TTree("statistics", "statistics");
-  //   TTree* t_scifi_hits = new TTree("scifi_hits", "scifi_hits");
-  //   uint planeCode, LHCbID;
-  //   float x0, z0, w, dxdy, dzdy, yMin, yMax;
-  //   float qop;
-  //   int n_tracks;
-  //   float state_x, state_y, state_z, state_tx, state_ty;
-
-  //   t_Forward_tracks->Branch("qop", &qop);
-  //   t_Forward_tracks->Branch("state_x", &state_x);
-  //   t_Forward_tracks->Branch("state_y", &state_y);
-  //   t_Forward_tracks->Branch("state_z", &state_z);
-  //   t_Forward_tracks->Branch("state_tx", &state_tx);
-  //   t_Forward_tracks->Branch("state_ty", &state_ty);
-  //   t_statistics->Branch("n_tracks", &n_tracks);
-  //   t_scifi_hits->Branch("planeCode", &planeCode);
-  //   t_scifi_hits->Branch("LHCbID", &LHCbID);
-  //   t_scifi_hits->Branch("x0", &x0);
-  //   t_scifi_hits->Branch("z0", &z0);
-  //   t_scifi_hits->Branch("w", &w);
-  //   t_scifi_hits->Branch("dxdy", &dxdy);
-  //   t_scifi_hits->Branch("dzdy", &dzdy);
-  //   t_scifi_hits->Branch("yMin", &yMin);
-  //   t_scifi_hits->Branch("yMax", &yMax);
-  // #endif
-
-  // // to do: set from configuration
-  // const float magnet_polarity = -1.f;
-
-  //   for (uint i_event = 0; i_event < number_of_events; ++i_event) {
-
-  //     // Velo consolidated types
-  //     const Velo::Consolidated::Tracks velo_tracks {
-  //       (uint*) host_velo_tracks_atomics, (uint*) host_velo_track_hit_number, i_event, number_of_events};
-  //     const uint event_tracks_offset = velo_tracks.tracks_offset(i_event);
-  //     const Velo::Consolidated::States host_velo_states_event {(char*) host_velo_states,
-  //                                                              velo_tracks.total_number_of_tracks};
-
-  //     // UT consolidated types
-  //     UT::Consolidated::Tracks ut_tracks {(uint*) host_atomics_ut,
-  //                                         (uint*) host_ut_track_hit_number,
-  //                                         (float*) host_ut_qop,
-  //                                         (uint*) host_ut_track_velo_indices,
-  //                                         i_event,
-  //                                         number_of_events};
-  //     const int n_veloUT_tracks_event = ut_tracks.number_of_tracks(i_event);
-
-  //     // SciFi non-consolidated types
-  //     int* n_forward_tracks = host_scifi_n_tracks + i_event;
-  //     SciFi::TrackHits* scifi_tracks_event = host_scifi_tracks + i_event * SciFi::Constants::max_tracks;
-  //  MiniState* scifi_states_event;
-
-  //     const uint total_number_of_hits = host_scifi_hit_count[number_of_events *
-  //     SciFi::Constants::n_mat_groups_and_mats]; SciFi::HitCount scifi_hit_count {(uint32_t*) host_scifi_hit_count,
-  //     i_event};
-
-  //     const SciFi::SciFiGeometry scifi_geometry(host_scifi_geometry);
-
-  //     SciFi::Hits scifi_hits(
-  //       (uint*) host_scifi_hits,
-  //       total_number_of_hits,
-  //       &scifi_geometry,
-  //       reinterpret_cast<const float*>(host_inv_clus_res.data()));
-
-  // #ifdef WITH_ROOT
-  //     // store hit variables in tree
-  //     for (size_t zone = 0; zone < SciFi::Constants::n_zones; zone++) {
-  //       const auto zone_offset = scifi_hit_count.zone_offset(zone);
-  //       for (size_t hit = 0; hit < scifi_hit_count.zone_number_of_hits(zone); hit++) {
-  //         const auto hit_offset = zone_offset + hit;
-
-  //         planeCode = scifi_hits.planeCode(hit_offset);
-  //         LHCbID = scifi_hits.LHCbID(hit_offset);
-  //         x0 = scifi_hits.x0[hit_offset];
-  //         z0 = scifi_hits.z0[hit_offset];
-  //         w = scifi_hits.w(hit_offset);
-  //         dxdy = scifi_hits.dxdy(hit_offset);
-  //         dzdy = scifi_hits.dzdy(hit_offset);
-  //         yMin = scifi_hits.yMin(hit_offset);
-  //         yMax = scifi_hits.yMax(hit_offset);
-  //         t_scifi_hits->Fill();
-  //       }
-  //     }
-  // #endif
-
-  //     // initialize TMVA vars
-  //     SciFi::Tracking::TMVA tmva1;
-  //     SciFi::Tracking::TMVA1_Init(tmva1);
-  //     SciFi::Tracking::TMVA tmva2;
-  //     SciFi::Tracking::TMVA2_Init(tmva2);
-
-  //     SciFi::Tracking::Arrays constArrays;
-
-  //     PrForwardWrapper(
-  //       scifi_hits,
-  //       scifi_hit_count,
-  //       host_velo_states_event,
-  //       event_tracks_offset,
-  //       ut_tracks,
-  //       n_veloUT_tracks_event,
-  //       &tmva1,
-  //       &tmva2,
-  //       &constArrays,
-  //       magnet_polarity,
-  //       scifi_tracks_event,
-  //       (uint*) n_forward_tracks);
-
-  // #ifdef WITH_ROOT
-  //     // store qop in tree
-  //     for (int i_track = 0; i_track < *n_forward_tracks; ++i_track) {
-  //       qop = scifi_tracks_event[i_track].qop;
-  //       state_x = scifi_tracks_event[i_track].state.x;
-  //       state_y = scifi_tracks_event[i_track].state.y;
-  //       state_z = scifi_tracks_event[i_track].state.z;
-  //       state_tx = scifi_tracks_event[i_track].state.tx;
-  //       state_ty = scifi_tracks_event[i_track].state.ty;
-  //       t_Forward_tracks->Fill();
-  //     }
-  //     n_tracks = n_forward_tracks[i_event];
-  //     t_statistics->Fill();
-  // #endif
-  //   }
-
-  // #ifdef WITH_ROOT
-  //   f->Write();
-  //   f->Close();
-  // #endif
-
-  return 0;
-}
diff --git a/x86/global_event_cut/CMakeLists.txt b/x86/global_event_cut/CMakeLists.txt
index 0973aa29a48226ab1c3a4b132f74ea9cd8335a1c..31ae99b937555015d378a14724ab918576742b98 100644
--- a/x86/global_event_cut/CMakeLists.txt
+++ b/x86/global_event_cut/CMakeLists.txt
@@ -7,7 +7,6 @@ include_directories(${CMAKE_SOURCE_DIR}/stream/setup/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/velo/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/UT/common/include)
-include_directories(${CMAKE_SOURCE_DIR}/cuda/SciFi/PrForward/include)
 include_directories(${CMAKE_SOURCE_DIR}/cuda/global_event_cut/include)
 
 add_library(CpuGEC STATIC