TrackFindingAlg.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "src/TrackFindingAlg.h"
#include "Acts/Propagator/PropagatorOptions.hpp"
#include "src/detail/FitterHelperFunctions.h"
 
// Athena
#include "AsgTools/ToolStore.h"
#include "AthenaMonitoringKernel/Monitored.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkTrackSummary/TrackSummary.h"
#include "InDetPrepRawData/PixelClusterCollection.h"
#include "InDetPrepRawData/SCT_ClusterCollection.h"
#include "HGTD_PrepRawData/HGTD_ClusterCollection.h"
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "InDetRIO_OnTrack/PixelClusterOnTrack.h"
#include "InDetRIO_OnTrack/SCT_ClusterOnTrack.h"
#include "HGTD_RIO_OnTrack/HGTD_ClusterOnTrack.h"
 
// ACTS
#include "Acts/Definitions/Units.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/TrackFinding/MeasurementSelector.hpp"
#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Utilities/TrackHelpers.hpp"
#include "Acts/TrackFitting/MbfSmoother.hpp"
 
// ActsTrk
#include "ActsEvent/TrackContainer.h"
#include "ActsGeometry/ATLASMagneticFieldWrapper.h"
#include "ActsGeometryInterfaces/ActsGeometryContext.h"
#include "ActsInterop/Logger.h"
#include "ActsInterop/TableUtils.h"
#include "src/detail/AtlasMeasurementSelector.h"
#include "src/detail/OnTrackCalibrator.h"
 
// STL
#include <sstream>
#include <functional>
#include <tuple>
#include <utility>
#include <algorithm>
 
namespace ActsTrk
{
  struct TrackFindingAlg::CKF_pimpl : public detail::CKF_config
  {
  };
 
  TrackFindingAlg::CKF_pimpl &TrackFindingAlg::trackFinder() { return *m_trackFinder; }
  const TrackFindingAlg::CKF_pimpl &TrackFindingAlg::trackFinder() const { return *m_trackFinder; }
 
  TrackFindingAlg::TrackFindingAlg(const std::string &name,
                                   ISvcLocator *pSvcLocator)
      : AthReentrantAlgorithm(name, pSvcLocator)
  {
  }
 
  TrackFindingAlg::~TrackFindingAlg() = default;
 
  // === initialize ==========================================================
 
  StatusCode TrackFindingAlg::initialize()
  {
    ATH_MSG_INFO("Initializing " << name() << " ... ");
    ATH_MSG_DEBUG("Properties Summary:");
    ATH_MSG_DEBUG("   " << m_maxPropagationStep);
    ATH_MSG_DEBUG("   " << m_skipDuplicateSeeds);
    ATH_MSG_DEBUG("   " << m_refitSeeds);
    ATH_MSG_DEBUG("   " << m_etaBins);
    ATH_MSG_DEBUG("   " << m_chi2CutOff);
    ATH_MSG_DEBUG("   " << m_chi2OutlierCutOff);
    ATH_MSG_DEBUG("   " << m_numMeasurementsCutOff);
    ATH_MSG_DEBUG("   " << m_ptMinMeasurements);
    ATH_MSG_DEBUG("   " << m_absEtaMaxMeasurements);
    ATH_MSG_DEBUG("   " << m_doBranchStopper);
    ATH_MSG_DEBUG("   " << m_addPixelStripCounts);
    ATH_MSG_DEBUG("   " << m_doTwoWay);
    ATH_MSG_DEBUG("   " << m_reverseSearch);
    ATH_MSG_DEBUG("   " << m_phiMin);
    ATH_MSG_DEBUG("   " << m_phiMax);
    ATH_MSG_DEBUG("   " << m_etaMin);
    ATH_MSG_DEBUG("   " << m_etaMax);
    ATH_MSG_DEBUG("   " << m_absEtaMin);
    ATH_MSG_DEBUG("   " << m_absEtaMax);
    ATH_MSG_DEBUG("   " << m_ptMin);
    ATH_MSG_DEBUG("   " << m_ptMax);
    ATH_MSG_DEBUG("   " << m_z0Min);
    ATH_MSG_DEBUG("   " << m_z0Max);
    ATH_MSG_DEBUG("   " << m_minMeasurements);
    ATH_MSG_DEBUG("   " << m_minPixelHits);
    ATH_MSG_DEBUG("   " << m_minStripHits);
    ATH_MSG_DEBUG("   " << m_maxHoles);
    ATH_MSG_DEBUG("   " << m_maxPixelHoles);
    ATH_MSG_DEBUG("   " << m_maxStripHoles);
    ATH_MSG_DEBUG("   " << m_maxOutliers);
    ATH_MSG_DEBUG("   " << m_maxPixelOutliers);
    ATH_MSG_DEBUG("   " << m_maxStripOutliers);
    ATH_MSG_DEBUG("   " << m_maxSharedHits);
    ATH_MSG_DEBUG("   " << m_maxChi2);
    ATH_MSG_DEBUG("   " << m_statEtaBins);
    ATH_MSG_DEBUG("   " << m_seedLabels);
    ATH_MSG_DEBUG("   " << m_dumpAllStatEtaBins);
    ATH_MSG_DEBUG("   " << m_branchStopperPtMinFactor);
    ATH_MSG_DEBUG("   " << m_branchStopperAbsEtaMaxExtra);
 
    // Read and Write handles
    ATH_CHECK(m_seedContainerKeys.initialize());
    ATH_CHECK(m_uncalibratedMeasurementContainerKeys.initialize());
    ATH_CHECK(m_detectorElementToGeometryIdMapKey.initialize());
    ATH_CHECK(m_estimatedTrackParametersKeys.initialize());
    ATH_CHECK(m_trackContainerKey.initialize());
    ATH_CHECK(m_tracksBackendHandlesHelper.initialize(ActsTrk::prefixFromTrackContainerName(m_trackContainerKey.key())));
 
    if (m_estimatedTrackParametersKeys.size() != m_seedLabels.size())
    {
      ATH_MSG_FATAL("There are " << m_seedLabels.size() << " SeedLabels, but " << m_estimatedTrackParametersKeys.size() << " EstimatedTrackParametersKeys");
      return StatusCode::FAILURE;
    }
 
    if (m_seedContainerKeys.size() != m_estimatedTrackParametersKeys.size())
    {
      ATH_MSG_FATAL("There are " << m_estimatedTrackParametersKeys.size() << " EstimatedTrackParametersKeys, but " << m_seedContainerKeys.size() << " SeedContainerKeys");
      return StatusCode::FAILURE;
    }
 
    ATH_CHECK(m_monTool.retrieve(EnableTool{not m_monTool.empty()}));
    ATH_CHECK(m_trackingGeometryTool.retrieve());
    ATH_CHECK(m_extrapolationTool.retrieve());
    ATH_CHECK(m_trackStatePrinter.retrieve(EnableTool{not m_trackStatePrinter.empty()}));
    ATH_CHECK(m_fitterTool.retrieve());
    ATH_CHECK(m_pixelCalibTool.retrieve(EnableTool{not m_pixelCalibTool.empty()}));
    ATH_CHECK(m_stripCalibTool.retrieve(EnableTool{not m_stripCalibTool.empty()}));
    ATH_CHECK(m_hgtdCalibTool.retrieve(EnableTool{not m_hgtdCalibTool.empty()}));
 
    m_logger = makeActsAthenaLogger(this, "Acts");
 
    auto magneticField = std::make_unique<ATLASMagneticFieldWrapper>();
    auto trackingGeometry = m_trackingGeometryTool->trackingGeometry();
 
    detail::Stepper stepper(std::move(magneticField));
    detail::Navigator::Config cfg{trackingGeometry};
    cfg.resolvePassive = false;
    cfg.resolveMaterial = true;
    cfg.resolveSensitive = true;
    detail::Navigator navigator(cfg, logger().cloneWithSuffix("Navigator"));
    detail::Propagator propagator(std::move(stepper), std::move(navigator), logger().cloneWithSuffix("Prop"));
 
    // Using the CKF propagator as extrapolator
    detail::Extrapolator extrapolator = propagator;
 
    std::vector<double> etaBins;
    // m_etaBins (from flags.Tracking.ActiveConfig.etaBins) includes a dummy first and last bin, which we ignore
    if (m_etaBins.size() > 2) {
      etaBins.assign(m_etaBins.begin() + 1, m_etaBins.end() - 1);
    }
    Acts::MeasurementSelectorCuts measurementSelectorCuts{etaBins};
 
    if (!m_chi2CutOff.empty())
      measurementSelectorCuts.chi2CutOff = m_chi2CutOff;
    if (!m_chi2OutlierCutOff.empty() && m_chi2OutlierCutOff.size() == m_chi2CutOff.size())
      measurementSelectorCuts.chi2CutOffOutlier = m_chi2OutlierCutOff;
    if (!m_numMeasurementsCutOff.empty())
      measurementSelectorCuts.numMeasurementsCutOff = m_numMeasurementsCutOff;
 
    Acts::MeasurementSelector::Config measurementSelectorCfg{{Acts::GeometryIdentifier(), std::move(measurementSelectorCuts)}};
    Acts::MeasurementSelector measurementSelector(measurementSelectorCfg);
 
    std::vector<double> absEtaEdges;
    absEtaEdges.reserve(etaBins.size() + 2);
    if (etaBins.empty())
    {
      absEtaEdges.push_back(0.0);
      absEtaEdges.push_back(std::numeric_limits<double>::infinity());
    }
    else
    {
      absEtaEdges.push_back(m_absEtaMin);
      absEtaEdges.insert(absEtaEdges.end(), etaBins.begin(), etaBins.end());
      absEtaEdges.push_back(m_absEtaMax);
    }
 
    auto setCut = [](auto &cfgVal, const auto &cuts, size_t ind) -> void
    {
      if (cuts.empty())
        return;
      cfgVal = (ind < cuts.size()) ? cuts[ind] : cuts[cuts.size() - 1];
    };
 
    Acts::TrackSelector::EtaBinnedConfig trackSelectorCfg{std::move(absEtaEdges)};
    if (etaBins.empty())
    {
      assert(trackSelectorCfg.cutSets.size() == 1);
      trackSelectorCfg.cutSets[0].absEtaMin = m_absEtaMin;
      trackSelectorCfg.cutSets[0].absEtaMax = m_absEtaMax;
    }
    size_t cutIndex = 0;
    for (auto &cfg : trackSelectorCfg.cutSets)
    {
      setCut(cfg.phiMin, m_phiMin, cutIndex);
      setCut(cfg.phiMax, m_phiMax, cutIndex);
      setCut(cfg.etaMin, m_etaMin, cutIndex);
      setCut(cfg.etaMax, m_etaMax, cutIndex);
      setCut(cfg.ptMin, m_ptMin, cutIndex);
      setCut(cfg.ptMax, m_ptMax, cutIndex);
      setCut(cfg.loc0Min, m_d0Min, cutIndex);
      setCut(cfg.loc0Max, m_d0Max, cutIndex);
      setCut(cfg.loc1Min, m_z0Min, cutIndex);
      setCut(cfg.loc1Max, m_z0Max, cutIndex);
      setCut(cfg.minMeasurements, m_minMeasurements, cutIndex);
      setCut(cfg.maxHoles, m_maxHoles, cutIndex);
      setCut(cfg.maxOutliers, m_maxOutliers, cutIndex);
      setCut(cfg.maxSharedHits, m_maxSharedHits, cutIndex);
      setCut(cfg.maxChi2, m_maxChi2, cutIndex);
      ++cutIndex;
    }
 
    ATH_MSG_INFO(trackSelectorCfg);
 
    if (!m_useDefaultMeasurementSelector.value()) {
       // initializer measurement selector and connect it to the delegates of the track finder optins
       ATH_CHECK( initializeMeasurementSelector());
    }
    else if (!m_chi2OutlierCutOff.empty()) {
       ATH_MSG_DEBUG("chi2OutlierCutOff set but not supported when using the default measurement selector.");
    }
 
    detail::CKF_config ckfConfig{
        std::move(extrapolator),
        {std::move(propagator), logger().cloneWithSuffix("CKF")},
        measurementSelector,
        {},
        trackSelectorCfg};
 
    m_trackFinder = std::make_unique<CKF_pimpl>(std::move(ckfConfig));
 
    trackFinder().ckfExtensions.updater.connect<&ActsTrk::detail::FitterHelperFunctions::gainMatrixUpdate<detail::RecoTrackStateContainer>>();
    trackFinder().ckfExtensions.measurementSelector.connect<&Acts::MeasurementSelector::select<detail::RecoTrackStateContainer>>(&trackFinder().measurementSelector);
    initStatTables();
 
    return StatusCode::SUCCESS;
  }
 
