From 4d2d6e0b5c37da202d4c6e82c341fd628d40ce2a Mon Sep 17 00:00:00 2001
From: Yizhou Zhang <zhangyz@ihep.ac.cn>
Date: Mon, 2 Dec 2024 14:10:09 +0000
Subject: [PATCH] add ACTS rec for sitrack
---
Reconstruction/CMakeLists.txt | 1 +
Reconstruction/RecActsTracking/CMakeLists.txt | 37 +
.../options/RecActsTracking.py | 79 ++
.../RecActsTracking/src/RecActsTracking.cpp | 928 ++++++++++++++++++
.../RecActsTracking/src/RecActsTracking.h | 288 ++++++
.../RecActsTracking/src/utils/CKFhelper.hpp | 388 ++++++++
.../src/utils/GeometryContainers.hpp | 355 +++++++
.../src/utils/MagneticField.hpp | 121 +++
.../RecActsTracking/src/utils/Options.cpp | 166 ++++
.../RecActsTracking/src/utils/Options.hpp | 159 +++
.../src/utils/SimSpacePoint.hpp | 243 +++++
.../src/utils/TGeoDetector.hpp | 596 +++++++++++
12 files changed, 3361 insertions(+)
create mode 100644 Reconstruction/RecActsTracking/CMakeLists.txt
create mode 100644 Reconstruction/RecActsTracking/options/RecActsTracking.py
create mode 100644 Reconstruction/RecActsTracking/src/RecActsTracking.cpp
create mode 100644 Reconstruction/RecActsTracking/src/RecActsTracking.h
create mode 100644 Reconstruction/RecActsTracking/src/utils/CKFhelper.hpp
create mode 100644 Reconstruction/RecActsTracking/src/utils/GeometryContainers.hpp
create mode 100644 Reconstruction/RecActsTracking/src/utils/MagneticField.hpp
create mode 100644 Reconstruction/RecActsTracking/src/utils/Options.cpp
create mode 100644 Reconstruction/RecActsTracking/src/utils/Options.hpp
create mode 100644 Reconstruction/RecActsTracking/src/utils/SimSpacePoint.hpp
create mode 100644 Reconstruction/RecActsTracking/src/utils/TGeoDetector.hpp
diff --git a/Reconstruction/CMakeLists.txt b/Reconstruction/CMakeLists.txt
index db8d7c2c..4ce55556 100644
--- a/Reconstruction/CMakeLists.txt
+++ b/Reconstruction/CMakeLists.txt
@@ -7,3 +7,4 @@ add_subdirectory(RecTrkGlobal)
add_subdirectory(RecDCHGenfit)
add_subdirectory(RecAssociationMaker)
add_subdirectory(ParticleID)
+add_subdirectory(RecActsTracking)
\ No newline at end of file
diff --git a/Reconstruction/RecActsTracking/CMakeLists.txt b/Reconstruction/RecActsTracking/CMakeLists.txt
new file mode 100644
index 00000000..d0c8c9d9
--- /dev/null
+++ b/Reconstruction/RecActsTracking/CMakeLists.txt
@@ -0,0 +1,37 @@
+find_package(Acts COMPONENTS
+ Core PluginFpeMonitoring PluginGeant4 PluginJson
+ PluginTGeo PluginDD4hep PluginEDM4hep Fatras)
+
+if(NOT Acts_FOUND)
+ message("Acts package not found. RecActsTracking module requires Acts.")
+ return()
+endif()
+
+gaudi_add_module(RecActsTracking
+ SOURCES
+ src/RecActsTracking.cpp
+ LINK DetInterface
+ k4FWCore::k4FWCore
+ Gaudi::GaudiAlgLib Gaudi::GaudiKernel
+ ${LCIO_LIBRARIES}
+ ${DD4hep_COMPONENT_LIBRARIES}
+ EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+ ${podio_LIBRARIES} podio::podioRootIO
+ GearSvc ${GEAR_LIBRARIES}
+ ActsCore ActsPluginFpeMonitoring ActsPluginGeant4
+ ActsPluginJson ActsPluginTGeo ActsPluginDD4hep
+ ActsPluginEDM4hep ActsFatras
+)
+
+target_include_directories(RecActsTracking PUBLIC $ENV{ACTS}/include)
+
+target_include_directories(RecActsTracking PUBLIC
+ ${LCIO_INCLUDE_DIRS}
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+ $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+install(TARGETS RecActsTracking
+ EXPORT CEPCSWTargets
+ RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+ LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+ COMPONENT dev)
diff --git a/Reconstruction/RecActsTracking/options/RecActsTracking.py b/Reconstruction/RecActsTracking/options/RecActsTracking.py
new file mode 100644
index 00000000..c570bad2
--- /dev/null
+++ b/Reconstruction/RecActsTracking/options/RecActsTracking.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+import os
+
+from Gaudi.Configuration import *
+
+NTupleSvc().Output = ["MyTuples DATAFILE='test.root' OPT='NEW' TYP='ROOT'"]
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc", input="rec00.root")
+
+from Configurables import PodioInput
+podioinput = PodioInput("PodioReader", collections=[
+ "MCParticle",
+ "VXDCollection",
+ "SITCollection",
+ "VXDTrackerHits",
+ "SITTrackerHits",
+])
+
+##############################################################################
+# Geometry Svc
+##############################################################################
+
+# geometry_option = "TDR_o1_v01/TDR_o1_v01-onlyTracker.xml"
+geometry_option = "TDR_o1_v01/TDR_o1_v01.xml"
+
+if not os.getenv("DETCRDROOT"):
+ print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import DetGeomSvc
+geosvc = DetGeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+from Configurables import GearSvc
+gearsvc = GearSvc("GearSvc")
+
+from Configurables import TrackSystemSvc
+tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
+
+##############################################################################
+# CEPCSWData
+##############################################################################
+cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
+
+##############################################################################
+# RecActsTracking
+##############################################################################
+
+from Configurables import RecActsTracking
+actstracking = RecActsTracking("RecActsTracking")
+actstracking.TGeoFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr24.12.0/SiTrack-tgeo.root")
+actstracking.TGeoConfigFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr24.12.0/SiTrack-tgeo-config.json")
+actstracking.MaterialMapFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr24.12.0/SiTrack-material-maps.json")
+
+##############################################################################
+# output
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("out")
+out.filename = "ACTS_rec00.root"
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [podioinput, actstracking, out],
+ EvtSel = 'NONE',
+ EvtMax = 5,
+ ExtSvc = [dsvc, geosvc, evtseeder, gearsvc],
+ # OutputLevel=DEBUG
+)
diff --git a/Reconstruction/RecActsTracking/src/RecActsTracking.cpp b/Reconstruction/RecActsTracking/src/RecActsTracking.cpp
new file mode 100644
index 00000000..cd228ee7
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/RecActsTracking.cpp
@@ -0,0 +1,928 @@
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+#include <sstream>
+
+// dependence
+#include "RecActsTracking.h"
+
+#include "GearSvc/IGearSvc.h"
+
+// csv parser
+// #include "csv2/writer.hpp"
+
+// MC
+#include "CLHEP/Units/SystemOfUnits.h"
+
+
+using namespace Acts::UnitLiterals;
+
+DECLARE_COMPONENT(RecActsTracking)
+
+RecActsTracking::RecActsTracking(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc)
+{
+ // input collections
+ // declareProperty("VXDTrackerHits", _inVTXTrackHdl, "Handle of the Input VTX TrackerHits collection");
+ // declareProperty("SITTrackerHits", _inSITTrackHdl, "Handle of the Input SIT TrackerHits collection");
+
+ // declareProperty("VXDCollection", _inVTXColHdl, "Handle of the Input VTX Sim Hit collection");
+ // declareProperty("SITCollection", _inSITColHdl, "Handle of the Input SIT Sim Hit collection");
+
+ // declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
+
+ // Output Collections
+ // declareProperty("OutputTracks", _outColHdl, "Handle of the ACTS SiTrack output collection");
+}
+
+StatusCode RecActsTracking::initialize()
+{
+ chronoStatSvc = service<IChronoStatSvc>("ChronoStatSvc");
+ _nEvt = 0;
+
+ if (!std::filesystem::exists(TGeo_path.value())) {
+ error() << "CEPC TGeo file: " << TGeo_path.value() << " does not exist!" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ if (!std::filesystem::exists(TGeo_config_path.value())) {
+ error() << "CEPC TGeo config file: " << TGeo_config_path.value() << " does not exist!" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ if (!std::filesystem::exists(MaterialMap_path.value())) {
+ error() << "CEPC Material map file: " << MaterialMap_path.value() << " does not exist!" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ TGeo_ROOTFilePath = TGeo_path.value();
+ TGeoConfig_jFilePath = TGeo_config_path.value();
+ MaterialMap_jFilePath = MaterialMap_path.value();
+
+ chronoStatSvc->chronoStart("read geometry");
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if (!m_geosvc) {
+ error() << "Failed to find GeomSvc." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ if(m_geosvc){
+ const dd4hep::Direction& field = m_geosvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
+ m_field = field.z()/dd4hep::tesla;
+ }
+
+ m_vtx_surfaces = m_geosvc->getSurfaceMap("VXD");
+ debug() << "Surface map size: " << m_vtx_surfaces->size() << endmsg;
+
+ m_sit_surfaces = m_geosvc->getSurfaceMap("SIT");
+ debug() << "Surface map size: " << m_sit_surfaces->size() << endmsg;
+
+ vxd_decoder = m_geosvc->getDecoder("VXDCollection");
+ if(!vxd_decoder){
+ return StatusCode::FAILURE;
+ }
+
+ sit_decoder = m_geosvc->getDecoder("SITCollection");
+ if(!sit_decoder){
+ return StatusCode::FAILURE;
+ }
+
+ info() << "ACTS Tracking Geometry initialized successfully!" << endmsg;
+ // initialize tgeo detector
+ auto logger = Acts::getDefaultLogger("TGeoDetector", Acts::Logging::INFO);
+ trackingGeometry = buildTGeoDetector(geoContext, detElementStore, TGeo_ROOTFilePath, TGeoConfig_jFilePath, MaterialMap_jFilePath, *logger);
+
+ info() << "Seeding tools initialized successfully!" << endmsg;
+
+ // configure the acts tools
+ seed_cfg.seedFinderOptions.bFieldInZ = 3_T;
+ seed_cfg.seedFinderConfig.deltaRMin = 4_mm;
+ seed_cfg.seedFinderConfig.deltaRMax = 25_mm;
+ seed_cfg.seedFinderConfig.rMax = 40_mm;
+ seed_cfg.seedFinderConfig.rMin = 8_mm;
+ seed_cfg.seedFinderConfig.impactMax = 4_mm;
+ seed_cfg.seedFinderConfig.useVariableMiddleSPRange = false;
+ seed_cfg.seedFinderConfig.rMinMiddle = 15_mm; // range for middle spacepoint (TDR)
+ seed_cfg.seedFinderConfig.rMaxMiddle = 30_mm; // range for middle spacepoint (TDR)
+
+ // initialize the acts tools
+ seed_cfg.seedFilterConfig = seed_cfg.seedFilterConfig.toInternalUnits();
+ seed_cfg.seedFinderConfig.seedFilter =
+ std::make_unique<Acts::SeedFilter<SimSpacePoint>>(seed_cfg.seedFilterConfig);
+ seed_cfg.seedFinderConfig =
+ seed_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
+ seed_cfg.seedFinderOptions =
+ seed_cfg.seedFinderOptions.toInternalUnits().calculateDerivedQuantities(seed_cfg.seedFinderConfig);
+ seed_cfg.gridConfig = seed_cfg.gridConfig.toInternalUnits();
+ seed_cfg.gridOptions = seed_cfg.gridOptions.toInternalUnits();
+
+ if (std::isnan(seed_cfg.seedFinderConfig.deltaRMaxTopSP)) {
+ seed_cfg.seedFinderConfig.deltaRMaxTopSP = seed_cfg.seedFinderConfig.deltaRMax;}
+ if (std::isnan(seed_cfg.seedFinderConfig.deltaRMinTopSP)) {
+ seed_cfg.seedFinderConfig.deltaRMinTopSP = seed_cfg.seedFinderConfig.deltaRMin;}
+ if (std::isnan(seed_cfg.seedFinderConfig.deltaRMaxBottomSP)) {
+ seed_cfg.seedFinderConfig.deltaRMaxBottomSP = seed_cfg.seedFinderConfig.deltaRMax;}
+ if (std::isnan(seed_cfg.seedFinderConfig.deltaRMinBottomSP)) {
+ seed_cfg.seedFinderConfig.deltaRMinBottomSP = seed_cfg.seedFinderConfig.deltaRMin;}
+
+ m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
+ seed_cfg.numPhiNeighbors, seed_cfg.zBinNeighborsBottom);
+ m_topBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
+ seed_cfg.numPhiNeighbors, seed_cfg.zBinNeighborsTop);
+
+ seed_cfg.seedFinderConfig.seedFilter =
+ std::make_unique<Acts::SeedFilter<SimSpacePoint>>(seed_cfg.seedFilterConfig);
+ m_seedFinder =
+ Acts::SeedFinder<SimSpacePoint, Acts::CylindricalSpacePointGrid<SimSpacePoint>>(seed_cfg.seedFinderConfig);
+
+ // initialize the ckf
+ findTracks = makeTrackFinderFunction(trackingGeometry, magneticField);
+
+ info() << "CKF Track Finder initialized successfully!" << endmsg;
+ chronoStatSvc->chronoStop("read geometry");
+
+ return GaudiAlgorithm::initialize();
+}
+
+StatusCode RecActsTracking::execute()
+{
+ auto trkCol = _outColHdl.createAndPut();
+
+ bool MCsuccess = true;
+ const edm4hep::MCParticleCollection* mcCols = nullptr;
+ mcCols = _inMCColHdl.get();
+ if(mcCols)
+ {
+ for(auto particle : *mcCols)
+ {
+ if (particle.getPDG() != 13)
+ {
+ MCsuccess = false;
+ continue;
+ }
+ break;
+ }
+ }
+
+ if (!MCsuccess)
+ {
+ _nEvt++;
+ return StatusCode::SUCCESS;
+ }
+
+ SpacePointPtrs.clear();
+ sourceLinks.clear();
+ measurements.clear();
+ initialParameters.clear();
+ Selected_Seeds.clear();
+
+ chronoStatSvc->chronoStart("read input hits");
+ int successVTX = InitialiseVTX();
+
+ if (successVTX == 0)
+ {
+ _nEvt++;
+ return StatusCode::SUCCESS;
+ }
+
+ int successSIT = InitialiseSIT();
+ if (successSIT == 0)
+ {
+ _nEvt++;
+ return StatusCode::SUCCESS;
+ }
+
+ chronoStatSvc->chronoStop("read input hits");
+ // info() << "Generated " << SpacePointPtrs.size() << " spacepoints for event " << _nEvt << "!" << endmsg;
+ // info() << "Generated " << measurements.size() << " measurements for event " << _nEvt << "!" << endmsg;
+
+ // --------------------------------------------
+ // Seeding
+ // --------------------------------------------
+ chronoStatSvc->chronoStart("seeding");
+
+ // construct the seeding tools
+ // covariance tool, extracts covariances per spacepoint as required
+ auto extractGlobalQuantities = [=](const SimSpacePoint& sp, float, float, float)
+ {
+ Acts::Vector3 position{sp.x(), sp.y(), sp.z()};
+ Acts::Vector2 covariance{sp.varianceR(), sp.varianceZ()};
+ return std::make_tuple(position, covariance, sp.t());
+ };
+
+ // extent used to store r range for middle spacepoint
+ Acts::Extent rRangeSPExtent;
+
+ // construct the seeding tool
+ Acts::CylindricalSpacePointGrid<SimSpacePoint> grid =
+ Acts::CylindricalSpacePointGridCreator::createGrid<SimSpacePoint>(seed_cfg.gridConfig, seed_cfg.gridOptions);
+ Acts::CylindricalSpacePointGridCreator::fillGrid(
+ seed_cfg.seedFinderConfig, seed_cfg.seedFinderOptions, grid,
+ SpacePointPtrs.begin(), SpacePointPtrs.end(), extractGlobalQuantities, rRangeSPExtent);
+
+ std::array<std::vector<std::size_t>, 2ul> navigation;
+ navigation[1ul] = seed_cfg.seedFinderConfig.zBinsCustomLooping;
+
+ auto spacePointsGrouping = Acts::CylindricalBinnedGroup<SimSpacePoint>(
+ std::move(grid), *m_bottomBinFinder, *m_topBinFinder, std::move(navigation));
+
+ // safely clamp double to float
+ float up = Acts::clampValue<float>(
+ std::floor(rRangeSPExtent.max(Acts::binR) / 2) * 2);
+
+ /// variable middle SP radial region of interest
+ const Acts::Range1D<float> rMiddleSPRange(
+ std::floor(rRangeSPExtent.min(Acts::binR) / 2) * 2 +
+ seed_cfg.seedFinderConfig.deltaRMiddleMinSPRange,
+ up - seed_cfg.seedFinderConfig.deltaRMiddleMaxSPRange);
+
+ // run the seeding
+ static thread_local SimSeedContainer seeds;
+ seeds.clear();
+ static thread_local decltype(m_seedFinder)::SeedingState state;
+ state.spacePointData.resize(
+ SpacePointPtrs.size(),
+ seed_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo);
+
+ // use double stripe measurement
+ if (seed_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo)
+ {
+ for (std::size_t grid_glob_bin(0); grid_glob_bin < spacePointsGrouping.grid().size(); ++grid_glob_bin)
+ {
+ const auto& collection = spacePointsGrouping.grid().at(grid_glob_bin);
+ for (const auto& sp : collection)
+ {
+ std::size_t index = sp->index();
+ const float topHalfStripLength =
+ seed_cfg.seedFinderConfig.getTopHalfStripLength(sp->sp());
+ const float bottomHalfStripLength =
+ seed_cfg.seedFinderConfig.getBottomHalfStripLength(sp->sp());
+ const Acts::Vector3 topStripDirection =
+ seed_cfg.seedFinderConfig.getTopStripDirection(sp->sp());
+ const Acts::Vector3 bottomStripDirection =
+ seed_cfg.seedFinderConfig.getBottomStripDirection(sp->sp());
+
+ state.spacePointData.setTopStripVector(
+ index, topHalfStripLength * topStripDirection);
+ state.spacePointData.setBottomStripVector(
+ index, bottomHalfStripLength * bottomStripDirection);
+ state.spacePointData.setStripCenterDistance(
+ index, seed_cfg.seedFinderConfig.getStripCenterDistance(sp->sp()));
+ state.spacePointData.setTopStripCenterPosition(
+ index, seed_cfg.seedFinderConfig.getTopStripCenterPosition(sp->sp()));
+ }
+ }
+ }
+
+ for (const auto [bottom, middle, top] : spacePointsGrouping)
+ {
+ m_seedFinder.createSeedsForGroup(
+ seed_cfg.seedFinderOptions, state, spacePointsGrouping.grid(),
+ std::back_inserter(seeds), bottom, middle, top, rMiddleSPRange);
+ }
+
+ // int seed_counter = 0;
+ // for (const auto& seed : seeds)
+ // {
+ // for (const auto& sp : seed.sp())
+ // {
+ // info() << "found seed #" << seed_counter << ": x:" << sp->x() << " y: " << sp->y() << " z: " << sp->z() << endmsg;
+ // }
+ // seed_counter++;
+ // }
+
+ chronoStatSvc->chronoStop("seeding");
+ debug() << "Found " << seeds.size() << " seeds for event " << _nEvt << "!" << endmsg;
+
+ // --------------------------------------------
+ // track estimation
+ // --------------------------------------------
+ chronoStatSvc->chronoStart("track_param");
+
+ IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
+
+ for (std::size_t iseed = 0; iseed < seeds.size(); ++iseed)
+ {
+ const auto& seed = seeds[iseed];
+ // Get the bottom space point and its reference surface
+ const auto bottomSP = seed.sp().front();
+ const auto& sourceLink = bottomSP->sourceLinks()[0];
+ const Acts::Surface* surface = surfaceAccessor(sourceLink);
+ if (surface == nullptr) {
+ debug() << "Surface from source link is not found in the tracking geometry: iseed " << iseed << endmsg;
+ continue;
+ }
+
+ auto optParams = Acts::estimateTrackParamsFromSeed(
+ geoContext, seed.sp().begin(), seed.sp().end(), *surface, acts_field_value, bFieldMin);
+
+ if (!optParams.has_value()) {
+ debug() << "Estimation of track parameters for seed " << iseed << " failed." << endmsg;
+ continue;
+ }
+
+ const auto& params = optParams.value();
+ Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
+ for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
+ double sigma = initialSigmas[i];
+ sigma += initialSimgaQoverPCoefficients[i] * params[Acts::eBoundQOverP];
+ double var = sigma * sigma;
+ if (i == Acts::eBoundTime && !bottomSP->t().has_value()) { var *= noTimeVarInflation; }
+ var *= initialVarInflation[i];
+
+ cov(i, i) = var;
+ }
+ initialParameters.emplace_back(surface->getSharedPtr(), params, cov, particleHypothesis);
+ Selected_Seeds.push_back(seed);
+ }
+
+ chronoStatSvc->chronoStop("track_param");
+ debug() << "Found " << initialParameters.size() << " tracks for event " << _nEvt << "!" << endmsg;
+
+ // --------------------------------------------
+ // CKF track finding
+ // --------------------------------------------
+ chronoStatSvc->chronoStart("ckf_findTracks");
+
+ // Construct a perigee surface as the target surface
+ auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
+ PassThroughCalibrator pcalibrator;
+ MeasurementCalibratorAdapter calibrator(pcalibrator, measurements);
+ Acts::GainMatrixUpdater kfUpdater;
+ Acts::MeasurementSelector::Config measurementSelectorCfg =
+ {
+ {Acts::GeometryIdentifier(), {{}, {30}, {1u}}},
+ };
+
+ MeasurementSelector measSel{ Acts::MeasurementSelector(measurementSelectorCfg) };
+ using Extensions = Acts::CombinatorialKalmanFilterExtensions<Acts::VectorMultiTrajectory>;
+ BranchStopper branchStopper(trackSelectorCfg);
+
+ // Construct the CKF
+ Extensions extensions;
+ extensions.calibrator.connect<&MeasurementCalibratorAdapter::calibrate>(&calibrator);
+ extensions.updater.connect<&Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(&kfUpdater);
+ extensions.measurementSelector.connect<&MeasurementSelector::select>(&measSel);
+ extensions.branchStopper.connect<&BranchStopper::operator()>(&branchStopper);
+
+ IndexSourceLinkAccessor slAccessor;
+ slAccessor.container = &sourceLinks;
+ Acts::SourceLinkAccessorDelegate<IndexSourceLinkAccessor::Iterator> slAccessorDelegate;
+ slAccessorDelegate.connect<&IndexSourceLinkAccessor::range>(&slAccessor);
+
+ Acts::PropagatorPlainOptions firstPropOptions;
+ firstPropOptions.maxSteps = maxSteps;
+ firstPropOptions.direction = Acts::Direction::Forward;
+
+ Acts::PropagatorPlainOptions secondPropOptions;
+ secondPropOptions.maxSteps = maxSteps;
+ secondPropOptions.direction = firstPropOptions.direction.invert();
+
+ // Set the CombinatorialKalmanFilter options
+ TrackFinderOptions firstOptions(
+ geoContext, magFieldContext, calibContext,
+ slAccessorDelegate, extensions, firstPropOptions);
+
+ TrackFinderOptions secondOptions(
+ geoContext, magFieldContext, calibContext,
+ slAccessorDelegate, extensions, secondPropOptions);
+ secondOptions.targetSurface = pSurface.get();
+
+ Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator> extrapolator(
+ Acts::EigenStepper<>(magneticField), Acts::Navigator({trackingGeometry}));
+
+ Acts::PropagatorOptions<Acts::ActionList<Acts::MaterialInteractor>, Acts::AbortList<Acts::EndOfWorldReached>>
+ extrapolationOptions(geoContext, magFieldContext);
+
+ auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
+ auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
+ auto trackContainerTemp = std::make_shared<Acts::VectorTrackContainer>();
+ auto trackStateContainerTemp = std::make_shared<Acts::VectorMultiTrajectory>();
+ TrackContainer tracks(trackContainer, trackStateContainer);
+ TrackContainer tracksTemp(trackContainerTemp, trackStateContainerTemp);
+
+ tracks.addColumn<unsigned int>("trackGroup");
+ tracksTemp.addColumn<unsigned int>("trackGroup");
+ Acts::ProxyAccessor<unsigned int> seedNumber("trackGroup");
+
+ unsigned int nSeed = 0;
+ // A map indicating whether a seed has been discovered already
+ std::unordered_map<SeedIdentifier, bool> discoveredSeeds;
+ auto addTrack = [&](const TrackProxy& track)
+ {
+ ++m_nFoundTracks;
+ // flag seeds which are covered by the track
+ visitSeedIdentifiers(track, [&](const SeedIdentifier& seedIdentifier)
+ {
+ if (auto it = discoveredSeeds.find(seedIdentifier); it != discoveredSeeds.end())
+ {
+ it->second = true;
+ }
+ });
+
+ if (m_trackSelector.has_value() && !m_trackSelector->isValidTrack(track)) { return; }
+
+ ++m_nSelectedTracks;
+ auto destProxy = tracks.makeTrack();
+ // make sure we copy track states!
+ destProxy.copyFrom(track, true);
+ };
+
+ for (const auto& seed : Selected_Seeds) {
+ SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
+ discoveredSeeds.emplace(seedIdentifier, false);
+ }
+
+ for (std::size_t iSeed = 0; iSeed < initialParameters.size(); ++iSeed)
+ {
+ m_nTotalSeeds++;
+ const auto& seed = Selected_Seeds[iSeed];
+
+ SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
+ // check if the seed has been discovered already
+ if (auto it = discoveredSeeds.find(seedIdentifier); it != discoveredSeeds.end() && it->second)
+ {
+ m_nDeduplicatedSeeds++;
+ continue;
+ }
+
+ /// Whether to stick on the seed measurements during track finding.
+ // measSel.setSeed(seed);
+
+ // Clear trackContainerTemp and trackStateContainerTemp
+ tracksTemp.clear();
+
+ const Acts::BoundTrackParameters& firstInitialParameters = initialParameters.at(iSeed);
+ auto firstResult = (*findTracks)(firstInitialParameters, firstOptions, tracksTemp);
+
+ nSeed++;
+ if (!firstResult.ok())
+ {
+ m_nFailedSeeds++;
+ continue;
+ }
+
+ auto& firstTracksForSeed = firstResult.value();
+ for (auto& firstTrack : firstTracksForSeed)
+ {
+ auto trackCandidate = tracksTemp.makeTrack();
+ trackCandidate.copyFrom(firstTrack, true);
+
+ auto firstSmoothingResult = Acts::smoothTrack(geoContext, trackCandidate);
+
+ if (!firstSmoothingResult.ok())
+ {
+ m_nFailedSmoothing++;
+ continue;
+ }
+
+ seedNumber(trackCandidate) = nSeed - 1;
+
+ auto firstExtrapolationResult = Acts::extrapolateTrackToReferenceSurface(
+ trackCandidate, *pSurface, extrapolator, extrapolationOptions, extrapolationStrategy);
+
+ if (!firstExtrapolationResult.ok())
+ {
+ m_nFailedExtrapolation++;
+ continue;
+ }
+
+ addTrack(trackCandidate);
+ }
+ }
+
+ chronoStatSvc->chronoStop("ckf_findTracks");
+ debug() << "CKF found " << tracks.size() << " tracks for event " << _nEvt << "!" << endmsg;
+
+ m_memoryStatistics.local().hist += tracks.trackStateContainer().statistics().hist;
+ auto constTrackStateContainer = std::make_shared<Acts::ConstVectorMultiTrajectory>(std::move(*trackStateContainer));
+ auto constTrackContainer = std::make_shared<Acts::ConstVectorTrackContainer>(std::move(*trackContainer));
+ ConstTrackContainer constTracks{constTrackContainer, constTrackStateContainer};
+
+ chronoStatSvc->chronoStart("writeout tracks");
+
+ if (constTracks.size() == 0) {
+ chronoStatSvc->chronoStop("writeout tracks");
+ _nEvt++;
+ return StatusCode::SUCCESS;
+ }
+
+ for (const auto& cur_track : constTracks)
+ {
+ auto writeout_track = trkCol->create();
+ int nVTXHit = 0, nFTDHit = 0, nSITHit = 0;
+ int getFirstHit = 0;
+
+ writeout_track.setChi2(cur_track.chi2());
+ writeout_track.setNdf(cur_track.nDoF());
+ writeout_track.setDEdx(cur_track.absoluteMomentum() / Acts::UnitConstants::GeV);
+ writeout_track.setDEdxError(cur_track.qOverP());
+ for (auto trackState : cur_track.trackStates()){
+ if (trackState.hasUncalibratedSourceLink()){
+ auto cur_measurement_sl = trackState.getUncalibratedSourceLink();
+ const auto& MeasSourceLink = cur_measurement_sl.get<IndexSourceLink>();
+ auto cur_measurement = measurements[MeasSourceLink.index()];
+ auto cur_measurement_gid = MeasSourceLink.geometryId();
+ if (std::find(VXD_inner_volume_ids.begin(), VXD_inner_volume_ids.end(), cur_measurement_gid.volume()) != VXD_inner_volume_ids.end()){
+ nVTXHit++;
+ } else if (std::find(SIT_acts_volume_ids.begin(), SIT_acts_volume_ids.end(), cur_measurement_gid.volume()) != SIT_acts_volume_ids.end()){
+ nSITHit++;
+ }
+
+ writeout_track.addToTrackerHits(MeasSourceLink.getTrackerHit());
+
+ if (!getFirstHit){
+ const auto& par = std::get<1>(cur_measurement).parameters();
+ const Acts::Surface* surface = surfaceAccessor(cur_measurement_sl);
+ auto acts_global_postion = surface->localToGlobal(geoContext, par, globalFakeMom);
+ writeout_track.setRadiusOfInnermostHit(
+ std::sqrt(acts_global_postion[0] * acts_global_postion[0] +
+ acts_global_postion[1] * acts_global_postion[1] +
+ acts_global_postion[2] * acts_global_postion[2])
+ );
+ getFirstHit = 1;
+ }
+ }
+ }
+
+ // SubdetectorHitNumbers: VXD = 0, FTD = 1, SIT = 2
+ writeout_track.addToSubdetectorHitNumbers(nVTXHit);
+ writeout_track.addToSubdetectorHitNumbers(nFTDHit);
+ writeout_track.addToSubdetectorHitNumbers(nSITHit);
+
+ std::array<float, 21> writeout_covMatrix;
+ auto cur_track_covariance = cur_track.covariance();
+ for (int i = 0; i < 6; i++) {
+ for (int j = 0; j < 6-i; j++) {
+ writeout_covMatrix[i * 6 + j] = cur_track_covariance(writeout_indices[i], writeout_indices[j]);
+ }
+ }
+ // location: At{Other|IP|FirstHit|LastHit|Calorimeter|Vertex}|LastLocation
+ // TrackState: location, d0, phi, omega, z0, tanLambda, time, referencePoint, covMatrix
+ edm4hep::TrackState writeout_trackState{
+ 0, // location: AtOther
+ cur_track.loc0() / Acts::UnitConstants::mm, // d0
+ cur_track.phi(), // phi
+ cur_track.qOverP() * sin(cur_track.theta()) * _FCT * m_field, // omega = qop * sin(theta) * _FCT * bf
+ cur_track.loc1() / Acts::UnitConstants::mm, // z0
+ 1 / tan(cur_track.theta()), // tanLambda = 1 / tan(theta)
+ cur_track.time(), // time
+ ::edm4hep::Vector3f(0, 0, 0), // referencePoint
+ writeout_covMatrix
+ };
+
+ debug() << "Found track with momentum " << cur_track.absoluteMomentum() / Acts::UnitConstants::GeV << " !" << endmsg;
+ writeout_track.addToTrackStates(writeout_trackState);
+ }
+
+ chronoStatSvc->chronoStop("writeout tracks");
+
+ _nEvt++;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode RecActsTracking::finalize()
+{
+ debug() << "finalize RecActsTracking" << endmsg;
+ info() << "Total number of events processed: " << _nEvt << endmsg;
+ info() << "Total number of **TotalSeeds** processed: " << m_nTotalSeeds << endmsg;
+ info() << "Total number of **FoundTracks** processed: " << m_nFoundTracks << endmsg;
+ info() << "Total number of **SelectedTracks** processed: " << m_nSelectedTracks << endmsg;
+
+ return GaudiAlgorithm::finalize();
+}
+
+int RecActsTracking::InitialiseVTX()
+{
+ int success = 1;
+
+ const edm4hep::TrackerHitCollection* hitVTXCol = nullptr;
+ const edm4hep::SimTrackerHitCollection* SimhitVTXCol = nullptr;
+
+ try {
+ hitVTXCol = _inVTXTrackHdl.get();
+ } catch (GaudiException& e) {
+ debug() << "Collection " << _inVTXTrackHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ success = 0;
+ }
+
+ try {
+ SimhitVTXCol = _inVTXColHdl.get();
+ } catch (GaudiException& e) {
+ debug() << "Sim Collection " << _inVTXColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ success = 0;
+ }
+
+ if(hitVTXCol && SimhitVTXCol)
+ {
+ int nelem = hitVTXCol->size();
+ debug() << "Number of VTX hits = " << nelem << endmsg;
+ if ((nelem < 3) or (nelem > 10)) { success = 0; }
+
+ // std::string truth_file = "obj/vtx/truth/event" + std::to_string(_nEvt) + ".csv";
+ // std::ofstream truth_stream(truth_file);
+ // csv2::Writer<csv2::delimiter<','>> truth_writer(truth_stream);
+ // std::vector<std::string> truth_header = {"layer", "x", "y", "z"};
+ // truth_writer.write_row(truth_header);
+
+ // std::string converted_file = "obj/vtx/converted/event" + std::to_string(_nEvt) + ".csv";
+ // std::ofstream converted_stream(converted_file);
+ // csv2::Writer<csv2::delimiter<','>> converted_writer(converted_stream);
+ // std::vector<std::string> converted_header = {"layer", "x", "y", "z"};
+ // converted_writer.write_row(converted_header);
+
+ for (int ielem = 0; ielem < nelem; ++ielem)
+ {
+ auto hit = hitVTXCol->at(ielem);
+ auto simhit = SimhitVTXCol->at(ielem);
+
+ auto simcellid = simhit.getCellID();
+ // system:5,side:-2,layer:9,module:8,sensor:32:8
+ uint64_t m_layer = vxd_decoder->get(simcellid, "layer");
+ uint64_t m_module = vxd_decoder->get(simcellid, "module");
+ uint64_t m_sensor = vxd_decoder->get(simcellid, "sensor");
+ double acts_x = simhit.getPosition()[0];
+ double acts_y = simhit.getPosition()[1];
+ double acts_z = simhit.getPosition()[2];
+
+ double momentum_x = simhit.getMomentum()[0];
+ double momentum_y = simhit.getMomentum()[1];
+ double momentum_z = simhit.getMomentum()[2];
+
+ dd4hep::rec::ISurface* surface = nullptr;
+ auto it = m_vtx_surfaces->find(simcellid);
+ if (it != m_vtx_surfaces->end()) {
+ surface = it->second;
+ if (!surface) {
+ fatal() << "found surface for VTX cell id " << simcellid << ", but NULL" << endmsg;
+ return 0;
+ }
+ }
+ else {
+ fatal() << "not found surface for VTX cell id " << simcellid << endmsg;
+ return 0;
+ }
+
+ // dd4hep::rec::Vector3D oldPos(simhit.getPosition()[0]*dd4hep::mm/CLHEP::mm, simhit.getPosition()[1]*dd4hep::mm/CLHEP::mm, simhit.getPosition()[2]*dd4hep::mm/CLHEP::mm);
+ // dd4hep::rec::Vector2D localPoint = surface->globalToLocal(oldPos);
+
+ // if (m_layer < current_layer){
+ // info() << "ring hits happend in layer " << m_layer << " before layer " << current_layer << ", at event " << _nEvt << endmsg;
+ // success = 0;
+ // break;
+ // }
+ // current_layer = m_layer;
+
+ if (m_layer <= 3){
+ // set acts geometry identifier
+ // VXD_inner_volume_ids{16, 17, 18, 19};
+ uint64_t acts_volume = VXD_inner_volume_ids[m_layer];
+ uint64_t acts_boundary = 0;
+ uint64_t acts_layer = 2;
+ uint64_t acts_approach = 0;
+ uint64_t acts_sensitive = m_module + 1;
+
+ Acts::GeometryIdentifier moduleGeoId;
+ moduleGeoId.setVolume(acts_volume);
+ moduleGeoId.setBoundary(acts_boundary);
+ moduleGeoId.setLayer(acts_layer);
+ moduleGeoId.setApproach(acts_approach);
+ moduleGeoId.setSensitive(acts_sensitive);
+
+ // create and store the source link
+ uint32_t measurementIdx = measurements.size();
+ IndexSourceLink sourceLink{moduleGeoId, measurementIdx, hit};
+ sourceLinks.insert(sourceLinks.end(), sourceLink);
+ Acts::SourceLink sl{sourceLink};
+ boost::container::static_vector<Acts::SourceLink, 2> slinks;
+ slinks.emplace_back(sl);
+
+ // get the surface of the hit
+ IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
+ const Acts::Surface* acts_surface = surfaceAccessor(sl);
+
+ // get the local position of the hit
+ const Acts::Vector3 globalPos{acts_x, acts_y, acts_z};
+ const Acts::Vector3 globalmom{momentum_x, momentum_y, momentum_z};
+ auto acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, onSurfaceTolerance);
+ if (!acts_local_postion.ok()){
+ info() << "Error: failed to get local position for VTX hit " << simcellid << endmsg;
+ acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, 100*onSurfaceTolerance);
+ }
+ const std::array<Acts::BoundIndices, 2> indices{Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc1};
+ const Acts::Vector2 par{acts_local_postion.value()[0], acts_local_postion.value()[1]};
+
+ // *** debug ***
+ debug() << "VXD measurements global position(x,y,z): " << simhit.getPosition()[0] << ", " << simhit.getPosition()[1] << ", " << simhit.getPosition()[2]
+ << "; local position(loc0, loc1): "<< acts_local_postion.value()[0] << ", " << acts_local_postion.value()[1] << endmsg;
+ auto acts_global_postion = acts_surface->localToGlobal(geoContext, par, globalFakeMom);
+ debug() << "debug surface at: x:" << acts_global_postion[0] << ", y:" << acts_global_postion[1] << ", z:" << acts_global_postion[2] << endmsg;
+
+ // SimSpacePoint *hitExt = new SimSpacePoint(hit, simhit, slinks);
+ SimSpacePoint *hitExt = new SimSpacePoint(hit, acts_global_postion[0], acts_global_postion[1], acts_global_postion[2], 0.002, slinks);
+ // debug() << "debug hitExt at: x:" << hitExt->x() << ", y:" << hitExt->y() << ", z:" << hitExt->z() << endmsg;
+ SpacePointPtrs.push_back(hitExt);
+
+ // create and store the measurement
+ Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Identity();
+ cov(0, 0) = std::max<double>(std::abs(acts_global_postion[0] - simhit.getPosition()[0]), eps);
+ cov(1, 1) = std::max<double>(std::abs(acts_global_postion[1] - simhit.getPosition()[1]), eps);
+ measurements.emplace_back(Acts::Measurement<Acts::BoundIndices, 2>(sl, indices, par, cov));
+
+ // std::vector<std::string> truth_col = {std::to_string(m_layer*2 + m_module), std::to_string(simhit.getPosition()[0]), std::to_string(simhit.getPosition()[1]), std::to_string(simhit.getPosition()[2])};
+ // truth_writer.write_row(truth_col);
+ // std::vector<std::string> converted_col = {std::to_string(m_layer*2 + m_module), std::to_string(acts_global_postion[0]), std::to_string(acts_global_postion[1]), std::to_string(acts_global_postion[2])};
+ // converted_writer.write_row(converted_col);
+
+ } else {
+ // set acts geometry identifier
+ // VXD_outer_volume_id = 20;
+ uint64_t acts_volume = VXD_outer_volume_id;
+ uint64_t acts_boundary = 0;
+ uint64_t acts_layer = 2;
+ uint64_t acts_approach = 0;
+ uint64_t acts_sensitive = (m_layer == 5) ? m_module*2 + 1 : m_module*2 + 2;
+
+ Acts::GeometryIdentifier moduleGeoId;
+ moduleGeoId.setVolume(acts_volume);
+ moduleGeoId.setBoundary(acts_boundary);
+ moduleGeoId.setLayer(acts_layer);
+ moduleGeoId.setApproach(acts_approach);
+ moduleGeoId.setSensitive(acts_sensitive);
+
+ // create and store the source link
+ uint32_t measurementIdx = measurements.size();
+ IndexSourceLink sourceLink{moduleGeoId, measurementIdx, hit};
+ sourceLinks.insert(sourceLinks.end(), sourceLink);
+ Acts::SourceLink sl{sourceLink};
+
+ // get the local position of the hit
+ IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
+ const Acts::Surface* acts_surface = surfaceAccessor(sl);
+ const Acts::Vector3 globalPos{acts_x, acts_y, acts_z};
+ const Acts::Vector3 globalmom{momentum_x, momentum_y, momentum_z};
+ auto acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, onSurfaceTolerance);
+ if (!acts_local_postion.ok()){
+ info() << "Error: failed to get local position for VTX hit " << simcellid << endmsg;
+ acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, 100*onSurfaceTolerance);
+ }
+ const std::array<Acts::BoundIndices, 2> indices{Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc1};
+ const Acts::Vector2 par{acts_local_postion.value()[0], acts_local_postion.value()[1]};
+
+ // *** debug ***
+ debug() << "VXD measurements global position(x,y,z): " << simhit.getPosition()[0] << ", " << simhit.getPosition()[1] << ", " << simhit.getPosition()[2]
+ << "; local position(loc0, loc1): "<< acts_local_postion.value()[0] << ", " << acts_local_postion.value()[1] << endmsg;
+ auto acts_global_postion = acts_surface->localToGlobal(geoContext, par, globalFakeMom);
+ debug() << "debug surface at: x:" << acts_global_postion[0] << ", y:" << acts_global_postion[1] << ", z:" << acts_global_postion[2] << endmsg;
+
+ // create and store the measurement
+ Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Identity();
+ cov(0, 0) = std::max<double>(std::abs(acts_global_postion[0] - simhit.getPosition()[0]), eps);
+ cov(1, 1) = std::max<double>(std::abs(acts_global_postion[1] - simhit.getPosition()[1]), eps);
+ measurements.emplace_back(Acts::Measurement<Acts::BoundIndices, 2>(sl, indices, par, cov));
+
+ // std::vector<std::string> truth_col = {std::to_string(m_layer+4), std::to_string(simhit.getPosition()[0]), std::to_string(simhit.getPosition()[1]), std::to_string(simhit.getPosition()[2])};
+ // truth_writer.write_row(truth_col);
+ // std::vector<std::string> converted_col = {std::to_string(m_layer+4), std::to_string(acts_global_postion[0]), std::to_string(acts_global_postion[1]), std::to_string(acts_global_postion[2])};
+ // converted_writer.write_row(converted_col);
+ }
+ }
+ } else { success = 0; }
+
+ return success;
+}
+
+int RecActsTracking::InitialiseSIT()
+{
+ int success = 1;
+
+ const edm4hep::TrackerHitCollection* hitSITCol = nullptr;
+ const edm4hep::SimTrackerHitCollection* SimhitSITCol = nullptr;
+
+ double min_z = 0;
+ try {
+ hitSITCol = _inSITTrackHdl.get();
+ } catch (GaudiException& e) {
+ debug() << "Collection " << _inSITTrackHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ success = 0;
+ }
+
+ try {
+ SimhitSITCol = _inSITColHdl.get();
+ } catch (GaudiException& e) {
+ debug() << "Sim Collection " << _inSITColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ success = 0;
+ }
+
+ if(hitSITCol && SimhitSITCol)
+ {
+ int nelem = hitSITCol->size();
+ debug() << "Number of SIT hits = " << nelem << endmsg;
+ // SpacePointPtrs.resize(nelem);
+
+ // std::string truth_file = "obj/sit/truth/event" + std::to_string(_nEvt) + ".csv";
+ // std::ofstream truth_stream(truth_file);
+ // csv2::Writer<csv2::delimiter<','>> truth_writer(truth_stream);
+ // std::vector<std::string> truth_header = {"layer", "x", "y", "z"};
+ // truth_writer.write_row(truth_header);
+
+ // std::string converted_file = "obj/sit/converted/event" + std::to_string(_nEvt) + ".csv";
+ // std::ofstream converted_stream(converted_file);
+ // csv2::Writer<csv2::delimiter<','>> converted_writer(converted_stream);
+ // std::vector<std::string> converted_header = {"layer", "x", "y", "z"};
+ // converted_writer.write_row(converted_header);
+
+ for (int ielem = 0; ielem < nelem; ++ielem)
+ {
+ auto hit = hitSITCol->at(ielem);
+ auto simhit = SimhitSITCol->at(ielem);
+
+ auto simcellid = simhit.getCellID();
+ // <id>system:5,side:-2,layer:9,stave:8,module:8,sensor:5,y:-11,z:-11</id>
+ uint64_t m_layer = sit_decoder->get(simcellid, "layer");
+ uint64_t m_stave = sit_decoder->get(simcellid, "stave");
+ uint64_t m_module = sit_decoder->get(simcellid, "module");
+ uint64_t m_sensor = sit_decoder->get(simcellid, "sensor");
+ double acts_x = simhit.getPosition()[0];
+ double acts_y = simhit.getPosition()[1];
+ double acts_z = simhit.getPosition()[2];
+ double momentum_x = simhit.getMomentum()[0];
+ double momentum_y = simhit.getMomentum()[1];
+ double momentum_z = simhit.getMomentum()[2];
+
+ dd4hep::rec::ISurface* surface = nullptr;
+ auto it = m_sit_surfaces->find(simcellid);
+ if (it != m_sit_surfaces->end()) {
+ surface = it->second;
+ if (!surface) {
+ fatal() << "found surface for SIT cell id " << simcellid << ", but NULL" << endmsg;
+ return 0;
+ }
+ }
+ else {
+ fatal() << "not found surface for SIT cell id " << simcellid << endmsg;
+ return 0;
+ }
+
+ if (ielem == 0) {
+ min_z = hitSITCol->at(ielem).getPosition()[2];
+ } else {
+ if (std::abs(acts_z - min_z) > 100) {
+ continue;
+ }
+ }
+
+ // set acts geometry identifier
+ uint64_t acts_volume = SIT_acts_volume_ids[m_layer];
+ uint64_t acts_boundary = 0;
+ uint64_t acts_layer = 2;
+ uint64_t acts_approach = 0;
+ uint64_t acts_sensitive = m_stave*SIT_module_nums[m_layer]*SIT_sensor_nums + m_module*SIT_sensor_nums + m_sensor + 1;
+
+ Acts::GeometryIdentifier moduleGeoId;
+ moduleGeoId.setVolume(acts_volume);
+ moduleGeoId.setBoundary(acts_boundary);
+ moduleGeoId.setLayer(acts_layer);
+ moduleGeoId.setApproach(acts_approach);
+ moduleGeoId.setSensitive(acts_sensitive);
+
+ // create and store the source link
+ uint32_t measurementIdx = measurements.size();
+ IndexSourceLink sourceLink{moduleGeoId, measurementIdx, hit};
+ sourceLinks.insert(sourceLinks.end(), sourceLink);
+ Acts::SourceLink sl{sourceLink};
+
+ // get the local position of the hit
+ IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
+ const Acts::Surface* acts_surface = surfaceAccessor(sl);
+ const Acts::Vector3 globalPos{acts_x, acts_y, acts_z};
+ const Acts::Vector3 globalmom{momentum_x, momentum_y, momentum_z};
+ auto acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, onSurfaceTolerance);
+ if (!acts_local_postion.ok()){
+ info() << "Error: failed to get local position for SIT hit " << simcellid << endmsg;
+ acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, 100*onSurfaceTolerance);
+ }
+ const std::array<Acts::BoundIndices, 2> indices{Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc1};
+ const Acts::Vector2 par{acts_local_postion.value()[0], acts_local_postion.value()[1]};
+
+ // *** debug ***
+ debug() << "SIT measurements global position(x,y,z): " << simhit.getPosition()[0] << ", " << simhit.getPosition()[1] << ", " << simhit.getPosition()[2]
+ << "; local position(loc0, loc1): "<< acts_local_postion.value()[0] << ", " << acts_local_postion.value()[1] << endmsg;
+ auto acts_global_postion = acts_surface->localToGlobal(geoContext, par, globalFakeMom);
+ debug() << "debug surface at: x:" << acts_global_postion[0] << ", y:" << acts_global_postion[1] << ", z:" << acts_global_postion[2] << endmsg;
+
+ // create and store the measurement
+ Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Identity();
+ cov(0, 0) = std::max<double>(std::abs(acts_global_postion[0] - simhit.getPosition()[0]), eps);
+ cov(1, 1) = std::max<double>(std::abs(acts_global_postion[1] - simhit.getPosition()[1]), eps);
+ measurements.emplace_back(Acts::Measurement<Acts::BoundIndices, 2>(sl, indices, par, cov));
+
+ // std::vector<std::string> truth_col = {std::to_string(m_layer), std::to_string(simhit.getPosition()[0]), std::to_string(simhit.getPosition()[1]), std::to_string(simhit.getPosition()[2])};
+ // truth_writer.write_row(truth_col);
+ // std::vector<std::string> converted_col = {std::to_string(m_layer), std::to_string(acts_global_postion[0]), std::to_string(acts_global_postion[1]), std::to_string(acts_global_postion[2])};
+ // converted_writer.write_row(converted_col);
+ }
+ } else { success = 0; }
+
+ return success;
+}
\ No newline at end of file
diff --git a/Reconstruction/RecActsTracking/src/RecActsTracking.h b/Reconstruction/RecActsTracking/src/RecActsTracking.h
new file mode 100644
index 00000000..be7c889a
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/RecActsTracking.h
@@ -0,0 +1,288 @@
+#ifndef RecActsTracking_H
+#define RecActsTracking_H
+
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+#include <sstream>
+#include <filesystem>
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "DD4hep/Detector.h"
+#include "DDRec/DetectorData.h"
+#include "DDRec/ISurface.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Vector3D.h"
+
+#include "UTIL/ILDConf.h"
+#include "GaudiKernel/NTuple.h"
+#include "DetInterface/IGeomSvc.h"
+
+// gear
+#include "GearSvc/IGearSvc.h"
+#include <gear/GEAR.h>
+#include <gear/GearMgr.h>
+#include <gear/GearParameters.h>
+#include <gear/VXDLayerLayout.h>
+#include <gear/VXDParameters.h>
+#include "gear/FTDLayerLayout.h"
+#include "gear/FTDParameters.h"
+#include <gear/BField.h>
+
+// edm4hep
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/MutableTrack.h"
+// #include "edm4hep/TrackerHit.h"
+// #include "edm4hep/SimTrackerHit.h"
+#include "edm4hep/TrackState.h"
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+// acts
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Direction.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/ProxyAccessor.hpp"
+#include "Acts/EventData/SpacePointData.hpp"
+#include <Acts/EventData/Measurement.hpp>
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/EventData/TrackContainer.hpp"
+#include "Acts/EventData/VectorMultiTrajectory.hpp"
+#include "Acts/EventData/VectorTrackContainer.hpp"
+#include "Acts/EventData/ParticleHypothesis.hpp"
+
+#include "Acts/Propagator/AbortList.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/MaterialInteractor.hpp"
+#include "Acts/Propagator/Navigator.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Propagator/StandardAborters.hpp"
+
+#include "Acts/Geometry/Extent.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include <Acts/Geometry/GeometryContext.hpp>
+
+#include "Acts/Seeding/BinnedGroup.hpp"
+#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
+#include "Acts/Seeding/InternalSpacePoint.hpp"
+#include "Acts/Seeding/SeedFilter.hpp"
+#include "Acts/Seeding/SeedFilterConfig.hpp"
+#include "Acts/Seeding/SeedFinder.hpp"
+#include "Acts/Seeding/SeedFinderConfig.hpp"
+#include "Acts/Seeding/SpacePointGrid.hpp"
+
+#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
+#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
+#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
+#include "Acts/TrackFitting/KalmanFitter.hpp"
+
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+
+#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/Delegate.hpp"
+#include "Acts/Utilities/Grid.hpp"
+#include "Acts/Utilities/GridBinFinder.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/Enumerate.hpp"
+#include "Acts/Utilities/TrackHelpers.hpp"
+
+// local include
+#include "utils/TGeoDetector.hpp"
+#include "utils/CKFhelper.hpp"
+#include "utils/MagneticField.hpp"
+
+namespace gear{
+ class GearMgr;
+}
+
+namespace dd4hep {
+ namespace DDSegmentation {
+ class BitFieldCoder;
+ }
+ namespace rec{
+ class ISurface;
+ }
+}
+
+// Seeding: Construct track seeds from space points.
+// config for seed finding
+struct SeedingConfig
+{
+ /// Input space point collections.
+ ///
+ /// We allow multiple space point collections to allow different parts of
+ /// the detector to use different algorithms for space point construction,
+ /// e.g. single-hit space points for pixel-like detectors or double-hit
+ /// space points for strip-like detectors.
+ std::vector<std::string> inputSpacePoints;
+
+ /// Output track seed collection.
+ std::string outputSeeds;
+
+
+ Acts::SeedFilterConfig seedFilterConfig;
+ Acts::SeedFinderConfig<SimSpacePoint> seedFinderConfig;
+ Acts::CylindricalSpacePointGridConfig gridConfig;
+ Acts::CylindricalSpacePointGridOptions gridOptions;
+ Acts::SeedFinderOptions seedFinderOptions;
+
+ // allow for different values of rMax in gridConfig and seedFinderConfig
+ bool allowSeparateRMax = false;
+
+ // vector containing the map of z bins in the top and bottom layers
+ std::vector<std::pair<int, int>> zBinNeighborsTop;
+ std::vector<std::pair<int, int>> zBinNeighborsBottom;
+
+ // number of phiBin neighbors at each side of the current bin that will be
+ // used to search for SPs
+ int numPhiNeighbors = 1;
+};
+
+class RecActsTracking : public GaudiAlgorithm
+{
+
+ public :
+
+ RecActsTracking(const std::string& name, ISvcLocator* svcLoc);
+
+ virtual StatusCode initialize();
+
+ virtual StatusCode execute();
+
+ virtual StatusCode finalize();
+
+ private :
+
+ // Input collections
+ DataHandle<edm4hep::TrackerHitCollection> _inVTXTrackHdl{"VXDTrackerHits", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> _inSITTrackHdl{"SITTrackerHits", Gaudi::DataHandle::Reader, this};
+
+ DataHandle<edm4hep::SimTrackerHitCollection> _inVTXColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> _inSITColHdl{"SITCollection", Gaudi::DataHandle::Reader, this};
+
+ DataHandle<edm4hep::MCParticleCollection> _inMCColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
+
+ // Output collections
+ DataHandle<edm4hep::TrackCollection> _outColHdl{"ACTSSiTracks", Gaudi::DataHandle::Writer, this};
+
+ // properties
+ Gaudi::Property<std::string> TGeo_path{this, "TGeoFile"};
+ Gaudi::Property<std::string> TGeo_config_path{this, "TGeoConfigFile"};
+ Gaudi::Property<std::string> MaterialMap_path{this, "MaterialMapFile"};
+ Gaudi::Property<double> m_field{this, "Field", 3.0}; // tesla
+ Gaudi::Property<double> onSurfaceTolerance{this, "onSurfaceTolerance", 1e-2}; // mm
+
+ SmartIF<IGeomSvc> m_geosvc;
+ SmartIF<IChronoStatSvc> chronoStatSvc;
+ dd4hep::DDSegmentation::BitFieldCoder *vxd_decoder;
+ dd4hep::DDSegmentation::BitFieldCoder *sit_decoder;
+ const dd4hep::rec::SurfaceMap* m_vtx_surfaces;
+ const dd4hep::rec::SurfaceMap* m_sit_surfaces;
+
+ // configs to build acts geometry
+ Acts::GeometryContext geoContext;
+ Acts::MagneticFieldContext magFieldContext;
+ Acts::CalibrationContext calibContext;
+ std::string TGeo_ROOTFilePath;
+ std::string TGeoConfig_jFilePath;
+ std::string MaterialMap_jFilePath;
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
+ std::vector<std::shared_ptr<const Acts::TGeoDetectorElement>> detElementStore;
+
+ // Store the space points & measurements & sourcelinks **per event**
+ // * only store the VXD tracker hits to the SpacePointPtrs
+ std::vector<const SimSpacePoint*> SpacePointPtrs;
+ IndexSourceLinkContainer sourceLinks;
+ std::vector<::Acts::BoundVariantMeasurement> measurements;
+ std::vector<::Acts::BoundTrackParameters> initialParameters;
+ SimSeedContainer Selected_Seeds;
+
+ // seed finder
+ SeedingConfig seed_cfg;
+ std::unique_ptr<const Acts::GridBinFinder<2ul>> m_bottomBinFinder;
+ std::unique_ptr<const Acts::GridBinFinder<2ul>> m_topBinFinder;
+ Acts::SeedFinder<SimSpacePoint, Acts::CylindricalSpacePointGrid<SimSpacePoint>> m_seedFinder;
+
+ int InitialiseVTX();
+ int InitialiseSIT();
+ const dd4hep::rec::ISurface* getISurface(edm4hep::TrackerHit* hit);
+ const Acts::GeometryIdentifier getVTXGid(uint64_t cellid);
+ const Acts::GeometryIdentifier getSITGid(uint64_t cellid);
+
+ // utils
+ int _nEvt;
+ double eps = 1.0e-05;
+ const double _FCT = 2.99792458E-4;
+ Acts::Vector3 acts_field_value = Acts::Vector3(0., 0., 3*_FCT); // tesla
+ // Acts::Vector3 acts_field_value = Acts::Vector3(0., 0., 3); // tesla
+ std::shared_ptr<const Acts::MagneticFieldProvider> magneticField
+ = std::make_shared<Acts::ConstantBField>(acts_field_value);
+ double bFieldMin = 0.3 * _FCT * Acts::UnitConstants::T; // tesla
+ // double bFieldMin = 0.1 * Acts::UnitConstants::T; // tesla
+ const Acts::Vector3 globalFakeMom{1.e6, 1.e6, 1.e6};
+ const std::array<Acts::BoundIndices, 6> writeout_indices{
+ Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundPhi,
+ Acts::BoundIndices::eBoundQOverP, Acts::BoundIndices::eBoundLoc1,
+ Acts::BoundIndices::eBoundTheta, Acts::BoundIndices::eBoundTime};
+
+ // param estimate configuration
+ double noTimeVarInflation = 100.;
+ std::array<double, 6> initialSigmas = {
+ 25 * Acts::UnitConstants::um,
+ 100 * Acts::UnitConstants::um,
+ 0.02 * Acts::UnitConstants::degree,
+ 0.02 * Acts::UnitConstants::degree,
+ 0.1 * Acts::UnitConstants::e / Acts::UnitConstants::GeV,
+ 1400 * Acts::UnitConstants::s};
+ std::array<double, 6> initialSimgaQoverPCoefficients = {
+ 0 * Acts::UnitConstants::mm / (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
+ 0 * Acts::UnitConstants::mm / (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
+ 0 * Acts::UnitConstants::degree / (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
+ 0 * Acts::UnitConstants::degree / (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
+ 0.1,
+ 0 * Acts::UnitConstants::ns / (Acts::UnitConstants::e * Acts::UnitConstants::GeV)};
+ std::array<double, 6> initialVarInflation = {10., 10., 10., 10., 10., 10.};
+ Acts::ParticleHypothesis particleHypothesis = Acts::ParticleHypothesis::muon();
+ // Acts::ParticleHypothesis particleHypothesis = Acts::ParticleHypothesis::chargedGeantino();
+
+ // gid convert configuration
+ std::vector<uint64_t> VXD_inner_volume_ids{16, 17, 18, 19};
+ uint64_t VXD_outer_volume_id = 20;
+ std::vector<uint64_t> SIT_acts_volume_ids{22, 24, 26};
+ std::vector<uint64_t> SIT_module_nums{7, 10, 14};
+ uint64_t SIT_sensor_nums = 14;
+
+ // CKF configuration
+ // Acts::MeasurementSelector::Config measurementSelectorCfg;
+ std::shared_ptr<TrackFinderFunction> findTracks;
+ std::optional<Acts::TrackSelector> m_trackSelector;
+ std::optional<std::variant<Acts::TrackSelector::Config, Acts::TrackSelector::EtaBinnedConfig>> trackSelectorCfg = std::nullopt;
+ mutable std::atomic<std::size_t> m_nTotalSeeds{0};
+ mutable std::atomic<std::size_t> m_nDeduplicatedSeeds{0};
+ mutable std::atomic<std::size_t> m_nFailedSeeds{0};
+ mutable std::atomic<std::size_t> m_nFailedSmoothing{0};
+ mutable std::atomic<std::size_t> m_nFailedExtrapolation{0};
+ mutable std::atomic<std::size_t> m_nFoundTracks{0};
+ mutable std::atomic<std::size_t> m_nSelectedTracks{0};
+ mutable std::atomic<std::size_t> m_nStoppedBranches{0};
+
+ mutable tbb::combinable<Acts::VectorMultiTrajectory::Statistics> m_memoryStatistics{[]() {
+ auto mtj = std::make_shared<Acts::VectorMultiTrajectory>();
+ return mtj->statistics();
+ }};
+ bool twoWay = false; /// Run finding in two directions
+ Acts::TrackExtrapolationStrategy extrapolationStrategy = Acts::TrackExtrapolationStrategy::firstOrLast; /// Extrapolation strategy
+ unsigned int maxSteps = 100000;
+};
+
+#endif // RecActsTracking_H
diff --git a/Reconstruction/RecActsTracking/src/utils/CKFhelper.hpp b/Reconstruction/RecActsTracking/src/utils/CKFhelper.hpp
new file mode 100644
index 00000000..71dae438
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/CKFhelper.hpp
@@ -0,0 +1,388 @@
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Direction.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+
+#include "Acts/EventData/ProxyAccessor.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/SourceLink.hpp"
+#include "Acts/EventData/TrackContainer.hpp"
+#include "Acts/EventData/TrackProxy.hpp"
+#include <Acts/EventData/Measurement.hpp>
+#include "Acts/EventData/VectorMultiTrajectory.hpp"
+#include "Acts/EventData/VectorTrackContainer.hpp"
+
+#include "ActsFatras/Digitization/Segmentizer.hpp"
+
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+
+#include "Acts/Propagator/AbortList.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/MaterialInteractor.hpp"
+#include "Acts/Propagator/Navigator.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Propagator/StandardAborters.hpp"
+
+
+#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
+#include "Acts/TrackFinding/MeasurementSelector.hpp"
+#include "Acts/TrackFinding/TrackSelector.hpp"
+#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
+#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
+#include "Acts/TrackFitting/KalmanFitter.hpp"
+
+#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/Result.hpp"
+#include "Acts/Utilities/TrackHelpers.hpp"
+#include "Acts/Utilities/Delegate.hpp"
+#include "Acts/Utilities/Enumerate.hpp"
+#include "Acts/Utilities/TrackHelpers.hpp"
+#include "Acts/Utilities/CalibrationContext.hpp"
+
+#include <atomic>
+#include <cstddef>
+#include <functional>
+#include <limits>
+#include <memory>
+#include <optional>
+#include <string>
+#include <variant>
+#include <vector>
+#include <cmath>
+#include <ostream>
+#include <stdexcept>
+#include <system_error>
+#include <unordered_map>
+#include <utility>
+
+#include <tbb/combinable.h>
+#include <boost/functional/hash.hpp>
+
+#include "SimSpacePoint.hpp"
+
+namespace Acts
+{
+ class MagneticFieldProvider;
+ class TrackingGeometry;
+}
+
+using Updater = Acts::GainMatrixUpdater;
+using Smoother = Acts::GainMatrixSmoother;
+using Stepper = Acts::EigenStepper<>;
+using Navigator = Acts::Navigator;
+using Propagator = Acts::Propagator<Stepper, Navigator>;
+using CKF = Acts::CombinatorialKalmanFilter<Propagator, Acts::VectorMultiTrajectory>;
+
+// track container types
+using TrackContainer = Acts::TrackContainer<Acts::VectorTrackContainer, Acts::VectorMultiTrajectory, std::shared_ptr>;
+using ConstTrackContainer = Acts::TrackContainer<Acts::ConstVectorTrackContainer, Acts::ConstVectorMultiTrajectory, std::shared_ptr>;
+using TrackParameters = ::Acts::BoundTrackParameters;
+using TrackParametersContainer = std::vector<TrackParameters>;
+using TrackIndexType = TrackContainer::IndexType;
+using TrackProxy = TrackContainer::TrackProxy;
+using ConstTrackProxy = ConstTrackContainer::ConstTrackProxy;
+
+// track finder types
+using TrackFinderOptions = Acts::CombinatorialKalmanFilterOptions<IndexSourceLinkAccessor::Iterator, Acts::VectorMultiTrajectory>;
+using TrackFinderResult = Acts::Result<std::vector<TrackContainer::TrackProxy>>;
+
+// measurement types
+using Measurement = ::Acts::BoundVariantMeasurement;
+using MeasurementContainer = std::vector<Measurement>;
+
+class TrackFinderFunction
+{
+public:
+ virtual ~TrackFinderFunction() = default;
+ virtual TrackFinderResult operator()(const TrackParameters&, const TrackFinderOptions&, TrackContainer&) const = 0;
+};
+
+static std::shared_ptr<TrackFinderFunction> makeTrackFinderFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ std::shared_ptr<const Acts::MagneticFieldProvider> magneticField,
+ const Acts::Logger& logger);
+
+struct TrackFinderFunctionImpl : public TrackFinderFunction {
+ CKF trackFinder;
+ TrackFinderFunctionImpl(CKF&& f) : trackFinder(std::move(f)) {}
+ TrackFinderResult operator()( const TrackParameters& initialParameters,
+ const TrackFinderOptions& options,
+ TrackContainer& tracks) const override
+ {
+ return trackFinder.findTracks(initialParameters, options, tracks);
+ };
+};
+
+std::shared_ptr<TrackFinderFunction> makeTrackFinderFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ std::shared_ptr<const Acts::MagneticFieldProvider> magneticField)
+{
+ Stepper stepper(magneticField);
+ Navigator::Config cfg{trackingGeometry};
+ cfg.resolvePassive = false;
+ cfg.resolveMaterial = true;
+ cfg.resolveSensitive = true;
+ Navigator navigator(cfg);
+ Propagator propagator(std::move(stepper), std::move(navigator));
+ CKF trackFinder(std::move(propagator));
+
+ return std::make_shared<TrackFinderFunctionImpl>(std::move(trackFinder));
+}
+
+struct Cluster
+{
+ using Cell = ActsFatras::Segmentizer::ChannelSegment;
+ std::size_t sizeLoc0 = 0;
+ std::size_t sizeLoc1 = 0;
+ std::vector<Cell> channels;
+};
+
+using ClusterContainer = std::vector<Cluster>;
+
+/// Abstract base class for measurement-based calibration
+class MeasurementCalibrator
+{
+ public:
+ virtual void calibrate(
+ const MeasurementContainer& measurements, const ClusterContainer* clusters,
+ const Acts::GeometryContext& gctx, const Acts::CalibrationContext& cctx,
+ const Acts::SourceLink& sourceLink,
+ Acts::VectorMultiTrajectory::TrackStateProxy& trackState) const = 0;
+ virtual ~MeasurementCalibrator() = default;
+ virtual bool needsClusters() const { return false; }
+};
+
+// Calibrator to convert an index source link to a measurement as-is
+class PassThroughCalibrator : public MeasurementCalibrator
+{
+ public:
+ /// Find the measurement corresponding to the source link.
+ ///
+ /// @tparam parameters_t Track parameters type
+ /// @param gctx The geometry context (unused)
+ /// @param trackState The track state to calibrate
+ void calibrate(
+ const MeasurementContainer& measurements,
+ const ClusterContainer* clusters, const Acts::GeometryContext& gctx,
+ const Acts::CalibrationContext& cctx, const Acts::SourceLink& sourceLink,
+ Acts::VectorMultiTrajectory::TrackStateProxy& trackState) const override;
+};
+
+// Adapter class that wraps a MeasurementCalibrator to conform to the
+// core ACTS calibration interface
+class MeasurementCalibratorAdapter
+{
+ public:
+ MeasurementCalibratorAdapter(const MeasurementCalibrator& calibrator,
+ const MeasurementContainer& measurements,
+ const ClusterContainer* clusters = nullptr);
+ MeasurementCalibratorAdapter() = delete;
+
+ void calibrate(const Acts::GeometryContext& gctx,
+ const Acts::CalibrationContext& cctx,
+ const Acts::SourceLink& sourceLink,
+ Acts::VectorMultiTrajectory::TrackStateProxy trackState) const;
+
+ private:
+ const MeasurementCalibrator& m_calibrator;
+ const MeasurementContainer& m_measurements;
+ const ClusterContainer* m_clusters;
+};
+
+void PassThroughCalibrator::calibrate(
+ const MeasurementContainer& measurements, const ClusterContainer* /*clusters*/,
+ const Acts::GeometryContext& /*gctx*/, const Acts::CalibrationContext& /*cctx*/,
+ const Acts::SourceLink& sourceLink, Acts::VectorMultiTrajectory::TrackStateProxy& trackState) const
+{
+ trackState.setUncalibratedSourceLink(sourceLink);
+ const IndexSourceLink& idxSourceLink = sourceLink.get<IndexSourceLink>();
+
+ assert((idxSourceLink.index() < measurements.size()) &&
+ "Source link index is outside the container bounds");
+
+ std::visit(
+ [&trackState](const auto& meas) { trackState.setCalibrated(meas); },
+ measurements[idxSourceLink.index()]);
+}
+
+MeasurementCalibratorAdapter::MeasurementCalibratorAdapter(
+ const MeasurementCalibrator& calibrator, const MeasurementContainer& measurements,
+ const ClusterContainer* clusters) : m_calibrator{calibrator}, m_measurements{measurements}, m_clusters{clusters} {}
+
+void MeasurementCalibratorAdapter::calibrate(
+ const Acts::GeometryContext& gctx, const Acts::CalibrationContext& cctx, const Acts::SourceLink& sourceLink,
+ Acts::VectorMultiTrajectory::TrackStateProxy trackState) const
+{
+ return m_calibrator.calibrate(m_measurements, m_clusters, gctx, cctx, sourceLink, trackState);
+}
+
+// Specialize std::hash for SeedIdentifier
+// This is required to use SeedIdentifier as a key in an `std::unordered_map`.
+template <class T, std::size_t N>
+struct std::hash<std::array<T, N>>
+{
+ std::size_t operator()(const std::array<T, N>& array) const
+ {
+ std::hash<T> hasher;
+ std::size_t result = 0;
+ for (auto&& element : array) { boost::hash_combine(result, hasher(element)); }
+ return result;
+ }
+};
+
+// Measurement selector for seed
+class MeasurementSelector
+{
+ public:
+ using Traj = Acts::VectorMultiTrajectory;
+ explicit MeasurementSelector(Acts::MeasurementSelector selector)
+ : m_selector(std::move(selector)) {}
+
+ void setSeed(const std::optional<SimSeed>& seed) { m_seed = seed; }
+
+ Acts::Result<std::pair<std::vector<Traj::TrackStateProxy>::iterator,
+ std::vector<Traj::TrackStateProxy>::iterator>>
+ select(std::vector<Traj::TrackStateProxy>& candidates,
+ bool& isOutlier, const Acts::Logger& logger) const
+ {
+ if (m_seed.has_value())
+ {
+ std::vector<Traj::TrackStateProxy> newCandidates;
+ for (const auto& candidate : candidates)
+ {
+ if (isSeedCandidate(candidate)) { newCandidates.push_back(candidate); }
+ }
+
+ if (!newCandidates.empty()) { candidates = std::move(newCandidates); }
+ }
+
+ return m_selector.select<Acts::VectorMultiTrajectory>(candidates, isOutlier, logger);
+ }
+
+ private:
+ Acts::MeasurementSelector m_selector;
+ std::optional<SimSeed> m_seed;
+
+ bool isSeedCandidate(const Traj::TrackStateProxy& candidate) const
+ {
+ assert(candidate.hasUncalibratedSourceLink());
+ const Acts::SourceLink& sourceLink = candidate.getUncalibratedSourceLink();
+ for (const auto& sp : m_seed->sp())
+ {
+ for (const auto& sl : sp->sourceLinks())
+ {
+ if (sourceLink.get<IndexSourceLink>() == sl.get<IndexSourceLink>()) { return true; }
+ }
+ }
+ return false;
+ }
+
+}; // class MeasurementSelector
+
+/// Source link indices of the bottom, middle, top measurements.
+/// * In case of strip seeds only the first source link of the pair is used.
+using SeedIdentifier = std::array<Index, 3>;
+
+/// Build a seed identifier from a seed.
+///
+/// @param seed The seed to build the identifier from.
+/// @return The seed identifier.
+SeedIdentifier makeSeedIdentifier(const SimSeed& seed)
+{
+ SeedIdentifier result;
+ for (const auto& [i, sp] : Acts::enumerate(seed.sp()))
+ {
+ const Acts::SourceLink& firstSourceLink = sp->sourceLinks().front();
+ result.at(i) = firstSourceLink.get<IndexSourceLink>().index();
+ }
+ return result;
+}
+
+/// Visit all possible seed identifiers of a track.
+///
+/// @param track The track to visit the seed identifiers of.
+/// @param visitor The visitor to call for each seed identifier.
+template <typename Visitor>
+void visitSeedIdentifiers(const TrackProxy& track, Visitor visitor)
+{
+ // first we collect the source link indices of the track states
+ std::vector<Index> sourceLinkIndices;
+ sourceLinkIndices.reserve(track.nMeasurements());
+ for (const auto& trackState : track.trackStatesReversed())
+ {
+ if (!trackState.hasUncalibratedSourceLink()) { continue; }
+ const Acts::SourceLink& sourceLink = trackState.getUncalibratedSourceLink();
+ sourceLinkIndices.push_back(sourceLink.get<IndexSourceLink>().index());
+ }
+
+ // then we iterate over all possible triplets and form seed identifiers
+ for (std::size_t i = 0; i < sourceLinkIndices.size(); ++i)
+ {
+ for (std::size_t j = i + 1; j < sourceLinkIndices.size(); ++j)
+ {
+ for (std::size_t k = j + 1; k < sourceLinkIndices.size(); ++k)
+ {
+ // Putting them into reverse order (k, j, i) to compensate for the `trackStatesReversed` above.
+ visitor({sourceLinkIndices.at(k), sourceLinkIndices.at(j), sourceLinkIndices.at(i)});
+ }
+ }
+ }
+}
+
+class BranchStopper
+{
+ public:
+ using Config = std::optional<std::variant<Acts::TrackSelector::Config, Acts::TrackSelector::EtaBinnedConfig>>;
+ using BranchStopperResult = Acts::CombinatorialKalmanFilterBranchStopperResult;
+
+ mutable std::atomic<std::size_t> m_nStoppedBranches{0};
+ explicit BranchStopper(const Config& config) : m_config(config) {}
+
+ BranchStopperResult operator()( const Acts::CombinatorialKalmanFilterTipState& tipState,
+ Acts::VectorMultiTrajectory::TrackStateProxy& trackState ) const
+ {
+ if (!m_config.has_value()) { return BranchStopperResult::Continue; }
+
+ const Acts::TrackSelector::Config* singleConfig = std::visit
+ (
+ [&](const auto& config) -> const Acts::TrackSelector::Config*
+ {
+ using T = std::decay_t<decltype(config)>;
+ if constexpr (std::is_same_v<T, Acts::TrackSelector::Config>) { return &config; }
+ else if constexpr (std::is_same_v<T, Acts::TrackSelector::EtaBinnedConfig>)
+ {
+ double theta = trackState.parameters()[Acts::eBoundTheta];
+ double eta = -std::log(std::tan(0.5 * theta));
+ return config.hasCuts(eta) ? &config.getCuts(eta) : nullptr;
+ }
+ }, *m_config
+ );
+
+ if (singleConfig == nullptr)
+ {
+ ++m_nStoppedBranches;
+ return BranchStopperResult::StopAndDrop;
+ }
+
+ bool enoughMeasurements =
+ tipState.nMeasurements >= singleConfig->minMeasurements;
+ bool tooManyHoles = tipState.nHoles > singleConfig->maxHoles;
+ bool tooManyOutliers = tipState.nOutliers > singleConfig->maxOutliers;
+
+ if (tooManyHoles || tooManyOutliers) {
+ ++m_nStoppedBranches;
+ return enoughMeasurements ? BranchStopperResult::StopAndKeep
+ : BranchStopperResult::StopAndDrop;
+ }
+
+ return BranchStopperResult::Continue;
+ }
+
+ private:
+ Config m_config;
+};
\ No newline at end of file
diff --git a/Reconstruction/RecActsTracking/src/utils/GeometryContainers.hpp b/Reconstruction/RecActsTracking/src/utils/GeometryContainers.hpp
new file mode 100644
index 00000000..6d4b6dcc
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/GeometryContainers.hpp
@@ -0,0 +1,355 @@
+// STL
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <iostream>
+#include <utility>
+#include <iterator>
+
+// boost
+#include <boost/container/flat_map.hpp>
+#include <boost/container/flat_set.hpp>
+
+// acts
+#include "Acts/EventData/SourceLink.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+
+template <typename Iterator>
+class Range
+{
+ public:
+ Range(Iterator b, Iterator e) : m_begin(b), m_end(e) {}
+ Range(Range&&) = default;
+ Range(const Range&) = default;
+ ~Range() = default;
+ Range& operator=(Range&&) = default;
+ Range& operator=(const Range&) = default;
+
+ Iterator begin() const { return m_begin; }
+ Iterator end() const { return m_end; }
+ bool empty() const { return m_begin == m_end; }
+ std::size_t size() const { return std::distance(m_begin, m_end); }
+
+ private:
+ Iterator m_begin;
+ Iterator m_end;
+};
+
+template <typename Iterator>
+Range<Iterator> makeRange(Iterator begin, Iterator end)
+{ return Range<Iterator>(begin, end); }
+
+template <typename Iterator>
+Range<Iterator> makeRange(std::pair<Iterator, Iterator> range)
+{ return Range<Iterator>(range.first, range.second); }
+
+
+/// Proxy for iterating over groups of elements within a container.
+///
+/// @note Each group will contain at least one element.
+///
+/// Consecutive elements with the same key (as defined by the KeyGetter) are
+/// placed in one group. The proxy should always be used as part of a
+/// range-based for loop. In combination with structured bindings to reduce the
+/// boilerplate, the group iteration can be written as
+///
+/// for (auto&& [key, elements] : GroupBy<...>(...)) {
+/// // do something with just the key
+/// ...
+///
+/// // iterate over the group elements
+/// for (const auto& element : elements) {
+/// ...
+/// }
+/// }
+///
+template <typename Iterator, typename KeyGetter>
+class GroupBy {
+ public:
+ /// The key type that identifies elements within a group.
+ using Key = std::decay_t<decltype(KeyGetter()(*Iterator()))>;
+ /// A Group is an iterator range with the associated key.
+ using Group = std::pair<Key, Range<Iterator>>;
+ /// Iterator type representing the end of the groups.
+ ///
+ /// The end iterator will not be dereferenced in C++17 range-based loops. It
+ /// can thus be a simpler type without the overhead of the full group iterator
+ /// below.
+ using GroupEndIterator = Iterator;
+ /// Iterator type representing a group of elements.
+ class GroupIterator {
+ public:
+ using iterator_category = std::input_iterator_tag;
+ using value_type = Group;
+ using difference_type = std::ptrdiff_t;
+ using pointer = Group*;
+ using reference = Group&;
+
+ constexpr GroupIterator(const GroupBy& groupBy, Iterator groupBegin)
+ : m_groupBy(groupBy),
+ m_groupBegin(groupBegin),
+ m_groupEnd(groupBy.findEndOfGroup(groupBegin)) {}
+ /// Pre-increment operator to advance to the next group.
+ constexpr GroupIterator& operator++() {
+ // make the current end the new group beginning
+ std::swap(m_groupBegin, m_groupEnd);
+ // find the end of the next group starting from the new beginning
+ m_groupEnd = m_groupBy.findEndOfGroup(m_groupBegin);
+ return *this;
+ }
+ /// Post-increment operator to advance to the next group.
+ constexpr GroupIterator operator++(int) {
+ GroupIterator retval = *this;
+ ++(*this);
+ return retval;
+ }
+ /// Dereference operator that returns the pointed-to group of elements.
+ constexpr Group operator*() const {
+ const Key key = (m_groupBegin != m_groupEnd)
+ ? m_groupBy.m_keyGetter(*m_groupBegin)
+ : Key();
+ return {key, makeRange(m_groupBegin, m_groupEnd)};
+ }
+
+ private:
+ const GroupBy& m_groupBy;
+ Iterator m_groupBegin;
+ Iterator m_groupEnd;
+
+ friend constexpr bool operator==(const GroupIterator& lhs,
+ const GroupEndIterator& rhs) {
+ return lhs.m_groupBegin == rhs;
+ }
+ friend constexpr bool operator!=(const GroupIterator& lhs,
+ const GroupEndIterator& rhs) {
+ return !(lhs == rhs);
+ }
+ };
+
+ /// Construct the group-by proxy for an iterator range.
+ constexpr GroupBy(Iterator begin, Iterator end,
+ KeyGetter keyGetter = KeyGetter())
+ : m_begin(begin), m_end(end), m_keyGetter(std::move(keyGetter)) {}
+ constexpr GroupIterator begin() const {
+ return GroupIterator(*this, m_begin);
+ }
+ constexpr GroupEndIterator end() const { return m_end; }
+ constexpr bool empty() const { return m_begin == m_end; }
+
+ private:
+ Iterator m_begin;
+ Iterator m_end;
+ KeyGetter m_keyGetter;
+
+ /// Find the end of the group that starts at the given position.
+ ///
+ /// This uses a linear search from the start position and thus has linear
+ /// complexity in the group size. It does not assume any ordering of the
+ /// underlying container and is a cache-friendly access pattern.
+ constexpr Iterator findEndOfGroup(Iterator start) const {
+ // check for end so we can safely dereference the start iterator.
+ if (start == m_end) {
+ return start;
+ }
+ // search the first element that does not share a key with the start.
+ return std::find_if_not(std::next(start), m_end,
+ [this, start](const auto& x) {
+ return m_keyGetter(x) == m_keyGetter(*start);
+ });
+ }
+};
+
+/// Construct the group-by proxy for a container.
+template <typename Container, typename KeyGetter>
+auto makeGroupBy(const Container& container, KeyGetter keyGetter)
+ -> GroupBy<decltype(std::begin(container)), KeyGetter> {
+ return {std::begin(container), std::end(container), std::move(keyGetter)};
+}
+
+// extract the geometry identifier from a variety of types
+struct GeometryIdGetter
+{
+ // explicit geometry identifier are just forwarded
+ constexpr Acts::GeometryIdentifier operator()(
+ Acts::GeometryIdentifier geometryId) const { return geometryId; }
+
+ // encoded geometry ids are converted back to geometry identifiers.
+ constexpr Acts::GeometryIdentifier operator()(
+ Acts::GeometryIdentifier::Value encoded) const { return Acts::GeometryIdentifier(encoded); }
+
+ // support elements in map-like structures.
+ template <typename T>
+ constexpr Acts::GeometryIdentifier operator()(
+ const std::pair<Acts::GeometryIdentifier, T>& mapItem) const { return mapItem.first; }
+
+ // support elements that implement `.geometryId()`.
+ template <typename T>
+ inline auto operator()(const T& thing) const ->
+ decltype(thing.geometryId(), Acts::GeometryIdentifier()) { return thing.geometryId(); }
+
+ // support reference_wrappers around such types as well
+ template <typename T>
+ inline auto operator()(std::reference_wrapper<T> thing) const ->
+ decltype(thing.get().geometryId(), Acts::GeometryIdentifier()) { return thing.get().geometryId(); }
+};
+
+struct CompareGeometryId
+{
+ // indicate that comparisons between keys and full objects are allowed.
+ using is_transparent = void;
+ // compare two elements using the automatic key extraction.
+ template <typename Left, typename Right>
+ constexpr bool operator()(Left&& lhs, Right&& rhs) const
+ { return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs); }
+};
+
+/// Store elements that know their detector geometry id, e.g. simulation hits.
+///
+/// @tparam T type to be stored, must be compatible with `CompareGeometryId`
+///
+/// The container stores an arbitrary number of elements for any geometry
+/// id. Elements can be retrieved via the geometry id; elements can be selected
+/// for a specific geometry id or for a larger range, e.g. a volume or a layer
+/// within the geometry hierarchy using the helper functions below. Elements can
+/// also be accessed by index that uniquely identifies each element regardless
+/// of geometry id.
+template <typename T>
+using GeometryIdMultiset = boost::container::flat_multiset<T, CompareGeometryId>;
+
+/// Store elements indexed by an geometry id.
+///
+/// @tparam T type to be stored
+///
+/// The behaviour is the same as for the `GeometryIdMultiset` except that the
+/// stored elements do not know their geometry id themself. When iterating
+/// the iterator elements behave as for the `std::map`, i.e.
+///
+/// for (const auto& entry: elements) {
+/// auto id = entry.first; // geometry id
+/// const auto& el = entry.second; // stored element
+/// }
+///
+template <typename T>
+using GeometryIdMultimap = GeometryIdMultiset<std::pair<Acts::GeometryIdentifier, T>>;
+
+/// Select all elements within the given volume.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator> selectVolume(
+ const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier::Value volume)
+{
+ auto cmp = Acts::GeometryIdentifier().setVolume(volume);
+ auto beg = std::lower_bound(container.begin(), container.end(), cmp, CompareGeometryId{});
+ // WARNING overflows to volume==0 if the input volume is the last one
+ cmp = Acts::GeometryIdentifier().setVolume(volume + 1u);
+ // optimize search by using the lower bound as start point. also handles
+ // volume overflows since the geo id would be located before the start of
+ // the upper edge search window.
+ auto end = std::lower_bound(beg, container.end(), cmp, CompareGeometryId{});
+ return makeRange(beg, end);
+}
+
+/// Select all elements within the given volume.
+template <typename T>
+inline auto selectVolume(const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier id)
+{
+ return selectVolume(container, id.volume());
+}
+
+/// Select all elements within the given layer.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator> selectLayer(
+ const GeometryIdMultiset<T>& container,
+ Acts::GeometryIdentifier::Value volume,
+ Acts::GeometryIdentifier::Value layer)
+{
+ auto cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer);
+ auto beg = std::lower_bound(container.begin(), container.end(), cmp, CompareGeometryId{});
+ // WARNING resets to layer==0 if the input layer is the last one
+ cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer + 1u);
+ // optimize search by using the lower bound as start point. also handles
+ // volume overflows since the geo id would be located before the start of
+ // the upper edge search window.
+ auto end = std::lower_bound(beg, container.end(), cmp, CompareGeometryId{});
+ return makeRange(beg, end);
+}
+
+// Select all elements within the given layer.
+template <typename T>
+inline auto selectLayer(const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier id)
+{
+ return selectLayer(container, id.volume(), id.layer());
+}
+
+/// Select all elements for the given module / sensitive surface.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator> selectModule(
+const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier geoId)
+{
+ // module is the lowest level and defines a single geometry id value
+ return makeRange(container.equal_range(geoId));
+}
+
+/// Select all elements for the given module / sensitive surface.
+template <typename T>
+inline auto selectModule(
+ const GeometryIdMultiset<T>& container,
+ Acts::GeometryIdentifier::Value volume,
+ Acts::GeometryIdentifier::Value layer,
+ Acts::GeometryIdentifier::Value module)
+{
+ return selectModule(container,
+ Acts::GeometryIdentifier().setVolume(volume).setLayer(layer).setSensitive(module));
+}
+
+/// Select all elements for the lowest non-zero identifier component.
+///
+/// Zero values of lower components are interpreted as wildcard search patterns
+/// that select all element at the given geometry hierarchy and below. This only
+/// applies to the lower components and not to intermediate zeros.
+///
+/// Examples:
+/// - volume=2,layer=0,module=3 -> select all elements in the module
+/// - volume=1,layer=2,module=0 -> select all elements in the layer
+/// - volume=3,layer=0,module=0 -> select all elements in the volume
+///
+/// @note An identifier with all components set to zero selects the whole input
+/// container.
+/// @note Boundary and approach surfaces do not really fit into the geometry
+/// hierarchy and must be set to zero for the selection. If they are set on an
+/// input identifier, the behaviour of this search method is undefined.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator>
+selectLowestNonZeroGeometryObject(const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier geoId)
+{
+ assert((geoId.boundary() == 0u) && "Boundary component must be zero");
+ assert((geoId.approach() == 0u) && "Approach component must be zero");
+
+ if (geoId.sensitive() != 0u) { return selectModule(container, geoId); }
+ else if (geoId.layer() != 0u) { return selectLayer(container, geoId); }
+ else if (geoId.volume() != 0u) { return selectVolume(container, geoId); }
+ else { return makeRange(container.begin(), container.end()); }
+}
+
+/// Iterate over groups of elements belonging to each module/ sensitive surface.
+template <typename T>
+inline GroupBy<typename GeometryIdMultiset<T>::const_iterator, GeometryIdGetter>
+groupByModule(const GeometryIdMultiset<T>& container)
+{
+ return makeGroupBy(container, GeometryIdGetter());
+}
+
+/// The accessor for the GeometryIdMultiset container
+///
+/// It wraps up a few lookup methods to be used in the Combinatorial Kalman
+/// Filter
+template <typename T>
+struct GeometryIdMultisetAccessor {
+using Container = GeometryIdMultiset<T>;
+using Key = Acts::GeometryIdentifier;
+using Value = typename GeometryIdMultiset<T>::value_type;
+using Iterator = typename GeometryIdMultiset<T>::const_iterator;
+
+// pointer to the container
+const Container* container = nullptr;
+};
\ No newline at end of file
diff --git a/Reconstruction/RecActsTracking/src/utils/MagneticField.hpp b/Reconstruction/RecActsTracking/src/utils/MagneticField.hpp
new file mode 100644
index 00000000..94ef4152
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/MagneticField.hpp
@@ -0,0 +1,121 @@
+// STL
+#include <memory>
+#include <variant>
+#include <vector>
+
+// acts
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
+#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
+#include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/MagneticField/NullBField.hpp"
+#include "Acts/Utilities/Grid.hpp"
+#include "Acts/Utilities/Result.hpp"
+#include "Acts/Utilities/detail/Axis.hpp"
+#include "Acts/Utilities/detail/AxisFwd.hpp"
+#include "Acts/Utilities/detail/grid_helper.hpp"
+
+/// The ScalableBField-specific magnetic field context.
+struct ScalableBFieldContext
+{ Acts::ActsScalar scalor = 1.; };
+
+/// A constant magnetic field that is scaled depending on the event context.
+class ScalableBField final : public Acts::MagneticFieldProvider
+{
+ public:
+ struct Cache
+ {
+ Acts::ActsScalar scalor = 1.;
+ /// @brief constructor with context
+ Cache(const Acts::MagneticFieldContext& mctx)
+ { scalor = mctx.get<const ScalableBFieldContext>().scalor; }
+ };
+
+ /// @brief construct constant magnetic field from field vector
+ ///
+ /// @param [in] B magnetic field vector in global coordinate system
+ explicit ScalableBField(Acts::Vector3 B) : m_BField(std::move(B)) {}
+
+ /// @brief construct constant magnetic field from components
+ ///
+ /// @param [in] Bx magnetic field component in global x-direction
+ /// @param [in] By magnetic field component in global y-direction
+ /// @param [in] Bz magnetic field component in global z-direction
+ ScalableBField(Acts::ActsScalar Bx = 0, Acts::ActsScalar By = 0, Acts::ActsScalar Bz = 0)
+ : m_BField(Bx, By, Bz) {}
+
+ /// @brief retrieve magnetic field value
+ ///
+ /// @param [in] position global position
+ /// @param [in] cache Cache object (is ignored)
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ Acts::Result<Acts::Vector3> getField(
+ const Acts::Vector3& /*position*/,
+ MagneticFieldProvider::Cache& gCache) const override
+ {
+ Cache& cache = gCache.as<Cache>();
+ return Acts::Result<Acts::Vector3>::success(m_BField * cache.scalor);
+ }
+
+ /// @brief retrieve magnetic field value & its gradient
+ ///
+ /// @param [in] position global position
+ /// @param [out] derivative gradient of magnetic field vector as (3x3)
+ /// matrix
+ /// @param [in] cache Cache object (is ignored)
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ /// @note currently the derivative is not calculated
+ /// @todo return derivative
+ Acts::Result<Acts::Vector3> getFieldGradient(
+ const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/,
+ MagneticFieldProvider::Cache& gCache) const override
+ {
+ Cache& cache = gCache.as<Cache>();
+ return Acts::Result<Acts::Vector3>::success(m_BField * cache.scalor);
+ }
+
+ Acts::MagneticFieldProvider::Cache makeCache(
+ const Acts::MagneticFieldContext& mctx) const override
+ {
+ return Acts::MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
+ }
+
+ /// @brief check whether given 3D position is inside look-up domain
+ ///
+ /// @param [in] position global 3D position
+ /// @return @c true if position is inside the defined look-up grid,
+ /// otherwise @c false
+ /// @note The method will always return true for the constant B-Field
+ bool isInside(const Acts::Vector3& /*position*/) const { return true; }
+
+ /// @brief update magnetic field vector from components
+ ///
+ /// @param [in] Bx magnetic field component in global x-direction
+ /// @param [in] By magnetic field component in global y-direction
+ /// @param [in] Bz magnetic field component in global z-direction
+ void setField(double Bx, double By, double Bz) { m_BField << Bx, By, Bz; }
+
+ /// @brief update magnetic field vector
+ ///
+ /// @param [in] B magnetic field vector in global coordinate system
+ void setField(const Acts::Vector3& B) { m_BField = B; }
+
+ private:
+ /// magnetic field vector
+ Acts::Vector3 m_BField;
+}; // ScalableBField
+
+using InterpolatedMagneticField2 = Acts::InterpolatedBFieldMap<
+ Acts::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>;
+
+using InterpolatedMagneticField3 = Acts::InterpolatedBFieldMap<
+ Acts::Grid<Acts::Vector3, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
diff --git a/Reconstruction/RecActsTracking/src/utils/Options.cpp b/Reconstruction/RecActsTracking/src/utils/Options.cpp
new file mode 100644
index 00000000..0ce24a87
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/Options.cpp
@@ -0,0 +1,166 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 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 "Options.hpp"
+
+#include <algorithm>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <utility>
+
+namespace {
+constexpr char s_separator = ':';
+}
+
+// interval
+
+std::istream& Options::operator>>(
+ std::istream& is, Options::Interval& interval) {
+ std::string buf;
+ is >> buf;
+
+ // default to an unbounded interval
+ interval.lower.reset();
+ interval.upper.reset();
+
+ // find the limit separator
+ auto pos = buf.find_first_of(s_separator);
+ // no separator -> invalid input -> unbounded interval
+ if (pos == std::string::npos) {
+ return is;
+ }
+
+ // if it exists, parse limit before separator
+ if (0 < pos) {
+ auto lowerStr = buf.substr(0, pos);
+ interval.lower = std::stod(lowerStr);
+ }
+ // if it exists, parse limit after separator
+ if ((pos + 1) < buf.size()) {
+ auto upperStr = buf.substr(pos + 1);
+ interval.upper = std::stod(upperStr);
+ }
+
+ return is;
+}
+
+std::ostream& Options::operator<<(
+ std::ostream& os, const Options::Interval& interval) {
+ if (!interval.lower.has_value() && !interval.upper.has_value()) {
+ os << "unbounded";
+ } else {
+ if (interval.lower.has_value()) {
+ os << interval.lower.value();
+ }
+ os << s_separator;
+ if (interval.upper.has_value()) {
+ os << interval.upper.value();
+ }
+ }
+ return os;
+}
+
+// helper functions to parse and print multiple values
+
+namespace {
+
+template <typename value_t, typename converter_t>
+void parseVariable(std::istream& is, std::vector<value_t>& values,
+ converter_t&& convert) {
+ values.clear();
+
+ std::string buf;
+ is >> buf;
+ std::string bufValue;
+ std::string::size_type pos = 0;
+ std::string::size_type end = std::string::npos;
+ do {
+ end = buf.find_first_of(s_separator, pos);
+ if (end == std::string::npos) {
+ // last element; take the rest of the buffer
+ bufValue = buf.substr(pos);
+ } else {
+ bufValue = buf.substr(pos, end - pos);
+ pos = end + 1u;
+ }
+ values.push_back(convert(bufValue));
+ } while (end != std::string::npos);
+}
+
+template <typename value_t, typename converter_t>
+void parseFixed(std::istream& is, std::size_t size, value_t* values,
+ converter_t&& convert) {
+ // reserve space for the expected number of values
+ std::vector<value_t> tmp(size, 0);
+ parseVariable(is, tmp, std::forward<converter_t>(convert));
+ if (tmp.size() < size) {
+ throw std::invalid_argument(
+ "Not enough values for fixed-size user option, expected " +
+ std::to_string(size) + " received " + std::to_string(tmp.size()));
+ }
+ if (size < tmp.size()) {
+ throw std::invalid_argument(
+ "Too many values for fixed-size user option, expected " +
+ std::to_string(size) + " received " + std::to_string(tmp.size()));
+ }
+ std::copy(tmp.begin(), tmp.end(), values);
+}
+
+template <typename value_t>
+void print(std::ostream& os, std::size_t size, const value_t* values) {
+ for (std::size_t i = 0; i < size; ++i) {
+ if (0u < i) {
+ os << s_separator;
+ }
+ os << values[i];
+ }
+}
+
+} // namespace
+
+// fixed and variable number of generic values
+
+void Options::detail::parseDoublesFixed(std::istream& is,
+ std::size_t size,
+ double* values) {
+ parseFixed(is, size, values,
+ [](const std::string& s) { return std::stod(s); });
+}
+
+void Options::detail::parseDoublesVariable(
+ std::istream& is, std::vector<double>& values) {
+ parseVariable(is, values, [](const std::string& s) { return std::stod(s); });
+}
+
+void Options::detail::printDoubles(std::ostream& os,
+ std::size_t size,
+ const double* values) {
+ print(os, size, values);
+}
+
+// fixed and variable number of integers
+
+void Options::detail::parseIntegersFixed(std::istream& is,
+ std::size_t size,
+ int* values) {
+ parseFixed(is, size, values,
+ [](const std::string& s) { return std::stoi(s); });
+}
+
+void Options::detail::parseIntegersVariable(
+ std::istream& is, std::vector<int>& values) {
+ parseVariable(is, values, [](const std::string& s) { return std::stoi(s); });
+}
+
+void Options::detail::printIntegers(std::ostream& os,
+ std::size_t size,
+ const int* values) {
+ print(os, size, values);
+}
diff --git a/Reconstruction/RecActsTracking/src/utils/Options.hpp b/Reconstruction/RecActsTracking/src/utils/Options.hpp
new file mode 100644
index 00000000..897233cb
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/Options.hpp
@@ -0,0 +1,159 @@
+#include <array>
+#include <cstddef>
+#include <iosfwd>
+#include <optional>
+#include <vector>
+
+namespace Options {
+
+/// @defgroup option-types Additional types for program options
+///
+/// All types are intended as utility type for the user options and not as a
+/// variable type for the configuration structs. They should only be used where
+/// a single option can not be represented by an existing primitive types.
+///
+/// They also must be distinct types and can not just be typedefs; otherwise we
+/// can not define the required operator{<<,>>} overloads in this namespace.
+///
+/// @{
+
+/// Half open [lower,upper) interval type for a single user option.
+///
+/// A missing limit represents an unbounded upper or lower limit. With just
+/// one defined limit the interval is just a lower/upper bound; with both
+/// limits undefined, the interval is unbounded everywhere and thus contains
+/// all possible values.
+struct Interval {
+ std::optional<double> lower;
+ std::optional<double> upper;
+};
+
+/// A fixed number of real values as one user option.
+///
+/// @note Implemented as a subclass so it is distinct from `std::array`
+/// and we can provide overloads in the same namespace.
+template <std::size_t kSize>
+class Reals : public std::array<double, kSize> {};
+
+/// An arbitrary number of revaluesal as one user option.
+///
+/// @note Making this a `std::vector<double>` typedef or subclass confuses
+/// program options, since `std::vector<double>` is interpreted as a `double`
+/// option that can be provided multiple times.
+struct VariableReals {
+ std::vector<double> values;
+};
+
+/// A fixed number of integers as one user option.
+///
+/// @note Implemented as a subclass so it is distinct from `std::array`
+/// and we can provide overloads in the same namespace.
+template <std::size_t kSize>
+class Integers : public std::array<int, kSize> {};
+
+/// An arbitrary number of integers as one user option.
+///
+/// @note Making this a `std::vector<int>` typedef or subclass confuses
+/// program options, since `std::vector<int>` is interpreted as an `int`
+/// option that can be provided multiple times.
+struct VariableIntegers {
+ std::vector<int> values;
+};
+
+/// @}
+
+/// Extract an interval from an input of the form 'lower:upper'.
+///
+/// An input of the form `lower:` or `:upper` sets just one of the limits. Any
+/// other input leads to an unbounded interval.
+///
+/// @note The more common range notation uses `lower-upper` but the `-`
+/// separator complicates the parsing of negative values.
+std::istream& operator>>(std::istream& is, Interval& interval);
+
+/// Print an interval as `lower:upper`.
+std::ostream& operator<<(std::ostream& os, const Interval& interval);
+
+namespace detail {
+void parseDoublesFixed(std::istream& is, std::size_t size, double* values);
+void parseDoublesVariable(std::istream& is, std::vector<double>& values);
+void printDoubles(std::ostream& os, std::size_t size, const double* values);
+} // namespace detail
+
+/// Extract a fixed number of doubles from an input of the form 'x:y:z'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+template <std::size_t kSize>
+inline std::istream& operator>>(std::istream& is, Reals<kSize>& values) {
+ detail::parseDoublesFixed(is, kSize, values.data());
+ return is;
+}
+
+/// Extract a variable number of doubles from an input of the form 'x:y:...'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+inline std::istream& operator>>(std::istream& is, VariableReals& values) {
+ detail::parseDoublesVariable(is, values.values);
+ return is;
+}
+
+/// Print a fixed number of doubles as `x:y:z`.
+template <std::size_t kSize>
+inline std::ostream& operator<<(std::ostream& os, const Reals<kSize>& values) {
+ detail::printDoubles(os, kSize, values.data());
+ return os;
+}
+
+/// Print a variable number of doubles as `x:y:z:...`.
+inline std::ostream& operator<<(std::ostream& os, const VariableReals& values) {
+ detail::printDoubles(os, values.values.size(), values.values.data());
+ return os;
+}
+
+namespace detail {
+void parseIntegersFixed(std::istream& is, std::size_t size, int* values);
+void parseIntegersVariable(std::istream& is, std::vector<int>& values);
+void printIntegers(std::ostream& os, std::size_t size, const int* values);
+} // namespace detail
+
+/// Extract a fixed number of integers from an input of the form 'x:y:z'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+template <std::size_t kSize>
+inline std::istream& operator>>(std::istream& is, Integers<kSize>& values) {
+ detail::parseIntegersFixed(is, kSize, values.data());
+ return is;
+}
+
+/// Extract a variable number of integers from an input of the form 'x:y:...'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+inline std::istream& operator>>(std::istream& is, VariableIntegers& values) {
+ detail::parseIntegersVariable(is, values.values);
+ return is;
+}
+
+/// Print a fixed number of integers as `x:y:z`.
+template <std::size_t kSize>
+inline std::ostream& operator<<(std::ostream& os,
+ const Integers<kSize>& values) {
+ detail::printIntegers(os, kSize, values.data());
+ return os;
+}
+
+/// Print a variable number of integers as `x:y:z:...`.
+inline std::ostream& operator<<(std::ostream& os,
+ const VariableIntegers& values) {
+ detail::printIntegers(os, values.values.size(), values.values.data());
+ return os;
+}
+
+} // namespace Options
diff --git a/Reconstruction/RecActsTracking/src/utils/SimSpacePoint.hpp b/Reconstruction/RecActsTracking/src/utils/SimSpacePoint.hpp
new file mode 100644
index 00000000..4c577ccd
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/SimSpacePoint.hpp
@@ -0,0 +1,243 @@
+// STL
+#include <cmath>
+#include <vector>
+#include <cstdint>
+
+// boost
+#include <boost/container/static_vector.hpp>
+
+// acts tools
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/EventData/SourceLink.hpp"
+#include "Acts/Seeding/Seed.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+
+// edm4hep
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/SimTrackerHit.h"
+
+// local
+#include "GeometryContainers.hpp"
+
+using Index = std::uint32_t;
+template <typename value_t>
+using IndexMultimap = boost::container::flat_multimap<Index, value_t>;
+
+template <typename value_t>
+inline boost::container::flat_multimap<value_t, Index> invertIndexMultimap(
+ const IndexMultimap<value_t>& multimap) {
+ using InverseMultimap = boost::container::flat_multimap<value_t, Index>;
+
+ // switch key-value without enforcing the new ordering (linear copy)
+ typename InverseMultimap::sequence_type unordered;
+ unordered.reserve(multimap.size());
+ for (auto&& [index, value] : multimap) {
+ // value is now the key and the index is now the value
+ unordered.emplace_back(value, index);
+ }
+
+ // adopting the unordered sequence will reestablish the correct order
+ InverseMultimap inverse;
+#if BOOST_VERSION < 107800
+ for (const auto& i : unordered) {
+ inverse.insert(i);
+ }
+#else
+ inverse.insert(unordered.begin(), unordered.end());
+#endif
+ return inverse;
+}
+
+/// Space point representation of a measurement suitable for track seeding.
+class SimSpacePoint {
+ using Scalar = Acts::ActsScalar;
+
+public:
+ /// Construct the space point from edm4hep::TrackerHit
+ ///
+ /// @param trackerhit tracker hit to construct the space point from
+ SimSpacePoint(const edm4hep::TrackerHit trackerhit,
+ float x, float y, float z, float t,
+ // const edm4hep::SimTrackerHit simtrackerhit,
+ // Acts::GeometryIdentifier geometryId,
+ boost::container::static_vector<Acts::SourceLink, 2> sourceLinks)
+ : m_x(x), m_y(y), m_z(z), m_t(t),
+ m_sourceLinks(std::move(sourceLinks))
+ // , m_geometryId(geometryId)
+ {
+ // m_trackerHit = trackerhit;
+ // m_simtrackerHit = simtrackerhit;
+ // m_x = simtrackerhit.getPosition()[0];
+ // m_y = simtrackerhit.getPosition()[1];
+ // m_z = simtrackerhit.getPosition()[2];
+ // m_t = simtrackerhit.getTime();
+ // m_cellid = simtrackerhit.getCellID();
+ m_rho = std::sqrt(m_x * m_x + m_y * m_y + m_z * m_z);
+ m_varianceRho = trackerhit.getCovMatrix()[0];
+ m_varianceZ = trackerhit.getCovMatrix()[5];
+ }
+
+ // edm4hep::TrackerHit getTrackerHit() { return m_trackerHit; }
+ // edm4hep::SimTrackerHit getSimTrackerHit() { return m_simtrackerHit; }
+
+ constexpr Scalar x() const { return m_x; }
+ constexpr Scalar y() const { return m_y; }
+ constexpr Scalar z() const { return m_z; }
+ constexpr std::optional<Scalar> t() const { return m_t; }
+ constexpr Scalar r() const { return m_rho; }
+ constexpr Scalar varianceR() const { return m_varianceRho; }
+ constexpr Scalar varianceZ() const { return m_varianceZ; }
+ constexpr std::optional<Scalar> varianceT() const { return m_varianceT; }
+
+ // constexpr std::uint64_t cellid() const { return m_cellid; }
+ // constexpr Acts::GeometryIdentifier geometryId() const { return m_geometryId; }
+
+ const boost::container::static_vector<Acts::SourceLink, 2>&
+ sourceLinks() const { return m_sourceLinks; }
+
+ constexpr float topHalfStripLength() const { return m_topHalfStripLength; }
+ constexpr float bottomHalfStripLength() const { return m_bottomHalfStripLength; }
+
+ Acts::Vector3 topStripDirection() const { return m_topStripDirection; }
+ Acts::Vector3 bottomStripDirection() const { return m_bottomStripDirection; }
+ Acts::Vector3 stripCenterDistance() const { return m_stripCenterDistance; }
+ Acts::Vector3 topStripCenterPosition() const { return m_topStripCenterPosition; }
+
+ constexpr bool validDoubleMeasurementDetails() const {
+ return m_validDoubleMeasurementDetails;
+ }
+
+private:
+
+ // edm4hep::TrackerHit m_trackerHit;
+ // edm4hep::SimTrackerHit m_simtrackerHit;
+ // std::uint64_t m_cellid;
+ // Acts::GeometryIdentifier m_geometryId;
+
+ // Global position
+ Scalar m_x;
+ Scalar m_y;
+ Scalar m_z;
+ std::optional<Scalar> m_t;
+ Scalar m_rho;
+ // Variance in rho/z of the global coordinates
+ Scalar m_varianceRho;
+ Scalar m_varianceZ;
+ std::optional<Scalar> m_varianceT;
+ // SourceLinks of the corresponding measurements. A Pixel (strip) SP has one
+ // (two) sourceLink(s).
+ boost::container::static_vector<Acts::SourceLink, 2> m_sourceLinks;
+
+ // half of the length of the top strip
+ float m_topHalfStripLength = 0;
+ // half of the length of the bottom strip
+ float m_bottomHalfStripLength = 0;
+ // direction of the top strip
+ Acts::Vector3 m_topStripDirection = {0, 0, 0};
+ // direction of the bottom strip
+ Acts::Vector3 m_bottomStripDirection = {0, 0, 0};
+ // distance between the center of the two strips
+ Acts::Vector3 m_stripCenterDistance = {0, 0, 0};
+ // position of the center of the bottom strip
+ Acts::Vector3 m_topStripCenterPosition = {0, 0, 0};
+ bool m_validDoubleMeasurementDetails = false;
+
+}; // SimSpacePoint
+
+inline bool operator==(const SimSpacePoint& lhs, const SimSpacePoint& rhs) {
+ // (TODO) use sourceLinks for comparison
+ return ((lhs.x() == rhs.x()) && (lhs.y() == rhs.y()) &&
+ (lhs.z() == rhs.z()) && (lhs.t() == rhs.t()));
+}
+
+/// Container of space points.
+using SimSpacePointContainer = std::vector<SimSpacePoint>;
+using SimSeed = Acts::Seed<SimSpacePoint>;
+using SimSeedContainer = std::vector<Acts::Seed<SimSpacePoint>>;
+
+
+// --------------------------
+// IndexSourceLink
+// --------------------------
+
+/// A source link that stores just an index.
+///
+/// This is intentionally kept as barebones as possible. The source link
+/// is just a reference and will be copied, moved around, etc. often.
+/// Keeping it small and separate from the actual, potentially large,
+/// measurement data should result in better overall performance.
+/// Using an index instead of e.g. a pointer, means source link and
+/// measurement are decoupled and the measurement representation can be
+/// easily changed without having to also change the source link.
+class IndexSourceLink final
+{
+ public:
+ /// Construct from geometry identifier and index.
+ IndexSourceLink(Acts::GeometryIdentifier gid, Index idx, edm4hep::TrackerHit trackhit)
+ : m_geometryId(gid), m_index(idx), m_trackhit(trackhit) {}
+
+ // Construct an invalid source link.
+ // Must be default constructible to satisfy SourceLinkConcept.
+ IndexSourceLink() = default;
+ IndexSourceLink(const IndexSourceLink&) = default;
+ IndexSourceLink(IndexSourceLink&&) = default;
+ IndexSourceLink& operator=(const IndexSourceLink&) = default;
+ IndexSourceLink& operator=(IndexSourceLink&&) = default;
+
+ /// Access the index.
+ Index index() const { return m_index; }
+ Acts::GeometryIdentifier geometryId() const { return m_geometryId; }
+ edm4hep::TrackerHit getTrackerHit() const { return m_trackhit; }
+
+ struct SurfaceAccessor
+ {
+ const Acts::TrackingGeometry& trackingGeometry;
+ const Acts::Surface* operator()(const Acts::SourceLink& sourceLink) const
+ {
+ const auto& indexSourceLink = sourceLink.get<IndexSourceLink>();
+ return trackingGeometry.findSurface(indexSourceLink.geometryId());
+ }
+ };
+
+ private:
+ Acts::GeometryIdentifier m_geometryId;
+ Index m_index = 0;
+ edm4hep::TrackerHit m_trackhit;
+
+ friend bool operator==(const IndexSourceLink& lhs, const IndexSourceLink& rhs)
+ {
+ return (lhs.geometryId() == rhs.geometryId()) && (lhs.m_index == rhs.m_index);
+ }
+
+ friend bool operator!=(const IndexSourceLink& lhs, const IndexSourceLink& rhs)
+ {
+ return !(lhs == rhs);
+ }
+}; // IndexSourceLink
+
+/// Container of index source links.
+///
+/// Since the source links provide a `.geometryId()` accessor,
+/// they can be stored in an ordered geometry container.
+using IndexSourceLinkContainer = GeometryIdMultiset<IndexSourceLink>;
+
+/// Accessor for the above source link container
+///
+/// It wraps up a few lookup methods to be used in the Combinatorial Kalman
+/// Filter
+struct IndexSourceLinkAccessor : GeometryIdMultisetAccessor<IndexSourceLink>
+{
+ using BaseIterator = GeometryIdMultisetAccessor<IndexSourceLink>::Iterator;
+
+ using Iterator = Acts::SourceLinkAdapterIterator<BaseIterator>;
+
+ // get the range of elements with requested geoId
+ std::pair<Iterator, Iterator> range(const Acts::Surface& surface) const
+ {
+ assert(container != nullptr);
+ auto [begin, end] = container->equal_range(surface.geometryId());
+ return {Iterator{begin}, Iterator{end}};
+ }
+};
\ No newline at end of file
diff --git a/Reconstruction/RecActsTracking/src/utils/TGeoDetector.hpp b/Reconstruction/RecActsTracking/src/utils/TGeoDetector.hpp
new file mode 100644
index 00000000..4163d23b
--- /dev/null
+++ b/Reconstruction/RecActsTracking/src/utils/TGeoDetector.hpp
@@ -0,0 +1,596 @@
+#include "Acts/Geometry/CylinderVolumeBuilder.hpp"
+#include "Acts/Geometry/CylinderVolumeHelper.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Geometry/ITrackingVolumeBuilder.hpp"
+#include "Acts/Geometry/LayerArrayCreator.hpp"
+#include "Acts/Geometry/LayerCreator.hpp"
+#include "Acts/Geometry/PassiveLayerBuilder.hpp"
+#include "Acts/Geometry/ProtoLayerHelper.hpp"
+#include "Acts/Geometry/SurfaceArrayCreator.hpp"
+#include "Acts/Geometry/SurfaceBinningMatcher.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Geometry/TrackingGeometryBuilder.hpp"
+#include "Acts/Geometry/TrackingVolumeArrayCreator.hpp"
+
+#include "Acts/Plugins/TGeo/TGeoCylinderDiscSplitter.hpp"
+#include "Acts/Plugins/TGeo/TGeoLayerBuilder.hpp"
+#include "Acts/Plugins/Json/JsonMaterialDecorator.hpp"
+#include "Acts/Plugins/Json/ActsJson.hpp"
+
+#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/Logger.hpp"
+
+#include <boost/program_options.hpp>
+#include <nlohmann/json.hpp>
+
+#include <cstddef>
+#include <map>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Acts
+{
+ class TGeoDetectorElement;
+ class TrackingGeometry;
+ class IMaterialDecorator;
+}
+
+/// ----------------------------
+/// pre definitions
+/// ----------------------------
+
+/// Half open [lower,upper) interval type for a single user option.
+///
+/// A missing limit represents an unbounded upper or lower limit. With just
+/// one defined limit the interval is just a lower/upper bound; with both
+/// limits undefined, the interval is unbounded everywhere and thus contains
+/// all possible values.
+struct Interval
+{
+ std::optional<double> lower;
+ std::optional<double> upper;
+};
+
+/// Extract an interval from an input of the form 'lower:upper'.
+///
+/// An input of the form `lower:` or `:upper` sets just one of the limits. Any
+/// other input leads to an unbounded interval.
+///
+/// @note The more common range notation uses `lower-upper` but the `-`
+/// separator complicates the parsing of negative values.
+std::istream& operator>>(std::istream& is, Interval& interval);
+
+/// Print an interval as `lower:upper`.
+std::ostream& operator<<(std::ostream& os, const Interval& interval);
+
+struct TGeoConfig {
+ Acts::Logging::Level surfaceLogLevel = Acts::Logging::WARNING;
+ Acts::Logging::Level layerLogLevel = Acts::Logging::WARNING;
+ Acts::Logging::Level volumeLogLevel = Acts::Logging::WARNING;
+
+ std::string fileName;
+ bool buildBeamPipe = false;
+ double beamPipeRadius{0};
+ double beamPipeHalflengthZ{0};
+ double beamPipeLayerThickness{0};
+ double beamPipeEnvelopeR{1.0};
+ double layerEnvelopeR{1.0};
+
+ double unitScalor = 1.0;
+
+ Acts::TGeoLayerBuilder::ElementFactory elementFactory =
+ Acts::TGeoLayerBuilder::defaultElementFactory;
+
+ /// Optional geometry identifier hook to be used during closure
+ std::shared_ptr<const Acts::GeometryIdentifierHook> geometryIdentifierHook =
+ std::make_shared<Acts::GeometryIdentifierHook>();
+
+ enum SubVolume : std::size_t { Negative = 0, Central, Positive };
+
+ template <typename T>
+ struct LayerTriplet
+ {
+ LayerTriplet() = default;
+
+ LayerTriplet(T value)
+ : negative{value}, central{value}, positive{value} {}
+
+ LayerTriplet(T _negative, T _central, T _positive)
+ : negative{_negative}, central{_central}, positive{_positive} {}
+
+ T negative;
+ T central;
+ T positive;
+
+ T& at(SubVolume i)
+ {
+ switch (i)
+ {
+ case Negative: return negative;
+ case Central: return central;
+ case Positive: return positive;
+ default: throw std::invalid_argument{"Unknown index"};
+ }
+ }
+
+ const T& at(SubVolume i) const
+ {
+ switch (i)
+ {
+ case Negative: return negative;
+ case Central: return central;
+ case Positive: return positive;
+ default: throw std::invalid_argument{"Unknown index"};
+ }
+ }
+ };
+
+ struct Volume {
+ std::string name;
+ LayerTriplet<bool> layers{false};
+ LayerTriplet<std::string> subVolumeName;
+ LayerTriplet<std::vector<std::string>> sensitiveNames;
+ LayerTriplet<std::string> sensitiveAxes;
+ LayerTriplet<Interval> rRange;
+ LayerTriplet<Interval> zRange;
+ LayerTriplet<double> splitTolR{0};
+ LayerTriplet<double> splitTolZ{0};
+ LayerTriplet<std::vector<std::pair<int, Acts::BinningType>>> binning0;
+ LayerTriplet<std::vector<std::pair<int, Acts::BinningType>>> binning1;
+
+ Interval binToleranceR;
+ Interval binTolerancePhi;
+ Interval binToleranceZ;
+
+ bool cylinderDiscSplit = false;
+ unsigned int cylinderNZSegments = 0;
+ unsigned int cylinderNPhiSegments = 0;
+ unsigned int discNRSegments = 0;
+ unsigned int discNPhiSegments = 0;
+
+ bool itkModuleSplit = false;
+ std::map<std::string, unsigned int> barrelMap;
+ std::map<std::string, std::vector<std::pair<double, double>>> discMap;
+ /// pairs of regular expressions to match sensor names and category keys
+ /// for either the barrelMap or the discMap
+ /// @TODO in principle vector<pair< > > would be good enough
+ std::map<std::string, std::string> splitPatterns;
+ };
+
+ std::vector<Volume> volumes;
+};
+
+/// ------------------------------------------------------------------------------------------------
+/// @note nlohmann::json must define functions *from_json* && *to_json* for json conversion
+/// ------------------------------------------------------------------------------------------------
+
+namespace Acts {
+/// Read & Write config for cylinder/disc module splitter
+void from_json(const nlohmann::json& j,
+ Acts::TGeoCylinderDiscSplitter::Config& cdc)
+{
+ /// Number of segments in phi for a disc
+ cdc.cylinderPhiSegments = j.at("geo-tgeo-cyl-nphi-segs");
+ /// Number of segments in r for a disk
+ cdc.cylinderLongitudinalSegments = j.at("geo-tgeo-cyl-nz-segs");
+ /// Number of segments in phi for a disc
+ cdc.discPhiSegments = j.at("geo-tgeo-disc-nphi-segs");
+ /// Number of segments in r for a disk
+ cdc.discRadialSegments = j.at("geo-tgeo-disc-nr-segs");
+}
+
+void to_json(nlohmann::json& j,
+ const Acts::TGeoCylinderDiscSplitter::Config& cdc)
+{
+ j = nlohmann::json{{"geo-tgeo-cyl-nphi-segs", cdc.cylinderPhiSegments},
+ {"geo-tgeo-cyl-nz-segs", cdc.cylinderLongitudinalSegments},
+ {"geo-tgeo-disc-nphi-segs", cdc.discPhiSegments},
+ {"geo-tgeo-disc-nr-segs", cdc.discRadialSegments}};
+}
+
+// enum specialization by nlohman library
+NLOHMANN_JSON_SERIALIZE_ENUM(Acts::BinningType,
+ {
+ {Acts::BinningType::equidistant, "equidistant"},
+ {Acts::BinningType::arbitrary, "arbitrary"},
+ })
+}
+
+/// Read & Write config for options interval
+void from_json(const nlohmann::json& j,
+ Interval& interval)
+{
+ interval.lower = j.at("lower");
+ interval.upper = j.at("upper");
+}
+
+void to_json(nlohmann::json& j,
+ const Interval& interval)
+{
+ j = nlohmann::json{{"lower", interval.lower.value_or(0)},
+ {"upper", interval.upper.value_or(0)}};
+}
+
+/// Read & Write layer configuration triplets
+template <typename T>
+void from_json(const nlohmann::json& j,
+ TGeoConfig::LayerTriplet<T>& ltr)
+{
+ ltr.negative = j.at("negative").get<T>();
+ ltr.central = j.at("central").get<T>();
+ ltr.positive = j.at("positive").get<T>();
+}
+
+template <typename T>
+void to_json(nlohmann::json& j,
+ const TGeoConfig::LayerTriplet<T>& ltr)
+{
+ j = nlohmann::json{{"negative", ltr.negative},
+ {"central", ltr.central},
+ {"positive", ltr.positive}};
+}
+
+/// Read & Write volume struct
+void from_json(const nlohmann::json& j, TGeoConfig::Volume& vol) {
+ // subdetector selection
+ vol.name = j.at("geo-tgeo-volume-name");
+
+ // configure surface autobinning
+ vol.binToleranceR = j.at("geo-tgeo-sfbin-r-tolerance");
+ vol.binToleranceZ = j.at("geo-tgeo-sfbin-z-tolerance");
+ vol.binTolerancePhi = j.at("geo-tgeo-sfbin-phi-tolerance");
+
+ // Fill layer triplets
+ vol.layers = j.at("geo-tgeo-volume-layers");
+ vol.subVolumeName = j.at("geo-tgeo-subvolume-names");
+ vol.sensitiveNames = j.at("geo-tgeo-sensitive-names");
+ vol.sensitiveAxes = j.at("geo-tgeo-sensitive-axes");
+ vol.rRange = j.at("geo-tgeo-layer-r-ranges");
+ vol.zRange = j.at("geo-tgeo-layer-z-ranges");
+ vol.splitTolR = j.at("geo-tgeo-layer-r-split");
+ vol.splitTolZ = j.at("geo-tgeo-layer-z-split");
+ // Set binning manually
+ vol.binning0 = j.at("geo-tgeo-binning0");
+ vol.binning1 = j.at("geo-tgeo-binning1");
+
+ vol.cylinderDiscSplit = j.at("geo-tgeo-cyl-disc-split");
+ if (vol.cylinderDiscSplit)
+ {
+ Acts::TGeoCylinderDiscSplitter::Config cdConfig = j.at("Splitters").at("CylinderDisk");
+ vol.cylinderNZSegments = cdConfig.cylinderLongitudinalSegments;
+ vol.cylinderNPhiSegments = cdConfig.cylinderPhiSegments;
+ vol.discNRSegments = cdConfig.discRadialSegments;
+ vol.discNPhiSegments = cdConfig.discPhiSegments;
+ }
+ vol.itkModuleSplit = false;
+}
+
+void to_json(nlohmann::json& j, const TGeoConfig::Volume& vol) {
+ j["geo-tgeo-volume-name"] = vol.name;
+
+ j["geo-tgeo-sfbin-r-tolerance"] = vol.binToleranceR;
+ j["geo-tgeo-sfbin-z-tolerance"] = vol.binToleranceZ;
+ j["geo-tgeo-sfbin-phi-tolerance"] = vol.binTolerancePhi;
+
+ j["geo-tgeo-volume-layers"] = vol.layers;
+ j["geo-tgeo-subvolume-names"] = vol.subVolumeName;
+ j["geo-tgeo-sensitive-names"] = vol.sensitiveNames;
+ j["geo-tgeo-sensitive-axes"] = vol.sensitiveAxes;
+ j["geo-tgeo-layer-r-ranges"] = vol.rRange;
+ j["geo-tgeo-layer-z-ranges"] = vol.zRange;
+ j["geo-tgeo-layer-r-split"] = vol.splitTolR;
+ j["geo-tgeo-layer-z-split"] = vol.splitTolZ;
+ j["geo-tgeo-binning0"] = vol.binning0;
+ j["geo-tgeo-binning1"] = vol.binning1;
+
+ j["geo-tgeo-cyl-disc-split"] = vol.cylinderDiscSplit;
+ j["geo-tgeo-itk-module-split"] = vol.itkModuleSplit;
+
+ Acts::TGeoCylinderDiscSplitter::Config cdConfig;
+ cdConfig.cylinderLongitudinalSegments = vol.cylinderNZSegments;
+ cdConfig.cylinderPhiSegments = vol.cylinderNPhiSegments;
+ cdConfig.discRadialSegments = vol.discNRSegments;
+ cdConfig.discPhiSegments = vol.discNPhiSegments;
+ j["Splitters"]["CylinderDisk"] = cdConfig;
+}
+
+/// @brief Function that constructs a set of layer builder configs from a central @c TGeoDetector config.
+///
+/// @param config The input config
+/// @return Vector of layer builder configs
+std::vector<Acts::TGeoLayerBuilder::Config> makeLayerBuilderConfigs(
+ const TGeoConfig& config, const Acts::Logger& logger)
+{
+ std::vector<Acts::TGeoLayerBuilder::Config> detLayerConfigs;
+
+ // iterate over all configured detector volumes
+ for (const auto& volume : config.volumes)
+ {
+ Acts::TGeoLayerBuilder::Config layerBuilderConfig;
+ layerBuilderConfig.configurationName = volume.name;
+ layerBuilderConfig.unit = config.unitScalor;
+ layerBuilderConfig.elementFactory = config.elementFactory;
+
+ // configure surface autobinning
+ std::vector<std::pair<double, double>> binTolerances(
+ static_cast<std::size_t>(Acts::binValues), {0., 0.});
+ binTolerances[Acts::binR] = {volume.binToleranceR.lower.value_or(0.),
+ volume.binToleranceR.upper.value_or(0.)};
+ binTolerances[Acts::binZ] = {volume.binToleranceZ.lower.value_or(0.),
+ volume.binToleranceZ.upper.value_or(0.)};
+ binTolerances[Acts::binPhi] = {volume.binTolerancePhi.lower.value_or(0.),
+ volume.binTolerancePhi.upper.value_or(0.)};
+
+ layerBuilderConfig.autoSurfaceBinning = true;
+ layerBuilderConfig.surfaceBinMatcher = Acts::SurfaceBinningMatcher(binTolerances);
+
+ // loop over the negative/central/positive layer configurations
+ for (auto ncp : { TGeoConfig::Negative,
+ TGeoConfig::Central,
+ TGeoConfig::Positive,} )
+ {
+ if (!volume.layers.at(ncp)) { continue; }
+
+ Acts::TGeoLayerBuilder::LayerConfig lConfig;
+ lConfig.volumeName = volume.subVolumeName.at(ncp);
+ lConfig.sensorNames = volume.sensitiveNames.at(ncp);
+ lConfig.localAxes = volume.sensitiveAxes.at(ncp);
+ lConfig.envelope = {config.layerEnvelopeR, config.layerEnvelopeR};
+
+ auto rR = volume.rRange.at(ncp);
+ auto rMin = rR.lower.value_or(0.);
+ auto rMax = rR.upper.value_or(std::numeric_limits<double>::max());
+ auto zR = volume.zRange.at(ncp);
+ auto zMin = zR.lower.value_or(-std::numeric_limits<double>::max());
+ auto zMax = zR.upper.value_or(std::numeric_limits<double>::max());
+ lConfig.parseRanges =
+ {
+ {Acts::binR, {rMin, rMax}},
+ {Acts::binZ, {zMin, zMax}},
+ };
+
+ // Fill the layer splitting parameters in r/z
+ auto str = volume.splitTolR.at(ncp);
+ auto stz = volume.splitTolZ.at(ncp);
+ if (0 < str) { lConfig.splitConfigs.emplace_back(Acts::binR, str); }
+ if (0 < stz) { lConfig.splitConfigs.emplace_back(Acts::binZ, stz); }
+ lConfig.binning0 = volume.binning0.at(ncp);
+ lConfig.binning1 = volume.binning1.at(ncp);
+
+ layerBuilderConfig.layerConfigurations[ncp].push_back(lConfig);
+ }
+
+ // Perform splitting of cylinders and discs
+ if (volume.cylinderDiscSplit)
+ {
+ Acts::TGeoCylinderDiscSplitter::Config cdsConfig;
+ cdsConfig.cylinderPhiSegments = volume.cylinderNPhiSegments;
+ cdsConfig.cylinderLongitudinalSegments = volume.cylinderNZSegments;
+ cdsConfig.discPhiSegments = volume.discNPhiSegments;
+ cdsConfig.discRadialSegments = volume.discNRSegments;
+ layerBuilderConfig.detectorElementSplitter =
+ std::make_shared<const Acts::TGeoCylinderDiscSplitter>(cdsConfig,
+ logger.clone("TGeoCylinderDiscSplitter", config.layerLogLevel));
+ }
+
+ detLayerConfigs.push_back(layerBuilderConfig);
+ }
+
+ return detLayerConfigs;
+}
+
+/// @brief Function to build the generic tracking geometry from a TGeo object.
+///
+/// @param TGeo_ROOTFilePath is the TGeo ROOT file path
+/// @param TGeoConfig_jFilePath is the TGeo configuration file path
+/// @param MaterialMap_jFilePath is the material map file path
+/// @param logger is the logger object
+/// @return a shared pointer to the tracking geometry
+std::shared_ptr<const Acts::TrackingGeometry> buildTGeoDetector(
+ const Acts::GeometryContext& context,
+ std::vector<std::shared_ptr<const Acts::TGeoDetectorElement>>& detElementStore,
+ const std::string& TGeo_ROOTFilePath,
+ const std::string& TGeoConfig_jFilePath,
+ const std::string& MaterialMap_jFilePath,
+ const Acts::Logger& logger)
+{
+ if (TGeo_ROOTFilePath.empty() | TGeoConfig_jFilePath.empty() |
+ MaterialMap_jFilePath.empty()) { return nullptr; }
+
+ TGeoConfig config;
+ nlohmann::json djson;
+ const Acts::MaterialMapJsonConverter::Config m_converter;
+
+ // read input files
+ config.fileName = TGeo_ROOTFilePath;
+ std::shared_ptr<const Acts::IMaterialDecorator> mdecorator = std::make_shared<const Acts::JsonMaterialDecorator>(m_converter, MaterialMap_jFilePath, Acts::Logging::Level::INFO);
+ std::ifstream infile(TGeoConfig_jFilePath, std::ifstream::in | std::ifstream::binary);
+
+ // ------------------------------------------
+ // read TGeo Layer Builder Configs File
+ // ------------------------------------------
+ infile >> djson;
+ config.unitScalor = djson["geo-tgeo-unit-scalor"];
+ config.buildBeamPipe = djson["geo-tgeo-build-beampipe"];
+ if (config.buildBeamPipe)
+ {
+ const auto beamPipeParameters = djson["geo-tgeo-beampipe-parameters"].get<std::array<double, 3>>();
+ config.beamPipeRadius = beamPipeParameters[0];
+ config.beamPipeHalflengthZ = beamPipeParameters[1];
+ config.beamPipeLayerThickness = beamPipeParameters[2];
+ }
+
+ // Fill nested volume configs
+ for (const auto& volume : djson["Volumes"])
+ {
+ auto& vol = config.volumes.emplace_back();
+ vol = volume;
+ }
+
+ // ------------------------------------------
+ // logger Configaration
+ // ------------------------------------------
+ Acts::SurfaceArrayCreator::Config sacConfig;
+ auto surfaceArrayCreator = std::make_shared<const Acts::SurfaceArrayCreator>(
+ sacConfig, logger.clone("SurfaceArrayCreator", config.surfaceLogLevel));
+ // configure the proto layer helper
+ Acts::ProtoLayerHelper::Config plhConfig;
+ auto protoLayerHelper = std::make_shared<const Acts::ProtoLayerHelper>(
+ plhConfig, logger.clone("ProtoLayerHelper", config.layerLogLevel));
+ // configure the layer creator that uses the surface array creator
+ Acts::LayerCreator::Config lcConfig;
+ lcConfig.surfaceArrayCreator = surfaceArrayCreator;
+ auto layerCreator = std::make_shared<const Acts::LayerCreator>(
+ lcConfig, logger.clone("LayerCreator", config.layerLogLevel));
+ // configure the layer array creator
+ Acts::LayerArrayCreator::Config lacConfig;
+ auto layerArrayCreator = std::make_shared<const Acts::LayerArrayCreator>(
+ lacConfig, logger.clone("LayerArrayCreator", config.layerLogLevel));
+ // tracking volume array creator
+ Acts::TrackingVolumeArrayCreator::Config tvacConfig;
+ auto tVolumeArrayCreator =
+ std::make_shared<const Acts::TrackingVolumeArrayCreator>( tvacConfig,
+ logger.clone("TrackingVolumeArrayCreator", config.volumeLogLevel));
+
+ // configure the cylinder volume helper
+ Acts::CylinderVolumeHelper::Config cvhConfig;
+ cvhConfig.layerArrayCreator = layerArrayCreator;
+ cvhConfig.trackingVolumeArrayCreator = tVolumeArrayCreator;
+ auto cylinderVolumeHelper =
+ std::make_shared<const Acts::CylinderVolumeHelper>(cvhConfig,
+ logger.clone("CylinderVolumeHelper", config.volumeLogLevel));
+
+ // ------------------------------------------
+ // logger Configaration
+ // ------------------------------------------
+ // list the volume builders
+ std::list<std::shared_ptr<const Acts::ITrackingVolumeBuilder>> volumeBuilders;
+
+ // Create a beam pipe if configured to do so
+ if (config.buildBeamPipe)
+ {
+ /// configure the beam pipe layer builder
+ Acts::PassiveLayerBuilder::Config bplConfig;
+ bplConfig.layerIdentification = "BeamPipe";
+ bplConfig.centralLayerRadii = {config.beamPipeRadius};
+ bplConfig.centralLayerHalflengthZ = {config.beamPipeHalflengthZ};
+ bplConfig.centralLayerThickness = {config.beamPipeLayerThickness};
+ auto beamPipeBuilder = std::make_shared<const Acts::PassiveLayerBuilder>(
+ bplConfig, logger.clone("BeamPipeLayerBuilder", config.layerLogLevel));
+ // create the volume for the beam pipe
+ Acts::CylinderVolumeBuilder::Config bpvConfig;
+ bpvConfig.trackingVolumeHelper = cylinderVolumeHelper;
+ bpvConfig.volumeName = "BeamPipe";
+ bpvConfig.layerBuilder = beamPipeBuilder;
+ bpvConfig.layerEnvelopeR = {config.beamPipeEnvelopeR, config.beamPipeEnvelopeR};
+ bpvConfig.buildToRadiusZero = true;
+ auto beamPipeVolumeBuilder =
+ std::make_shared<const Acts::CylinderVolumeBuilder>(bpvConfig,
+ logger.clone("BeamPipeVolumeBuilder", config.volumeLogLevel));
+ // add to the list of builders
+ volumeBuilders.push_back(beamPipeVolumeBuilder);
+ }
+
+ // Import the file from
+ TGeoManager::Import(config.fileName.c_str());
+ auto layerBuilderConfigs = makeLayerBuilderConfigs(config, logger);
+ // Remember the layer builders to collect the detector elements
+ std::vector<std::shared_ptr<const Acts::TGeoLayerBuilder>> tgLayerBuilders;
+ for (auto& lbc : layerBuilderConfigs)
+ {
+ std::shared_ptr<const Acts::LayerCreator> layerCreatorLB = nullptr;
+ if (lbc.autoSurfaceBinning)
+ {
+ // Configure surface array creator (optionally) per layer builder
+ // (in order to configure them to work appropriately)
+ Acts::SurfaceArrayCreator::Config sacConfigLB;
+ sacConfigLB.surfaceMatcher = lbc.surfaceBinMatcher;
+ auto surfaceArrayCreatorLB =
+ std::make_shared<const Acts::SurfaceArrayCreator>(sacConfigLB,
+ logger.clone(lbc.configurationName + "SurfaceArrayCreator", config.surfaceLogLevel));
+ // configure the layer creator that uses the surface array creator
+ Acts::LayerCreator::Config lcConfigLB;
+ lcConfigLB.surfaceArrayCreator = surfaceArrayCreatorLB;
+ layerCreatorLB = std::make_shared<const Acts::LayerCreator>(lcConfigLB,
+ logger.clone(lbc.configurationName + "LayerCreator", config.layerLogLevel));
+ }
+
+ // Configure the proto layer helper
+ Acts::ProtoLayerHelper::Config plhConfigLB;
+ auto protoLayerHelperLB = std::make_shared<const Acts::ProtoLayerHelper>( plhConfigLB,
+ logger.clone(lbc.configurationName + "ProtoLayerHelper", config.layerLogLevel));
+
+ lbc.layerCreator =
+ (layerCreatorLB != nullptr) ? layerCreatorLB : layerCreator;
+ lbc.protoLayerHelper =
+ (protoLayerHelperLB != nullptr) ? protoLayerHelperLB : protoLayerHelper;
+
+ auto layerBuilder = std::make_shared<const Acts::TGeoLayerBuilder>( lbc,
+ logger.clone(lbc.configurationName + "LayerBuilder", config.layerLogLevel));
+
+ // remember the layer builder
+ tgLayerBuilders.push_back(layerBuilder);
+
+ // build the pixel volume
+ Acts::CylinderVolumeBuilder::Config volumeConfig;
+ volumeConfig.trackingVolumeHelper = cylinderVolumeHelper;
+ volumeConfig.volumeName = lbc.configurationName;
+ volumeConfig.buildToRadiusZero = volumeBuilders.empty();
+ volumeConfig.layerEnvelopeR = {config.layerEnvelopeR, config.layerEnvelopeR};
+ auto ringLayoutConfiguration =
+ [&](const std::vector<Acts::TGeoLayerBuilder::LayerConfig>& lConfigs) -> void
+ {
+ for (const auto& lcfg : lConfigs)
+ {
+ for (const auto& scfg : lcfg.splitConfigs)
+ {
+ if (scfg.first == Acts::binR && scfg.second > 0.)
+ {
+ volumeConfig.ringTolerance = std::max(volumeConfig.ringTolerance, scfg.second);
+ volumeConfig.checkRingLayout = true;
+ }
+ }
+ }
+ };
+ ringLayoutConfiguration(lbc.layerConfigurations[0]);
+ ringLayoutConfiguration(lbc.layerConfigurations[2]);
+ volumeConfig.layerBuilder = layerBuilder;
+ auto volumeBuilder = std::make_shared<const Acts::CylinderVolumeBuilder>( volumeConfig,
+ logger.clone(lbc.configurationName + "VolumeBuilder", config.volumeLogLevel));
+ // add to the list of builders
+ volumeBuilders.push_back(volumeBuilder);
+ }
+
+ //-------------------------------------------------------------------------------------
+ // create the tracking geometry
+ Acts::TrackingGeometryBuilder::Config tgConfig;
+ // Add the builders
+ tgConfig.materialDecorator = std::move(mdecorator);
+ tgConfig.geometryIdentifierHook = config.geometryIdentifierHook;
+
+ for (auto& vb : volumeBuilders)
+ {
+ tgConfig.trackingVolumeBuilders.push_back(
+ [=](const auto& gcontext, const auto& inner, const auto&)
+ { return vb->trackingVolume(gcontext, inner); });
+ }
+ // Add the helper
+ tgConfig.trackingVolumeHelper = cylinderVolumeHelper;
+ auto cylinderGeometryBuilder = std::make_shared<const Acts::TrackingGeometryBuilder>(
+ tgConfig, logger.clone("TrackerGeometryBuilder", config.volumeLogLevel));
+ // get the geometry
+ auto trackingGeometry = cylinderGeometryBuilder->trackingGeometry(context);
+ // collect the detector element store
+ for (auto& lBuilder : tgLayerBuilders)
+ {
+ auto detElements = lBuilder->detectorElements();
+ detElementStore.insert(detElementStore.begin(), detElements.begin(), detElements.end());
+ }
+
+ /// return the tracking geometry
+ return trackingGeometry;
+}
\ No newline at end of file
--
GitLab