SeedFinderOrthogonal.ipp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SeedFinderOrthogonal.ipp

// This file is part of the Acts project.
//
// Copyright (C) 2023 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "Acts/Geometry/Extent.hpp"
#include "Acts/Seeding/SeedFilter.hpp"
#include "Acts/Seeding/SeedFinder.hpp"
#include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp"
#include "Acts/Seeding/SeedFinderUtils.hpp"
#include "Acts/Utilities/BinningType.hpp"

#include <cmath>
#include <functional>
#include <numeric>
#include <type_traits>

namespace Acts {
template <typename external_spacepoint_t>
auto SeedFinderOrthogonal<external_spacepoint_t>::validTupleOrthoRangeLH(
    const internal_sp_t &low) const -> typename tree_t::range_t {
  float colMin = m_config.collisionRegionMin;
  float colMax = m_config.collisionRegionMax;
  float pL = low.phi();
  float rL = low.radius();
  float zL = low.z();

  typename tree_t::range_t res;

  /*
   * Cut: Ensure that we search only in φ_min ≤ φ ≤ φ_max, as defined by the
   * seeding configuration.
   */
  res[DimPhi].shrinkMin(m_config.phiMin);
  res[DimPhi].shrinkMax(m_config.phiMax);

  /*
   * Cut: Ensure that we search only in r ≤ r_max, as defined by the seeding
   * configuration.
   */
  res[DimR].shrinkMax(m_config.rMax);

  /*
   * Cut: Ensure that we search only in z_min ≤ z ≤ z_max, as defined by the
   * seeding configuration.
   */
  res[DimZ].shrinkMin(m_config.zMin);
  res[DimZ].shrinkMax(m_config.zMax);

  /*
   * Cut: Ensure that we search only in Δr_min ≤ r - r_L ≤ Δr_max, as defined
   * by the seeding configuration and the given lower spacepoint.
   */
  res[DimR].shrinkMin(rL + m_config.deltaRMinTopSP);
  res[DimR].shrinkMax(rL + m_config.deltaRMaxTopSP);

  /*
   * Cut: Now that we have constrained r, we can use that new r range to
   * further constrain z.
   */
  float zMax = (res[DimR].max() / rL) * (zL - colMin) + colMin;
  float zMin = colMax - (res[DimR].max() / rL) * (colMax - zL);

  /*
   * This cut only works if z_low is outside the collision region for z.
   */
  if (zL > colMin) {
    res[DimZ].shrinkMax(zMax);
  } else if (zL < colMax) {
    res[DimZ].shrinkMin(zMin);
  }

  /*
   * Cut: Shrink the z-range using the maximum cot(θ).
   */
  res[DimZ].shrinkMin(zL - m_config.cotThetaMax * (res[DimR].max() - rL));
  res[DimZ].shrinkMax(zL + m_config.cotThetaMax * (res[DimR].max() - rL));

  /*
   * Cut: Shrink the φ range, such that Δφ_min ≤ φ - φ_L ≤ Δφ_max
   */
  res[DimPhi].shrinkMin(pL - m_config.deltaPhiMax);
  res[DimPhi].shrinkMax(pL + m_config.deltaPhiMax);

  return res;
}

template <typename external_spacepoint_t>
auto SeedFinderOrthogonal<external_spacepoint_t>::validTupleOrthoRangeHL(
    const internal_sp_t &high) const -> typename tree_t::range_t {
  float pM = high.phi();
  float rM = high.radius();

  typename tree_t::range_t res;

  /*
   * Cut: Ensure that we search only in φ_min ≤ φ ≤ φ_max, as defined by the
   * seeding configuration.
   */
  res[DimPhi].shrinkMin(m_config.phiMin);
  res[DimPhi].shrinkMax(m_config.phiMax);

  /*
   * Cut: Ensure that we search only in r ≤ r_max, as defined by the seeding
   * configuration.
   */
  res[DimR].shrinkMax(m_config.rMax);

  /*
   * Cut: Ensure that we search only in z_min ≤ z ≤ z_max, as defined by the
   * seeding configuration.
   */
  res[DimZ].shrinkMin(m_config.zMin);
  res[DimZ].shrinkMax(m_config.zMax);

  /*
   * Cut: Ensure that we search only in Δr_min ≤ r_H - r ≤ Δr_max, as defined
   * by the seeding configuration and the given higher spacepoint.
   */
  res[DimR].shrinkMin(rM - m_config.deltaRMaxBottomSP);
  res[DimR].shrinkMax(rM - m_config.deltaRMinBottomSP);

  /*
   * Cut: Now that we have constrained r, we can use that new r range to
   * further constrain z.
   */
  float fracR = res[DimR].min() / rM;

  float zMin = (high.z() - m_config.collisionRegionMin) * fracR +
               m_config.collisionRegionMin;
  float zMax = (high.z() - m_config.collisionRegionMax) * fracR +
               m_config.collisionRegionMax;

  res[DimZ].shrinkMin(std::min(zMin, high.z()));
  res[DimZ].shrinkMax(std::max(zMax, high.z()));

  /*
   * Cut: Shrink the φ range, such that Δφ_min ≤ φ - φ_H ≤ Δφ_max
   */
  res[DimPhi].shrinkMin(pM - m_config.deltaPhiMax);
  res[DimPhi].shrinkMax(pM + m_config.deltaPhiMax);

  return res;
}

template <typename external_spacepoint_t>
bool SeedFinderOrthogonal<external_spacepoint_t>::validTuple(
    const SeedFinderOptions &options, const internal_sp_t &low,
    const internal_sp_t &high, bool isMiddleInverted) const {
  float rL = low.radius();
  float rH = high.radius();

  float zL = low.z();
  float zH = high.z();

  float deltaR = rH - rL;

  float deltaZ = (zH - zL);
  float cotTheta = deltaZ / deltaR;
  /*
   * Cut: Ensure that the origin of the dublet (the intersection of the line
   * between them with the z axis) lies within the collision region.
   */
  float zOrigin = zL - rL * cotTheta;
  if (zOrigin < m_config.collisionRegionMin ||
      zOrigin > m_config.collisionRegionMax) {
    return false;
  }

  // cut on the max curvature calculated from dublets
  // first transform the space point coordinates into a frame such that the
  // central space point SPm is in the origin of the frame and the x axis
  // points away from the interaction point in addition to a translation
  // transformation we also perform a rotation in order to keep the
  // curvature of the circle tangent to the x axis
  if (m_config.interactionPointCut) {
    float xVal = (high.x() - low.x()) * (low.x() / rL) +
                 (high.y() - low.y()) * (low.y() / rL);
    float yVal = (high.y() - low.y()) * (low.x() / rL) -
                 (high.x() - low.x()) * (low.y() / rL);

    const int sign = isMiddleInverted ? -1 : 1;

    if (std::abs(rL * yVal) > sign * m_config.impactMax * xVal) {
      // conformal transformation u=x/(x²+y²) v=y/(x²+y²) transform the
      // circle into straight lines in the u/v plane the line equation can
      // be described in terms of aCoef and bCoef, where v = aCoef * u +
      // bCoef
      float uT = xVal / (xVal * xVal + yVal * yVal);
      float vT = yVal / (xVal * xVal + yVal * yVal);
      // in the rotated frame the interaction point is positioned at x = -rM
      // and y ~= impactParam
      float uIP = -1. / rL;
      float vIP = m_config.impactMax / (rL * rL);
      if (sign * yVal > 0.) {
        vIP = -vIP;
      }
      // we can obtain aCoef as the slope dv/du of the linear function,
      // estimated using du and dv between the two SP bCoef is obtained by
      // inserting aCoef into the linear equation
      float aCoef = (vT - vIP) / (uT - uIP);
      float bCoef = vIP - aCoef * uIP;
      // the distance of the straight line from the origin (radius of the
      // circle) is related to aCoef and bCoef by d^2 = bCoef^2 / (1 +
      // aCoef^2) = 1 / (radius^2) and we can apply the cut on the curvature
      if ((bCoef * bCoef) > (1 + aCoef * aCoef) / options.minHelixDiameter2) {
        return false;
      }
    }
  }

  /*
   * Cut: Ensure that the forward angle (z / r) lies within reasonable bounds,
   * which is to say the absolute value must be smaller than the max cot(θ).
   */
  if (std::fabs(cotTheta) > m_config.cotThetaMax) {
    return false;
  }
  /*
   * Cut: Ensure that z-distance between SPs is within max and min values.
   */
  if (std::abs(deltaZ) > m_config.deltaZMax) {
    return false;
  }

  return true;
}

template <typename external_spacepoint_t>
SeedFinderOrthogonal<external_spacepoint_t>::SeedFinderOrthogonal(
    const SeedFinderOrthogonalConfig<external_spacepoint_t> &config)
    : m_config(config) {
  if (not config.isInInternalUnits) {
    throw std::runtime_error(
        "SeedFinderOrthogonalConfig not in ACTS internal units in "
        "SeedFinderOrthogonal");
  }
}

template <typename external_spacepoint_t>
void SeedFinderOrthogonal<external_spacepoint_t>::filterCandidates(
    const SeedFinderOptions &options, internal_sp_t &middle,
    std::vector<internal_sp_t *> &bottom, std::vector<internal_sp_t *> &top,
    SeedFilterState seedFilterState,
    CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>
        &candidates_collector,
    Acts::SpacePointData &spacePointData) const {
  float rM = middle.radius();
  float zM = middle.z();
  float varianceRM = middle.varianceR();
  float varianceZM = middle.varianceZ();

  // apply cut on the number of top SP if seedConfirmation is true
  if (m_config.seedConfirmation == true) {
    // check if middle SP is in the central or forward region
    SeedConfirmationRangeConfig seedConfRange =
        (zM > m_config.centralSeedConfirmationRange.zMaxSeedConf ||
         zM < m_config.centralSeedConfirmationRange.zMinSeedConf)
            ? m_config.forwardSeedConfirmationRange
            : m_config.centralSeedConfirmationRange;
    // set the minimum number of top SP depending on whether the middle SP is
    // in the central or forward region
    seedFilterState.nTopSeedConf = rM > seedConfRange.rMaxSeedConf
                                       ? seedConfRange.nTopForLargeR
                                       : seedConfRange.nTopForSmallR;
    // set max bottom radius for seed confirmation
    seedFilterState.rMaxSeedConf = seedConfRange.rMaxSeedConf;
    // continue if number of top SPs is smaller than minimum
    if (top.size() < seedFilterState.nTopSeedConf) {
      return;
    }
  }

  std::vector<const internal_sp_t *> top_valid;
  std::vector<float> curvatures;
  std::vector<float> impactParameters;

  // contains parameters required to calculate circle with linear equation
  // ...for bottom-middle
  std::vector<LinCircle> linCircleBottom;
  // ...for middle-top
  std::vector<LinCircle> linCircleTop;

  // transform coordinates
  transformCoordinates(spacePointData, bottom, middle, true, linCircleBottom);
  transformCoordinates(spacePointData, top, middle, false, linCircleTop);

  // sort: make index vector
  std::vector<std::size_t> sorted_bottoms(linCircleBottom.size());
  for (std::size_t i(0); i < sorted_bottoms.size(); ++i) {
    sorted_bottoms[i] = i;
  }

  std::vector<std::size_t> sorted_tops(linCircleTop.size());
  for (std::size_t i(0); i < sorted_tops.size(); ++i) {
    sorted_tops[i] = i;
  }

  std::sort(
      sorted_bottoms.begin(), sorted_bottoms.end(),
      [&linCircleBottom](const std::size_t a, const std::size_t b) -> bool {
        return linCircleBottom[a].cotTheta < linCircleBottom[b].cotTheta;
      });

  std::sort(sorted_tops.begin(), sorted_tops.end(),
            [&linCircleTop](const std::size_t a, const std::size_t b) -> bool {
              return linCircleTop[a].cotTheta < linCircleTop[b].cotTheta;
            });

  std::vector<float> tanMT;
  tanMT.reserve(top.size());

  size_t numTopSP = top.size();
  for (size_t t = 0; t < numTopSP; t++) {
    tanMT.push_back(
        std::atan2(top[t]->radius() - middle.radius(), top[t]->z() - zM));
  }

  size_t t0 = 0;

  for (const std::size_t b : sorted_bottoms) {
    // break if we reached the last top SP
    if (t0 == numTopSP) {
      break;
    }

    auto lb = linCircleBottom[b];

    // 1+(cot^2(theta)) = 1/sin^2(theta)
    float iSinTheta2 = (1. + lb.cotTheta * lb.cotTheta);
    float sigmaSquaredPtDependent = iSinTheta2 * options.sigmapT2perRadius;
    // calculate max scattering for min momentum at the seed's theta angle
    // scaling scatteringAngle^2 by sin^2(theta) to convert pT^2 to p^2
    // accurate would be taking 1/atan(thetaBottom)-1/atan(thetaTop) <
    // scattering
    // but to avoid trig functions we approximate cot by scaling by
    // 1/sin^4(theta)
    // resolving with pT to p scaling --> only divide by sin^2(theta)
    // max approximation error for allowed scattering angles of 0.04 rad at
    // eta=infinity: ~8.5%
    float scatteringInRegion2 = options.multipleScattering2 * iSinTheta2;

    // minimum number of compatible top SPs to trigger the filter for a certain
    // middle bottom pair if seedConfirmation is false we always ask for at
    // least one compatible top to trigger the filter
    size_t minCompatibleTopSPs = 2;
    if (!m_config.seedConfirmation or
        bottom[b]->radius() > seedFilterState.rMaxSeedConf) {
      minCompatibleTopSPs = 1;
    }
    if (m_config.seedConfirmation and seedFilterState.numQualitySeeds) {
      minCompatibleTopSPs++;
    }

    // clear all vectors used in each inner for loop
    top_valid.clear();
    curvatures.clear();
    impactParameters.clear();

    float tanLM = std::atan2(rM - bottom[b]->radius(), zM - bottom[b]->z());

    for (size_t index_t = t0; index_t < numTopSP; index_t++) {
      const std::size_t t = sorted_tops[index_t];
      auto lt = linCircleTop[t];

      if (std::abs(tanLM - tanMT[t]) > 0.005) {
        continue;
      }

      // add errors of spB-spM and spM-spT pairs and add the correlation term
      // for errors on spM
      float error2 = lt.Er + lb.Er +
                     2 * (lb.cotTheta * lt.cotTheta * varianceRM + varianceZM) *
                         lb.iDeltaR * lt.iDeltaR;

      float deltaCotTheta = lb.cotTheta - lt.cotTheta;
      float deltaCotTheta2 = deltaCotTheta * deltaCotTheta;

      // Apply a cut on the compatibility between the r-z slope of the two
      // seed segments. This is done by comparing the squared difference
      // between slopes, and comparing to the squared uncertainty in this
      // difference - we keep a seed if the difference is compatible within
      // the assumed uncertainties. The uncertainties get contribution from
      // the  space-point-related squared error (error2) and a scattering term
      // calculated assuming the minimum pt we expect to reconstruct
      // (scatteringInRegion2). This assumes gaussian error propagation which
      // allows just adding the two errors if they are uncorrelated (which is
      // fair for scattering and measurement uncertainties)
      if (deltaCotTheta2 > (error2 + scatteringInRegion2)) {
        // skip top SPs based on cotTheta sorting when producing triplets
        // break if cotTheta from bottom SP < cotTheta from top SP because
        // the SP are sorted by cotTheta
        if (deltaCotTheta < 0) {
          break;
        }
        t0 = index_t + 1;
        continue;
      }

      float dU = lt.U - lb.U;

      // A and B are evaluated as a function of the circumference parameters
      // x_0 and y_0
      float A = (lt.V - lb.V) / dU;
      float S2 = 1. + A * A;
      float B = lb.V - A * lb.U;
      float B2 = B * B;
      // sqrt(S2)/B = 2 * helixradius
      // calculated radius must not be smaller than minimum radius
      if (S2 < B2 * options.minHelixDiameter2) {
        continue;
      }

      // 1/helixradius: (B/sqrt(S2))*2 (we leave everything squared)
      // convert p(T) to p scaling by sin^2(theta) AND scale by 1/sin^4(theta)
      // from rad to deltaCotTheta
      float p2scatterSigma = B2 / S2 * sigmaSquaredPtDependent;
      if (!std::isinf(m_config.maxPtScattering)) {
        // if pT > maxPtScattering, calculate allowed scattering angle using
        // maxPtScattering instead of pt.
        if (B2 == 0 or options.pTPerHelixRadius * std::sqrt(S2 / B2) >
                           2. * m_config.maxPtScattering) {
          float pTscatterSigma =
              (m_config.highland / m_config.maxPtScattering) *
              m_config.sigmaScattering;
          p2scatterSigma = pTscatterSigma * pTscatterSigma * iSinTheta2;
        }
      }
      // if deltaTheta larger than allowed scattering for calculated pT, skip
      if (deltaCotTheta2 > (error2 + p2scatterSigma)) {
        if (deltaCotTheta < 0) {
          break;
        }
        t0 = index_t;
        continue;
      }

      // A and B allow calculation of impact params in U/V plane with linear
      // function
      // (in contrast to having to solve a quadratic function in x/y plane)
      float Im = std::abs((A - B * rM) * rM);

      if (Im <= m_config.impactMax) {
        top_valid.push_back(top[t]);
        // inverse diameter is signed depending on if the curvature is
        // positive/negative in phi
        curvatures.push_back(B / std::sqrt(S2));
        impactParameters.push_back(Im);
      }
    }

    // continue if number of top SPs is smaller than minimum required for filter
    if (top.size() < minCompatibleTopSPs) {
      continue;
    }

    seedFilterState.zOrigin = zM - rM * lb.cotTheta;

    m_config.seedFilter->filterSeeds_2SpFixed(
        spacePointData, *bottom[b], middle, top_valid, curvatures,
        impactParameters, seedFilterState, candidates_collector);

  }  // loop on bottoms
}

template <typename external_spacepoint_t>
template <typename output_container_t>
void SeedFinderOrthogonal<external_spacepoint_t>::processFromMiddleSP(
    const SeedFinderOptions &options, const tree_t &tree,
    output_container_t &out_cont, const typename tree_t::pair_t &middle_p,
    Acts::SpacePointData &spacePointData) const {
  using range_t = typename tree_t::range_t;
  internal_sp_t &middle = *middle_p.second;

  /*
   * Prepare four output vectors for seed candidates:
   *
   * bottom_lh_v denotes the candidates bottom seed points, assuming that the
   * track has monotonically _increasing_ z position. bottom_hl_v denotes the
   * candidate bottom points assuming that the track has monotonically
   * _decreasing_ z position. top_lh_v are the candidate top points for an
   * increasing z track, and top_hl_v are the candidate top points for a
   * decreasing z track.
   */
  std::vector<internal_sp_t *> bottom_lh_v, bottom_hl_v, top_lh_v, top_hl_v;

  /*
   * Storage for seed candidates
   */
  std::size_t max_num_quality_seeds_per_spm =
      m_config.seedFilter->getSeedFilterConfig().maxQualitySeedsPerSpMConf;
  std::size_t max_num_seeds_per_spm =
      m_config.seedFilter->getSeedFilterConfig().maxSeedsPerSpMConf;

  CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>
      candidates_collector;
  candidates_collector.setMaxElements(max_num_seeds_per_spm,
                                      max_num_quality_seeds_per_spm);

  /*
   * Calculate the search ranges for bottom and top candidates for this middle
   * space point.
   */
  range_t bottom_r = validTupleOrthoRangeHL(middle);
  range_t top_r = validTupleOrthoRangeLH(middle);

  /*
   * Calculate the value of cot(θ) for this middle spacepoint.
   */
  float myCotTheta =
      std::max(std::abs(middle.z() / middle.radius()), m_config.cotThetaMax);

  /*
   * Calculate the maximum Δr, given that we have already constrained our
   * search space.
   */
  float deltaRMaxTop = top_r[DimR].max() - middle.radius();
  float deltaRMaxBottom = middle.radius() - bottom_r[DimR].min();

  /*
   * Create the search range for the bottom spacepoint assuming a
   * monotonically increasing z track, by calculating the minimum z value from
   * the cot(θ), and by setting the maximum to the z position of the middle
   * spacepoint - if the z position is higher than the middle point, then it
   * would be a decreasing z track!
   */
  range_t bottom_lh_r = bottom_r;
  bottom_lh_r[DimZ].shrink(middle.z() - myCotTheta * deltaRMaxBottom,
                           middle.z());

  /*
   * Calculate the search ranges for the other four sets of points in a
   * similar fashion.
   */
  range_t top_lh_r = top_r;
  top_lh_r[DimZ].shrink(middle.z(), middle.z() + myCotTheta * deltaRMaxTop);

  range_t bottom_hl_r = bottom_r;
  bottom_hl_r[DimZ].shrink(middle.z(),
                           middle.z() + myCotTheta * deltaRMaxBottom);
  range_t top_hl_r = top_r;
  top_hl_r[DimZ].shrink(middle.z() - myCotTheta * deltaRMaxTop, middle.z());

  /*
   * Make sure the candidate vectors are clear, in case we've used them
   * before.
   */
  bottom_lh_v.clear();
  bottom_hl_v.clear();
  top_lh_v.clear();
  top_hl_v.clear();

  /*
   * Now, we will actually search for the spaces. Remembering that we combine
   * bottom and top candidates for increasing and decreasing tracks
   * separately, we will first check whether both the search ranges for
   * increasing tracks are not degenerate - if they are, we will never find
   * any seeds and we do not need to bother doing the search.
   */
  if (!bottom_lh_r.degenerate() && !top_lh_r.degenerate()) {
    /*
     * Search the trees for points that lie in the given search range.
     */
    tree.rangeSearchMapDiscard(top_lh_r,
                               [this, &options, &middle, &top_lh_v](
                                   const typename tree_t::coordinate_t &,
                                   const typename tree_t::value_t &top) {
                                 if (validTuple(options, *top, middle, true)) {
                                   top_lh_v.push_back(top);
                                 }
                               });
  }

  /*
   * Perform the same search for candidate bottom spacepoints, but for
   * monotonically decreasing z tracks.
   */
  if (!bottom_hl_r.degenerate() && !top_hl_r.degenerate()) {
    tree.rangeSearchMapDiscard(top_hl_r,
                               [this, &options, &middle, &top_hl_v](
                                   const typename tree_t::coordinate_t &,
                                   const typename tree_t::value_t &top) {
                                 if (validTuple(options, middle, *top, false)) {
                                   top_hl_v.push_back(top);
                                 }
                               });
  }

  // apply cut on the number of top SP if seedConfirmation is true
  SeedFilterState seedFilterState;
  bool search_bot_hl = true;
  bool search_bot_lh = true;
  if (m_config.seedConfirmation) {
    // check if middle SP is in the central or forward region
    SeedConfirmationRangeConfig seedConfRange =
        (middle.z() > m_config.centralSeedConfirmationRange.zMaxSeedConf ||
         middle.z() < m_config.centralSeedConfirmationRange.zMinSeedConf)
            ? m_config.forwardSeedConfirmationRange
            : m_config.centralSeedConfirmationRange;
    // set the minimum number of top SP depending on whether the middle SP is
    // in the central or forward region
    seedFilterState.nTopSeedConf = middle.radius() > seedConfRange.rMaxSeedConf
                                       ? seedConfRange.nTopForLargeR
                                       : seedConfRange.nTopForSmallR;
    // set max bottom radius for seed confirmation
    seedFilterState.rMaxSeedConf = seedConfRange.rMaxSeedConf;
    // continue if number of top SPs is smaller than minimum
    if (top_lh_v.size() < seedFilterState.nTopSeedConf) {
      search_bot_lh = false;
    }
    if (top_hl_v.size() < seedFilterState.nTopSeedConf) {
      search_bot_hl = false;
    }
  }

  /*
   * Next, we perform a search for bottom candidates in increasing z tracks,
   * which only makes sense if we found any bottom candidates.
   */
  if (!top_lh_v.empty() and search_bot_lh) {
    tree.rangeSearchMapDiscard(
        bottom_lh_r, [this, &options, &middle, &bottom_lh_v](
                         const typename tree_t::coordinate_t &,
                         const typename tree_t::value_t &bottom) {
          if (validTuple(options, *bottom, middle, false)) {
            bottom_lh_v.push_back(bottom);
          }
        });
  }

  /*
   * And repeat for the top spacepoints for decreasing z tracks!
   */
  if (!top_hl_v.empty() and search_bot_hl) {
    tree.rangeSearchMapDiscard(
        bottom_hl_r, [this, &options, &middle, &bottom_hl_v](
                         const typename tree_t::coordinate_t &,
                         const typename tree_t::value_t &bottom) {
          if (validTuple(options, middle, *bottom, true)) {
            bottom_hl_v.push_back(bottom);
          }
        });
  }

  /*
   * If we have candidates for increasing z tracks, we try to combine them.
   */
  if (!bottom_lh_v.empty() && !top_lh_v.empty()) {
    filterCandidates(options, middle, bottom_lh_v, top_lh_v, seedFilterState,
                     candidates_collector, spacePointData);
  }
  /*
   * Try to combine candidates for decreasing z tracks.
   */
  if (!bottom_hl_v.empty() && !top_hl_v.empty()) {
    filterCandidates(options, middle, bottom_hl_v, top_hl_v, seedFilterState,
                     candidates_collector, spacePointData);
  }
  /*
   * Run a seed filter, just like in other seeding algorithms.
   */
  if ((!bottom_lh_v.empty() && !top_lh_v.empty()) or
      (!bottom_hl_v.empty() && !top_hl_v.empty())) {
    m_config.seedFilter->filterSeeds_1SpFixed(
        spacePointData, candidates_collector, seedFilterState.numQualitySeeds,
        std::back_inserter(out_cont));
  }
}

template <typename external_spacepoint_t>
auto SeedFinderOrthogonal<external_spacepoint_t>::createTree(
    const std::vector<internal_sp_t *> &spacePoints) const -> tree_t {
  std::vector<typename tree_t::pair_t> points;

  /*
   * For every input point, we create a coordinate-pointer pair, which we then
   * linearly pass to the k-d tree constructor. That constructor will take
   * care of sorting the pairs and splitting the space.
   */
  for (internal_sp_t *sp : spacePoints) {
    typename tree_t::coordinate_t point;

    point[DimPhi] = sp->phi();
    point[DimR] = sp->radius();
    point[DimZ] = sp->z();

    points.emplace_back(point, sp);
  }

  return tree_t(std::move(points));
}

template <typename external_spacepoint_t>
template <typename input_container_t, typename output_container_t,
          typename callable_t>
void SeedFinderOrthogonal<external_spacepoint_t>::createSeeds(
    const Acts::SeedFinderOptions &options,
    const input_container_t &spacePoints, output_container_t &out_cont,
    callable_t &&extract_coordinates) const {
  if (not options.isInInternalUnits) {
    throw std::runtime_error(
        "SeedFinderOptions not in ACTS internal units in "
        "SeedFinderOrthogonal");
  }
  /*
   * The template parameters we accept are a little too generic, so we want to
   * run some basic checks to make sure the containers have the correct value
   * types.
   */
  static_assert(std::is_same_v<typename output_container_t::value_type,
                               Seed<external_spacepoint_t>>,
                "Output iterator container type must accept seeds.");
  static_assert(std::is_same_v<typename input_container_t::value_type,
                               const external_spacepoint_t *>,
                "Input container must contain external spacepoints.");

  /*
   * Sadly, for the time being, we will need to construct our internal space
   * points on the heap. This adds some additional overhead and work. Here we
   * take each external spacepoint, allocate a corresponding internal space
   * point, and save it in a vector.
   */
  Acts::Extent rRangeSPExtent;
  std::size_t counter = 0;
  std::vector<internal_sp_t *> internalSpacePoints;
  Acts::SpacePointData spacePointData;
  spacePointData.resize(spacePoints.size());

  for (const external_spacepoint_t *p : spacePoints) {
    auto [position, variance] = extract_coordinates(p);
    internalSpacePoints.push_back(new InternalSpacePoint<external_spacepoint_t>(
        counter++, *p, position, options.beamPos, variance));
    // store x,y,z values in extent
    rRangeSPExtent.extend(position);
  }
  // variable middle SP radial region of interest
  const Acts::Range1D<float> rMiddleSPRange(
      std::floor(rRangeSPExtent.min(Acts::binR) / 2) * 2 +
          m_config.deltaRMiddleMinSPRange,
      std::floor(rRangeSPExtent.max(Acts::binR) / 2) * 2 -
          m_config.deltaRMiddleMaxSPRange);

  /*
   * Construct the k-d tree from these points. Note that this not consume or
   * take ownership of the points.
   */
  tree_t tree = createTree(internalSpacePoints);
  /*
   * Run the seeding algorithm by iterating over all the points in the tree
   * and seeing what happens if we take them to be our middle spacepoint.
   */
  for (const typename tree_t::pair_t &middle_p : tree) {
    internal_sp_t &middle = *middle_p.second;
    auto rM = middle.radius();

    /*
     * Cut: Ensure that the middle spacepoint lies within a valid r-region for
     * middle points.
     */
    if (m_config.useVariableMiddleSPRange) {
      if (rM < rMiddleSPRange.min() || rM > rMiddleSPRange.max()) {
        continue;
      }
    } else {
      if (rM > m_config.rMaxMiddle || rM < m_config.rMinMiddle) {
        continue;
      }
    }

    // remove all middle SPs outside phi and z region of interest
    if (middle.z() < m_config.zOutermostLayers.first or
        middle.z() > m_config.zOutermostLayers.second) {
      continue;
    }
    float spPhi = middle.phi();
    if (spPhi > m_config.phiMax or spPhi < m_config.phiMin) {
      continue;
    }

    processFromMiddleSP(options, tree, out_cont, middle_p, spacePointData);
  }

  /*
   * Don't forget to get rid of all the spacepoints we just allocated!
   */
  for (const internal_sp_t *p : internalSpacePoints) {
    delete p;
  }
}

template <typename external_spacepoint_t>
template <typename input_container_t, typename callable_t>
std::vector<Seed<external_spacepoint_t>>
SeedFinderOrthogonal<external_spacepoint_t>::createSeeds(
    const Acts::SeedFinderOptions &options,
    const input_container_t &spacePoints,
    callable_t &&extract_coordinates) const {
  std::vector<seed_t> r;
  createSeeds(options, spacePoints, r,
              std::forward<callable_t>(extract_coordinates));
  return r;
}

}  // namespace Acts