  // === finalize ============================================================
 
  StatusCode TrackFindingAlg::finalize()
  {
    printStatTables();
    return StatusCode::SUCCESS;
  }
 
  // === execute =============================================================
 
  StatusCode TrackFindingAlg::execute(const EventContext &ctx) const
  {
    ATH_MSG_DEBUG("Executing " << name() << " ... ");
 
    auto timer = Monitored::Timer<std::chrono::milliseconds>("TIME_execute");
    auto mon_nTracks = Monitored::Scalar<int>("nTracks");
    auto mon = Monitored::Group(m_monTool, timer, mon_nTracks);
 
    // ================================================== //
    // ===================== INPUTS ===================== //
    // ================================================== //
 
    // SEED PARAMETERS
    std::vector<const ActsTrk::BoundTrackParametersContainer *> estimatedTrackParametersContainers;
    estimatedTrackParametersContainers.reserve(m_estimatedTrackParametersKeys.size());
    for (const auto &estimatedTrackParametersKey : m_estimatedTrackParametersKeys)
    {
      ATH_MSG_DEBUG("Reading input collection with key " << estimatedTrackParametersKey.key());
      SG::ReadHandle<ActsTrk::BoundTrackParametersContainer> estimatedTrackParametersHandle = SG::makeHandle(estimatedTrackParametersKey, ctx);
      ATH_CHECK(estimatedTrackParametersHandle.isValid());
      estimatedTrackParametersContainers.push_back(estimatedTrackParametersHandle.cptr());
      ATH_MSG_DEBUG("Retrieved " << estimatedTrackParametersContainers.back()->size() << " input elements from key " << estimatedTrackParametersKey.key());
    }
 
    // SEED TRIPLETS
    std::vector<const ActsTrk::SeedContainer *> seedContainers;
    seedContainers.reserve(m_seedContainerKeys.size());
    std::size_t total_seeds = 0;
    for (const auto &seedContainerKey : m_seedContainerKeys)
    {
      ATH_MSG_DEBUG("Reading input collection with key " << seedContainerKey.key());
      SG::ReadHandle<ActsTrk::SeedContainer> seedsHandle = SG::makeHandle(seedContainerKey, ctx);
      ATH_CHECK(seedsHandle.isValid());
      seedContainers.push_back(seedsHandle.cptr());
      ATH_MSG_DEBUG("Retrieved " << seedContainers.back()->size() << " input elements from key " << seedContainerKey.key());
      total_seeds += seedContainers.back()->size();
    }
 
    // MEASUREMENTS
    std::vector<const xAOD::UncalibratedMeasurementContainer *> uncalibratedMeasurementContainers;
    uncalibratedMeasurementContainers.reserve(m_uncalibratedMeasurementContainerKeys.size());
    for (const auto &uncalibratedMeasurementContainerKey : m_uncalibratedMeasurementContainerKeys)
    {
      ATH_MSG_DEBUG("Reading input collection with key " << uncalibratedMeasurementContainerKey.key());
      SG::ReadHandle<xAOD::UncalibratedMeasurementContainer> uncalibratedMeasurementContainerHandle = SG::makeHandle(uncalibratedMeasurementContainerKey, ctx);
      ATH_CHECK(uncalibratedMeasurementContainerHandle.isValid());
      uncalibratedMeasurementContainers.push_back(uncalibratedMeasurementContainerHandle.cptr());
      ATH_MSG_DEBUG("Retrieved " << uncalibratedMeasurementContainers.back()->size() << " input elements from key " << uncalibratedMeasurementContainerKey.key());
    }
 
    SG::ReadCondHandle<ActsTrk::DetectorElementToActsGeometryIdMap>
       detectorElementToGeometryIdMap{m_detectorElementToGeometryIdMapKey, ctx};
    ATH_CHECK(detectorElementToGeometryIdMap.isValid());
 
    detail::DuplicateSeedDetector duplicateSeedDetector(total_seeds, m_skipDuplicateSeeds);
    for (std::size_t icontainer = 0; icontainer < seedContainers.size(); ++icontainer)
    {
      if (!seedContainers[icontainer])
        continue;
      duplicateSeedDetector.addSeeds(icontainer, *seedContainers[icontainer]);
    }
 
    detail::TrackFindingMeasurements measurements(uncalibratedMeasurementContainers.size() /* number of measurement containers*/);
    const Acts::TrackingGeometry *
       acts_tracking_geometry = m_trackingGeometryTool->trackingGeometry().get();
    ATH_CHECK(acts_tracking_geometry != nullptr);
 
    for (std::size_t icontainer = 0; icontainer < uncalibratedMeasurementContainers.size(); ++icontainer) {
      ATH_MSG_DEBUG("Create " << uncalibratedMeasurementContainers[icontainer]->size() << " source links from measurements in " << m_uncalibratedMeasurementContainerKeys[icontainer].key());
      measurements.addMeasurements(icontainer,
                                   *uncalibratedMeasurementContainers[icontainer],
                                   **detectorElementToGeometryIdMap);
    }
 
    if (!m_trackStatePrinter.empty()) {
      m_trackStatePrinter->printMeasurements(ctx, uncalibratedMeasurementContainers, **detectorElementToGeometryIdMap, measurements.measurementOffsets());
    }
 
 
    // ================================================== //
    // ===================== COMPUTATION ================ //
    // ================================================== //
    ActsTrk::MutableTrackContainer tracksContainer;
    EventStats event_stat;
    event_stat.resize(m_stat.size());
 
    // Perform the track finding for all initial parameters.
    for (std::size_t icontainer = 0; icontainer < estimatedTrackParametersContainers.size(); ++icontainer)
    {
      if (estimatedTrackParametersContainers[icontainer]->empty())
        continue;
      ATH_CHECK(findTracks(ctx,
                           *acts_tracking_geometry,
                           **detectorElementToGeometryIdMap,
                           measurements,
                           duplicateSeedDetector,
                           *estimatedTrackParametersContainers[icontainer],
                           seedContainers[icontainer],
                           tracksContainer,
                           icontainer,
                           icontainer < m_seedLabels.size() ? m_seedLabels[icontainer].c_str() : m_seedContainerKeys[icontainer].key().c_str(),
                           event_stat));
    }
 
    ATH_MSG_DEBUG("    \\__ Created " << tracksContainer.size() << " tracks");
 
    mon_nTracks = tracksContainer.size();
 
    copyStats(event_stat);
 
    std::unique_ptr<ActsTrk::TrackContainer> constTracksContainer = m_tracksBackendHandlesHelper.moveToConst(std::move(tracksContainer), 
      m_trackingGeometryTool->getGeometryContext(ctx).context(), ctx);
    // ================================================== //
    // ===================== OUTPUTS ==================== //
    // ================================================== //
    auto trackContainerHandle = SG::makeHandle(m_trackContainerKey, ctx);
    ATH_MSG_DEBUG("    \\__ Tracks Container `" << m_trackContainerKey.key() << "` created ...");
 
    ATH_CHECK(trackContainerHandle.record(std::move(constTracksContainer)));
    if (!trackContainerHandle.isValid())
    {
      ATH_MSG_FATAL("Failed to write TrackContainer with key " << m_trackContainerKey.key());
      return StatusCode::FAILURE;
    }
 
    return StatusCode::SUCCESS;
  }
 
