diff --git a/.github/workflows/downstream-build.yaml b/.github/workflows/downstream-build.yaml index f46a6646..fddbec59 100644 --- a/.github/workflows/downstream-build.yaml +++ b/.github/workflows/downstream-build.yaml @@ -17,5 +17,5 @@ jobs: strategy: fail-fast: false steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: key4hep/key4hep-actions/downstream-build@main diff --git a/.github/workflows/mucoll-ci.yml b/.github/workflows/mucoll-ci.yml index c33ad962..7912a3ca 100644 --- a/.github/workflows/mucoll-ci.yml +++ b/.github/workflows/mucoll-ci.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest steps: - name: Free Disk Space (Ubuntu) - uses: endersonmenezes/free-disk-space@v2 + uses: endersonmenezes/free-disk-space@v3 with: remove_android: true remove_dotnet: true @@ -25,12 +25,14 @@ jobs: remove_swap: true remove_packages: "azure-cli google-cloud-cli microsoft-edge-stable google-chrome-stable firefox postgresql* temurin-* *llvm* mysql* dotnet-sdk-*" remove_packages_one_command: true - remove_folders: "/usr/share/swift /usr/share/miniconda /usr/share/az* /usr/share/glade* /usr/local/lib/node_modules /usr/local/share/chromium /usr/local/share/powershell" + remove_folders: "/usr/share/swift /usr/share/miniconda /usr/share/az* /usr/local/lib/node_modules /usr/local/share/chromium /usr/local/share/powershell /usr/local/julia /usr/local/aws-cli /usr/local/aws-sam-cli /usr/share/gradle" + rm_cmd: "rmz" # Use 'rmz' for faster deletion (default: 'rm') + rmz_version: "3.1.1" # Required when rm_cmd is 'rmz' testing: false - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - shell: bash run: echo "NOW=$(date +'%Y-%m-%d')" >> $GITHUB_ENV - - uses: actions/cache/restore@v4 + - uses: actions/cache/restore@v5 with: path: ~/.cache/ccache key: ccache-mucoll-image-${{ env.NOW }} @@ -59,8 +61,7 @@ jobs: -DCMAKE_INSTALL_PREFIX=$(pwd)/install \ -DCMAKE_CXX_FLAGS=" -fdiagnostics-color=always -Werror -Wno-error=deprecated-declarations" cmake --build build - # Tests need OpenDataDetector, which is not present in the mucoll image atm - # ctest --test-dir build -j$(nproc) --output-on-failure + ctest --test-dir build -j$(nproc) --output-on-failure -E "(CLD|ILD)" cmake --build build --target install ccache -s @@ -69,7 +70,7 @@ jobs: chmod +x ${GITHUB_WORKSPACE}/build.sh docker exec ci_container /bin/bash -c "/packages/build.sh" - - uses: actions/cache/save@v4 + - uses: actions/cache/save@v5 if: always() with: path: ~/.cache/ccache diff --git a/.gitignore b/.gitignore index d988e494..0c80bde2 100644 --- a/.gitignore +++ b/.gitignore @@ -251,3 +251,10 @@ test/gaudi_opts/testConverterConstants.py .claude/ .codex/ .aider.* + +# OS junk +.DS_Store + +# Detector/material files fetched by data/download_files.sh +/data/*.root +/data/*.json diff --git a/README.md b/README.md index 9e0baf9e..a223e264 100644 --- a/README.md +++ b/README.md @@ -1,25 +1,112 @@ # k4ActsTracking +[![Key4hep build](https://github.com/key4hep/k4ActsTracking/actions/workflows/key4hep-build.yaml/badge.svg)](https://github.com/key4hep/k4ActsTracking/actions/workflows/key4hep-build.yaml) +[![downstream-build](https://github.com/key4hep/k4ActsTracking/actions/workflows/downstream-build.yaml/badge.svg)](https://github.com/key4hep/k4ActsTracking/actions/workflows/downstream-build.yaml) +[![MuColl Image Build](https://github.com/key4hep/k4ActsTracking/actions/workflows/mucoll-ci.yml/badge.svg)](https://github.com/key4hep/k4ActsTracking/actions/workflows/mucoll-ci.yml) +[![pre-commit](https://github.com/key4hep/k4ActsTracking/actions/workflows/pre-commit.yml/badge.svg)](https://github.com/key4hep/k4ActsTracking/actions/workflows/pre-commit.yml) -This repository contains the necessary tools to use ACTS functionality in Key4hep - +`k4ActsTracking` provides [ACTS](https://acts.readthedocs.io/)-based track +reconstruction components for the [Key4hep](https://key4hep.github.io/key4hep-doc/) +software stack. It converts DD4hep tracking geometries into ACTS geometries and +exposes a set of Gaudi algorithms and services (geometry conversion, seeding, +combinatorial Kalman filter track finding, duplicate removal, track filtering +and truth matching) usable from `k4run` option files. ## Dependencies -* Acts +* [ACTS](https://github.com/acts-project/acts) (`Core`, `PluginDD4hep`, + `PluginJson`, `PluginRoot`) +* [DD4hep](https://github.com/AIDASoft/DD4hep) (`DDCore`, `DDRec`) +* [EDM4hep](https://github.com/key4hep/EDM4hep) +* [k4FWCore](https://github.com/key4hep/k4FWCore) +* [Gaudi](https://gitlab.cern.ch/gaudi/Gaudi) +* TBB -* DD4hep +All of these are provided by a Key4hep release. The easiest way to get a working +environment is to source a nightly or stable Key4hep stack, e.g.: -* k4FWCore +```sh +source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh +``` -## Installation +## Building +```sh +mkdir -p build install +cmake -B build -S . -GNinja \ + -DCMAKE_CXX_STANDARD=20 \ + -DCMAKE_INSTALL_PREFIX=$(pwd)/install \ + -DCMAKE_CXX_FLAGS=" -fdiagnostics-color=always -Werror -Wno-error=deprecated-declarations" +cmake --build build +cmake --build build --target install +``` +After installing, make the package visible to Gaudi/`k4run`: + +```sh +source ../install/setup.sh # if generated, otherwise set the paths below +# or, manually: +export LD_LIBRARY_PATH=$PWD/install/lib:$LD_LIBRARY_PATH +export PYTHONPATH=$PWD/install/python:$PYTHONPATH ``` -mkdir build install -cd build; -cmake .. -DCMAKE_INSTALL_PREFIX=../install -make install +## Components + +The Gaudi plugin module `k4ActsTrackingPlugins` provides, among others: + +* **`ActsGeoSvc`** — builds an ACTS tracking geometry from a DD4hep compact + file (consumed via `GeoSvc`). Can optionally dump the converted geometry to an + `.obj` file for visualization. +* **`CKFTrackingAlg`** — seeding plus combinatorial Kalman filter (CKF) track + finding on EDM4hep tracker hits. +* **`ACTSSeededCKFTrackingAlg`** — alternative seeded CKF tracking algorithm. +* **`ACTSDuplicateRemoval`** — removes duplicate tracks produced by the CKF. +* **`FilterTracksAlg`** — applies quality cuts to a track collection. +* **`TrackTruthAlg`** — associates reconstructed tracks with truth particles. +* **`ActsTestPropagator`** — propagates ACTS particle-gun tracks through the + converted geometry (useful for geometry validation). + +## Usage + +The algorithms are configured and run through `k4run` option files. See the +examples and the test option files for working configurations: + +* [`k4ActsTracking/examples/test_visualize_acts_geo.py`](k4ActsTracking/examples/test_visualize_acts_geo.py) + — load a compact file, convert it to ACTS geometry, optionally dump an `.obj` + and run the test propagator: + ```sh + k4run k4ActsTracking/examples/test_visualize_acts_geo.py \ + --compactFile .xml --test-propagation + ``` + +* [`test/options/MAIA_CKFTrackingAlg.py`](test/options/MAIA_CKFTrackingAlg.py), + [`test/options/CLD_CKFTracking.py`](test/options/CLD_CKFTracking.py), + [`test/options/ILD_CKFTracking.py`](test/options/ILD_CKFTracking.py) — full + digitization + CKF tracking chains for the MAIA, CLD and ILD detectors. The + shared helpers live in + [`test/options/_ckf_helpers.py`](test/options/_ckf_helpers.py). + + > **Note:** the parameters in these option files are tuned only for technical + > tests and are *not* a meaningful physics tracking configuration. + +## Tests + +Tests are built when `BUILD_TESTING` is on (default) and run with `ctest` from +the build directory: + +```sh +cd build +ctest --output-on-failure ``` + +The suite covers C++ unit tests plus end-to-end chains (`ddsim` simulation → +reconstruction → CKF tracking) for the MAIA, CLD and ILD geometries. The CLD and +ILD chains clone the corresponding upstream config repositories +([`CLDConfig`](https://github.com/key4hep/CLDConfig), +[`ILDConfig`](https://github.com/iLCSoft/ILDConfig)) on the fly, and geometries +are taken from `k4geo` (via `$k4geo_DIR`). + +## License + +Licensed under the Apache License, Version 2.0. See [LICENSE](LICENSE). diff --git a/k4ActsTracking/CMakeLists.txt b/k4ActsTracking/CMakeLists.txt index 9bb0f90a..4acc93b7 100644 --- a/k4ActsTracking/CMakeLists.txt +++ b/k4ActsTracking/CMakeLists.txt @@ -21,7 +21,9 @@ set(ACTS_LIBRARY_TARGETS "Acts::Core") gaudi_add_library(k4ActsTracking SOURCES src/ActsGaudiLogger.cpp + src/CellIDSelector.cxx LINK + DD4hep::DDCore k4FWCore::k4FWCore ${ACTS_LIBRARY_TARGETS} ) @@ -43,6 +45,7 @@ set(_plugin_sources src/components/TrackTruthAlg.cxx src/components/ActsTestPropagator.cpp src/components/DD4hepBlueprintConstruction.cpp + src/components/CKFTrackingAlg.cpp ) gaudi_add_module(k4ActsTrackingPlugins @@ -69,12 +72,4 @@ install(TARGETS k4ActsTrackingPlugins COMPONENT dev) set(GAUDI_GENCONF_DIR "genConfDir") - -function(set_test_env _testname) - set_property(TEST ${_testname} APPEND PROPERTY ENVIRONMENT - LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$:$:$:$:$:$ENV{LD_LIBRARY_PATH} - PYTHONPATH=${CMAKE_BINARY_DIR}/${CMAKE_PROJECT_NAME}/${GAUDI_GENCONF_DIR}:$/../python:$ENV{PYTHONPATH} - PATH=$/../bin:$ENV{PATH} - K4ACTSTRACKING=${CMAKE_CURRENT_LIST_DIR}/ - ) -endfunction() +set_property(TARGET k4ActsTrackingPlugins PROPERTY GAUDI_GENCONF TRUE) diff --git a/k4ActsTracking/include/k4ActsTracking/ACTSSeededCKFTrackingAlg.hxx b/k4ActsTracking/include/k4ActsTracking/ACTSSeededCKFTrackingAlg.hxx index 2916c729..1651b900 100644 --- a/k4ActsTracking/include/k4ActsTracking/ACTSSeededCKFTrackingAlg.hxx +++ b/k4ActsTracking/include/k4ActsTracking/ACTSSeededCKFTrackingAlg.hxx @@ -203,9 +203,6 @@ protected: Gaudi::Property m_numThreads{this, "NumThreads", 1, "Number of threads to use for internal multithreading."}; ///@} - // Thread-safe counter - // mutable Gaudi::Accumulators::Counter<> m_fitFails{this, "FitFails"}; - private: // Mutexes for threadsafe container filling of seeds and tracks mutable std::mutex m_seedMutex{}; @@ -302,7 +299,7 @@ std::vector ACTSSeededCKFTrackingAlg::findSeeds( seedTrack.addToTrackStates(seedTrackState); } - debug() << "Seed Paramemeters" << std::endl << paramseed << endmsg; + debug() << "Seed Parameters" << std::endl << paramseed << endmsg; } debug() << "Seeds found: " << std::endl << paramseeds.size() << endmsg; diff --git a/k4ActsTracking/include/k4ActsTracking/CKFTrackingAlg.hxx b/k4ActsTracking/include/k4ActsTracking/CKFTrackingAlg.hxx new file mode 100644 index 00000000..03817753 --- /dev/null +++ b/k4ActsTracking/include/k4ActsTracking/CKFTrackingAlg.hxx @@ -0,0 +1,51 @@ +/* + * Copyright (c) 2014-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#pragma once + +// ACTS +#include + +// Standard +#include + +namespace ACTSTracking { + + //! Build a diagonal initial track covariance matrix from per-parameter error estimates + /** + * \param p Absolute momentum [ACTS units]. + * \param errPos Local-position error [ACTS units]. + * \param errPhi Azimuthal angle error [rad]. + * \param errLambda Polar angle (lambda) error [rad]. + * \param errRelP Relative momentum error (dimensionless fraction). + * \param errTime Time error [ACTS units]. + * \return 5×5 diagonal bound covariance matrix. + */ + inline Acts::BoundMatrix makeInitialCovariance(double p, double errPos, double errPhi, double errLambda, + double errRelP, double errTime) { + Acts::BoundMatrix cov = Acts::BoundMatrix::Zero(); + cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = std::pow(errPos, 2); + cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = std::pow(errPos, 2); + cov(Acts::eBoundTime, Acts::eBoundTime) = std::pow(errTime, 2); + cov(Acts::eBoundPhi, Acts::eBoundPhi) = std::pow(errPhi, 2); + cov(Acts::eBoundTheta, Acts::eBoundTheta) = std::pow(errLambda, 2); + cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = std::pow(errRelP * p / (p * p), 2); + return cov; + } + +} // namespace ACTSTracking diff --git a/k4ActsTracking/include/k4ActsTracking/CellIDSelector.h b/k4ActsTracking/include/k4ActsTracking/CellIDSelector.h new file mode 100644 index 00000000..c5cfc6cd --- /dev/null +++ b/k4ActsTracking/include/k4ActsTracking/CellIDSelector.h @@ -0,0 +1,125 @@ +/* + * Copyright (c) 2014-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#ifndef K4ACTSTRACKING_CELLIDSELECTOR_H +#define K4ACTSTRACKING_CELLIDSELECTOR_H + +#include +#include + +#include +#include + +namespace k4ActsTracking { + + /// Selects detector hits by matching their DD4hep CellID against a set of + /// user-defined selection strings. + /// + /// Each selection string describes a conjunction (AND) of field constraints. + /// Multiple selection strings passed to the constructor are evaluated as a + /// disjunction (OR): a CellID is accepted if it satisfies **at least one** of + /// them. + /// + /// ## Selection string grammar + /// + /// A selection string is a comma-separated list of `:` + /// constraints, where `` must be a field name present in the DD4hep + /// encoding string supplied to the constructor. + /// + /// ### Value syntax + /// + /// | Syntax | Meaning | + /// |-----------------|----------------------------------------------| + /// | `field:N` | Field must equal the integer `N` | + /// | `field:N\|M\|…` | Field must equal one of `N`, `M`, … (OR) | + /// + /// Fields that appear in the encoding string but are absent from the + /// selection are unconstrained. Specifying the same field more than + /// once in a single selection string is not detected as an error: only the + /// last occurrence takes effect. + /// + /// ### Examples + /// + /// ``` + /// // Accept hits in system 8, layer 3 (any module, sensor, side, …) + /// CellIDSelector sel(enc, {"system:8,layer:3"}); + /// + /// // Accept hits in system 5 on layers 1, 4, or 5 OR + /// // in system 3 on layers 2, 6, or 8 OR + /// // any hit whose sensor field equals 42 + /// CellIDSelector sel(enc, {"system:5,layer:1|4|5", + /// "system:3,layer:2|6|8", + /// "sensor:42"}); + /// + /// // Empty selection list — rejects every CellID + /// CellIDSelector sel(enc, {}); + /// ``` + /// + /// ### Expansion to (mask, value) pairs + /// + /// Internally each selection string is expanded into one or more + /// `Selector{mask, value}` pairs via the Cartesian product of the + /// per-field value lists. A CellID passes the selector when + /// `(cellID & mask) == (value & mask)`. Fields absent from the selection + /// contribute neither to `mask` nor to `value`. + class CellIDSelector { + public: + /// A single bitmask/value pair produced from one element of the Cartesian + /// product expansion of a selection string. + struct Selector { + dd4hep::CellID mask; ///< Bits that are relevant for this selector + dd4hep::CellID value; ///< Expected bit pattern after masking + }; + + CellIDSelector() = default; + CellIDSelector(const CellIDSelector&) = default; + CellIDSelector& operator=(const CellIDSelector&) = default; + CellIDSelector(CellIDSelector&&) = default; + CellIDSelector& operator=(CellIDSelector&&) = default; + ~CellIDSelector() = default; + + /// Construct a selector from a DD4hep encoding string and a list of + /// selection strings. + /// + /// @param encodingString DD4hep BitFieldCoder descriptor, e.g. + /// `"system:8,side:-2,layer:5,module:7,sensor:10"` + /// @param selections List of selection strings (OR-ed together). + /// An empty list rejects every CellID. + /// @throws std::invalid_argument if a selection string is malformed or + /// references an unknown field. + CellIDSelector(const std::string& encodingString, const std::vector& selections); + + /// Return `true` if @p cellID satisfies at least one configured selector. + bool accept(const dd4hep::CellID cellID) const; + + /// Expand a single selection string into its `Selector` pairs. + /// + /// Useful for inspection and testing. The returned vector contains one + /// entry per element of the Cartesian product of all constrained field + /// value lists. + /// + /// @param selection A single selection string (see class-level grammar). + /// @throws std::invalid_argument if the string is malformed. + static std::vector getSelectionMasks(const std::string& selection, const dd4hep::BitFieldCoder& decoder); + + private: + std::vector m_selectors{}; + }; +} // namespace k4ActsTracking + +#endif // K4ACTSTRACKING_CELLIDSELECTOR_H diff --git a/k4ActsTracking/include/k4ActsTracking/IActsGeoSvc.h b/k4ActsTracking/include/k4ActsTracking/IActsGeoSvc.h index 250ed777..63a44626 100644 --- a/k4ActsTracking/include/k4ActsTracking/IActsGeoSvc.h +++ b/k4ActsTracking/include/k4ActsTracking/IActsGeoSvc.h @@ -22,8 +22,10 @@ #include -#include +#include + #include +#include #include namespace dd4hep { @@ -40,14 +42,15 @@ namespace Acts { class GAUDI_API IActsGeoSvc : virtual public IService { public: - using CellIDSurfaceMap = std::unordered_map; + using CellIDSurfaceMap = std::unordered_map; public: DeclareInterfaceID(IActsGeoSvc, 1, 0); - virtual std::shared_ptr trackingGeometry() const = 0; - virtual std::shared_ptr magneticField() const = 0; - virtual const CellIDSurfaceMap& cellIdToSurfaceMap() const = 0; + virtual std::shared_ptr trackingGeometry() const = 0; + virtual std::shared_ptr magneticField() const = 0; + virtual const CellIDSurfaceMap& cellIdToSurfaceMap() const = 0; + virtual std::string cellIDEncodingString() const = 0; virtual ~IActsGeoSvc() = default; }; diff --git a/k4ActsTracking/src/CellIDSelector.cxx b/k4ActsTracking/src/CellIDSelector.cxx new file mode 100644 index 00000000..ce6c7fbc --- /dev/null +++ b/k4ActsTracking/src/CellIDSelector.cxx @@ -0,0 +1,162 @@ +/* + * Copyright (c) 2014-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "k4ActsTracking/CellIDSelector.h" + +#include +#include +#include +#include +#include +#include + +namespace { + auto splitString(const std::string_view str, const char delim) { + namespace rv = std::ranges::views; + + return str | rv::split(delim) | + rv::transform([](auto&& subrange) { return std::string_view(subrange.begin(), subrange.end()); }); + } + + // Implement a very simple polyfill to maintain compatibility with c++23 + template auto to_vector(R&& range) { + std::vector container; + if constexpr (std::ranges::sized_range) { + container.reserve(std::ranges::size(range)); + } + std::ranges::copy(range, std::back_inserter(container)); + return container; + } + + std::pair> getFieldAndValues(std::string_view partialSelection) { + const auto fieldConfig = to_vector(splitString(partialSelection, ':')); + if (fieldConfig.size() != 2) { + throw std::invalid_argument(std::string(partialSelection) + "' is not a valid selection string"); + } + + const auto fieldName = std::string(fieldConfig[0]); + const auto fieldValue = fieldConfig[1]; + if (fieldValue.empty()) { + throw std::invalid_argument("'" + std::string(partialSelection) + "' has an empty value"); + } + + namespace rv = std::ranges::views; + auto layers = splitString(fieldValue, '|') | rv::transform([&](auto&& elem) { + int val; + if (std::from_chars(elem.data(), elem.data() + elem.size(), val).ec != std::errc{}) { + throw std::invalid_argument("'" + std::string(elem) + "' in '" + std::string(partialSelection) + + "' cannot be converted to an integer"); + } + return val; + }); + return {fieldName, to_vector(layers)}; + } + + struct FieldValue { + std::string name; + int value; + }; + + std::vector> cartesianProductFields(auto&& fieldsAndValues) { + // fieldsAndValues is a range with elements of shape pair>. Each element in the range represents a field and values for + // that field. What we want is each combination of {(field_i, value_i), + // (field_j, value_j), (...)} + std::vector> product{}; + + for (const auto& [field, values] : fieldsAndValues) { + std::vector thisFieldValues{}; + thisFieldValues.reserve(values.size()); + for (const int value : values) { + thisFieldValues.emplace_back(field, value); + } + + // If we are in the first field product is still empty and we simply store + // each individual value from this field there. + if (product.empty()) { + for (auto&& fv : thisFieldValues) { + product.emplace_back(std::vector{std::move(fv)}); + } + } else { + // Otherwise we create new products by attaching the new field / values + // to the existing ones and then simply replace the existing ones + std::vector> newProduct{}; + for (auto existing : product) { + for (const auto& fv : thisFieldValues) { + existing.emplace_back(fv); + newProduct.emplace_back(existing); + } + } + product = std::move(newProduct); + } + } + + return product; + } + +} // namespace + +namespace k4ActsTracking { + CellIDSelector::CellIDSelector(const std::string& encodingString, const std::vector& selections) { + dd4hep::BitFieldCoder decoder{encodingString}; + + for (const auto& selection : selections) { + for (auto&& sel : getSelectionMasks(selection, decoder)) { + m_selectors.emplace_back(std::move(sel)); + } + } + } + + std::vector CellIDSelector::getSelectionMasks(const std::string& selection, + const dd4hep::BitFieldCoder& decoder) { + if (selection.empty()) { + throw std::invalid_argument("selection string must not be empty"); + } + namespace rv = std::ranges::views; + auto fieldsAndValues = + splitString(selection, ',') | rv::transform([](auto&& part) { return getFieldAndValues(part); }); + + dd4hep::CellID mask = 0; + for (const auto& [field, _] : fieldsAndValues) { + mask |= decoder[field].mask(); + } + + std::vector selectors{}; + for (const auto& fieldsValues : cartesianProductFields(fieldsAndValues)) { + dd4hep::CellID value{}; + for (const auto& [name, val] : fieldsValues) { + decoder.set(value, std::string(name), val); + } + selectors.emplace_back(mask, value); + } + + return selectors; + } + + bool CellIDSelector::accept(const dd4hep::CellID cellID) const { + for (const auto& selector : m_selectors) { + // The passed CellID has to be equal to any of the configured values after + // masking irrelevant parts + if ((selector.value & selector.mask) == (cellID & selector.mask)) { + return true; + } + } + return false; + } +} // namespace k4ActsTracking diff --git a/k4ActsTracking/src/components/ActsGeoSvc.cpp b/k4ActsTracking/src/components/ActsGeoSvc.cpp index 0b1f7bce..97806db3 100644 --- a/k4ActsTracking/src/components/ActsGeoSvc.cpp +++ b/k4ActsTracking/src/components/ActsGeoSvc.cpp @@ -73,6 +73,9 @@ StatusCode ActsGeoSvc::initialize() { m_geoSvc = Gaudi::svcLocator()->service("GeoSvc"); K4_GAUDI_CHECK(m_geoSvc); + m_cellIDEncodingString = m_geoSvc->getDetector()->constantAsString(m_encodingStringConstant.value()); + debug() << "CellID encoding string: " << m_cellIDEncodingString << endmsg; + std::array magneticFieldVector = {0, 0, 0}; std::array position = {0, 0, 0}; m_geoSvc->getDetector()->field().magneticField(position.data(), magneticFieldVector.data()); @@ -149,7 +152,7 @@ StatusCode ActsGeoSvc::initialize() { } if (m_dumpVisualization.value()) { - info() << "Creating visualiztion" << endmsg; + info() << "Creating visualization" << endmsg; // Adjust the scale here to make it easier to import in blender Acts::ObjVisualization3D vis{4, 0.001}; m_trackingGeo->visualize(vis, gctxt); diff --git a/k4ActsTracking/src/components/ActsGeoSvc.h b/k4ActsTracking/src/components/ActsGeoSvc.h index 53a743e5..fcc62d46 100644 --- a/k4ActsTracking/src/components/ActsGeoSvc.h +++ b/k4ActsTracking/src/components/ActsGeoSvc.h @@ -63,8 +63,12 @@ class ActsGeoSvc : public extends { "Name of the 3D visualization file"}; Gaudi::Property m_dumpVisualization{this, "DumpVisualization", false, "Whether or not to create a 3D visualization dump"}; + Gaudi::Property m_encodingStringConstant{ + this, "EncodingStringVariable", "GlobalTrackerReadoutID", + "Name of the DD4hep constant holding the CellID encoding string."}; const CellIDSurfaceMap& cellIdToSurfaceMap() const override { return m_cellIDToSurface; } + std::string cellIDEncodingString() const override { return m_cellIDEncodingString; } private: using BlueprintBuilder = ActsPlugins::DD4hep::BlueprintBuilder; @@ -76,6 +80,7 @@ class ActsGeoSvc : public extends { std::shared_ptr m_magneticField{nullptr}; std::unordered_map m_cellIDToSurface{}; std::unordered_map m_bluePrintPopulationFuncs{}; + std::string m_cellIDEncodingString{}; }; inline std::shared_ptr ActsGeoSvc::trackingGeometry() const { return m_trackingGeo; } diff --git a/k4ActsTracking/src/components/CKFTrackingAlg.cpp b/k4ActsTracking/src/components/CKFTrackingAlg.cpp new file mode 100644 index 00000000..55b765ff --- /dev/null +++ b/k4ActsTracking/src/components/CKFTrackingAlg.cpp @@ -0,0 +1,656 @@ +/* + * Copyright (c) 2014-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +// k4ActsTracking +#include "k4ActsTracking/CKFTrackingAlg.hxx" +#include "k4ActsTracking/CellIDSelector.h" +#include "k4ActsTracking/Helpers.hxx" +#include "k4ActsTracking/IActsGeoSvc.h" +#include "k4ActsTracking/Measurement.hxx" +#include "k4ActsTracking/MeasurementCalibrator.hxx" +#include "k4ActsTracking/SeedSpacePoint.hxx" +#include "k4ActsTracking/SourceLink.hxx" +#include "k4ActsTracking/SpacePointContainer.hxx" + +// k4FWCore +#include +#include + +#include +// edm4hep +#include +#include +#include + +// Gaudi +#include +#include + +// DD4hep +#include + +// ACTS +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// TBB +#include +#include +#include +#include + +#include + +// Standard +#include +#include +#include +#include +#include +#include + +using namespace Acts::UnitLiterals; + +template <> struct fmt::formatter : fmt::ostream_formatter {}; + +/** + * @brief Seeded CKF tracking algorithm using ActsGeoSvc. + * + * Functionally equivalent to ACTSSeededCKFTrackingAlg but obtains the + * tracking geometry, magnetic field and surface–cellID mapping directly + * from IActsGeoSvc instead of inheriting from ACTSAlgBase. + * + * ACTS contexts (geometry, magnetic-field, calibration) are default- + * constructed as documented by the ACTS framework. + */ +struct CKFTrackingAlg final + : k4FWCore::MultiTransformer( + const edm4hep::TrackerHitPlaneCollection&, const edm4hep::TrackerHitSimTrackerHitLinkCollection&)> { + using TrackContainer = Acts::TrackContainer; + using TrackFinderOptions = Acts::CombinatorialKalmanFilterOptions; + + using SSPoint = typename Acts::SpacePointContainer< + ACTSTracking::SpacePointContainer>, + Acts::detail::RefHolder>::SpacePointProxyType; + + using SSPointGrid = Acts::CylindricalSpacePointGrid; + + using Stepper = Acts::EigenStepper<>; + using Navigator = Acts::Navigator; + using Propagator = Acts::Propagator; + using CKF = Acts::CombinatorialKalmanFilter; + + CKFTrackingAlg(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize() override; + + std::tuple operator()( + const edm4hep::TrackerHitPlaneCollection& trackerHitCollection, + const edm4hep::TrackerHitSimTrackerHitLinkCollection& trackerHitRelations) const override; + + // ----- private helpers --------------------------------------------------- +private: + std::vector findSeeds(const Acts::SeedFinder& finder, + const Acts::SeedFinderOptions& finderOpts, + const auto& spacePointGroup, const SSPointGrid& grid, + Acts::Range1D middleSpRange, std::size_t mutSpDataSize, + edm4hep::TrackCollection& seedCollection, + Acts::MagneticFieldProvider::Cache& magCache) const; + + StatusCode tracking(const std::vector& paramseeds, const CKF& trackFinder, + const TrackFinderOptions& ckfOptions, Acts::MagneticFieldProvider::Cache& magCache, + edm4hep::TrackCollection& trackCollection) const; + + // ----- Gaudi properties -------------------------------------------------- + + /// @name Run control + ///@{ + Gaudi::Property m_runCKF{this, "RunCKF", true, "Run tracking using CKF. False means stop at seeding."}; + Gaudi::Property m_propagateBackward{this, "PropagateBackward", false, "Extrapolates tracks towards beamline."}; + ///@} + + /// @name Seed-finding configuration + ///@{ + Gaudi::Property m_seedFinding_rMax{this, "SeedFinding_RMax", 150, "Maximum radius of hits to consider."}; + Gaudi::Property m_seedFinding_deltaRMin{this, "SeedFinding_DeltaRMin", 5, + "Minimum dR between hits in a seed."}; + Gaudi::Property m_seedFinding_deltaRMax{this, "SeedFinding_DeltaRMax", 80, + "Maximum dR between hits in a seed."}; + Gaudi::Property m_seedFinding_deltaRMinTop{this, "SeedFinding_DeltaRMinTop", 0.f, + "Minimum dR between reference and outer hits in a seed."}; + Gaudi::Property m_seedFinding_deltaRMaxTop{this, "SeedFinding_DeltaRMaxTop", 0.f, + "Maximum dR between reference and outer hits in a seed."}; + Gaudi::Property m_seedFinding_deltaRMinBottom{this, "SeedFinding_DeltaRMinBottom", 0.f, + "Minimum dR between reference and inner hits in a seed."}; + Gaudi::Property m_seedFinding_deltaRMaxBottom{this, "SeedFinding_DeltaRMaxBottom", 0.f, + "Maximum dR between reference and inner hits in a seed."}; + Gaudi::Property m_seedFinding_collisionRegion{this, "SeedFinding_CollisionRegion", 75.0, + "Half-size of the collision region along z."}; + Gaudi::Property m_seedFinding_zMax{this, "SeedFinding_ZMax", 600.0, "Maximum |z| of hits to consider."}; + Gaudi::Property m_seedFinding_sigmaScattering{this, "SeedFinding_SigmaScattering", 50.0, + "Number of sigmas allowed in scattering angle."}; + Gaudi::Property m_seedFinding_radLengthPerSeed{this, "SeedFinding_RadLengthPerSeed", 0.1, + "Average radiation length per seed."}; + Gaudi::Property m_seedFinding_minPt{this, "SeedFinding_MinPt", 500.0, "Minimum pT of tracks to seed [MeV]."}; + Gaudi::Property m_seedFinding_impactMax{this, "SeedFinding_ImpactMax", 3.0, + "Maximum d0 of tracks to seed [mm]."}; + + std::vector m_default_empty_vec; + Gaudi::Property> m_seedFinding_zBinEdges{this, "SeedFinding_zBinEdges", m_default_empty_vec, + "Custom z bin edges for the seeding grid."}; + Gaudi::Property m_zTopBinLen{this, "SeedFinding_zTopBinLen", 1, "Number of top bins along Z for seeding."}; + Gaudi::Property m_zBottomBinLen{this, "SeedFinding_zBottomBinLen", 1, + "Number of bottom bins along Z for seeding."}; + Gaudi::Property m_phiTopBinLen{this, "SeedFinding_phiTopBinLen", 1, "Number of top bins along phi for seeding."}; + Gaudi::Property m_phiBottomBinLen{this, "SeedFinding_phiBottomBinLen", 1, + "Number of bottom bins along phi for seeding."}; + ///@} + + /// @name Track-fit initial error estimates + ///@{ + Gaudi::Property m_initialTrackError_pos{this, "InitialTrackError_Pos", 10 * Acts::UnitConstants::um, + "Initial track error for local position."}; + Gaudi::Property m_initialTrackError_phi{this, "InitialTrackError_Phi", 1 * Acts::UnitConstants::degree, + "Initial track error for phi."}; + Gaudi::Property m_initialTrackError_relP{this, "InitialTrackError_RelP", 0.25, + "Initial track error for momentum (relative)."}; + Gaudi::Property m_initialTrackError_lambda{this, "InitialTrackError_Lambda", 1 * Acts::UnitConstants::degree, + "Initial track error for lambda."}; + Gaudi::Property m_initialTrackError_time{this, "InitialTrackError_Time", 100 * Acts::UnitConstants::ns, + "Initial track error for time."}; + Gaudi::Property m_CKF_chi2CutOff{this, "CKF_Chi2CutOff", 15, "Maximum local chi2 contribution."}; + Gaudi::Property m_CKF_numMeasurementsCutOff{this, "CKF_NumMeasurementsCutOff", 10, + "Maximum measurements on a single surface."}; + ///@} + + /// @name Seeding layer selection + ///@{ + Gaudi::Property> m_seedingSensorsCellIDs{ + this, "SeedingSensorsCellIDs", m_default_empty_vec, + "CellIDSelector selection strings for seed space-point filtering. " + "Each entry is a comma-separated list of field:value constraints (e.g. \"system:1,layer:2|3\"). " + "Multiple entries are OR-ed together. Omitted fields act as wildcards."}; + ///@} + + /// @name Multi-threading + ///@{ + Gaudi::Property m_numThreads{this, "NumThreads", 1, "Number of threads for internal TBB parallelism."}; + ///@} + + SmartIF m_actsGeoSvc; + + // Track finder (propagator) is geometry/field dependent only, so it is built + // once in initialize() and reused (read-only) across events and threads. + std::optional m_trackFinder{}; + + k4ActsTracking::CellIDSelector m_seedSelector{}; + + mutable std::mutex m_seedMutex{}; + mutable std::mutex m_trackMutex{}; +}; + +// --------------------------------------------------------------------------- +// Constructor +// --------------------------------------------------------------------------- + +CKFTrackingAlg::CKFTrackingAlg(const std::string& name, ISvcLocator* svcLoc) + : MultiTransformer(name, svcLoc, + {KeyValue("InputTrackerHitCollection", "TrackerHits"), + KeyValue("InputTrackerHitRelationCollection", "TrackerHitRelations")}, + {KeyValue("OutputSeedCollection", "SeedTracks"), KeyValue("OutputTrackCollection", "Tracks")}) {} + +StatusCode CKFTrackingAlg::initialize() { + m_actsGeoSvc = svcLoc()->service("ActsGeoSvc"); + K4_GAUDI_CHECK(m_actsGeoSvc); + + m_seedSelector = + k4ActsTracking::CellIDSelector(m_actsGeoSvc->cellIDEncodingString(), m_seedingSensorsCellIDs.value()); + + // Apply deltaR fallback defaults + if (m_seedFinding_deltaRMinTop == 0.f) + m_seedFinding_deltaRMinTop = m_seedFinding_deltaRMin; + if (m_seedFinding_deltaRMaxTop == 0.f) + m_seedFinding_deltaRMaxTop = m_seedFinding_deltaRMax; + if (m_seedFinding_deltaRMinBottom == 0.f) + m_seedFinding_deltaRMinBottom = m_seedFinding_deltaRMin; + if (m_seedFinding_deltaRMaxBottom == 0.f) + m_seedFinding_deltaRMaxBottom = m_seedFinding_deltaRMax; + + // The propagator and CKF only depend on the tracking geometry and magnetic + // field, both available here, so build them once instead of per event. + Navigator::Config navigatorCfg{m_actsGeoSvc->trackingGeometry()}; + navigatorCfg.resolvePassive = false; + navigatorCfg.resolveMaterial = true; + navigatorCfg.resolveSensitive = true; + + Stepper stepper(m_actsGeoSvc->magneticField()); + Navigator navigator(navigatorCfg); + Propagator propagator(std::move(stepper), std::move(navigator)); + m_trackFinder.emplace(std::move(propagator)); + + return StatusCode::SUCCESS; +} + +std::tuple CKFTrackingAlg::operator()( + const edm4hep::TrackerHitPlaneCollection& trackerHitCollection, + const edm4hep::TrackerHitSimTrackerHitLinkCollection& /*trackerHitRelations*/) const { + edm4hep::TrackCollection seedCollection; + edm4hep::TrackCollection trackCollection; + + // Default-construct ACTS contexts + const Acts::GeometryContext geoCtx = Acts::GeometryContext::dangerouslyDefaultConstruct(); + const Acts::MagneticFieldContext magCtx{}; + const Acts::CalibrationContext calCtx{}; + + const auto& cellIdToSurface = m_actsGeoSvc->cellIdToSurfaceMap(); + + dd4hep::DDSegmentation::BitFieldCoder decoder{m_actsGeoSvc->cellIDEncodingString()}; + + std::vector> sortedHits; + ACTSTracking::SourceLinkContainer sourceLinks; + ACTSTracking::MeasurementContainer measurements; + ACTSTracking::SeedSpacePointContainer spacePoints; + + sortedHits.reserve(trackerHitCollection.size()); + + for (const auto& hit : trackerHitCollection) { + verbose() << fmt::format("Adding hit {} with cell id {:x}", hit.id(), hit.getCellID()) << endmsg; + auto it = cellIdToSurface.find(hit.getCellID()); + if (it == cellIdToSurface.end()) { + warning() << "No surface found for cellID " << hit.getCellID() << ". skipping hit for tracking." << endmsg; + continue; + } + sortedHits.push_back({it->second->geometryId(), hit}); + } + debug() << "Working with " << sortedHits.size() << " hits." << endmsg; + + // Sort hits by geometry ID for efficient SourceLink multiset insertion + auto compare = [](const auto& a, const auto& b) { return a.first < b.first; }; + tbb::task_arena arena(m_numThreads.value()); + if (m_numThreads > 1) { + arena.execute([&] { tbb::parallel_sort(sortedHits.begin(), sortedHits.end(), compare); }); + } else { + std::sort(sortedHits.begin(), sortedHits.end(), compare); + } + + sourceLinks.reserve(sortedHits.size()); + + for (const auto& hitPair : sortedHits) { + const Acts::Surface* surface = m_actsGeoSvc->trackingGeometry()->findSurface(hitPair.first); + if (surface == nullptr) { + warning() << "Surface with geoID " << hitPair.first + << " not found in tracking geometry. Skipping hit for tracking." << endmsg; + continue; + } + + const edm4hep::Vector3d& edmGlobalPos = hitPair.second.getPosition(); + Acts::Vector3 globalPos = {edmGlobalPos.x, edmGlobalPos.y, edmGlobalPos.z}; + + verbose() << "Converting hit " << hitPair.second.id() << " to local position (pos = " << edmGlobalPos + << ") using surface with geoId " << hitPair.first << endmsg; + + Acts::Result lpResult = surface->globalToLocal(geoCtx, globalPos, {0, 0, 0}, 0.5_um); + if (!lpResult.ok()) { + warning() << "Global to local transformation did not succeed for hit. Skipping it in tracking." << endmsg; + continue; + } + + Acts::Vector2 loc = lpResult.value(); + + Acts::SquareMatrix2 localCov = Acts::SquareMatrix2::Zero(); + const edm4hep::TrackerHitPlane hitplane = hitPair.second; + localCov(0, 0) = std::pow(hitplane.getDu() * Acts::UnitConstants::mm, 2); + localCov(1, 1) = std::pow(hitplane.getDv() * Acts::UnitConstants::mm, 2); + + ACTSTracking::SourceLink sourceLink(surface->geometryId(), measurements.size(), hitPair.second); + Acts::SourceLink srcWrap{sourceLink}; + ACTSTracking::Measurement meas = + ACTSTracking::makeMeasurement(srcWrap, loc, localCov, Acts::eBoundLoc0, Acts::eBoundLoc1); + + measurements.push_back(meas); + sourceLinks.emplace_hint(sourceLinks.end(), sourceLink); + + // Create space point for seeding if this surface is selected + if (m_seedSelector.accept(hitPair.second.getCellID())) { + Acts::RotationMatrix3 rotLocalToGlobal = surface->referenceFrame(geoCtx, globalPos, {0, 0, 0}); + + // Jacobian from global (x,y,z) to (rho, z) + double x = globalPos[Acts::ePos0]; + double y = globalPos[Acts::ePos1]; + double scale = 2 / std::hypot(x, y); + Acts::Matrix<2, 3> jacXyzToRhoZ = Acts::Matrix<2, 3>::Zero(); + jacXyzToRhoZ(0, Acts::ePos0) = scale * x; + jacXyzToRhoZ(0, Acts::ePos1) = scale * y; + jacXyzToRhoZ(1, Acts::ePos2) = 1; + const auto jac = jacXyzToRhoZ * rotLocalToGlobal.block<3, 2>(Acts::ePos0, Acts::ePos0); + const auto var = (jac * localCov * jac.transpose()).diagonal(); + + spacePoints.push_back(ACTSTracking::SeedSpacePoint(globalPos, var[0], var[1], sourceLink)); + } + } + + debug() << fmt::format("Created {} sourceLinks and {} space points for seeding", sourceLinks.size(), + spacePoints.size()) + << endmsg; + + Acts::MagneticFieldProvider::Cache magCache = m_actsGeoSvc->magneticField()->makeCache(magCtx); + + static const Acts::Vector3 zeropos(0, 0, 0); + + Acts::SeedFinderConfig finderCfg; + finderCfg.rMax = m_seedFinding_rMax; + finderCfg.deltaRMin = m_seedFinding_deltaRMin; + finderCfg.deltaRMax = m_seedFinding_deltaRMax; + finderCfg.deltaRMinTopSP = m_seedFinding_deltaRMinTop; + finderCfg.deltaRMaxTopSP = m_seedFinding_deltaRMaxTop; + finderCfg.deltaRMinBottomSP = m_seedFinding_deltaRMinBottom; + finderCfg.deltaRMaxBottomSP = m_seedFinding_deltaRMaxBottom; + finderCfg.collisionRegionMin = -m_seedFinding_collisionRegion; + finderCfg.collisionRegionMax = m_seedFinding_collisionRegion; + finderCfg.zMin = -m_seedFinding_zMax; + finderCfg.zMax = m_seedFinding_zMax; + finderCfg.maxSeedsPerSpM = 1; + finderCfg.cotThetaMax = 7.40627; // ~2.7 η + finderCfg.sigmaScattering = m_seedFinding_sigmaScattering; + finderCfg.radLengthPerSeed = m_seedFinding_radLengthPerSeed; + finderCfg.minPt = m_seedFinding_minPt * Acts::UnitConstants::MeV; + finderCfg.impactMax = m_seedFinding_impactMax * Acts::UnitConstants::mm; + finderCfg.useVariableMiddleSPRange = true; + + Acts::SeedFilterConfig filterCfg; + filterCfg.maxSeedsPerSpM = finderCfg.maxSeedsPerSpM; + finderCfg.seedFilter = std::make_unique>(filterCfg); + finderCfg = finderCfg.calculateDerivedQuantities(); + + Acts::SeedFinderOptions finderOpts; + finderOpts.bFieldInZ = (*m_actsGeoSvc->magneticField()->getField(zeropos, magCache))[2]; + finderOpts.beamPos = {0, 0}; + finderOpts = finderOpts.calculateDerivedQuantities(finderCfg); + + Acts::CylindricalSpacePointGridConfig gridCfg; + gridCfg.cotThetaMax = finderCfg.cotThetaMax; + gridCfg.deltaRMax = finderCfg.deltaRMax; + gridCfg.minPt = finderCfg.minPt; + gridCfg.rMax = finderCfg.rMax; + gridCfg.zMax = finderCfg.zMax; + gridCfg.zMin = finderCfg.zMin; + gridCfg.impactMax = finderCfg.impactMax; + if (!m_seedFinding_zBinEdges.empty()) { + gridCfg.zBinEdges.resize(m_seedFinding_zBinEdges.size()); + for (std::size_t k = 0; k < m_seedFinding_zBinEdges.size(); k++) { + float pos = std::atof(m_seedFinding_zBinEdges[k].c_str()); + if (pos >= finderCfg.zMin && pos < finderCfg.zMax) { + gridCfg.zBinEdges[k] = pos; + } else { + warning() << "Wrong parameter SeedFinding_zBinEdges; default used" << endmsg; + gridCfg.zBinEdges.clear(); + break; + } + } + } + + Acts::CylindricalSpacePointGridOptions gridOpts; + gridOpts.bFieldInZ = (*m_actsGeoSvc->magneticField()->getField(zeropos, magCache))[2]; + + // Wrap space points for ACTS seed finder + std::vector spacePointPtrs(spacePoints.size(), nullptr); + std::transform(spacePoints.begin(), spacePoints.end(), spacePointPtrs.begin(), + [](const ACTSTracking::SeedSpacePoint& sp) { return &sp; }); + + Acts::SpacePointContainerConfig spConfig; + spConfig.useDetailedDoubleMeasurementInfo = finderCfg.useDetailedDoubleMeasurementInfo; + Acts::SpacePointContainerOptions spOptions; + spOptions.beamPos = {0., 0.}; + + ACTSTracking::SpacePointContainer container(spacePointPtrs); + Acts::SpacePointContainer spContainer(spConfig, spOptions, container); + + SSPointGrid grid = Acts::CylindricalSpacePointGridCreator::createGrid(gridCfg, gridOpts); + Acts::CylindricalSpacePointGridCreator::fillGrid(finderCfg, finderOpts, grid, spContainer); + + const Acts::GridBinFinder<3ul> bottomBinFinder(m_phiBottomBinLen.value(), m_zBottomBinLen.value(), 0); + const Acts::GridBinFinder<3ul> topBinFinder(m_phiTopBinLen.value(), m_zTopBinLen.value(), 0); + + Acts::SeedFinder finder(finderCfg); + + const CKF& trackFinder = *m_trackFinder; + + Acts::MeasurementSelector::Config measurementSelectorCfg = { + {Acts::GeometryIdentifier(), {{}, {m_CKF_chi2CutOff}, {(std::size_t)(m_CKF_numMeasurementsCutOff)}}}}; + + Acts::PropagatorPlainOptions pOptions{geoCtx, magCtx}; + pOptions.maxSteps = 10000; + if (m_propagateBackward) { + pOptions.direction = Acts::Direction::Backward(); + } + + Acts::GainMatrixUpdater kfUpdater; + Acts::MeasurementSelector measSel{measurementSelectorCfg}; + ACTSTracking::MeasurementCalibrator measCal{measurements}; + + ACTSTracking::SourceLinkAccessor slAccessor; + slAccessor.container = &sourceLinks; + + using TrackStateCreatorType = Acts::TrackStateCreator; + TrackStateCreatorType trackStateCreator; + trackStateCreator.sourceLinkAccessor.template connect<&ACTSTracking::SourceLinkAccessor::range>(&slAccessor); + trackStateCreator.calibrator.template connect<&ACTSTracking::MeasurementCalibrator::calibrate>(&measCal); + trackStateCreator.measurementSelector + .template connect<&Acts::MeasurementSelector::select>(&measSel); + + Acts::CombinatorialKalmanFilterExtensions extensions; + extensions.updater.connect<&Acts::GainMatrixUpdater::operator()>(&kfUpdater); + extensions.createTrackStates.template connect<&TrackStateCreatorType::createTrackStates>(&trackStateCreator); + + TrackFinderOptions ckfOptions = TrackFinderOptions(geoCtx, magCtx, calCtx, extensions, pOptions); + + float minRange = std::numeric_limits::max(); + float maxRange = std::numeric_limits::lowest(); + for (const auto& coll : grid) { + if (coll.empty()) + continue; + minRange = std::min(coll.front()->radius(), minRange); + maxRange = std::max(coll.back()->radius(), maxRange); + } + + auto spacePointsGrouping = Acts::CylindricalBinnedGroup(std::move(grid), bottomBinFinder, topBinFinder); + + const Acts::Range1D rMiddleSPRange(std::floor(minRange / 2) * 2 + finderCfg.deltaRMiddleMinSPRange, + std::floor(maxRange / 2) * 2 - finderCfg.deltaRMiddleMaxSPRange); + + using GroupIterator = decltype(spacePointsGrouping.begin()); + using GroupValue = std::decay_t())>; + std::vector spacePointGroups; + spacePointGroups.reserve(spacePointsGrouping.grid().size()); + for (auto spGroup : spacePointsGrouping) { + spacePointGroups.push_back(spGroup); + } + + auto parallelSeedingAndTracking = [&](const tbb::blocked_range& r) { + // The magnetic-field cache is mutated on every field lookup, so each + // parallel invocation needs its own cache rather than sharing one. + Acts::MagneticFieldProvider::Cache localMagCache = m_actsGeoSvc->magneticField()->makeCache(magCtx); + for (size_t i = r.begin(); i != r.end(); ++i) { + const auto paramseeds = findSeeds(finder, finderOpts, spacePointGroups[i], spacePointsGrouping.grid(), + rMiddleSPRange, spContainer.size(), seedCollection, localMagCache); + if (!m_runCKF) + continue; + if (!tracking(paramseeds, trackFinder, ckfOptions, localMagCache, trackCollection).isSuccess()) { + warning() << "Tracking failed for this event" << endmsg; + } + } + }; + + if (m_numThreads > 1) { + arena.execute( + [&] { tbb::parallel_for(tbb::blocked_range(0, spacePointGroups.size()), parallelSeedingAndTracking); }); + } else { + for (size_t i = 0; i < spacePointGroups.size(); ++i) { + parallelSeedingAndTracking(tbb::blocked_range(i, i + 1)); + } + } + + debug() << "Track Collection Size: " << trackCollection.size() << endmsg; + return std::make_tuple(std::move(seedCollection), std::move(trackCollection)); +} + +std::vector CKFTrackingAlg::findSeeds( + const Acts::SeedFinder& finder, const Acts::SeedFinderOptions& finderOpts, + const auto& spacePointGroup, const SSPointGrid& grid, Acts::Range1D middleSpRange, std::size_t mutSpDataSize, + edm4hep::TrackCollection& seedCollection, Acts::MagneticFieldProvider::Cache& magCache) const { + const Acts::GeometryContext geoCtx = Acts::GeometryContext::dangerouslyDefaultConstruct(); + + const auto& [bottom, middle, top] = spacePointGroup; + + std::vector> seeds; + std::vector paramseeds; + Acts::SeedFinder::SeedingState state; + state.spacePointMutableData.resize(mutSpDataSize); + + finder.createSeedsForGroup(finderOpts, state, grid, seeds, bottom, middle, top, middleSpRange); + + // Unwrap proxy types back to concrete SeedSpacePoint + std::vector> f_seeds; + f_seeds.reserve(seeds.size()); + for (const Acts::Seed& seed : seeds) { + const auto& sps = seed.sp(); + f_seeds.emplace_back(*sps[0]->externalSpacePoint(), *sps[1]->externalSpacePoint(), *sps[2]->externalSpacePoint()); + } + + for (const auto& seed : f_seeds) { + const ACTSTracking::SeedSpacePoint* bottomSP = seed.sp().front(); + const auto& sourceLink = bottomSP->sourceLink(); + const Acts::GeometryIdentifier& geoId = sourceLink.geometryId(); + const Acts::Surface* surface = m_actsGeoSvc->trackingGeometry()->findSurface(geoId); + if (surface == nullptr) { + warning() << "Surface with geoID " << geoId << " not found in tracking geometry" << endmsg; + continue; + } + + // Magnetic field at the seed position + const Acts::Vector3 seedPos(bottomSP->x(), bottomSP->y(), bottomSP->z()); + Acts::Result seedField = m_actsGeoSvc->magneticField()->getField(seedPos, magCache); + if (!seedField.ok()) { + throw std::runtime_error("Field lookup error: " + std::to_string(seedField.error().value())); + } + + Acts::Result optParams = + Acts::estimateTrackParamsFromSeed(geoCtx, seed.sp(), *surface, *seedField); + if (!optParams.ok()) { + debug() << "Failed estimation of track parameters for seed." << endmsg; + continue; + } + + const Acts::BoundVector& params = *optParams; + float p = std::abs(1.f / params[Acts::eBoundQOverP]); + + Acts::BoundMatrix cov = ACTSTracking::makeInitialCovariance(p, m_initialTrackError_pos, m_initialTrackError_phi, + m_initialTrackError_lambda, m_initialTrackError_relP, + m_initialTrackError_time); + + Acts::BoundTrackParameters paramseed(surface->getSharedPtr(), params, cov, Acts::ParticleHypothesis::pion()); + paramseeds.push_back(paramseed); + + // Compute seed track state before acquiring the lock + Acts::Vector3 globalPos = + surface->localToGlobal(geoCtx, {params[Acts::eBoundLoc0], params[Acts::eBoundLoc1]}, {0, 0, 0}); + Acts::Result hitField = m_actsGeoSvc->magneticField()->getField(globalPos, magCache); + if (!hitField.ok()) { + throw std::runtime_error("Field lookup error: " + std::to_string(hitField.error().value())); + } + + auto seedTrackState = ACTSTracking::ACTS2edm4hep_trackState(edm4hep::TrackState::AtFirstHit, paramseed, + (*hitField)[2] / Acts::UnitConstants::T); + + { + std::lock_guard lock(m_seedMutex); + auto seedTrack = seedCollection.create(); + for (const ACTSTracking::SeedSpacePoint* sp : seed.sp()) { + seedTrack.addToTrackerHits(sp->sourceLink().edm4hepHit()); + } + seedTrack.addToTrackStates(seedTrackState); + } + + debug() << "Seed Parameters" << std::endl << paramseed << endmsg; + } + + debug() << "Seeds found: " << paramseeds.size() << endmsg; + return paramseeds; +} + +StatusCode CKFTrackingAlg::tracking(const std::vector& paramseeds, const CKF& trackFinder, + const TrackFinderOptions& ckfOptions, Acts::MagneticFieldProvider::Cache& magCache, + edm4hep::TrackCollection& trackCollection) const { + const Acts::GeometryContext geoCtx = Acts::GeometryContext::dangerouslyDefaultConstruct(); + + debug() << "Starting CKF track finding with " << paramseeds.size() << " seeds." << endmsg; + + auto trackContainer = std::make_shared(); + auto trackStateContainer = std::make_shared(); + TrackContainer tracks(trackContainer, trackStateContainer); + + for (std::size_t iseed = 0; iseed < paramseeds.size(); ++iseed) { + tracks.clear(); + auto result = trackFinder.findTracks(paramseeds.at(iseed), ckfOptions, tracks); + if (result.ok()) { + const auto& fitOutput = result.value(); + for (const TrackContainer::TrackProxy& trackItem : fitOutput) { + auto trackTip = tracks.makeTrack(); + trackTip.copyFrom(trackItem); + auto smoothResult = Acts::smoothTrack(geoCtx, trackTip); + if (!smoothResult.ok()) { + warning() << "Track smoothing error: " << smoothResult.error() << endmsg; + continue; + } + + auto track = ACTSTracking::ACTS2edm4hep_track(trackTip, m_actsGeoSvc->magneticField(), magCache); + { + std::lock_guard lock{m_trackMutex}; + trackCollection.push_back(track); + } + } + } else { + warning() << "Track fit error: " << result.error() << endmsg; + } + } + + return StatusCode::SUCCESS; +} + +DECLARE_COMPONENT(CKFTrackingAlg); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6898adc9..17fe74ea 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,15 +17,30 @@ See the License for the specific language governing permissions and limitations under the License. ]] +add_subdirectory(unittests) + +# Force-load the freshly built Gaudi plugin module. +# +# When the sourced key4hep stack already ships this package (e.g. the nightly +# view contains an older libk4ActsTrackingPlugins.so), Gaudi's PluginService +# resolves plugin libraries relative to where libGaudiPluginService.so lives +# (the view) and ignores LD_LIBRARY_PATH/GAUDI_PLUGIN_PATH for the dlopen. The +# local build is then shadowed by the stale view copy, so components added on +# this branch (e.g. CKFTrackingAlg) are reported as "No factory available". +# LD_PRELOADing the locally built module makes its factories register first; +# PluginService keeps the first registration, so the local components win. set(test_environment "\ +LD_PRELOAD=$;\ LD_LIBRARY_PATH=\ ${PROJECT_BINARY_DIR}:\ ${PROJECT_BINARY_DIR}/${PROJECT_NAME}:\ -$ENV{LD_LIBRARY_PATH}; +$ENV{LD_LIBRARY_PATH};\ PYTHONPATH=\ +${CMAKE_CURRENT_SOURCE_DIR}/options:\ ${PROJECT_BINARY_DIR}/${PROJECT_NAME}/genConfDir:\ $ENV{PYTHONPATH};\ GAUDI_PLUGIN_PATH=\ +${PROJECT_BINARY_DIR}:\ ${PROJECT_BINARY_DIR}/${PROJECT_NAME}:\ $ENV{GAUDI_PLUGIN_PATH}" ) @@ -34,19 +49,92 @@ function(set_test_env _testname) set_tests_properties(${_testname} PROPERTIES ENVIRONMENT "${test_environment}") endfunction() +# add_clone_config(URL ) +# +# Derives the clone directory, test name (clone_), and fixture name () +# from the URL basename. Skips the clone if the directory already exists. +function(add_clone_config) + cmake_parse_arguments(ARG "" "URL" "" ${ARGN}) + get_filename_component(_dir "${ARG_URL}" NAME_WE) + set(_cmd "if [ ! -d ${_dir} ]; then git clone ${ARG_URL} --depth 1; fi") + add_test(NAME clone_${_dir} COMMAND bash -c "${_cmd}") + set_tests_properties(clone_${_dir} PROPERTIES FIXTURES_SETUP ${_dir}) +endfunction() + +# add_sim_test( +# DDSIM_EXTRA_ARGS ... +# [WORKING_DIRECTORY ] +# [FIXTURES_REQUIRED ...] [FIXTURES_SETUP ]) +# +# Runs ddsim with a fixed 10-event mu- particle gun (multiplicity 10, uniform +# distribution). Detector-specific args (--compactFile, --outputFile, +# --gun.momentumMin/Max, etc.) go in DDSIM_EXTRA_ARGS and are forwarded verbatim. +function(add_sim_test _name) + cmake_parse_arguments(ARG "" + "WORKING_DIRECTORY;FIXTURES_SETUP" + "EXTRA_ARGS;FIXTURES_REQUIRED" ${ARGN}) + if(NOT ARG_WORKING_DIRECTORY) + set(ARG_WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}") + endif() + add_test(NAME ${_name} + WORKING_DIRECTORY ${ARG_WORKING_DIRECTORY} + COMMAND ddsim --runType=batch -G -N=10 + --gun.multiplicity 10 + --gun.particle mu- + --gun.distribution uniform + ${ARG_EXTRA_ARGS} + ) + if(ARG_FIXTURES_REQUIRED) + set_tests_properties(${_name} PROPERTIES FIXTURES_REQUIRED "${ARG_FIXTURES_REQUIRED}") + endif() + if(ARG_FIXTURES_SETUP) + set_tests_properties(${_name} PROPERTIES FIXTURES_SETUP "${ARG_FIXTURES_SETUP}") + endif() +endfunction() + +# add_k4run_test( +# OPTIONS_FILE +# [WORKING_DIRECTORY ] +# [EXTRA_ARGS ...] +# [FIXTURES_REQUIRED ...] [FIXTURES_SETUP ] +# [USE_TEST_ENV]) +# +# USE_TEST_ENV applies the build-tree env (needed when the option file loads +# this project's Gaudi plugins). +function(add_k4run_test _name) + cmake_parse_arguments(ARG "USE_TEST_ENV" + "OPTIONS_FILE;WORKING_DIRECTORY;FIXTURES_SETUP" + "EXTRA_ARGS;FIXTURES_REQUIRED" ${ARGN}) + if(NOT ARG_WORKING_DIRECTORY) + set(ARG_WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}") + endif() + add_test(NAME ${_name} + WORKING_DIRECTORY ${ARG_WORKING_DIRECTORY} + COMMAND k4run ${ARG_OPTIONS_FILE} ${ARG_EXTRA_ARGS} + ) + if(ARG_USE_TEST_ENV) + set_test_env(${_name}) + endif() + if(ARG_FIXTURES_REQUIRED) + set_tests_properties(${_name} PROPERTIES FIXTURES_REQUIRED "${ARG_FIXTURES_REQUIRED}") + endif() + if(ARG_FIXTURES_SETUP) + set_tests_properties(${_name} PROPERTIES FIXTURES_SETUP "${ARG_FIXTURES_SETUP}") + endif() +endfunction() + # add_geometry_load_test( ) # # Add a test that loads the passed compact file and converts the geometry into # an ACTS geometry before dumping it to an .obj file function(add_geometry_load_test _compact_file) get_filename_component(_name "${_compact_file}" NAME_WE) - add_test(NAME load_geo_${_name} - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND k4run ${PROJECT_SOURCE_DIR}/k4ActsTracking/examples/test_visualize_acts_geo.py - --compactFile "${_compact_file}" - --test-propagation - ) - set_test_env(load_geo_${_name}) + add_k4run_test(load_geo_${_name} + OPTIONS_FILE ${PROJECT_SOURCE_DIR}/k4ActsTracking/examples/test_visualize_acts_geo.py + EXTRA_ARGS --compactFile "${_compact_file}" + --test-propagation + USE_TEST_ENV + ) endfunction() add_geometry_load_test($ENV{k4geo_DIR}/MuColl/MAIA/compact/MAIA_v0/MAIA_v0.xml) @@ -55,3 +143,120 @@ add_geometry_load_test($ENV{k4geo_DIR}/FCCee/ILD_FCCee/compact/ILD_FCCee_v01/ILD add_geometry_load_test($ENV{k4geo_DIR}/FCCee/ILD_FCCee/compact/ILD_FCCee_v02/ILD_FCCee_v02.xml) add_geometry_load_test($ENV{k4geo_DIR}/FCCee/CLD/compact/CLD_o2_v07/CLD_o2_v07.xml) add_geometry_load_test($ENV{k4geo_DIR}/FCCee/CLD/compact/CLD_o2_v08/CLD_o2_v08.xml) + + +# === MAIA chain (no upstream config, no standard reco baseline) === +set(_maia_compact $ENV{k4geo_DIR}/MuColl/MAIA/compact/MAIA_v0/MAIA_v0.xml) + +add_sim_test(sim_MAIA + EXTRA_ARGS + --compactFile=${_maia_compact} + --outputFile=particle_gun_MAIA_SIM.edm4hep.root + --gun.momentumMin 5000 --gun.momentumMax 500000 + --gun.thetaMin -2.96 --gun.thetaMax 2.96 # roughly +/-170 degrees + FIXTURES_SETUP sim_maia_fixture +) + +add_k4run_test(reco_MAIA_Gen3 + OPTIONS_FILE ${CMAKE_CURRENT_LIST_DIR}/options/MAIA_CKFTrackingAlg.py + EXTRA_ARGS --compactFile=${_maia_compact} + --IOSvc.Input particle_gun_MAIA_SIM.edm4hep.root + --IOSvc.Output particle_gun_MAIA_CKFTracks.edm4hep.root + USE_TEST_ENV + FIXTURES_REQUIRED sim_maia_fixture +) + + +# === CLD chain === +set(_cld_compact $ENV{k4geo_DIR}/FCCee/CLD/compact/CLD_o2_v08/CLD_o2_v08.xml) +set(_cld_dir ${CMAKE_CURRENT_BINARY_DIR}/CLDConfig/CLDConfig) + +add_clone_config(URL https://github.com/key4hep/CLDConfig.git) + +add_sim_test(sim_CLD_o2_v08 + WORKING_DIRECTORY ${_cld_dir} + EXTRA_ARGS + --compactFile=${_cld_compact} + --outputFile=particle_gun_CLD_o2_v08_SIM.edm4hep.root + --gun.momentumMin 1500 --gun.momentumMax 100000 + FIXTURES_REQUIRED CLDConfig + FIXTURES_SETUP sim_CLD_fixture +) + +add_k4run_test(std_reco_CLD_o2_v08 + WORKING_DIRECTORY ${_cld_dir} + OPTIONS_FILE CLDReconstruction.py + EXTRA_ARGS --compactFile ${_cld_compact} + --inputFiles=particle_gun_CLD_o2_v08_SIM.edm4hep.root + --trackingOnly + --outputBasename particle_gun_CLD_o2_v08 + FIXTURES_REQUIRED sim_CLD_fixture + FIXTURES_SETUP reco_CLD_fixture +) + +add_k4run_test(ckf_tracking_CLD_o2_v08 + WORKING_DIRECTORY ${_cld_dir} + OPTIONS_FILE ${CMAKE_CURRENT_LIST_DIR}/options/CLD_CKFTracking.py + EXTRA_ARGS --compactFile ${_cld_compact} + --IOSvc.Input particle_gun_CLD_o2_v08_REC.edm4hep.root + --IOSvc.Output particle_gun_CLD_o2_v08_CKFTracks.edm4hep.root + USE_TEST_ENV + FIXTURES_REQUIRED reco_CLD_fixture +) + + +# === ILD chain === +set(_ild_compact $ENV{k4geo_DIR}/FCCee/ILD_FCCee/compact/ILD_FCCee_v01/ILD_FCCee_v01.xml) +set(_ild_dir ${CMAKE_CURRENT_BINARY_DIR}/ILDConfig/StandardConfig/production) + +add_clone_config(URL https://github.com/iLCSoft/ILDConfig.git) + +add_sim_test(sim_ILD_FCCee_v01 + WORKING_DIRECTORY ${_ild_dir} + EXTRA_ARGS + --compactFile=${_ild_compact} + --outputFile=particle_gun_ILD_FCCee_v01_SIM.edm4hep.root + --gun.momentumMin 1500 --gun.momentumMax 1000000 + FIXTURES_REQUIRED ILDConfig + FIXTURES_SETUP sim_ILD_FCCee_v01_fixture +) + +add_k4run_test(std_reco_ILD_FCCee_v01 + WORKING_DIRECTORY ${_ild_dir} + OPTIONS_FILE ILDReconstruction.py + EXTRA_ARGS --compactFile ${_ild_compact} + --inputFiles particle_gun_ILD_FCCee_v01_SIM.edm4hep.root + --trackingOnly + --outputFileBase particle_gun_ILD_FCCee_v01 + FIXTURES_REQUIRED sim_ILD_FCCee_v01_fixture + FIXTURES_SETUP reco_ILD_fixture + USE_TEST_ENV +) + +add_k4run_test(ckf_tracking_ILD_FCCee_v01 + WORKING_DIRECTORY ${_ild_dir} + OPTIONS_FILE ${CMAKE_CURRENT_LIST_DIR}/options/ILD_CKFTracking.py + EXTRA_ARGS --compactFile ${_ild_compact} + --IOSvc.Input particle_gun_ILD_FCCee_v01_REC.edm4hep.root + --IOSvc.Output particle_gun_ILD_FCCee_v01_CKFTracks.edm4hep.root + USE_TEST_ENV + FIXTURES_REQUIRED reco_ILD_fixture +) + +# Explicitly add the ILD working directory to PYTHONPATH for ILD tests. +# +# This APPENDs a second PYTHONPATH= entry to the test ENVIRONMENT; CTest applies +# env entries in order, so a later duplicate of the same variable OVERRIDES the +# earlier one set by USE_TEST_ENV. We must therefore keep the local genConfDir +# (and options) ahead of the inherited $ENV{PYTHONPATH} here too — otherwise the +# CVMFS view's k4ActsTrackingPluginsConf (which lacks CKFTrackingAlg) shadows the +# local configurable and k4run fails with "Configurable CKFTrackingAlg not found". +set(_ild_pythonpath "\ +${_ild_dir}:\ +${CMAKE_CURRENT_SOURCE_DIR}/options:\ +${PROJECT_BINARY_DIR}/${PROJECT_NAME}/genConfDir:\ +$ENV{PYTHONPATH}") +set_property(TEST std_reco_ILD_FCCee_v01 APPEND PROPERTY + ENVIRONMENT "PYTHONPATH=${_ild_pythonpath}") +set_property(TEST ckf_tracking_ILD_FCCee_v01 APPEND PROPERTY + ENVIRONMENT "PYTHONPATH=${_ild_pythonpath}") diff --git a/test/options/CLD_CKFTracking.py b/test/options/CLD_CKFTracking.py new file mode 100644 index 00000000..ccb11128 --- /dev/null +++ b/test/options/CLD_CKFTracking.py @@ -0,0 +1,69 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import os +import sys + +from Gaudi.Configuration import INFO +from k4FWCore import ApplicationMgr, IOSvc + +sys.path.insert(0, os.path.dirname(__file__)) + +from _ckf_helpers import ( + make_ckf_tracking, + make_services, + make_hit_mergers, +) + +svcList = make_services() +iosvc = IOSvc("IOSvc") + +hit_merger, hit_rel_merger = make_hit_mergers( + hit_collections=[ + "VXDTrackerHits", + "VXDEndcapTrackerHits", + "ITrackerHits", + "OTrackerHits", + "ITrackerEndcapHits", + "OTrackerEndcapHits", + ], + relation_collections=[ + "VXDTrackerHitRelations", + "VXDEndcapTrackerHitRelations", + "InnerTrackerBarrelHitsRelations", + "InnerTrackerEndcapHitsRelations", + "OuterTrackerBarrelHitsRelations", + "OuterTrackerEndcapHitsRelations", + ], +) + +ckf_tracking = make_ckf_tracking( + hit_merger, + hit_rel_merger, + seeding_cellids=["system:1|2", "system:3|4,layer:1"], +) + +ApplicationMgr( + TopAlg=[hit_merger, hit_rel_merger, ckf_tracking], + ExtSvc=svcList, + OutputLevel=INFO, + EvtSel="NONE", + EvtMax=-1, +) diff --git a/test/options/ILD_CKFTracking.py b/test/options/ILD_CKFTracking.py new file mode 100644 index 00000000..d725a702 --- /dev/null +++ b/test/options/ILD_CKFTracking.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import os +import sys + +from Gaudi.Configuration import INFO +from k4FWCore import ApplicationMgr, IOSvc + +sys.path.insert(0, os.path.dirname(__file__)) + +from _ckf_helpers import ( + make_ckf_tracking, + make_services, + make_hit_mergers, +) + +svcList = make_services() +iosvc = IOSvc("IOSvc") + +hit_merger, hit_rel_merger = make_hit_mergers( + hit_collections=[ + "VertexBarrelTrackerHits", + "VertexEndcapTrackerHits", + "InnerTrackerBarrelHits", + "InnerTrackerEndcapHits", + "SETTrackerHits", + ], + relation_collections=[ + "VertexBarrelTrackerHitRelations", + "VertexEndcapTrackerHitRelations", + "InnerTrackerBarrelHitRelations", + "InnerTrackerEndcapHitRelations", + "SETTrackerHitRelations", + ], +) + +ckf_tracking = make_ckf_tracking( + hit_merger, + hit_rel_merger, + seeding_cellids=["system:1|2", "system:3|5,layer:1"], +) + +ApplicationMgr( + TopAlg=[hit_merger, hit_rel_merger, ckf_tracking], + ExtSvc=svcList, + OutputLevel=INFO, + EvtSel="NONE", + EvtMax=-1, +) diff --git a/test/options/MAIA_CKFTrackingAlg.py b/test/options/MAIA_CKFTrackingAlg.py new file mode 100644 index 00000000..15449090 --- /dev/null +++ b/test/options/MAIA_CKFTrackingAlg.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import os +import sys + +from Gaudi.Configuration import INFO +from Gaudi.Configurables import DDPlanarDigi +from k4FWCore import ApplicationMgr, IOSvc + +sys.path.insert(0, os.path.dirname(__file__)) + +from _ckf_helpers import ( + make_ckf_tracking, + make_services, + make_hit_mergers, +) + +svcList = make_services() +iosvc = IOSvc("IOSvc") + +algList = [] + +for name in ("VertexBarrel", "VertexEndcap"): + algList.append( + DDPlanarDigi( + f"{name}Digitizer", + CorrectTimesForPropagation=True, + IsStrip=False, + ResolutionT=[0.03], + ResolutionU=[0.005], + ResolutionV=[0.005], + SubDetectorName="Vertex", + TimeWindowMax=[0.15], + TimeWindowMin=[-0.09], + UseTimeWindow=True, + SimTrackHitCollectionName=[f"{name}Collection"], + SimTrkHitRelCollection=[f"{name}HitsRelations"], + TrackerHitCollectionName=[f"{name}Hits"], + ) + ) + +for name in ( + "InnerTrackerBarrel", + "InnerTrackerEndcap", + "OuterTrackerBarrel", + "OuterTrackerEndcap", +): + algList.append( + DDPlanarDigi( + f"{name}Digitizer", + CorrectTimesForPropagation=True, + IsStrip="InnerBarrel" in name or "OuterEndcap" in name, + ResolutionT=[0.06], + ResolutionU=[0.007], + ResolutionV=[0.09], + SubDetectorName=name.replace("Barrel", "s").replace("Endcap", "s"), + TimeWindowMax=[0.3], + TimeWindowMin=[-0.18], + UseTimeWindow=True, + SimTrackHitCollectionName=[f"{name}Collection"], + SimTrkHitRelCollection=[f"{name}HitsRelations"], + TrackerHitCollectionName=[f"{name}Hits"], + ) + ) + +hit_merger, hit_rel_merger = make_hit_mergers( + hit_collections=[ + "VertexBarrelHits", + "VertexEndcapHits", + "InnerTrackerBarrelHits", + "InnerTrackerEndcapHits", + "OuterTrackerBarrelHits", + "OuterTrackerEndcapHits", + ], + relation_collections=[ + "VertexBarrelHitsRelations", + "VertexEndcapHitsRelations", + "InnerTrackerBarrelHitsRelations", + "InnerTrackerEndcapHitsRelations", + "OuterTrackerBarrelHitsRelations", + "OuterTrackerEndcapHitsRelations", + ], +) +algList.extend([hit_merger, hit_rel_merger]) + +ckf_tracking = make_ckf_tracking( + hit_merger, + hit_rel_merger, + seeding_cellids=["system:1", "system:2,layer:1|2|3"], +) +algList.append(ckf_tracking) + +ApplicationMgr( + TopAlg=algList, ExtSvc=svcList, OutputLevel=INFO, EvtSel="NONE", EvtMax=-1 +) diff --git a/test/options/_ckf_helpers.py b/test/options/_ckf_helpers.py new file mode 100644 index 00000000..0bf487f9 --- /dev/null +++ b/test/options/_ckf_helpers.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import os +import sys + +from Gaudi.Configuration import INFO +from Gaudi.Configurables import ( + ActsGeoSvc, + CKFTrackingAlg, + CollectionMerger, + EventDataSvc, + GeoSvc, +) +from k4FWCore.parseArgs import parser + + +def _get_compact_file(): + parser.add_argument( + "--compactFile", + help="The geometry compact file to use for reconstruction", + type=str, + ) + args = parser.parse_known_args()[0] + return args.compactFile + + +def make_services(): + """Configure all the necessary services (including getting the geometry from + the command line geometry)""" + compact_file = _get_compact_file() + return [ + GeoSvc("GeoSvc", detectors=[compact_file], EnableGeant4Geo=False), + ActsGeoSvc("ActsGeoSvc"), + EventDataSvc("EventDataSvc"), + ] + + +def make_hit_mergers( + hit_collections, + relation_collections, + hits_output="AllTrackerHits", + relations_output="AllTrackerHitRelations", +): + """Configure the hit (and relation) mergers""" + hit_merger = CollectionMerger( + "MergeHits", + InputCollections=hit_collections, + OutputCollection=hits_output, + ) + hit_rel_merger = CollectionMerger( + "MergeHitRelations", + InputCollections=relation_collections, + OutputCollection=relations_output, + ) + return hit_merger, hit_rel_merger + + +def make_ckf_tracking( + hit_merger, + hit_rel_merger, + seeding_cellids, +): + """Configure the CKFTrackingAlg such that it does **some** tracking + + NOTE: This is really just an example on how to configure it for some + technical tests! This will not work for any meaningful tracking! + """ + + return CKFTrackingAlg( + "CKFTracking", + RunCKF=True, + CKF_Chi2CutOff=10, + SeedFinding_RMax=150, + SeedFinding_MinPt=500, + SeedFinding_ImpactMax=3, + CKF_NumMeasurementsCutOff=1, + SeedFinding_SigmaScattering=50, + SeedFinding_CollisionRegion=6, + SeedFinding_RadLengthPerSeed=0.1, + SeedingSensorsCellIDs=seeding_cellids, + OutputTrackCollection="CKFTracks", + OutputSeedCollection="CKFTrackSeeds", + InputTrackerHitCollection=hit_merger.OutputCollection, + InputTrackerHitRelationCollection=hit_rel_merger.OutputCollection, + OutputLevel=INFO, + ) diff --git a/test/unittests/CMakeLists.txt b/test/unittests/CMakeLists.txt new file mode 100644 index 00000000..a283d78e --- /dev/null +++ b/test/unittests/CMakeLists.txt @@ -0,0 +1,77 @@ +#[[ +Copyright (c) 2014-2024 Key4hep-Project. + +This file is part of Key4hep. +See https://key4hep.github.io/key4hep-doc/ for further info. + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +]] + +set(CATCH2_MIN_VERSION 3.5.0) + +set(USE_EXTERNAL_CATCH2 AUTO CACHE STRING "Link against an external Catch2 v3 static library, otherwise build it locally") +set_property(CACHE USE_EXTERNAL_CATCH2 PROPERTY STRINGS AUTO ON OFF) + +if(USE_EXTERNAL_CATCH2) + if (USE_EXTERNAL_CATCH2 STREQUAL AUTO) + find_package(Catch2 ${CATCH2_MIN_VERSION}) + else() + find_package(Catch2 ${CATCH2_MIN_VERSION} REQUIRED) + endif() +endif() + +if(NOT Catch2_FOUND) + message(STATUS "Fetching local copy of Catch2 library for unit-tests...") + # Build Catch2 with the default flags, to avoid generating warnings when we + # build it + set(CXX_FLAGS_CMAKE_USED ${CMAKE_CXX_FLAGS}) + set(CMAKE_CXX_FLAGS ${CXX_FLAGS_CMAKE_DEFAULTS}) + Include(FetchContent) + FetchContent_Declare( + Catch2 + GIT_REPOSITORY https://github.com/catchorg/Catch2.git + GIT_TAG v${CATCH2_MIN_VERSION} + ) + FetchContent_MakeAvailable(Catch2) + set(CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras ${CMAKE_MODULE_PATH}) + + # Disable clang-tidy on external contents + set_target_properties(Catch2 PROPERTIES CXX_CLANG_TIDY "") + + # Hack around the fact, that the include directories are not declared as + # SYSTEM for the targets defined this way. Otherwise warnings can still occur + # in Catch2 code when templates are evaluated (which happens quite a bit) + get_target_property(CATCH2_IF_INC_DIRS Catch2 INTERFACE_INCLUDE_DIRECTORIES) + set_target_properties(Catch2 PROPERTIES INTERFACE_SYSTEM_INCLUDE_DIRECTORIES "${CATCH2_IF_INC_DIRS}") + + # Reset the flags + set(CMAKE_CXX_FLAGS ${CXX_FLAGS_CMAKE_USED}) +else() + # Same hack as above for the system-installed Catch2: mark its include + # directories as SYSTEM to suppress warnings from Catch2 headers + get_target_property(CATCH2_IF_INC_DIRS Catch2::Catch2 INTERFACE_INCLUDE_DIRECTORIES) + set_target_properties(Catch2::Catch2 PROPERTIES INTERFACE_SYSTEM_INCLUDE_DIRECTORIES "${CATCH2_IF_INC_DIRS}") +endif() + +add_executable(unittest_k4ActsTracking CellIDSelectionTests.cpp) +target_link_libraries(unittest_k4ActsTracking PUBLIC k4ActsTracking PRIVATE Catch2::Catch2WithMain) +# Catch2 TEST_CASE/SECTION macros expand __COUNTER__ at the call site in user +# code, so marking Catch2 headers as system includes is not sufficient to +# suppress this warning, present with Clang 22 +target_compile_options(unittest_k4ActsTracking PRIVATE -Wno-c2y-extensions) + +include(Catch) +catch_discover_tests(unittest_k4ActsTracking + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR} + TEST_PREFIX "UT_" # make it possible to filter easily with -R ^UT +) diff --git a/test/unittests/CellIDSelectionTests.cpp b/test/unittests/CellIDSelectionTests.cpp new file mode 100644 index 00000000..4b539fc0 --- /dev/null +++ b/test/unittests/CellIDSelectionTests.cpp @@ -0,0 +1,229 @@ +/* + * Copyright (c) 2014-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "k4ActsTracking/CellIDSelector.h" + +#include +#include + +#include +#include "catch2/catch_test_macros.hpp" +#include "catch2/generators/catch_generators.hpp" +#include "catch2/generators/catch_generators_adapters.hpp" +#include "catch2/generators/catch_generators_random.hpp" +#include "catch2/matchers/catch_matchers_range_equals.hpp" + +#include +#include + +using namespace k4ActsTracking; + +TEST_CASE("CellIDSelector::accept single selection") { + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + const auto encoder = dd4hep::BitFieldCoder(encodingString); + const auto selector = CellIDSelector{encodingString, {"system:8,layer:3"}}; + + dd4hep::CellID cellID{0}; + encoder.set(cellID, "system", 8); + encoder.set(cellID, "layer", 3); + encoder.set(cellID, "side", 0); + encoder.set(cellID, "module", 5); + encoder.set(cellID, "sensor", 42); + + REQUIRE(selector.accept(cellID)); + + encoder.set(cellID, "sensor", 123); + encoder.set(cellID, "module", 12); + REQUIRE(selector.accept(cellID)); + + encoder.set(cellID, "system", 5); + REQUIRE_FALSE(selector.accept(cellID)); + + encoder.set(cellID, "system", 8); + encoder.set(cellID, "layer", 16); + REQUIRE_FALSE(selector.accept(cellID)); +} + +TEST_CASE("CellIDSelector::accept multiple selections") { + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + const auto encoder = dd4hep::BitFieldCoder(encodingString); + + const auto selector = CellIDSelector{encodingString, {"system:5,layer:1|4|5", "system:3,layer:2|6|8", "sensor:42"}}; + + dd4hep::CellID cellID{0}; + encoder.set(cellID, "side", -1); + encoder.set(cellID, "module", 5); + + encoder.set(cellID, "system", 8); + encoder.set(cellID, "layer", 4); + // Layer would match but system doesn't + REQUIRE_FALSE(selector.accept(cellID)); + + encoder.set(cellID, "sensor", 42); + encoder.set(cellID, "layer", 3); + // Nothing matches except sensor:42 from the last selection + REQUIRE(selector.accept(cellID)); + + encoder.set(cellID, "sensor", 123); + encoder.set(cellID, "system", 5); + // system:5 would match, but layers are off + REQUIRE_FALSE(selector.accept(cellID)); + encoder.set(cellID, "layer", 5); + REQUIRE(selector.accept(cellID)); + + encoder.set(cellID, "system", 3); + REQUIRE_FALSE(selector.accept(cellID)); + encoder.set(cellID, "layer", 8); + REQUIRE(selector.accept(cellID)); +} + +TEST_CASE("CellIDSelector::accept empty selection") { + const auto selector = CellIDSelector{"system:8,side:-2,layer:5,module:7,sensor:10", {}}; + const auto cellID = GENERATE(take(100, random(0UL, std::numeric_limits::max()))); + REQUIRE_FALSE(selector.accept(cellID)); +} + +TEST_CASE("CellIDSelector::getSelectionMasks") { + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + + constexpr dd4hep::CellID systemMask = (0x0001ULL << 8) - 1; + constexpr dd4hep::CellID layerMask = ((0x0001ULL << 5) - 1) << (8 + 2); + + constexpr auto allValues = [](const auto& sel) { + namespace rv = std::ranges::views; + return rv::transform(sel, [](const auto& elem) { return elem.value; }); + }; + + auto sel = CellIDSelector::getSelectionMasks("system:8", encodingString); + REQUIRE(sel.size() == 1); + REQUIRE(sel[0].mask == systemMask); + REQUIRE(sel[0].value == 8); + + sel = CellIDSelector::getSelectionMasks("layer:2", encodingString); + REQUIRE(sel.size() == 1); + REQUIRE(sel[0].mask == layerMask); + REQUIRE(sel[0].value == ((2 << (8 + 2)))); + + sel = CellIDSelector::getSelectionMasks("system:4,layer:3", encodingString); + REQUIRE(sel.size() == 1); + REQUIRE(sel[0].mask == (systemMask | layerMask)); + REQUIRE(sel[0].value == (4 + (3 << (8 + 2)))); + + using Catch::Matchers::UnorderedRangeEquals; + + sel = CellIDSelector::getSelectionMasks("system:3|5", encodingString); + REQUIRE(sel.size() == 2); + for (const auto s : sel) { + REQUIRE(s.mask == systemMask); + } + REQUIRE_THAT(allValues(sel), UnorderedRangeEquals(std::vector{3, 5})); + + sel = CellIDSelector::getSelectionMasks("system:3|5,layer:1|12|10", encodingString); + REQUIRE(sel.size() == 6); + for (const auto s : sel) { + REQUIRE(s.mask == (systemMask | layerMask)); + } + REQUIRE_THAT(allValues(sel), UnorderedRangeEquals(std::vector{ + 3 + (1 << (8 + 2)), + 5 + (1 << (8 + 2)), + 3 + (10 << (8 + 2)), + 5 + (10 << (8 + 2)), + 3 + (12 << (8 + 2)), + 5 + (12 << (8 + 2)), + })); +} + +TEST_CASE("CellIDSelector failure modes: malformed selection strings") { + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + + // No colon separator: "field:value" format is required + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system5"}), std::invalid_argument); + + // Multiple colons make the split ambiguous (field name contains ":") + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:3:5"}), std::invalid_argument); + + // Non-integer value that cannot be parsed + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:abc"}), std::invalid_argument); + + // Empty value after colon + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:"}), std::invalid_argument); + + // Empty string as selection is not valid + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {""}), std::invalid_argument); + + // Pipe-separated list with a non-integer entry + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:1|abc|3"}), std::invalid_argument); + + // Pipe-separated list with an empty entry + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:1|"}), std::invalid_argument); + + // Multi-field: valid first field, non-integer value in second + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:3,layer:abc"}), std::invalid_argument); + + // Multi-field: valid first field, missing colon in second + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:3,layer3"}), std::invalid_argument); + + // Multi-field: trailing comma produces an empty token for the second field + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {"system:3,"}), std::invalid_argument); + + // Multi-field: leading comma produces an empty token for the first field + REQUIRE_THROWS_AS(CellIDSelector(encodingString, {",layer:3"}), std::invalid_argument); +} + +TEST_CASE("CellIDSelector failure modes: unknown field name") { + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + + // Field name not present in the encoding string + REQUIRE_THROWS(CellIDSelector(encodingString, {"nonexistent:5"})); + + // Unknown field among valid ones in the same selection + REQUIRE_THROWS(CellIDSelector(encodingString, {"system:3,typo:1"})); +} + +TEST_CASE("CellIDSelector failure modes: value overflow throws") { + // system is 8 bits wide (values 0-255). DD4hep's BitFieldCoder::set is + // range-checked: specifying 300 throws rather than truncating silently. + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + REQUIRE_THROWS(CellIDSelector(encodingString, {"system:300"})); + + // Same for a value in a pipe list that is in range — succeeds. + REQUIRE_NOTHROW(CellIDSelector(encodingString, {"system:255"})); + REQUIRE_NOTHROW(CellIDSelector(encodingString, {"system:0|255"})); + + // Out-of-range value inside a pipe list also throws. + REQUIRE_THROWS(CellIDSelector(encodingString, {"system:1|300"})); +} + +TEST_CASE("CellIDSelector failure modes: duplicate field uses last value") { + // "system:3,system:5" contains system twice. The Cartesian product expansion + // produces one combination {(system,3),(system,5)} and the second set() call + // overwrites the first, so the effective selector is system==5 only. + const std::string encodingString = "system:8,side:-2,layer:5,module:7,sensor:10"; + const auto encoder = dd4hep::BitFieldCoder(encodingString); + const auto selector = CellIDSelector{encodingString, {"system:3,system:5"}}; + + dd4hep::CellID cellID{0}; + encoder.set(cellID, "system", 5); + // system:5 matches because the second value wins + REQUIRE(selector.accept(cellID)); + + encoder.set(cellID, "system", 3); + // system:3 does NOT match even though it appears in the selection string + REQUIRE_FALSE(selector.accept(cellID)); +}