  // === findTracks ==========================================================
 
  StatusCode
  TrackFindingAlg::findTracks(const EventContext &ctx,
                              const Acts::TrackingGeometry &trackingGeometry,
                              const ActsTrk::DetectorElementToActsGeometryIdMap &detectorElementToGeoId,
                              const detail::TrackFindingMeasurements &measurements,
                              detail::DuplicateSeedDetector &duplicateSeedDetector,
                              const ActsTrk::BoundTrackParametersContainer &estimatedTrackParameters,
                              const ActsTrk::SeedContainer *seeds,
                              ActsTrk::MutableTrackContainer &tracksContainer,
                              size_t typeIndex,
                              const char *seedType,
                              EventStats &event_stat) const
  {
    ATH_MSG_DEBUG(name() << "::" << __FUNCTION__);
 
    if (seeds && seeds->size() != estimatedTrackParameters.size())
    {
      // should be the same, but we can cope if not
      ATH_MSG_WARNING("Have " << seeds->size() << " " << seedType << " seeds, but " << estimatedTrackParameters.size() << "estimated track parameters");
    }
 
    // Construct a perigee surface as the target surface
    auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3::Zero());
 
    Acts::GeometryContext tgContext = m_trackingGeometryTool->getGeometryContext(ctx).context();
    Acts::MagneticFieldContext mfContext = m_extrapolationTool->getMagneticFieldContext(ctx);
    // CalibrationContext converter not implemented yet.
    Acts::CalibrationContext calContext = Acts::CalibrationContext();
 
    using AtlUncalibSourceLinkAccessor = detail::UncalibSourceLinkAccessor;
 
    AtlUncalibSourceLinkAccessor slAccessor(measurements.measurementRanges());
    Acts::SourceLinkAccessorDelegate<detail::UncalibSourceLinkAccessor::Iterator> slAccessorDelegate;
    slAccessorDelegate.connect<&detail::UncalibSourceLinkAccessor::range>(&slAccessor);
 
    Acts::PropagatorPlainOptions plainOptions{tgContext, mfContext};
    Acts::PropagatorPlainOptions plainSecondOptions{tgContext, mfContext};
 
    const bool reverseSearch = (typeIndex < m_reverseSearch.size() && m_reverseSearch[typeIndex]);
    plainOptions.maxSteps = m_maxPropagationStep;
    plainOptions.direction = reverseSearch ? Acts::Direction::Backward : Acts::Direction::Forward;
    plainSecondOptions.maxSteps = m_maxPropagationStep;
    plainSecondOptions.direction = plainOptions.direction.invert();
 
    // Set the CombinatorialKalmanFilter options
    using TrackFinderOptions = Acts::CombinatorialKalmanFilterOptions<detail::UncalibSourceLinkAccessor::Iterator, detail::RecoTrackContainer>;
    TrackFinderOptions options(tgContext,
                               mfContext,
                               calContext,
                               slAccessorDelegate,
                               trackFinder().ckfExtensions,
                               plainOptions,
                               pSurface.get());
    if (reverseSearch) options.targetSurface = pSurface.get();
    if (!m_useDefaultMeasurementSelector.value()) {
       m_measurementSelector->connect( &options.trackStateCandidateCreator );
    }
    std::optional<TrackFinderOptions> secondOptions;
    if (m_doTwoWay) {
      secondOptions.emplace(tgContext,
                            mfContext,
                            calContext,
                            slAccessorDelegate,
                            trackFinder().ckfExtensions,
                            plainSecondOptions,
                            pSurface.get());
      if (!reverseSearch) secondOptions->targetSurface = pSurface.get();
      if (!m_useDefaultMeasurementSelector.value()) {
         m_measurementSelector->connect( &secondOptions->trackStateCandidateCreator);
      }
      secondOptions->skipPrePropagationUpdate = true;
    }
 
    // ActsTrk::MutableTrackContainer tracksContainerTemp;
    Acts::VectorTrackContainer trackBackend;
    Acts::VectorMultiTrajectory trackStateBackend;
    detail::RecoTrackContainer tracksContainerTemp(trackBackend, trackStateBackend);
 
    if (m_addPixelStripCounts) {
      addPixelStripCounts(tracksContainerTemp);
    }
 
    // Measurement calibration
    // N.B. OnTrackCalibrator expects disabled tool handles when no calibration is requested.
    // Therefore, passing them without checking if they are enabled is safe.
 
    auto calibrator = detail::OnTrackCalibrator<detail::RecoTrackStateContainer>(trackingGeometry,
                                                                             detectorElementToGeoId,
                                                                             m_pixelCalibTool,
                                                                             m_stripCalibTool,
                                                                             m_hgtdCalibTool);
 
    if (m_useDefaultMeasurementSelector.value()) {
      // for default measurement selector need connect calibrator
      options.extensions.calibrator.connect<&detail::OnTrackCalibrator<detail::RecoTrackStateContainer>::calibrate>(&calibrator);
      if (m_doTwoWay) {
    secondOptions->extensions.calibrator.connect<&detail::OnTrackCalibrator<detail::RecoTrackStateContainer>::calibrate>(&calibrator);
      }
    }
 
    const auto &trackSelectorCfg = trackFinder().trackSelector.config();
    auto getCuts = [&trackSelectorCfg](double eta) -> const Acts::TrackSelector::Config & {
      // return the last bin for |eta|>=4 or nan
      return (!(std::abs(eta) < trackSelectorCfg.absEtaEdges.back())) ? trackSelectorCfg.cutSets.back()
             : (std::abs(eta) < trackSelectorCfg.absEtaEdges.front()) ? trackSelectorCfg.cutSets.front()
                                                                      : trackSelectorCfg.getCuts(eta);
    };
 
    std::size_t category_i = 0;
    const auto measurementContainerOffsets = measurements.measurementContainerOffsets();
 
    using BranchStopperResult = Acts::CombinatorialKalmanFilterBranchStopperResult;
    auto stopBranch = [&](const detail::RecoTrackContainer::TrackProxy &track,
                          const detail::RecoTrackContainer::TrackStateProxy &trackState) -> BranchStopperResult {
 
      if (m_addPixelStripCounts) {
        updatePixelStripCounts(track, trackState.typeFlags(), measurementType(trackState));
        checkPixelStripCounts(track);
      }
 
      if (!m_trackStatePrinter.empty()) {
        m_trackStatePrinter->printTrackState(tgContext, trackState, measurementContainerOffsets, true);
      }
 
      if (!m_doBranchStopper)
        return BranchStopperResult::Continue;
 
      const auto &parameters = trackState.hasFiltered() ? trackState.filtered() : trackState.predicted();
      double eta = -std::log(std::tan(0.5 * parameters[Acts::eBoundTheta]));
      const auto &cutSet = getCuts(eta);
 
      if (typeIndex < m_ptMinMeasurements.size() &&
          !(track.nMeasurements() < m_ptMinMeasurements[typeIndex])) {
        double pT = std::sin(parameters[Acts::eBoundTheta]) / parameters[Acts::eBoundQOverP];
        if (std::abs(pT) < cutSet.ptMin * m_branchStopperPtMinFactor) {
          ++event_stat[category_i][kNStoppedTracksMinPt];
          ATH_MSG_DEBUG("CkfBranchStopper: drop branch with q*pT="
                        << pT << " after "
                        << track.nMeasurements() << " measurements");
          return BranchStopperResult::StopAndDrop;
        }
      }
 
      if (typeIndex < m_absEtaMaxMeasurements.size() &&
          !(track.nMeasurements() < m_absEtaMaxMeasurements[typeIndex]) &&
          !(std::abs(eta) < trackSelectorCfg.absEtaEdges.back() + m_branchStopperAbsEtaMaxExtra)) {
        ++event_stat[category_i][kNStoppedTracksMaxEta];
        ATH_MSG_DEBUG("CkfBranchStopper: drop branch with eta="
                      << eta << " after "
                      << track.nMeasurements() << " measurements");
        return BranchStopperResult::StopAndDrop;
      }
 
      bool enoughMeasurements = (track.nMeasurements() >= cutSet.minMeasurements);
      bool tooManyHoles = (track.nHoles() > cutSet.maxHoles);
      bool tooManyOutliers = (track.nOutliers() > cutSet.maxOutliers);
 
      if (m_addPixelStripCounts) {
        auto [enoughMeasurementsPS, tooManyHolesPS, tooManyOutliersPS] = selectPixelStripCounts(track, eta);
        enoughMeasurements = enoughMeasurements && enoughMeasurementsPS;
        tooManyHoles = tooManyHoles || tooManyHolesPS;
        tooManyOutliers = tooManyOutliers || tooManyOutliersPS;
      }
 
      if (!(tooManyHoles || tooManyOutliers))
        return BranchStopperResult::Continue;
 
      if (!enoughMeasurements)
        ++event_stat[category_i][kNStoppedTracksMaxHoles];
      if (m_addPixelStripCounts) {
        ATH_MSG_DEBUG("CkfBranchStopper: stop and "
                      << (enoughMeasurements ? "keep" : "drop")
                      << " branch with nHoles=" << track.nHoles()
                      << " (" << s_branchState.nPixelHoles(track)
                      << " pixel+" << s_branchState.nStripHoles(track)
                      << " strip), nOutliers=" << track.nOutliers()
                      << " (" << s_branchState.nPixelOutliers(track)
                      << "+" << s_branchState.nStripOutliers(track)
                      << "), nMeasurements=" << track.nMeasurements()
                      << " (" << s_branchState.nPixelHits(track)
                      << "+" << s_branchState.nStripHits(track)
                      << ")");
      } else {
        ATH_MSG_DEBUG("CkfBranchStopper: stop and "
                      << (enoughMeasurements ? "keep" : "drop")
                      << " branch with nHoles=" << track.nHoles()
                      << ", nOutliers=" << track.nOutliers()
                      << ", nMeasurements=" << track.nMeasurements());
      }
 
      return enoughMeasurements ? BranchStopperResult::StopAndKeep
                                : BranchStopperResult::StopAndDrop;
    };
 
    options.extensions.branchStopper.connect(stopBranch);
    if (m_doTwoWay) {
      // If we extend the track inwards we can utilize the CKF propagation
      // as extrapolator and therefore do not need to stop branches.
      // This only works for outside-in extensions.
      if (reverseSearch)
        secondOptions->extensions.branchStopper.connect(stopBranch);
    }
 
    Acts::PropagatorOptions<detail::Stepper::Options, detail::Navigator::Options,
                            Acts::ActorList<Acts::MaterialInteractor>>
    extrapolationOptions(tgContext, mfContext);
 
    Acts::TrackExtrapolationStrategy extrapolationStrategy =
        Acts::TrackExtrapolationStrategy::first;
 
    // Perform the track finding for all initial parameters
    ATH_MSG_DEBUG("Invoke track finding with " << estimatedTrackParameters.size() << ' ' << seedType << " seeds.");
 
    std::size_t nPrinted = 0;
    auto printSeed = [&](std::size_t iseed, const Acts::BoundTrackParameters &seedParameters, bool isKF = false)
    {
      if (m_trackStatePrinter.empty() || !seeds)
        return;
      if (!nPrinted++)
      {
        ATH_MSG_INFO("CKF results for " << estimatedTrackParameters.size() << ' ' << seedType << " seeds:");
      }
      m_trackStatePrinter->printSeed(tgContext, *(*seeds)[iseed], seedParameters, measurementContainerOffsets, iseed, isKF);
    };
 
    // Loop over the track finding results for all initial parameters
    for (std::size_t iseed = 0; iseed < estimatedTrackParameters.size(); ++iseed)
    {
      category_i = typeIndex * (m_statEtaBins.size() + 1);
      tracksContainerTemp.clear();
 
      if (!estimatedTrackParameters[iseed])
      {
        ATH_MSG_WARNING("No " << seedType << " seed " << iseed);
        ++event_stat[category_i][kNoTrackParam];
        continue;
      }
 
      const Acts::BoundTrackParameters *initialParameters = estimatedTrackParameters[iseed];
      printSeed(iseed, *initialParameters);
 
      double etaInitial = -std::log(std::tan(0.5 * initialParameters->theta()));
      category_i = getStatCategory(typeIndex, etaInitial);
      ++event_stat[category_i][kNTotalSeeds];
 
      if (duplicateSeedDetector.isDuplicate(typeIndex, iseed))
      {
        ATH_MSG_DEBUG("skip " << seedType << " seed " << iseed << " - already found");
        ++event_stat[category_i][kNDuplicateSeeds];
        continue;
      }
 
      // Get the Acts tracks, given this seed
      // Result here contains a vector of TrackProxy objects
      ++event_stat[category_i][kNUsedSeeds];
 
      std::unique_ptr<Acts::BoundTrackParameters> seedParameters;
      const bool refitSeeds = (typeIndex < m_refitSeeds.size() && m_refitSeeds[typeIndex]);
      if (refitSeeds)
      {
        // Perform KF before CKF
        const auto fittedSeedCollection = m_fitterTool->fit(ctx, *(*seeds)[iseed], *initialParameters,
                                                            tgContext, mfContext, calContext,
                                                            detectorElementToGeoId);
        if (not fittedSeedCollection)
        {
          ATH_MSG_WARNING("KF Fitted Track is nullptr");
        }
        else if (fittedSeedCollection->size() != 1)
        {
          ATH_MSG_WARNING("KF produced " << fittedSeedCollection->size() << " tracks but should produce 1!");
        }
        else
        {
          // Check pTmin requirement
          const auto fittedSeed = fittedSeedCollection->getTrack(0);
 
          double etaSeed = -std::log(std::tan(0.5 * fittedSeed.parameters()[Acts::eBoundTheta]));
          const auto &cutSet = getCuts(etaSeed);
          if (fittedSeed.transverseMomentum() < cutSet.ptMin)
          {
            ATH_MSG_VERBOSE("min pt requirement not satisfied after param refinement: pt min is " << cutSet.ptMin << " but Refined params have pt of " << fittedSeed.transverseMomentum());
            ++event_stat[category_i][kNRejectedRefinedSeeds];
            continue;
          }
 
          seedParameters.reset(new Acts::BoundTrackParameters(fittedSeed.referenceSurface().getSharedPtr(),
                                                              fittedSeed.parameters(),
                                                              initialParameters->covariance(),
                                                              fittedSeed.particleHypothesis()));
          printSeed(iseed, *seedParameters, true);
 
          // Pass the refined params to the CKF
          initialParameters = seedParameters.get();
        }
      }
 
      auto result = trackFinder().ckf.findTracks(*initialParameters, options, tracksContainerTemp);
 
      // The result for this seed
      if (not result.ok()) {
        ATH_MSG_WARNING("Track finding failed for " << seedType << " seed " << iseed << " with error" << result.error());
        continue;
      }
      auto &tracksForSeed = result.value();
 
      size_t ntracks = 0;
 
      // lambda to collect together all the things we do with a viable track.
      auto addTrack = [&](detail::RecoTrackContainerProxy &track) {
        // if the the perigeeSurface was not hit (in particular the case for the inside-out pass,
        // the track has no reference surface and the extrapolation to the perigee has not been done
        // yet.
        if (!track.hasReferenceSurface()) {
           auto extrapolationResult = Acts::extrapolateTrackToReferenceSurface(
                   track, *pSurface, trackFinder().extrapolator, extrapolationOptions,
                   extrapolationStrategy, logger());
           if (!extrapolationResult.ok()) {
              ATH_MSG_WARNING("Extrapolation for seed "
                              << iseed << " and " << track.index()
                              << " failed with error " << extrapolationResult.error()
                              << " dropping track candidate.");
              return;
           }
        }
 
        Acts::trimTrack(track, true, true, true);
        Acts::calculateTrackQuantities(track);
        if (m_addPixelStripCounts) {
          initPixelStripCounts(track);
          for (const auto& trackState : track.trackStatesReversed()) {
            updatePixelStripCounts(track, trackState.typeFlags(), measurementType(trackState));
          }
          checkPixelStripCounts(track);
        }
 
        if (!m_trackStatePrinter.empty()) {
          m_trackStatePrinter->printTrack(tgContext, tracksContainerTemp, track, measurementContainerOffsets);
        }
 
        ++ntracks;
        ++event_stat[category_i][kNOutputTracks];
 
        // copy selected tracks into output tracksContainer
        auto selectPixelStripCountsFinal = [this](const detail::RecoTrackContainer::TrackProxy &track) {
          if (!m_addPixelStripCounts) return true;
          double eta = -std::log(std::tan(0.5 * track.theta()));
          auto [enoughMeasurementsPS, tooManyHolesPS, tooManyOutliersPS] = selectPixelStripCounts(track, eta);
          return enoughMeasurementsPS && !tooManyHolesPS && !tooManyOutliersPS;
        };
        if (trackFinder().trackSelector.isValidTrack(track) &&
            selectPixelStripCountsFinal(track)) {
          auto destProxy = tracksContainer.getTrack(tracksContainer.addTrack());
          destProxy.copyFrom(track, true);  // make sure we copy track states!
          ++event_stat[category_i][kNSelectedTracks];
 
          // Fill the track infos into the duplicate seed detector
          if (m_skipDuplicateSeeds) {
            storeSeedInfo(tracksContainerTemp, track, duplicateSeedDetector);
          }
        } else {
          ATH_MSG_DEBUG("Track " << ntracks << " from " << seedType << " seed " << iseed << " failed track selection");
        }
      };
 
      std::size_t nfirst = 0;
      for (auto &firstTrack : tracksForSeed) {
        std::size_t nsecond = 0;
 
        auto smoothingResult = Acts::smoothTrack(tgContext, firstTrack, logger(), Acts::MbfSmoother());
        if (!smoothingResult.ok()) {
          ATH_MSG_DEBUG("Smoothing for seed "
                     << iseed << " and first track " << firstTrack.index()
                     << " failed with error " << smoothingResult.error());
          continue;
        }
 
        if (m_doTwoWay) {
          std::optional<detail::RecoTrackStateContainerProxy> firstMeasurement;
          for (auto st : firstTrack.trackStatesReversed()) {
            bool isMeasurement = st.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag);
            bool isOutlier = st.typeFlags().test(Acts::TrackStateFlag::OutlierFlag);
            // We are excluding non measurement states and outlier here. Those can
            // decrease resolution because only the smoothing corrected the very
            // first prediction as filtering is not possible.
            if (isMeasurement && !isOutlier)
              firstMeasurement = st;
          }
 
          if (firstMeasurement.has_value()) {
            Acts::BoundTrackParameters secondInitialParameters = firstTrack.createParametersFromState(*firstMeasurement);
 
            if (!secondInitialParameters.referenceSurface().insideBounds(secondInitialParameters.localPosition())) {  // #3751
              ATH_MSG_DEBUG("Smoothing of first pass fit produced out-of-bounds parameters relative to the surface, '"
                            << secondInitialParameters.referenceSurface().name()
                            << "'. Skipping second pass for " << seedType << " seed " << iseed << " track " << nfirst);
            } else {
              auto rootBranch = tracksContainerTemp.makeTrack();
              rootBranch.copyFrom(firstTrack, false);  // #3534
              if (m_addPixelStripCounts)
                copyPixelStripCounts(rootBranch, firstTrack);
              auto secondResult = trackFinder().ckf.findTracks(secondInitialParameters, *secondOptions, tracksContainerTemp, rootBranch);
 
              if (not secondResult.ok()) {
                ATH_MSG_WARNING("Second track finding failed for " << seedType << " seed " << iseed << " track " << nfirst << " with error" << secondResult.error());
              } else {
                // store the original previous state to restore it later
                auto originalFirstMeasurementPrevious = firstMeasurement->previous();
 
                auto &secondTracksForSeed = secondResult.value();
                for (auto &secondTrack : secondTracksForSeed) {
                  secondTrack.reverseTrackStates(true);
 
                  firstMeasurement->previous() = secondTrack.outermostTrackState().index();
                  secondTrack.tipIndex() = firstTrack.tipIndex();
 
                  if (reverseSearch) {
                    // smooth the full track
                    auto secondSmoothingResult = Acts::smoothTrack(tgContext, secondTrack, logger());
                    if (!secondSmoothingResult.ok()) {
                      ATH_MSG_WARNING("Second smoothing for seed " << iseed << " and track " << secondTrack.index() << " failed with error " << secondSmoothingResult.error());
                      continue;
                    }
 
                    secondTrack.reverseTrackStates(true);
                  }
 
                  addTrack(secondTrack);
 
                  ++nsecond;
                }
 
                // restore the original previous state for the first track
                firstMeasurement->previous() = originalFirstMeasurementPrevious;
              }
            }
          }
        }
        if (nsecond == 0) {
          if (m_doTwoWay) {
            ATH_MSG_DEBUG("No viable result from second track finding for " << seedType << " seed " << iseed << " track " << nfirst);
            ++event_stat[category_i][kNoSecond];
          }
 
          addTrack(firstTrack);
        }
        nfirst++;
      }
      if (ntracks == 0) {
        ATH_MSG_DEBUG("Track finding found no track candidates for " << seedType << " seed " << iseed);
        ++event_stat[category_i][kNoTrack];
      } else if (ntracks >= 2) {
        ++event_stat[category_i][kMultipleBranches];
      }
      if (!m_trackStatePrinter.empty())
        std::cout << std::flush;
    }
 
    ATH_MSG_DEBUG("Completed " << seedType << " track finding with " << computeStatSum(typeIndex, kNOutputTracks, event_stat) << " track candidates.");
 
    return StatusCode::SUCCESS;
  }
 
  void 
  TrackFindingAlg::storeSeedInfo(const detail::RecoTrackContainer &tracksContainer,
                                 const detail::RecoTrackContainerProxy &track,
                                 detail::DuplicateSeedDetector &duplicateSeedDetector) const {
 
      const auto lastMeasurementIndex = track.tipIndex();
      duplicateSeedDetector.newTrajectory();
 
      tracksContainer.trackStateContainer().visitBackwards(
          lastMeasurementIndex,
          [&duplicateSeedDetector](const detail::RecoTrackStateContainer::ConstTrackStateProxy &state) -> void
          {
            // Check there is a source link
            if (not state.hasUncalibratedSourceLink())
              return;
 
            // Fill the duplicate selector
            auto sl = state.getUncalibratedSourceLink().template get<ATLASUncalibSourceLink>();
            duplicateSeedDetector.addMeasurement(sl);
          }); // end visitBackwards
  }
 
  xAOD::UncalibMeasType
  TrackFindingAlg::measurementType (const detail::RecoTrackContainer::TrackStateProxy &trackState)
  {
    if (trackState.hasReferenceSurface()) {
      if (const auto *actsDetElem = dynamic_cast<const ActsDetectorElement *>(trackState.referenceSurface().associatedDetectorElement())) {
        if (const auto *detElem = dynamic_cast<const InDetDD::SiDetectorElement *>(actsDetElem->upstreamDetectorElement())) {
          if (detElem->isPixel()) {
            return xAOD::UncalibMeasType::PixelClusterType;
          } else if (detElem->isSCT()) {
            return xAOD::UncalibMeasType::StripClusterType;
          }
        }
      }
    }
    return xAOD::UncalibMeasType::Other;
  }
 
  void
  TrackFindingAlg::addPixelStripCounts(detail::RecoTrackContainer& tracksContainer)
  {
    tracksContainer.addColumn<unsigned int>("nPixelHits");
    tracksContainer.addColumn<unsigned int>("nStripHits");
    tracksContainer.addColumn<unsigned int>("nPixelHoles");
    tracksContainer.addColumn<unsigned int>("nStripHoles");
    tracksContainer.addColumn<unsigned int>("nPixelOutliers");
    tracksContainer.addColumn<unsigned int>("nStripOutliers");
  }
 
  void
  TrackFindingAlg::initPixelStripCounts(const detail::RecoTrackContainer::TrackProxy &track)
  {
    s_branchState.nPixelHits(track) = 0;
    s_branchState.nStripHits(track) = 0;
    s_branchState.nPixelHoles(track) = 0;
    s_branchState.nStripHoles(track) = 0;
    s_branchState.nPixelOutliers(track) = 0;
    s_branchState.nStripOutliers(track) = 0;
  }
 
  void
  TrackFindingAlg::updatePixelStripCounts(const detail::RecoTrackContainer::TrackProxy &track,
                                          Acts::ConstTrackStateType typeFlags,
                                          xAOD::UncalibMeasType detType)
  {
    if (detType == xAOD::UncalibMeasType::PixelClusterType) {
      if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
        s_branchState.nPixelHoles(track)++;
      } else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
        s_branchState.nPixelOutliers(track)++;
      } else if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
        s_branchState.nPixelHits(track)++;
      }
    } else if (detType == xAOD::UncalibMeasType::StripClusterType) {
      if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
        s_branchState.nStripHoles(track)++;
      } else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
        s_branchState.nStripOutliers(track)++;
      } else if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
        s_branchState.nStripHits(track)++;
      }
    }
  }
 
  void
  TrackFindingAlg::copyPixelStripCounts(const detail::RecoTrackContainer::TrackProxy &track,
                                        const detail::RecoTrackContainer::TrackProxy &other)
  {
    s_branchState.nPixelHits(track) = s_branchState.nPixelHits(other);
    s_branchState.nStripHits(track) = s_branchState.nStripHits(other);
    s_branchState.nPixelHoles(track) = s_branchState.nPixelHoles(other);
    s_branchState.nStripHoles(track) = s_branchState.nStripHoles(other);
    s_branchState.nPixelOutliers(track) = s_branchState.nPixelOutliers(other);
    s_branchState.nStripOutliers(track) = s_branchState.nStripOutliers(other);
  }
 
  void
  TrackFindingAlg::checkPixelStripCounts(const detail::RecoTrackContainer::TrackProxy &track) const
  {
    // This check will fail if there are other types (HGTD, MS?) of hits, holes, or outliers.
    // The check can be removed when it is no longer appropriate.
    if (track.nMeasurements() != s_branchState.nPixelHits(track) + s_branchState.nStripHits(track))
      ATH_MSG_WARNING("mismatched hit count: total (" << track.nMeasurements()
                      << ") != pixel (" << s_branchState.nPixelHits(track)
                      << ") + strip (" << s_branchState.nStripHits(track) << ")");
    if (track.nHoles() < s_branchState.nPixelHoles(track) + s_branchState.nStripHoles(track))  // allow extra HGTD holes
      ATH_MSG_WARNING("mismatched hole count: total (" << track.nHoles()
                      << ") < pixel (" << s_branchState.nPixelHoles(track)
                      << ") + strip (" << s_branchState.nStripHoles(track) << ")");
    if (track.nOutliers() != s_branchState.nPixelOutliers(track) + s_branchState.nStripOutliers(track))
      ATH_MSG_WARNING("mismatched outlier count: total (" << track.nOutliers()
                      << ") != pixel (" << s_branchState.nPixelOutliers(track)
                      << ") + strip (" << s_branchState.nStripOutliers(track) << ")");
  };
 
  std::array<bool, 3>
  TrackFindingAlg::selectPixelStripCounts(const detail::RecoTrackContainer::TrackProxy &track, double eta) const
  {
    bool enoughMeasurements = true, tooManyHoles = false, tooManyOutliers = false;
    const auto &trackSelectorCfg = trackFinder().trackSelector.config();
    std::size_t etaBin = (std::abs(eta) < trackSelectorCfg.absEtaEdges.front())   ? 0
                         : (std::abs(eta) >= trackSelectorCfg.absEtaEdges.back()) ? trackSelectorCfg.absEtaEdges.size() - 1
                                                                                  : trackSelectorCfg.binIndex(eta);
    auto cutMin = [etaBin](std::size_t val, const std::vector<std::size_t> &cutSet) {
      return !cutSet.empty() && (val < (etaBin < cutSet.size() ? cutSet[etaBin] : cutSet.back()));
    };
    auto cutMax = [etaBin](std::size_t val, const std::vector<std::size_t> &cutSet) {
      return !cutSet.empty() && (val > (etaBin < cutSet.size() ? cutSet[etaBin] : cutSet.back()));
    };
    enoughMeasurements = enoughMeasurements && !cutMin(s_branchState.nPixelHits(track), m_minPixelHits);
    enoughMeasurements = enoughMeasurements && !cutMin(s_branchState.nStripHits(track), m_minStripHits);
    tooManyHoles = tooManyHoles || cutMax(s_branchState.nPixelHoles(track), m_maxPixelHoles);
    tooManyHoles = tooManyHoles || cutMax(s_branchState.nStripHoles(track), m_maxStripHoles);
    tooManyOutliers = tooManyOutliers || cutMax(s_branchState.nPixelOutliers(track), m_maxPixelOutliers);
    tooManyOutliers = tooManyOutliers || cutMax(s_branchState.nStripOutliers(track), m_maxStripOutliers);
    return {enoughMeasurements, tooManyHoles, tooManyOutliers};
  }
 
  // === Statistics printout =================================================
 
  void TrackFindingAlg::initStatTables()
  {
    if (!m_statEtaBins.empty())
    {
      m_useAbsEtaForStat = (m_statEtaBins[0] > 0.);
      float last_eta = m_statEtaBins[0];
      for (float eta : m_statEtaBins)
      {
        if (eta < last_eta)
        {
          ATH_MSG_FATAL("Eta bins for statistics counter not in ascending order.");
        }
        last_eta = eta;
      }
    }
    m_stat.resize(nSeedCollections() * seedCollectionStride());
  }
 
  // copy statistics
  void TrackFindingAlg::copyStats(const EventStats &event_stat) const
  {
    std::lock_guard<std::mutex> lock(m_mutex);
    std::size_t category_i = 0;
    for (const std::array<unsigned int, kNStat> &src_stat : event_stat)
    {
      std::array<std::size_t, kNStat> &dest_stat = m_stat[category_i++];
      for (std::size_t i = 0; i < src_stat.size(); ++i)
      {
        assert(i < dest_stat.size());
        dest_stat[i] += src_stat[i];
      }
    }
  }
 
  // print statistics
  void TrackFindingAlg::printStatTables() const
  {
    if (msgLvl(MSG::INFO))
    {
      std::vector<std::string> stat_labels =
          TableUtils::makeLabelVector(kNStat,
                                      {
                                          std::make_pair(kNTotalSeeds, "Input seeds"),
                                          std::make_pair(kNoTrackParam, "No track parameters"),
                                          std::make_pair(kNUsedSeeds, "Used   seeds"),
                                          std::make_pair(kNoTrack, "Cannot find track"),
                                          std::make_pair(kNDuplicateSeeds, "Duplicate seeds"),
                                          std::make_pair(kNRejectedRefinedSeeds, "Rejected refined parameters"),
                                          std::make_pair(kNOutputTracks, "CKF tracks"),
                                          std::make_pair(kNSelectedTracks, "selected tracks"),
                                          std::make_pair(kNStoppedTracksMaxHoles, "Stopped tracks reaching max holes"),
                                          std::make_pair(kMultipleBranches, "Seeds with more than one branch"),
                                          std::make_pair(kNoSecond, "Tracks failing second CKF"),
                                          std::make_pair(kNStoppedTracksMinPt, "Stopped tracks below pT cut"),
                                          std::make_pair(kNStoppedTracksMaxEta, "Stopped tracks above max eta")
                                      });
      assert(stat_labels.size() == kNStat);
      std::vector<std::string> categories;
      categories.reserve(m_seedLabels.size() + 1);
      categories.insert(categories.end(), m_seedLabels.begin(), m_seedLabels.end());
      categories.push_back("ALL");
 
      std::vector<std::string> eta_labels;
      eta_labels.reserve(m_statEtaBins.size() + 2);
      for (std::size_t eta_bin_i = 0; eta_bin_i < m_statEtaBins.size() + 2; ++eta_bin_i)
      {
        eta_labels.push_back(TableUtils::makeEtaBinLabel(m_statEtaBins, eta_bin_i, m_useAbsEtaForStat));
      }
 
      // vector used as 3D array stat[ eta_bin ][ stat_i ][ seed_type]
      // stat_i = [0, kNStat)
      // eta_bin = [0, m_statEtaBins.size()+2 ); eta_bin == m_statEtaBinsSize()+1 means sum of all etaBins
      // seed_type = [0, nSeedCollections()+1)  seed_type == nSeedCollections() means sum of all seed collections
      std::vector<std::size_t> stat =
          TableUtils::createCounterArrayWithProjections<std::size_t>(nSeedCollections(),
                                                                     m_statEtaBins.size() + 1,
                                                                     m_stat);
 
      // the extra columns and rows for the projections are addeded internally:
      std::size_t stat_stride =
          TableUtils::counterStride(nSeedCollections(),
                                    m_statEtaBins.size() + 1,
                                    kNStat);
      std::size_t eta_stride =
          TableUtils::subCategoryStride(nSeedCollections(),
                                        m_statEtaBins.size() + 1,
                                        kNStat);
      std::stringstream table_out;
 
      if (m_dumpAllStatEtaBins.value())
      {
        // dump for each counter a table with one row per eta bin
        std::size_t max_label_width = TableUtils::maxLabelWidth(stat_labels) + TableUtils::maxLabelWidth(eta_labels);
        for (std::size_t stat_i = 0; stat_i < kNStat; ++stat_i)
        {
          std::size_t dest_idx_offset = stat_i * stat_stride;
          table_out << makeTable(stat, dest_idx_offset, eta_stride,
                                 eta_labels,
                                 categories)
                           .columnWidth(10)
                           // only dump the footer for the last eta bin i.e. total
                           .dumpHeader(stat_i == 0)
                           .dumpFooter(stat_i + 1 == kNStat)
                           .separateLastRow(true) // separate the sum of all eta bins
                           .minLabelWidth(max_label_width)
                           .labelPrefix(stat_labels.at(stat_i));
        }
      }
      else
      {
        // dump one table with one row per counter showing the total eta range
        for (std::size_t eta_bin_i = (m_dumpAllStatEtaBins.value() ? 0 : m_statEtaBins.size() + 1);
             eta_bin_i < m_statEtaBins.size() + 2;
             ++eta_bin_i)
        {
          std::size_t dest_idx_offset = eta_bin_i * eta_stride;
          table_out << makeTable(stat, dest_idx_offset, stat_stride,
                                 stat_labels,
                                 categories,
                                 eta_labels.at(eta_bin_i))
                           .columnWidth(10)
                           // only dump the footer for the last eta bin i.e. total
                           .dumpFooter(!m_dumpAllStatEtaBins.value() || eta_bin_i == m_statEtaBins.size() + 1);
        }
      }
      ATH_MSG_INFO("statistics:\n"
                   << table_out.str());
      table_out.str("");
 
      // define retios first element numerator, second element denominator
      // each element contains a vector of counter and a multiplier e.g. +- 1
      // ratios are computed as  (sum_i stat[stat_i] *  multiplier_i ) / (sum_j stat[stat_j] *  multiplier_j )
      auto [ratio_labels, ratio_def] =
          TableUtils::splitRatioDefinitionsAndLabels({TableUtils::makeRatioDefinition("failed / seeds ",
                                                                                      std::vector<TableUtils::SummandDefinition>{
                                                                                          TableUtils::defineSummand(kNTotalSeeds, 1),
                                                                                          TableUtils::defineSummand(kNUsedSeeds, -1),
                                                                                          TableUtils::defineSummand(kNDuplicateSeeds, -1),
                                                                                          // no track counted  as used but want to include it as failed
                                                                                          TableUtils::defineSummand(kNoTrack, 1),
                                                                                      }, // failed seeds i.e. seeds which are not duplicates but did not produce a track
                                                                                      std::vector<TableUtils::SummandDefinition>{TableUtils::defineSummand(kNTotalSeeds, 1)}),
                                                      TableUtils::defineSimpleRatio("duplication / seeds", kNDuplicateSeeds, kNTotalSeeds),
                                                      TableUtils::defineSimpleRatio("Rejected refined params / seeds", kNRejectedRefinedSeeds, kNTotalSeeds),
                                                      TableUtils::defineSimpleRatio("selected / CKF tracks", kNSelectedTracks, kNOutputTracks),
                                                      TableUtils::defineSimpleRatio("selected tracks / used seeds", kNSelectedTracks, kNUsedSeeds),
                                                      TableUtils::defineSimpleRatio("branched tracks / used seeds", kMultipleBranches, kNUsedSeeds),
                                                      TableUtils::defineSimpleRatio("no 2nd CKF / CKF tracks", kNoSecond, kNOutputTracks)});
 
      std::vector<float> ratio = TableUtils::computeRatios(ratio_def,
                                                           nSeedCollections() + 1,
                                                           m_statEtaBins.size() + 2,
                                                           stat);
 
      // the extra columns and rows for the projections are _not_ added internally
      std::size_t ratio_stride = TableUtils::ratioStride(nSeedCollections() + 1,
                                                         m_statEtaBins.size() + 2,
                                                         ratio_def);
      std::size_t ratio_eta_stride = TableUtils::subCategoryStride(nSeedCollections() + 1,
                                                                   m_statEtaBins.size() + 2,
                                                                   ratio_def);
 
      std::size_t max_label_width = TableUtils::maxLabelWidth(ratio_labels) + TableUtils::maxLabelWidth(eta_labels);
      if (m_dumpAllStatEtaBins.value())
      {
        // show for each ratio a table with one row per eta bin
        for (std::size_t ratio_i = 0; ratio_i < ratio_labels.size(); ++ratio_i)
        {
          table_out << makeTable(ratio,
                                 ratio_i * ratio_stride,
                                 ratio_eta_stride,
                                 eta_labels,
                                 categories)
                           .columnWidth(10)
                           // only dump the footer for the last eta bin i.e. total
                           .dumpHeader(ratio_i == 0)
                           .dumpFooter(ratio_i + 1 == ratio_labels.size())
                           .separateLastRow(true) // separate the sum of las
                           .minLabelWidth(max_label_width)
                           .labelPrefix(ratio_labels.at(ratio_i));
        }
      }
      else
      {
        // dump one table with one row per ratio showing  the total eta range
        table_out << makeTable(ratio,
                               (m_statEtaBins.size() + 1) * ratio_eta_stride + 0 * ratio_stride,
                               ratio_stride,
                               ratio_labels,
                               categories)
                         .columnWidth(10)
                         // only dump the footer for the last eta bin i.e. total
                         .minLabelWidth(max_label_width)
                         .dumpFooter(false);
 
        // also dump a table for final tracks over seeds (ratio_i==3) showing one row per eta bin
        eta_labels.erase(eta_labels.end() - 1); // drop last line of table which shows again all eta bins summed.
        constexpr std::size_t ratio_i = 3;
        table_out << makeTable(ratio,
                               ratio_i * ratio_stride,
                               ratio_eta_stride,
                               eta_labels,
                               categories)
                         .columnWidth(10)
                         .dumpHeader(false)
                         // only dump the footer for the last eta bin i.e. total
                         .dumpFooter(!m_dumpAllStatEtaBins.value() || ratio_i + 1 == ratio_labels.size())
                         .separateLastRow(false)
                         .minLabelWidth(max_label_width)
                         .labelPrefix(ratio_labels.at(ratio_i));
      }
 
      ATH_MSG_INFO("Ratios:\n"
                   << table_out.str());
    }
  }
 
  inline std::size_t TrackFindingAlg::getStatCategory(std::size_t seed_collection, float eta) const
  {
    std::vector<float>::const_iterator bin_iter = std::upper_bound(m_statEtaBins.begin(),
                                                                   m_statEtaBins.end(),
                                                                   m_useAbsEtaForStat ? std::abs(eta) : eta);
    std::size_t category_i = seed_collection * seedCollectionStride() + static_cast<std::size_t>(bin_iter - m_statEtaBins.begin());
    assert(category_i < m_stat.size());
    return category_i;
  }
 
  inline std::size_t TrackFindingAlg::computeStatSum(std::size_t seed_collection, EStat counter_i, const EventStats &stat) const
  {
    std::size_t out = 0u;
    for (std::size_t category_i = seed_collection * seedCollectionStride() + static_cast<std::size_t>(counter_i);
         category_i < (seed_collection + 1) * seedCollectionStride();
         ++category_i)
    {
      assert(category_i < stat.size());
      out += stat[category_i][counter_i];
    }
    return out;
  }
 
 
  StatusCode TrackFindingAlg::initializeMeasurementSelector() {
    std::vector<std::pair<float, float> > chi2CutOffOutlier ;
    chi2CutOffOutlier .reserve( m_chi2CutOff.size() );
    if (!m_chi2OutlierCutOff.empty()) {
       if (m_chi2CutOff.size() !=  m_chi2OutlierCutOff.size()) {
          ATH_MSG_ERROR("Outlier chi2 cut off provided but number of elements does not agree with"
                        " chi2 cut off for measurements which however is required: "
                        << m_chi2CutOff.size() << " != " <<  m_chi2OutlierCutOff.size());
          return StatusCode::FAILURE;
       }
    }
    unsigned int idx=0;
    for (const auto &elm : m_chi2CutOff) {
       chi2CutOffOutlier.push_back( std::make_pair(static_cast<float>(elm),
                                                   idx < m_chi2OutlierCutOff.size()
                                                   ? static_cast<float>(m_chi2OutlierCutOff[idx])
                                                   : std::numeric_limits<float>::max()) );
       ++idx;
    }
    std::vector<float> etaBinsf;
    if (m_etaBins.size() > 2) {
      etaBinsf.assign(m_etaBins.begin() + 1, m_etaBins.end() - 1);
    }
 
    m_measurementSelector = ActsTrk::detail::getMeasurementSelector(m_pixelCalibTool.isEnabled() ? &(*m_pixelCalibTool) : nullptr,
                                    etaBinsf,
                                    chi2CutOffOutlier,
                                    m_numMeasurementsCutOff.value());
 
    return m_measurementSelector ? StatusCode::SUCCESS : StatusCode::FAILURE;
  }
 
} // namespace