diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml index d3aa6f32c99..97f66117320 100644 --- a/.github/workflows/checks.yml +++ b/.github/workflows/checks.yml @@ -150,3 +150,13 @@ jobs: - name: Check run: > CI/check_fpe_masks.py --token ${{ secrets.GITHUB_TOKEN }} + unused_files: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: actions/setup-python@v4 + with: + python-version: '3.10' + - name: Check + run: > + CI/check_unused_files.py diff --git a/.github/workflows/iwyu.yml b/.github/workflows/iwyu.yml index de0f6c2c006..9ed9323f27e 100644 --- a/.github/workflows/iwyu.yml +++ b/.github/workflows/iwyu.yml @@ -55,7 +55,7 @@ jobs: -DACTS_BUILD_EXAMPLES_EDM4HEP=ON - name: Run IWYU - run: python3 iwyu-install/bin/iwyu_tool.py -p acts-build/ -j2 -- --mapping_file=acts/CI/iwyu/mapping.imp | tee iwyu-output.txt + run: python3 iwyu-install/bin/iwyu_tool.py -p acts-build/ -j2 -- -Xiwyu --mapping_file=acts/CI/iwyu/mapping.imp | tee iwyu-output.txt - name: Filter IWYU output run: python3 acts/CI/iwyu/filter.py acts/CI/iwyu/filter.yaml iwyu-output.txt iwyu-filtered.txt - name: Apply IWYU diff --git a/CI/check_unused_files.py b/CI/check_unused_files.py new file mode 100755 index 00000000000..0de9e8c5c12 --- /dev/null +++ b/CI/check_unused_files.py @@ -0,0 +1,228 @@ +#!/usr/bin/env python3 + +from pathlib import Path +import os +import sys +import subprocess + + +def file_can_be_removed(searchstring, scope): + cmd = "grep -IR '" + searchstring + "' " + " ".join(scope) + + p = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE) + output, _ = p.communicate() + return output == b"" + + +def count_files(path=".", exclude_dirs=(), exclude_files=()): + count = 0 + for root, dirs, files in os.walk(path): + dirs[:] = [d for d in dirs if d not in exclude_dirs] + files[:] = [f for f in files if f not in exclude_files] + count += len(files) + + return count + + +def main(): + print("\033[32mINFO\033[0m Start check_unused_files.py ...") + exclude_dirs = ( + "Scripts", + "thirdparty", + "CI", + "git", + "cmake", + ".git", + ".github", + ".", + ".idea", + ) + exclude_files = ( + "acts_logo_colored.svg", + ".gitignore", + "README.md", + "CMakeLists.txt", + "DetUtils.h", + "CommandLineArguments.h", + # Filename not completed in source + "vertexing_event_mu20_beamspot.csv", + "vertexing_event_mu20_tracks.csv", + "vertexing_event_mu20_vertices_AMVF.csv", + # TODO Move the following files to a better place? + "Magfield.ipynb", + "SolenoidField.ipynb", + # TODO Add README next to the following files? + "default-input-config-generic.json", + "geoSelection-openDataDetector.json", + "alignment-geo-contextualDetector.json", + ) + + suffix_header = ( + ".hpp", + ".cuh", + ) + suffix_source = ( + ".ipp", + ".cpp", + ".cu", + ) + suffix_image = ( + ".png", + ".svg", + ".jpg", + ".gif", + ) + suffix_python = (".py",) + suffix_doc = ( + ".md", + ".rst", + ) + suffix_other = ( + "", + ".csv", + ".css", + ".gdml", + ".hepmc3", + ".in", + ".ipynb", + ".json", + ".j2", + ".onnx", + ".root", + ".toml", + ".txt", + ".yml", + ) + suffix_allowed = ( + suffix_header + + suffix_source + + suffix_image + + suffix_python + + suffix_doc + + suffix_other + ) + + exit = 0 + + dirs_base = next(os.walk("."))[1] + dirs_base[:] = [d for d in dirs_base if d not in exclude_dirs] + dirs_base_docs = ("docs",) + dirs_base_code = [d for d in dirs_base if d not in dirs_base_docs] + + # Collectors + wrong_extension = () + unused_files = () + + # walk over all files + for root, dirs, files in os.walk("."): + dirs[:] = [d for d in dirs if d not in exclude_dirs] + files[:] = [f for f in files if f not in exclude_files] + + # Skip base-directory + if str(Path(root)) == ".": + continue + + # Print progress + if root[2:] in dirs_base: + processed_files = 0 + current_base_dir = root + number_files = count_files(root, exclude_dirs, exclude_files) + # print empty to start a new line + print("") + + # Skip "white-paper-figures" + # TODO Find a more elegant way + if str(root).find("white_papers/figures") != -1: + processed_files += count_files(root, exclude_dirs, exclude_files) + continue + + # Skip "DD4hep-tests" since their cmake looks a bit different + # TODO Find a more elegant way + if str(root).find("Tests/UnitTests/Plugins/DD4hep") != -1: + processed_files += count_files(root, exclude_dirs, exclude_files) + continue + + root = Path(root) + for filename in files: + processed_files += 1 + # get the full path of the file + filepath = root / filename + + # Check for wrong extensions + if filepath.suffix not in suffix_allowed: + wrong_extension += (str(filepath),) + + # Check header files and remove + if filepath.suffix in suffix_header + suffix_source: + if file_can_be_removed(filepath.stem, dirs_base_code): + unused_files += (str(filepath),) + remove_cmd = "rm " + str(filepath) + os.system(remove_cmd) + + # TODO Find test to check python files + if filepath.suffix in suffix_python: + continue + + # Check documentation files (weak tests) + # TODO find more reliable test for this + if filepath.suffix in suffix_doc: + if file_can_be_removed(filepath.stem, dirs_base_docs): + unused_files += (str(filepath),) + remove_cmd = "rm " + str(filepath) + os.system(remove_cmd) + + # Check and print other files + if filepath.suffix in suffix_image + suffix_other: + if file_can_be_removed(filename, dirs_base): + unused_files += (str(filepath),) + remove_cmd = "rm " + str(filepath) + os.system(remove_cmd) + + # Print the progress in place + progress = int(20 * processed_files / number_files) + sys.stdout.write("\r") + sys.stdout.write( + "Checked [%-20s] %d/%d files in %s" + % ("=" * progress, processed_files, number_files, current_base_dir) + ) + sys.stdout.flush() + + if len(wrong_extension) != 0: + print( + "\n\n\033[31mERROR\033[0m " + + f"The following {len(wrong_extension)} files have an unsupported extension:\n\n" + + "\033[31m" + + "\n".join(wrong_extension) + + "\033[0m" + + "\nCheck if you can change the format to one of the following:\n" + + "\n".join(suffix_allowed) + + "\nIf you really need that specific extension, add it to the list above.\n" + ) + + exit += 1 + + if len(unused_files) != 0: + print( + "\n\n\033[31mERROR\033[0m " + + f"The following {len(unused_files)} files seem to be unused:\n" + + "\033[31m" + + "\n".join(unused_files) + + "\033[0m" + + "\nYou have 3 options:" + + "\n\t- Remove them" + + "\n\t- Use them (check proper include)" + + "\n\t- Modify the ignore list of this check\n" + ) + + exit += 1 + + if exit == 0: + print( + "\n\n\033[32mINFO\033[0m Finished check_unused_files.py without any errors." + ) + + return exit + + +if "__main__" == __name__: + sys.exit(main()) diff --git a/CI/iwyu/README.md b/CI/iwyu/README.md new file mode 100644 index 00000000000..6204d806f20 --- /dev/null +++ b/CI/iwyu/README.md @@ -0,0 +1,16 @@ +# IWYU = Include what you use + +This tool finds unused includes and suggestes additional includes and declarations. + +It is not very stable at the moment and takes a few hours to complete within ACTS therefor it only runs once a week and can be triggered manually. + +There is also not sufficient filtering offered by IWYU at the moment. For that reason there is a specific filtering script for now which tries to get rid of unwanted changes. + +offline resources +- [GitHub Actions Workflow](../../.github/workflows/iwyu.yml) +- [Custom filter script](./filter.py) + +online resources + +- https://include-what-you-use.org/ +- https://github.com/include-what-you-use/include-what-you-use diff --git a/CI/iwyu/filter.yaml b/CI/iwyu/filter.yaml index 8571a2bceb1..07d3aaf643b 100644 --- a/CI/iwyu/filter.yaml +++ b/CI/iwyu/filter.yaml @@ -1,5 +1,5 @@ remove_lines: - # ignore unexiting std + # ignore if not existing in std - "^#include " - "^#include ": "#include " - "^#include ": "#include " - "^#include ": "#include " - # don use ipp + # don't use ipp - "^#include ([<\"].*)\\.ipp": "#include \\1.hpp" ignore_files: diff --git a/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root b/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root index fac29773401..9f9d071a76a 100644 Binary files a/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root and b/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root index 56555864be8..505e1dc7daa 100644 Binary files a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_orthogonal_hist.root b/CI/physmon/reference/performance_amvf_orthogonal_hist.root index b3a6f0f2279..d9e301ffc4d 100644 Binary files a/CI/physmon/reference/performance_amvf_orthogonal_hist.root and b/CI/physmon/reference/performance_amvf_orthogonal_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_seeded_hist.root b/CI/physmon/reference/performance_amvf_seeded_hist.root index 49b4e17dbda..599dc1d7944 100644 Binary files a/CI/physmon/reference/performance_amvf_seeded_hist.root and b/CI/physmon/reference/performance_amvf_seeded_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_truth_estimated_hist.root b/CI/physmon/reference/performance_amvf_truth_estimated_hist.root index 101c8c5fdc9..2f96311fff2 100644 Binary files a/CI/physmon/reference/performance_amvf_truth_estimated_hist.root and b/CI/physmon/reference/performance_amvf_truth_estimated_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_truth_smeared_hist.root b/CI/physmon/reference/performance_amvf_truth_smeared_hist.root index e4e9cc4723b..b5afc584b60 100644 Binary files a/CI/physmon/reference/performance_amvf_truth_smeared_hist.root and b/CI/physmon/reference/performance_amvf_truth_smeared_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_ttbar_hist.root b/CI/physmon/reference/performance_amvf_ttbar_hist.root index 60504fc9b18..0d33046b84a 100644 Binary files a/CI/physmon/reference/performance_amvf_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_ttbar_hist.root differ diff --git a/CI/physmon/reference/performance_ivf_orthogonal_hist.root b/CI/physmon/reference/performance_ivf_orthogonal_hist.root index 5ddabf0c0f1..16a242b7c5b 100644 Binary files a/CI/physmon/reference/performance_ivf_orthogonal_hist.root and b/CI/physmon/reference/performance_ivf_orthogonal_hist.root differ diff --git a/CI/physmon/reference/performance_ivf_seeded_hist.root b/CI/physmon/reference/performance_ivf_seeded_hist.root index ff5a7052db4..1dfaabb0dcf 100644 Binary files a/CI/physmon/reference/performance_ivf_seeded_hist.root and b/CI/physmon/reference/performance_ivf_seeded_hist.root differ diff --git a/CI/physmon/reference/performance_ivf_truth_estimated_hist.root b/CI/physmon/reference/performance_ivf_truth_estimated_hist.root index 303ce9ab961..5e44bb77359 100644 Binary files a/CI/physmon/reference/performance_ivf_truth_estimated_hist.root and b/CI/physmon/reference/performance_ivf_truth_estimated_hist.root differ diff --git a/CI/physmon/reference/performance_ivf_truth_smeared_hist.root b/CI/physmon/reference/performance_ivf_truth_smeared_hist.root index 7e2b13fd5dd..48c2e94a2de 100644 Binary files a/CI/physmon/reference/performance_ivf_truth_smeared_hist.root and b/CI/physmon/reference/performance_ivf_truth_smeared_hist.root differ diff --git a/CI/physmon/summary.py b/CI/physmon/summary.py index ae20ddfce4f..01aead939da 100755 --- a/CI/physmon/summary.py +++ b/CI/physmon/summary.py @@ -7,7 +7,7 @@ import os import csv -HERALD_URL = "https://herald.dokku.paulgessinger.com/view/{repo}/runs/{run_id}/artifacts/{artifact_name}/{path}" +HERALD_URL = "https://acts-herald.app.cern.ch/view/{repo}/runs/{run_id}/artifacts/{artifact_name}/{path}" IS_CI = "GITHUB_ACTIONS" in os.environ diff --git a/CI/physmon/vertexing_config.yml b/CI/physmon/vertexing_config.yml index 8ed875b686e..53c76240a10 100644 --- a/CI/physmon/vertexing_config.yml +++ b/CI/physmon/vertexing_config.yml @@ -84,6 +84,11 @@ histograms: nbins: 100 min: 0.999 max: 1 + + "trk_weight": + nbins: 100 + min: -0.01 + max: 1.01 extra_histograms: - expression: df["nRecoVtx"] / df["nTrueVtx"] diff --git a/Core/include/Acts/Detector/Portal.hpp b/Core/include/Acts/Detector/Portal.hpp index bde88aa523f..a38677cb7a3 100644 --- a/Core/include/Acts/Detector/Portal.hpp +++ b/Core/include/Acts/Detector/Portal.hpp @@ -41,14 +41,13 @@ struct NavigationState; /// The surface can carry material to allow mapping onto /// portal positions if required. /// -class Portal : public std::enable_shared_from_this { - protected: +class Portal { + public: /// Constructor from surface w/o portal links /// /// @param surface is the representing surface Portal(std::shared_ptr surface); - public: /// The volume links forward/backward with respect to the surface normal using DetectorVolumeUpdaters = std::array; @@ -60,34 +59,7 @@ class Portal : public std::enable_shared_from_this { /// Declare the DetectorVolume friend for portal setting friend class DetectorVolume; - /// Factory for producing memory managed instances of Portal. - static std::shared_ptr makeShared( - std::shared_ptr surface); - - /// Retrieve a @c std::shared_ptr for this surface (non-const version) - /// - /// @note Will error if this was not created through the @c makeShared factory - /// since it needs access to the original reference. In C++14 this is - /// undefined behavior (but most likely implemented as a @c bad_weak_ptr - /// exception), in C++17 it is defined as that exception. - /// @note Only call this if you need shared ownership of this object. - /// - /// @return The shared pointer - std::shared_ptr getSharedPtr(); - - /// Retrieve a @c std::shared_ptr for this surface (const version) - /// - /// @note Will error if this was not created through the @c makeShared factory - /// since it needs access to the original reference. In C++14 this is - /// undefined behavior, but most likely implemented as a @c bad_weak_ptr - /// exception, in C++17 it is defined as that exception. - /// @note Only call this if you need shared ownership of this object. - /// - /// @return The shared pointer - std::shared_ptr getSharedPtr() const; - Portal() = delete; - virtual ~Portal() = default; /// Const access to the surface representation const RegularSurface& surface() const; @@ -110,18 +82,22 @@ class Portal : public std::enable_shared_from_this { /// Fuse with another portal, this one is kept /// - /// @param other is the portal that will be fused + /// @param aPortal is the first portal to fuse + /// @param bPortal is the second portal to fuse /// - /// @note this will move the portal links from the other - /// into this volume, it will throw an exception if the + /// @note this will combine the portal links from the both + /// portals into a new one, it will throw an exception if the /// portals are not fusable /// /// @note if one portal carries material, it will be kept, /// however, if both portals carry material, an exception /// will be thrown and the portals are not fusable /// - /// @note that other will be overwritten to point to this - void fuse(std::shared_ptr& other) noexcept(false); + /// @note Both input portals become invalid, in that their update + /// delegates and attached volumes are reset + static std::shared_ptr fuse( + std::shared_ptr& aPortal, + std::shared_ptr& bPortal) noexcept(false); /// Update the volume link /// diff --git a/Core/include/Acts/Geometry/TrackingVolume.hpp b/Core/include/Acts/Geometry/TrackingVolume.hpp index df88e9f065c..62ef13e74dd 100644 --- a/Core/include/Acts/Geometry/TrackingVolume.hpp +++ b/Core/include/Acts/Geometry/TrackingVolume.hpp @@ -69,17 +69,19 @@ using LayerArray = BinnedArray; using LayerVector = std::vector; /// Intersection with @c Layer -using LayerIntersection = ObjectIntersection; +using LayerIntersection = std::pair; /// Multi-intersection with @c Layer -using LayerMultiIntersection = ObjectMultiIntersection; +using LayerMultiIntersection = + std::pair; /// BoundarySurface of a volume using BoundarySurface = BoundarySurfaceT; /// Intersection with a @c BoundarySurface -using BoundaryIntersection = ObjectIntersection; +using BoundaryIntersection = + std::pair; /// Multi-intersection with a @c BoundarySurface using BoundaryMultiIntersection = - ObjectMultiIntersection; + std::pair; /// @class TrackingVolume /// diff --git a/Core/include/Acts/Material/MaterialCollector.hpp b/Core/include/Acts/Material/MaterialCollector.hpp deleted file mode 100644 index f7b523549b4..00000000000 --- a/Core/include/Acts/Material/MaterialCollector.hpp +++ /dev/null @@ -1,145 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018-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/. - -#pragma once - -#include "Acts/Material/ISurfaceMaterial.hpp" -#include "Acts/Material/Material.hpp" -#include "Acts/Material/MaterialSlab.hpp" -#include "Acts/Surfaces/Surface.hpp" - -#include - -namespace Acts { - -/// The information to be writtern out per hit surface -struct MaterialHit { - const Surface* surface = nullptr; - Vector3 position; - Vector3 direction; - Material material; - double pathLength; -}; - -/// A Material Collector struct -struct MaterialCollector { - /// In the detailed collection mode the material - /// per surface is collected, otherwise only the total - /// pathlength in X0 or L0 are recorded - bool detailedCollection = false; - - /// Simple result struct to be returned - /// - /// Result of the material collection process - /// It collects the overall X0 and L0 path lengths, - /// and optionally a detailed per-material breakdown - struct this_result { - std::vector collected; - double materialInX0 = 0.; - double materialInL0 = 0.; - }; - - using result_type = this_result; - - /// Collector action for the ActionList of the Propagator - /// It checks if the state has a current surface, - /// in which case the action is performed: - /// - it records the surface given the configuration - /// - /// @tparam propagator_state_t is the type of Propagator state - /// @tparam stepper_t Type of the stepper of the propagation - /// @tparam navigator_t Type of the navigator of the propagation - /// - /// @param state is the mutable propagator state object - /// @param stepper The stepper in use - /// @param navigator The navigator in use - /// @param result is the result object to be filled - /// @param logger a logger instance - template - void operator()(propagator_state_t& state, const stepper_t& stepper, - const navigator_t& navigator, result_type& result, - const Logger& logger) const { - if (navigator.currentSurface(state.navigation)) { - if (navigator.currentSurface(state.navigation) == - navigator.targetSurface(state.navigation) && - !navigator.targetReached(state.navigation)) { - return; - } - - ACTS_VERBOSE("Material check on surface " - << navigator.currentSurface(state.navigation)->geometryId()); - - if (navigator.currentSurface(state.navigation)->surfaceMaterial()) { - // get the material propertices and only continue - const MaterialSlab* mProperties = - navigator.currentSurface(state.navigation) - ->surfaceMaterial() - ->material(stepper.position(state.stepping)); - if (mProperties) { - // pre/post/full update - double prepofu = 1.; - if (navigator.startSurface(state.navigation) == - navigator.currentSurface(state.navigation)) { - ACTS_VERBOSE("Update on start surface: post-update mode."); - prepofu = navigator.currentSurface(state.navigation) - ->surfaceMaterial() - ->factor(state.options.direction, - MaterialUpdateStage::PostUpdate); - } else if (navigator.targetSurface(state.navigation) == - navigator.currentSurface(state.navigation)) { - ACTS_VERBOSE("Update on target surface: pre-update mode"); - prepofu = navigator.currentSurface(state.navigation) - ->surfaceMaterial() - ->factor(state.options.direction, - MaterialUpdateStage::PreUpdate); - } else { - ACTS_VERBOSE("Update while pass through: full mode."); - } - - // the pre/post factor has been applied - // now check if there's still something to do - if (prepofu == 0.) { - ACTS_VERBOSE("Pre/Post factor set material to zero."); - return; - } - // more debugging output to the screen - ACTS_VERBOSE("Material properties found for this surface."); - - // the path correction from the surface intersection - double pCorrection = - prepofu * navigator.currentSurface(state.navigation) - ->pathCorrection(stepper.position(state.stepping), - stepper.direction(state.stepping)); - // the full material - result.materialInX0 += pCorrection * mProperties->thicknessInX0(); - result.materialInL0 += pCorrection * mProperties->thicknessInL0(); - - ACTS_VERBOSE("t/X0 (t/L0) increased to " - << result.materialInX0 << " (" << result.materialInL0 - << " )"); - - // if configured, record the individual material hits - if (detailedCollection) { - // create for recording - MaterialHit mHit; - mHit.surface = navigator.currentSurface(state.navigation); - mHit.position = stepper.position(state.stepping); - mHit.direction = stepper.direction(state.stepping); - // get the material & path length - mHit.material = mProperties->material(); - mHit.pathLength = pCorrection * mProperties->thickness(); - // save if in the result - result.collected.push_back(mHit); - } - } - } - } - } -}; // namespace Acts -} // namespace Acts diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp index 2e9cb6e206f..a49176d0297 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp @@ -662,7 +662,7 @@ class MultiEigenStepperLoop template void updateStepSize(State& state, const object_intersection_t& oIntersection, Direction direction, bool release = true) const { - const Surface& surface = *oIntersection.representation(); + const Surface& surface = *oIntersection.object(); for (auto& component : state.components) { auto intersection = surface.intersect( diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp index 46126777d01..158f0a7f185 100644 --- a/Core/include/Acts/Propagator/Navigator.hpp +++ b/Core/include/Acts/Propagator/Navigator.hpp @@ -525,7 +525,7 @@ class Navigator { state.navigation.lastHierarchySurfaceReached = false; // Update volume information // get the attached volume information - auto boundary = state.navigation.navBoundary().object(); + auto boundary = state.navigation.navBoundary().second; state.navigation.currentVolume = boundary->attachedVolume( state.geoContext, stepper.position(state.stepping), stepper.direction(state.stepping), state.options.direction); @@ -568,6 +568,19 @@ class Navigator { } private: + const SurfaceIntersection& candidateIntersection( + const NavigationSurfaces& surfaces, std::size_t index) const { + return surfaces.at(index); + } + const SurfaceIntersection& candidateIntersection( + const NavigationLayers& surfaces, std::size_t index) const { + return surfaces.at(index).first; + } + const SurfaceIntersection& candidateIntersection( + const NavigationBoundaries& surfaces, std::size_t index) const { + return surfaces.at(index).first; + } + /// @brief Status call for test surfaces (surfaces, layers, boundaries) /// /// If there are surfaces to be handled, check if the current @@ -592,9 +605,9 @@ class Navigator { if (navSurfaces.empty() || navIndex == navSurfaces.size()) { return false; } - const auto& intersection = navSurfaces.at(navIndex); + const auto& intersection = candidateIntersection(navSurfaces, navIndex); // Take the current surface - auto surface = intersection.representation(); + const auto* surface = intersection.object(); // Check if we are at a surface // If we are on the surface pointed at by the index, we can make // it the current one to pass it to the other actors @@ -844,9 +857,9 @@ class Navigator { // loop over the available navigation layer candidates while (state.navigation.navLayerIndex != state.navigation.navLayers.size()) { - const auto& intersection = state.navigation.navLayer(); + const auto& intersection = state.navigation.navLayer().first; // The layer surface - const auto* layerSurface = intersection.representation(); + const auto* layerSurface = intersection.object(); // We are on the layer if (state.navigation.currentSurface == layerSurface) { ACTS_VERBOSE(volInfo(state) << "We are on a layer, resolve Surfaces."); @@ -976,7 +989,7 @@ class Navigator { os << state.navigation.navBoundaries.size(); os << " boundary candidates found at path(s): "; for (auto& bc : state.navigation.navBoundaries) { - os << bc.pathLength() << " "; + os << bc.first.pathLength() << " "; } logger().log(Logging::VERBOSE, os.str()); } @@ -984,7 +997,8 @@ class Navigator { state.navigation.navBoundaryIndex = 0; if (!state.navigation.navBoundaries.empty()) { // Set to the first and return to the stepper - stepper.updateStepSize(state.stepping, state.navigation.navBoundary(), + stepper.updateStepSize(state.stepping, + state.navigation.navBoundary().first, state.options.direction, true); ACTS_VERBOSE(volInfo(state) << "Navigation stepSize updated to " << stepper.outputStepSize(state.stepping)); @@ -1001,9 +1015,9 @@ class Navigator { // Loop over the boundary surface while (state.navigation.navBoundaryIndex != state.navigation.navBoundaries.size()) { - const auto& intersection = state.navigation.navBoundary(); + const auto& intersection = state.navigation.navBoundary().first; // That is the current boundary surface - const auto* boundarySurface = intersection.representation(); + const auto* boundarySurface = intersection.object(); // Step towards the boundary surfrace auto boundaryStatus = stepper.updateSurfaceStatus( state.stepping, *boundarySurface, intersection.index(), @@ -1133,8 +1147,8 @@ class Navigator { const Layer* cLayer = nullptr) const { // get the layer and layer surface auto layerSurface = cLayer ? state.navigation.startSurface - : state.navigation.navLayer().representation(); - auto navLayer = cLayer ? cLayer : state.navigation.navLayer().object(); + : state.navigation.navLayer().first.object(); + auto navLayer = cLayer ? cLayer : state.navigation.navLayer().second; // are we on the start layer bool onStart = (navLayer == state.navigation.startLayer); auto startSurface = onStart ? state.navigation.startSurface : layerSurface; @@ -1248,7 +1262,7 @@ class Navigator { os << state.navigation.navLayers.size(); os << " layer candidates found at path(s): "; for (auto& lc : state.navigation.navLayers) { - os << lc.pathLength() << " "; + os << lc.first.pathLength() << " "; } logger().log(Logging::VERBOSE, os.str()); } @@ -1256,10 +1270,11 @@ class Navigator { state.navigation.navLayerIndex = 0; // Setting the step size towards first if (state.navigation.startLayer && - state.navigation.navLayer().object() != state.navigation.startLayer) { + state.navigation.navLayer().second != state.navigation.startLayer) { ACTS_VERBOSE(volInfo(state) << "Target at layer."); // The stepper updates the step size ( single / multi component) - stepper.updateStepSize(state.stepping, state.navigation.navLayer(), + stepper.updateStepSize(state.stepping, + state.navigation.navLayer().first, state.options.direction, true); ACTS_VERBOSE(volInfo(state) << "Navigation stepSize updated to " << stepper.outputStepSize(state.stepping)); diff --git a/Core/include/Acts/Seeding/GNN_DataStorage.hpp b/Core/include/Acts/Seeding/GNN_DataStorage.hpp index fb630b43a3f..9b0ca47ba76 100644 --- a/Core/include/Acts/Seeding/GNN_DataStorage.hpp +++ b/Core/include/Acts/Seeding/GNN_DataStorage.hpp @@ -180,10 +180,6 @@ class TrigFTF_GNN_DataStorage { float max_tau = 100.0; // can't do this bit yet as dont have cluster width if (useClusterWidth) { - // const Trk::SpacePoint* osp = sp.offlineSpacePoint(); - // const InDet::PixelCluster* pCL = dynamic_cast(osp->clusterList().first); - // float cluster_width = pCL->width().widthPhiRZ().y(); float cluster_width = 1; // temporary while cluster width not available min_tau = 6.7 * (cluster_width - 0.2); max_tau = @@ -194,10 +190,6 @@ class TrigFTF_GNN_DataStorage { sp, min_tau, max_tau)); // adding ftf member to nodes } else { if (useClusterWidth) { - // const Trk::SpacePoint* osp = sp.offlineSpacePoint(); - // const InDet::PixelCluster* pCL = dynamic_cast(osp->clusterList().first); - // float cluster_width = pCL->width().widthPhiRZ().y(); float cluster_width = 1; // temporary while cluster width not available if (cluster_width > 0.2) { return -3; diff --git a/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp b/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp index 5c388c83504..a4a4fe52af2 100644 --- a/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp +++ b/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp @@ -229,8 +229,8 @@ class TrigFTF_GNN_TrackingFilter { const float weight_x = 0.5; const float weight_y = 0.5; - const float maxDChi2_x = 60.0; // 35.0; - const float maxDChi2_y = 60.0; // 31.0; + const float maxDChi2_x = 60.0; // was 35.0; + const float maxDChi2_y = 60.0; // was 31.0; const float add_hit = 14.0; @@ -306,7 +306,6 @@ class TrigFTF_GNN_TrackingFilter { Cy[1][1] = ts.m_Cy[1][1]; // chi2 test - float resid_x = mx - X[0]; float resid_y = my - Y[0]; @@ -317,7 +316,7 @@ class TrigFTF_GNN_TrackingFilter { int type = getLayerType(pS->m_n1->m_sp_FTF.combined_ID); - if (type == 0) { // barrel TO-DO: split into barrel Pixel and barrel SCT + if (type == 0) { // barrel TODO: split into barrel Pixel and barrel SCT sigma_rz = sigma_y * sigma_y; } else { sigma_rz = sigma_y * ts.m_Y[1]; @@ -338,7 +337,6 @@ class TrigFTF_GNN_TrackingFilter { ts.m_J += add_hit - dchi2_x * weight_x - dchi2_y * weight_y; // state update - float Kx[3] = {Dx * Cx[0][0], Dx * Cx[0][1], Dx * Cx[0][2]}; float Ky[2] = {Dy * Cy[0][0], Dy * Cy[0][1]}; @@ -360,6 +358,7 @@ class TrigFTF_GNN_TrackingFilter { ts.m_Cy[i][j] = Cy[i][j] - Ky[i] * CHy[j]; } } + ts.m_refX = refX; ts.m_refY = refY; return true; @@ -371,7 +370,7 @@ class TrigFTF_GNN_TrackingFilter { }); // iterator to vector member with this id int index = std::distance(m_geo.begin(), iterator); - return m_geo.at(index).m_type; + return m_geo.at(index).m_type; // needs to be 0, 2, or -2 } const std::vector& m_geo; diff --git a/Core/include/Acts/Seeding/SeedConfirmationRangeConfig.hpp b/Core/include/Acts/Seeding/SeedConfirmationRangeConfig.hpp index 45be4a4f98a..e901fcfa0e4 100644 --- a/Core/include/Acts/Seeding/SeedConfirmationRangeConfig.hpp +++ b/Core/include/Acts/Seeding/SeedConfirmationRangeConfig.hpp @@ -15,36 +15,50 @@ namespace Acts { -// defaults experimental cuts to no operation in both seeding algorithms +/// defaults experimental cuts to no operation in both seeding algorithms inline bool noopExperimentCuts(float /*bottomRadius*/, float /*cotTheta*/) { return true; } -/// @brief contains parameters for seed confirmation +/// @brief Contains parameters for quality seed confirmation +/// @note Requirements on the number of compatible space-points and impact parameters can be defined +/// for different (r, z) regions of the detector (e.g. forward or central +/// region) by SeedConfirmationRange. Seeds are classified as "high-quality" +/// seeds and normal quality seeds. Normal quality seeds are only selected if +/// no other "high-quality" seed has been found for that inner-middle doublet. +/// For optimization reasons, the algorithm only calls the seed confirmation +/// for a certain inner-middle doublet, in case a configurable minimum number +/// of inner-middle-outer triplets have been found. struct SeedConfirmationRangeConfig { - // z minimum and maximum of middle component of the seed used to define the - // region of the detector for seed confirmation + /// Minimum z position of middle component of the seed used to + /// split the region of the detector for seed confirmation float zMinSeedConf = std::numeric_limits::lowest(); // Acts::UnitConstants::mm + + /// Maximum z position of middle component of the seed used to + /// split the region of the detector for seed confirmation float zMaxSeedConf = std::numeric_limits::max(); // Acts::UnitConstants::mm - // radius of bottom component of seed that is used to define the number of - // compatible top required + + /// Radius position of inner seed component that is used to + /// split the region of the detector for seed confirmation float rMaxSeedConf = std::numeric_limits::max(); // Acts::UnitConstants::mm - // number of compatible top SPs of seed if bottom radius is larger than - // rMaxSeedConf + /// Minimum number of compatible outer space-points required in quality seed + /// confirmation if inner space-points radius is larger than rMaxSeedConf std::size_t nTopForLargeR = 0; - // number of compatible top SPs of seed if bottom radius is smaller than - // rMaxSeedConf + /// Minimum number of compatible outer space-points required in quality seed + /// confirmation if inner space-points radius is smaller than rMaxSeedConf std::size_t nTopForSmallR = 0; - // minimum radius for bottom SP in seed confirmation + /// Minimum radius for inner seed component required in quality seed + /// confirmation float seedConfMinBottomRadius = 60. * Acts::UnitConstants::mm; - // maximum zOrigin in seed confirmation + /// Maximum longitudinal impact parameter of seed required in quality seed + /// confirmation float seedConfMaxZOrigin = 150. * Acts::UnitConstants::mm; - // minimum impact parameter for seed confirmation + /// Minimum impact parameter of seed required in quality seed confirmation float minImpactSeedConf = 1. * Acts::UnitConstants::mm; }; diff --git a/Core/include/Acts/Seeding/SeedFilterConfig.hpp b/Core/include/Acts/Seeding/SeedFilterConfig.hpp index bd700b5a4c8..879ad372652 100644 --- a/Core/include/Acts/Seeding/SeedFilterConfig.hpp +++ b/Core/include/Acts/Seeding/SeedFilterConfig.hpp @@ -16,55 +16,73 @@ #include namespace Acts { + +/// @brief Structure that holds configuration parameters for the seed filter algorithm struct SeedFilterConfig { - // the allowed delta between two inverted seed radii for them to be considered - // compatible. + /// Allowed difference in curvature (inverted seed radii) between two + /// compatible seeds float deltaInvHelixDiameter = 0.00003 * 1. / Acts::UnitConstants::mm; - // the transverse impact parameters (d0) is multiplied by this factor and - // subtracted from weight + /// Minimum distance between compatible outer space-points to be considered. + /// This is used to avoid counting space-points from the same layer + float deltaRMin = 5. * Acts::UnitConstants::mm; + /// Seed weight/score is increased by this value if a compatible seed has been + /// found. This is the c1 factor in the seed score calculation (w = c1 * Nt - + /// c2 * d0 - c3 * z0) + float compatSeedWeight = 200.; + /// The transverse impact parameters (d0) is multiplied by this factor and + /// subtracted from weight. This is the c2 factor in the seed score + /// calculation (w = c1 * Nt - c2 * d0 - c3 * z0) float impactWeightFactor = 1.; - // the logitudinal impact parameters (z0) is multiplied by this factor and - // subtracted from weight + /// The logitudinal impact parameters (z0) is multiplied by this factor and + /// subtracted from weight. This is the c3 factor in the seed score + /// calculation (w = c1 * Nt - c2 * d0 - c3 * z0) float zOriginWeightFactor = 1.; - // seed weight increased by this value if a compatible seed has been found. - float compatSeedWeight = 200.; - // minimum distance between compatible seeds to be considered for weight boost - float deltaRMin = 5. * Acts::UnitConstants::mm; - // in dense environments many seeds may be found per middle space point. - // only seeds with the highest weight will be kept if this limit is reached. + /// Maximum number (minus one) of accepted seeds per middle space-point + /// In dense environments many seeds may be found per middle space-point + /// Only seeds with the highest weight will be kept if this limit is reached unsigned int maxSeedsPerSpM = 10; - // how often do you want to increase the weight of a seed for finding a - // compatible seed? + /// Maximum limit to number of compatible space-point used in score + /// calculation. We increase by c1 the weight calculation for each compatible + /// space-point until we reach compatSeedLimit std::size_t compatSeedLimit = 2; - // Tool to apply experiment specific cuts on collected middle space points - // increment in seed weight if the number of compatible seeds is larger than - // numSeedIncrement, this is used in case of high occupancy scenarios if we - // want to increase the weight of the seed by seedWeightIncrement when the - // number of compatible seeds is higher than a certain value + /// Increment in seed weight if the number of compatible seeds is larger than + /// numSeedIncrement, this is used in case of high occupancy scenarios if we + /// want to increase the weight of the seed by seedWeightIncrement when the + /// number of compatible seeds is higher than a certain value float seedWeightIncrement = 0; float numSeedIncrement = std::numeric_limits::infinity(); - // seedConfirmation enables seed confirmation cuts - keep seeds if they have - // specific values of impact parameter, z-origin and number of compatible - // seeds inside a pre-defined range that also depends on the region of the - // detector (i.e. forward or central region) defined by SeedConfirmationRange + /// Seeding parameters used for quality seed confirmation + + /// Enable quality seed confirmation, this is different than default seeding + /// confirmation because it can also be defined for different (r, z) regions + /// of the detector (e.g. forward or central region) by SeedConfirmationRange. + /// Seeds are classified as "high-quality" seeds and normal quality seeds. + /// Normal quality seeds are only selected if no other "high-quality" seed + /// has been found for that inner-middle doublet. bool seedConfirmation = false; - // contains parameters for central seed confirmation + /// Contains parameters for central seed confirmation SeedConfirmationRangeConfig centralSeedConfirmationRange; - // contains parameters for forward seed confirmation + /// Contains parameters for forward seed confirmation SeedConfirmationRangeConfig forwardSeedConfirmationRange; - // maximum number of lower quality seeds in seed confirmation + /// If seedConfirmation is true we classify seeds as "high-quality" seeds. + /// Seeds that are not confirmed as "high-quality" are only selected if no + /// other "high-quality" seed has been found for that inner-middle doublet + /// Maximum number of normal seeds (not classified as "high-quality" seeds) + /// in seed confirmation std::size_t maxSeedsPerSpMConf = std::numeric_limits::max(); - // maximum number of quality seeds for each middle-bottom SP-dublet in seed - // confirmation if the limit is reached we check if there is a lower quality - // seed to be replaced + /// Maximum number of "high-quality" seeds for each inner-middle SP-dublet in + /// seed confirmation. If the limit is reached we check if there is a normal + /// quality seed to be replaced std::size_t maxQualitySeedsPerSpMConf = std::numeric_limits::max(); - // use deltaR between top and middle SP instead of top radius to search for - // compatible SPs + /// Other parameters + + /// Use deltaR between top and middle SP instead of top radius to search for + /// compatible SPs bool useDeltaRorTopRadius = false; bool isInInternalUnits = false; diff --git a/Core/include/Acts/Seeding/SeedFinderConfig.hpp b/Core/include/Acts/Seeding/SeedFinderConfig.hpp index aaf50d6fcbd..b4645241b48 100644 --- a/Core/include/Acts/Seeding/SeedFinderConfig.hpp +++ b/Core/include/Acts/Seeding/SeedFinderConfig.hpp @@ -23,155 +23,185 @@ namespace Acts { template class SeedFilter; +/// @brief Structure that holds configuration parameters for the seed finder algorithm template struct SeedFinderConfig { std::shared_ptr> seedFilter; - // Seed Cuts - // lower cutoff for seeds - float minPt = 400. * Acts::UnitConstants::MeV; - // cot of maximum theta angle - // equivalent to 3 eta (pseudorapidity) - float cotThetaMax = 10.01788; - // minimum distance in r between two measurements within one seed - float deltaRMin = 5 * Acts::UnitConstants::mm; - // maximum distance in r between two measurements within one seed - float deltaRMax = 270 * Acts::UnitConstants::mm; - // minimum distance in r between middle and top SP - float deltaRMinTopSP = std::numeric_limits::quiet_NaN(); - // maximum distance in r between middle and top SP - float deltaRMaxTopSP = std::numeric_limits::quiet_NaN(); - // minimum distance in r between middle and bottom SP - float deltaRMinBottomSP = std::numeric_limits::quiet_NaN(); - // maximum distance in r between middle and bottom SP - float deltaRMaxBottomSP = std::numeric_limits::quiet_NaN(); - // radial bin size for filling space point grid + /// Seeding parameters used in the space-point grid creation and bin finding + + /// Geometry Settings + Detector ROI + /// (r, z, phi) range for limiting location of all measurements and grid + /// creation + float phiMin = -M_PI; + float phiMax = M_PI; + float zMin = -2800 * Acts::UnitConstants::mm; + float zMax = 2800 * Acts::UnitConstants::mm; + float rMax = 600 * Acts::UnitConstants::mm; + /// WARNING: if rMin is smaller than impactMax, the bin size will be 2*pi, + /// which will make seeding very slow! + float rMin = 33 * Acts::UnitConstants::mm; + + /// Vector containing the z-bin edges for non equidistant binning in z + std::vector zBinEdges; + + /// Number of z bins to skip during the search for middle space-points. This + /// is useful for removing the outermost bins and speeding up seeding. Should + /// be used in conjunction with zBinsCustomLooping (skipZMiddleBinSearch + /// determines the first N bins in zBinsCustomLooping to be avoided). + std::size_t skipZMiddleBinSearch = 0; + /// Order of z bins to loop over when searching for SPs + std::vector zBinsCustomLooping = {}; + + /// Radial bin size used in space-point grid float binSizeR = 1. * Acts::UnitConstants::mm; - // radial range for middle SP - // variable range based on SP radius + /// Seeding parameters used to define the region of interest for middle + /// space-point + + /// Radial range for middle space-point + /// The range can be defined manually with (rMinMiddle, rMaxMiddle). If + /// useVariableMiddleSPRange is set to false and the vector rRangeMiddleSP is + /// empty, we use (rMinMiddle, rMaxMiddle) to cut the middle space-points + float rMinMiddle = 60.f * Acts::UnitConstants::mm; + float rMaxMiddle = 120.f * Acts::UnitConstants::mm; + /// If useVariableMiddleSPRange is set to false, the vector rRangeMiddleSP can + /// be used to define a fixed r range for each z bin: {{rMin, rMax}, ...} bool useVariableMiddleSPRange = false; + /// Range defined in vector for each z bin + std::vector> rRangeMiddleSP; + /// If useVariableMiddleSPRange is true, the radial range will be calculated + /// based on the maximum and minimum r values of the space-points in the event + /// and a deltaR (deltaRMiddleMinSPRange, deltaRMiddleMaxSPRange) float deltaRMiddleMinSPRange = 10. * Acts::UnitConstants::mm; float deltaRMiddleMaxSPRange = 10. * Acts::UnitConstants::mm; - // range defined in vector for each z region - std::vector> rRangeMiddleSP; - // range defined by rMinMiddle and rMaxMiddle - float rMinMiddle = 60.f * Acts::UnitConstants::mm; - float rMaxMiddle = 120.f * Acts::UnitConstants::mm; - // z of last layers to avoid iterations + /// Vector containing minimum and maximum z boundaries for cutting middle + /// space-points std::pair zOutermostLayers{-2700 * Acts::UnitConstants::mm, 2700 * Acts::UnitConstants::mm}; - // cut to the maximum value of delta z between SPs - float deltaZMax = - std::numeric_limits::infinity() * Acts::UnitConstants::mm; + /// Seeding parameters used to define the cuts on space-point doublets - // enable cut on the compatibility between interaction point and SPs - bool interactionPointCut = false; + /// Minimum radial distance between two doublet components (prefer + /// deltaRMinTopSP and deltaRMinBottomSP to set separate values for outer and + /// inner space-points) + float deltaRMin = 5 * Acts::UnitConstants::mm; + /// Maximum radial distance between two doublet components (prefer + /// deltaRMaxTopSP and deltaRMacBottomSP to set separate values for outer and + /// inner space-points) + float deltaRMax = 270 * Acts::UnitConstants::mm; + /// Minimum radial distance between middle-outer doublet components + float deltaRMinTopSP = std::numeric_limits::quiet_NaN(); + /// Maximum radial distance between middle-outer doublet components + float deltaRMaxTopSP = std::numeric_limits::quiet_NaN(); + /// Minimum radial distance between inner-middle doublet components + float deltaRMinBottomSP = std::numeric_limits::quiet_NaN(); + /// Maximum radial distance between inner-middle doublet components + float deltaRMaxBottomSP = std::numeric_limits::quiet_NaN(); - // non equidistant binning in z - std::vector zBinEdges; + /// Maximum value of z-distance between space-points in doublet + float deltaZMax = + std::numeric_limits::infinity() * Acts::UnitConstants::mm; - // seed confirmation - bool seedConfirmation = false; - // parameters for central seed confirmation - SeedConfirmationRangeConfig centralSeedConfirmationRange; - // parameters for forward seed confirmation - SeedConfirmationRangeConfig forwardSeedConfirmationRange; + /// Maximum allowed cotTheta between two space-points in doublet, used to + /// check if forward angle is within bounds + float cotThetaMax = 10.01788; // equivalent to eta = 3 (pseudorapidity) - // FIXME: this is not used yet - // float upperPtResolutionPerSeed = 20* Acts::GeV; + /// Limiting location of collision region in z-axis used to check if doublet + /// origin is within reasonable bounds + float collisionRegionMin = -150 * Acts::UnitConstants::mm; + float collisionRegionMax = +150 * Acts::UnitConstants::mm; - // the delta for inverse helix radius up to which compared seeds - // are considered to have a compatible radius. delta of inverse radius - // leads to this value being the cutoff. unit is 1/mm. default value - // of 0.00003 leads to all helices with radius>33m to be considered compatible + /// Enable cut on the compatibility between interaction point and doublet, + /// this is an useful approximation to speed up the seeding + bool interactionPointCut = false; - // impact parameter - float impactMax = 20. * Acts::UnitConstants::mm; + /// Seeding parameters used to define the cuts on space-point triplets - // how many sigmas of scattering angle should be considered? + /// Minimum transverse momentum (pT) used to check the r-z slope compatibility + /// of triplets with maximum multiple scattering effect (produced by the + /// minimum allowed pT particle) + a certain uncertainty term. Check the + /// documentation for more information + /// https://acts.readthedocs.io/en/latest/core/reconstruction/pattern_recognition/seeding.html + float minPt = 400. * Acts::UnitConstants::MeV; + /// Number of sigmas of scattering angle to be considered in the minimum pT + /// scattering term float sigmaScattering = 5; - // Upper pt limit for scattering calculation + /// Term that accounts for the thickness of scattering medium in radiation + /// lengths in the Lynch & Dahl correction to the Highland equation default is + /// 5% + /// TODO: necessary to make amount of material dependent on detector region? + float radLengthPerSeed = 0.05; + /// Maximum transverse momentum for scattering calculation float maxPtScattering = 10 * Acts::UnitConstants::GeV; - - // for how many seeds can one SpacePoint be the middle SpacePoint? - unsigned int maxSeedsPerSpM = 5; - - // tolerance parameter used to check the compatibility of SPs coordinates in - // xyz - float toleranceParam = 1.1 * Acts::UnitConstants::mm; - - // Parameter which can loosen the tolerance of the track seed to form a - // helix. This is useful for e.g. misaligned seeding. + /// Maximum value of impact parameter estimation of the seed candidates + float impactMax = 20. * Acts::UnitConstants::mm; + /// Parameter which can loosen the tolerance of the track seed to form a + /// helix. This is useful for e.g. misaligned seeding. float helixCutTolerance = 1.; - // Geometry Settings - // Detector ROI - // limiting location of collision region in z - float collisionRegionMin = -150 * Acts::UnitConstants::mm; - float collisionRegionMax = +150 * Acts::UnitConstants::mm; - float phiMin = -M_PI; - float phiMax = M_PI; - // limiting location of measurements - float zMin = -2800 * Acts::UnitConstants::mm; - float zMax = 2800 * Acts::UnitConstants::mm; - float rMax = 600 * Acts::UnitConstants::mm; - // WARNING: if rMin is smaller than impactMax, the bin size will be 2*pi, - // which will make seeding very slow! - float rMin = 33 * Acts::UnitConstants::mm; + /// Seeding parameters used for quality seed confirmation - // Order of z bins to loop over when searching for SPs - std::vector zBinsCustomLooping = {}; - // Number of Z bins to skip the search for middle SPs - std::size_t skipZMiddleBinSearch = 0; + /// Enable quality seed confirmation, this is different than default seeding + /// confirmation because it can also be defined for different (r, z) regions + /// of the detector (e.g. forward or central region) by SeedConfirmationRange. + /// Seeds are classified as "high-quality" seeds and normal quality seeds. + /// Normal quality seeds are only selected if no other "high-quality" seeds + /// has been found for that inner-middle doublet. + bool seedConfirmation = false; + /// Contains parameters for central seed confirmation + SeedConfirmationRangeConfig centralSeedConfirmationRange; + /// Contains parameters for forward seed confirmation + SeedConfirmationRangeConfig forwardSeedConfirmationRange; + /// Maximum number (minus one) of accepted seeds per middle space-point + unsigned int maxSeedsPerSpM = 5; - // average radiation lengths of material on the length of a seed. used for - // scattering. - // default is 5% - // TODO: necessary to make amount of material dependent on detector region? - float radLengthPerSeed = 0.05; + /// Other parameters - // alignment uncertainties, used for uncertainties in the - // non-measurement-plane of the modules - // which otherwise would be 0 - // will be added to spacepoint measurement uncertainties (and therefore also - // multiplied by sigmaError) - // FIXME: call align1 and align2 + /// Alignment uncertainties, used for uncertainties in the + /// non-measurement-plane of the modules + /// which otherwise would be 0 + /// will be added to spacepoint measurement uncertainties (and therefore also + /// multiplied by sigmaError) + /// FIXME: call align1 and align2 float zAlign = 0 * Acts::UnitConstants::mm; float rAlign = 0 * Acts::UnitConstants::mm; - // used for measurement (+alignment) uncertainties. - // find seeds within 5sigma error ellipse + /// used for measurement (+alignment) uncertainties. + /// find seeds within 5sigma error ellipse float sigmaError = 5; - // derived values, set on SeedFinder construction + /// derived values, set on SeedFinder construction float highland = 0; float maxScatteringAngle2 = 0; - // only for Cuda plugin + /// only for Cuda plugin int maxBlockSize = 1024; int nTrplPerSpBLimit = 100; int nAvgTrplPerSpBLimit = 2; - // Delegates for accessors to detailed information on double measurement that - // produced the space point. - // This is mainly referring to space points produced when combining - // measurement from strips on back-to-back modules. - // Enables setting of the following delegates. + /// Delegates for accessors to detailed information on double measurement that + /// produced the space point. + /// This is mainly referring to space points produced when combining + /// measurement from strips on back-to-back modules. + /// Enables setting of the following delegates. bool useDetailedDoubleMeasurementInfo = false; - // Returns half of the length of the top strip. + /// Returns half of the length of the top strip. Delegate getTopHalfStripLength; - // Returns half of the length of the bottom strip. + /// Returns half of the length of the bottom strip. Delegate getBottomHalfStripLength; - // Returns direction of the top strip. + /// Returns direction of the top strip. Delegate getTopStripDirection; - // Returns direction of the bottom strip. + /// Returns direction of the bottom strip. Delegate getBottomStripDirection; - // Returns distance between the centers of the two strips. + /// Returns distance between the centers of the two strips. Delegate getStripCenterDistance; - // Returns position of the center of the top strip. + /// Returns position of the center of the top strip. Delegate getTopStripCenterPosition; + /// Tolerance parameter used to check the compatibility of space-point + /// coordinates in xyz. This is only used in a detector specific check for + /// strip modules + float toleranceParam = 1.1 * Acts::UnitConstants::mm; // Delegate to apply experiment specific cuts Delegate experimentCuts{ diff --git a/Core/include/Acts/Seeding/SeedFinderFTF.hpp b/Core/include/Acts/Seeding/SeedFinderFTF.hpp index b03d7306082..5136e86614e 100644 --- a/Core/include/Acts/Seeding/SeedFinderFTF.hpp +++ b/Core/include/Acts/Seeding/SeedFinderFTF.hpp @@ -63,22 +63,16 @@ class SeedFinderFTF { void loadSpacePoints( const std::vector> &FTF_SP_vect); + // inner + template void createSeeds( - const Acts::RoiDescriptor &roi, - const Acts::TrigFTF_GNN_Geometry &gnngeo); - - // create seeeds function - template - void createSeeds_old(const Acts::SeedFinderOptions &options, - const input_container_t &spacePoints, - output_container_t &out_cont, - callable_t &&extract_coordinates) const; - - template - std::vector createSeeds_old(const Acts::SeedFinderOptions &options, - const input_container_t &spacePoints, - callable_t &&extract_coordinates) const; + const Acts::RoiDescriptor & /*roi*/, + const Acts::TrigFTF_GNN_Geometry & /*gnngeo*/, + output_container_t & /*out_cont*/); + // outer + std::vector createSeeds( + const Acts::RoiDescriptor & /*roi*/, + const Acts::TrigFTF_GNN_Geometry & /*gnngeo*/); private: enum Dim { DimPhi = 0, DimR = 1, DimZ = 2 }; diff --git a/Core/include/Acts/Seeding/SeedFinderFTF.ipp b/Core/include/Acts/Seeding/SeedFinderFTF.ipp index 38696c7ade1..002f3f8cc09 100644 --- a/Core/include/Acts/Seeding/SeedFinderFTF.ipp +++ b/Core/include/Acts/Seeding/SeedFinderFTF.ipp @@ -57,7 +57,6 @@ void SeedFinderFTF::loadSpacePoints( m_storage->addSpacePoint(FTF_sp, (m_config.m_useClusterWidth > 0)); } - m_config.m_nMaxPhiSlice = 1; m_config.m_phiSliceWidth = 2 * M_PI / m_config.m_nMaxPhiSlice; m_storage->sortByPhi(); @@ -70,25 +69,13 @@ void SeedFinderFTF::runGNN_TrackFinder( std::vector>& vTracks, const Acts::RoiDescriptor& roi, const Acts::TrigFTF_GNN_Geometry& gnngeo) { - // long term move these to ftf finder config, then m_config. to access them - const int MaxEdges = 2000000; - - const float cut_dphi_max = 0.012; - const float cut_dcurv_max = 0.001; - const float cut_tau_ratio_max = 0.007; - const float min_z0 = -2800; // roiDescriptor->zedMinus(); //blank for now, - // get from config eventually - const float max_z0 = 2800; // roiDescriptor->zedPlus(); - - const float maxOuterRadius = 550.0; - const float cut_zMinU = - min_z0 + - maxOuterRadius * roi.dzdrMinus(); // dzdr can only find =0 in athena - const float cut_zMaxU = max_z0 + maxOuterRadius * roi.dzdrPlus(); - - float m_minR_squ = 1; // set earlier - float m_maxCurv = 1; - + const float min_z0 = roi.zedMinus(); + const float max_z0 = roi.zedPlus(); + const float cut_zMinU = min_z0 + m_config.maxOuterRadius * roi.dzdrMinus(); + const float cut_zMaxU = max_z0 + m_config.maxOuterRadius * roi.dzdrPlus(); + float m_minR_squ = m_config.m_tripletPtMin * m_config.m_tripletPtMin / + std::pow(m_config.ptCoeff, 2); // from athena + float m_maxCurv = m_config.ptCoeff / m_config.m_tripletPtMin; const float maxKappa_high_eta = 0.8 / m_minR_squ; const float maxKappa_low_eta = 0.6 / m_minR_squ; @@ -100,7 +87,7 @@ void SeedFinderFTF::runGNN_TrackFinder( std::vector> edgeStorage; - edgeStorage.reserve(MaxEdges); + edgeStorage.reserve(m_config.MaxEdges); int nEdges = 0; @@ -256,7 +243,7 @@ void SeedFinderFTF::runGNN_TrackFinder( continue; } - float zouter = z0 + maxOuterRadius * tau; + float zouter = z0 + m_config.maxOuterRadius * tau; if (zouter < cut_zMinU || zouter > cut_zMaxU) { continue; @@ -300,8 +287,9 @@ void SeedFinderFTF::runGNN_TrackFinder( float tau2 = edgeStorage.at(n2_in_idx).m_p[0]; float tau_ratio = tau2 * uat_1 - 1.0f; - if (std::fabs(tau_ratio) > cut_tau_ratio_max) { // bad - // match + if (std::fabs(tau_ratio) > + m_config.cut_tau_ratio_max) { // bad + // match continue; } isGood = true; // good match found @@ -316,7 +304,7 @@ void SeedFinderFTF::runGNN_TrackFinder( float dPhi2 = std::asin(curv * r2); float dPhi1 = std::asin(curv * r1); - if (nEdges < MaxEdges) { + if (nEdges < m_config.MaxEdges) { edgeStorage.emplace_back(n1, n2, exp_eta, curv, phi1 + dPhi1, phi2 + dPhi2); @@ -395,11 +383,11 @@ void SeedFinderFTF::runGNN_TrackFinder( float tau2 = pNS->m_p[0]; float tau_ratio = tau2 * uat_1 - 1.0f; - if (tau_ratio < -cut_tau_ratio_max) { + if (tau_ratio < -m_config.cut_tau_ratio_max) { last_out = out_idx; continue; } - if (tau_ratio > cut_tau_ratio_max) { + if (tau_ratio > m_config.cut_tau_ratio_max) { break; } @@ -411,14 +399,14 @@ void SeedFinderFTF::runGNN_TrackFinder( dPhi -= 2 * M_PI; } - if (dPhi < -cut_dphi_max || dPhi > cut_dphi_max) { + if (dPhi < -m_config.cut_dphi_max || dPhi > m_config.cut_dphi_max) { continue; } float curv2 = pNS->m_p[1]; float dcurv = curv2 - curv1; - if (dcurv < -cut_dcurv_max || dcurv > cut_dcurv_max) { + if (dcurv < -m_config.cut_dcurv_max || dcurv > m_config.cut_dcurv_max) { continue; } @@ -497,7 +485,7 @@ void SeedFinderFTF::runGNN_TrackFinder( std::vector*> vSeeds; - vSeeds.reserve(MaxEdges / 2); + vSeeds.reserve(m_config.MaxEdges / 2); for (int edgeIndex = 0; edgeIndex < nEdges; edgeIndex++) { Acts::TrigFTF_GNN_Edge* pS = @@ -662,9 +650,11 @@ void SeedFinderFTF::runGNN_TrackFinder( } template +template void SeedFinderFTF::createSeeds( const Acts::RoiDescriptor& roi, - const Acts::TrigFTF_GNN_Geometry& gnngeo) { + const Acts::TrigFTF_GNN_Geometry& gnngeo, + output_container_t& out_cont) { std::vector> vTracks; // make empty vector @@ -705,27 +695,29 @@ void SeedFinderFTF::createSeeds( } } vTracks.clear(); -} -// // still to be developed -template -template -void SeedFinderFTF::createSeeds_old( - const Acts::SeedFinderOptions& /*options*/, - const input_container_t& /*spacePoints*/, output_container_t& /*out_cont*/, - callable_t&& /*extract_coordinates*/) const {} + for (auto& triplet : m_triplets) { + const external_spacepoint_t* S1 = + triplet.s1().SP; // triplet-> FTF_SP-> simspacepoint + const external_spacepoint_t* S2 = triplet.s2().SP; + const external_spacepoint_t* S3 = triplet.s3().SP; + + // input to seed + float Vertex = 0; + float Quality = triplet.Q(); + // make a new seed, add to vector of seeds + out_cont.emplace_back(*S1, *S2, *S3, Vertex, Quality); + } +} +// outer called in alg template -template std::vector> -SeedFinderFTF::createSeeds_old( - const Acts::SeedFinderOptions& options, - const input_container_t& spacePoints, - callable_t&& extract_coordinates) const { +SeedFinderFTF::createSeeds( + const Acts::RoiDescriptor& roi, + const Acts::TrigFTF_GNN_Geometry& gnngeo) { std::vector r; - createSeeds_old(options, spacePoints, r, - std::forward(extract_coordinates)); + createSeeds(roi, gnngeo, r); return r; } diff --git a/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp b/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp index 420935ecb82..7103a9a1a3d 100644 --- a/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp +++ b/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp @@ -46,22 +46,32 @@ struct SeedFinderFTFConfig { // helix. This is useful for e.g. misaligned seeding. float helixCutTolerance = 1.; - float m_phiSliceWidth{}; - float m_nMaxPhiSlice{}; - bool m_useClusterWidth = false; - std::string fastrack_input_file; + float m_phiSliceWidth{}; // initialised in loadSpacePoints function + float m_nMaxPhiSlice = 53; // used to calculate phi slices + bool m_useClusterWidth = + false; // bool for use of cluster width in loadSpacePoints function + std::string fastrack_input_file; // input file for fastrack object std::vector m_layerGeometry; - // for run function - // m_settings: + // for runGNN_TrackFinder bool m_LRTmode = true; // eventually want to set from full chain - bool m_useEtaBinning = true; - bool m_doubletFilterRZ = true; - float m_minDeltaRadius = 5.0; // eventually set in config or to equivalent - // acts 2.0 but increasing to test loops - // float m_maxDeltaRadius = 270.0 ; - float m_tripletD0Max = 4.0; // m_settings - unsigned int m_maxTripletBufferLength = 3; + bool m_useEtaBinning = + true; // bool to use eta binning from geometry structure + bool m_doubletFilterRZ = true; // bool applies new Z cuts on doublets + float m_minDeltaRadius = 2.0; // min dr for doublet + float m_tripletD0Max = 4.0; // D0 cut for triplets + unsigned int m_maxTripletBufferLength = + 3; // maximum number of space points per triplet + int MaxEdges = 2000000; // max number of GNN edges/doublets + float cut_dphi_max = 0.012; // phi cut for triplets + float cut_dcurv_max = 0.001; // curv cut for triplets + float cut_tau_ratio_max = 0.007; // tau cut for doublets and triplets + float maxOuterRadius = 550.0; // used to calculate Z cut on doublets + float m_PtMin = 1000.0; + float m_tripletPtMinFrac = 0.3; + float m_tripletPtMin = m_PtMin * m_tripletPtMinFrac; // Limit on triplet pt + double ptCoeff = + 0.29997 * 1.9972 / 2.0; // ~0.3*B/2 - assumes nominal field of 2*T // ROI: bool containsPhi() { diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp b/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp index f735c467c4c..706cb1e90cc 100644 --- a/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp @@ -20,93 +20,124 @@ namespace Acts { template class SeedFilter; +/// @brief Structure that holds configuration parameters for the orthogonal seed finder algorithm template struct SeedFinderOrthogonalConfig { std::shared_ptr> seedFilter; - // Seed Cuts - // lower cutoff for seeds - float minPt = 400. * Acts::UnitConstants::MeV; - // cot of maximum theta angle - // equivalent to 2.7 eta (pseudorapidity) - float cotThetaMax = 7.40627; - // minimum distance in r between middle and top SP in one seed - float deltaRMinTopSP = 5 * Acts::UnitConstants::mm; - // maximum distance in r between middle and top SP in one seed - float deltaRMaxTopSP = 270 * Acts::UnitConstants::mm; - // minimum distance in r between middle and bottom SP in one seed - float deltaRMinBottomSP = 5 * Acts::UnitConstants::mm; - // maximum distance in r between middle and bottom SP in one seed - float deltaRMaxBottomSP = 270 * Acts::UnitConstants::mm; - - // impact parameter - float impactMax = 20. * Acts::UnitConstants::mm; - - // how many sigmas of scattering angle should be considered? - float sigmaScattering = 5; - // Upper pt limit for scattering calculation - float maxPtScattering = 10 * Acts::UnitConstants::GeV; + /// Seeding parameters for geometry settings and detector ROI - // for how many seeds can one SpacePoint be the middle SpacePoint? - unsigned int maxSeedsPerSpM = 5; - - // Geometry Settings - // Detector ROI - // limiting location of collision region in z - float collisionRegionMin = -150 * Acts::UnitConstants::mm; - float collisionRegionMax = +150 * Acts::UnitConstants::mm; + // Limiting location of all measurements float phiMin = -M_PI; float phiMax = M_PI; - // limiting location of measurements + /// limiting location of measurements float zMin = -2800 * Acts::UnitConstants::mm; float zMax = 2800 * Acts::UnitConstants::mm; float rMax = 600 * Acts::UnitConstants::mm; - // WARNING: if rMin is smaller than impactMax, the bin size will be 2*pi, - // which will make seeding very slow! + /// @warning If rMin is smaller than impactMax, the bin size will be 2*pi, + /// which will make seeding very slow! float rMin = 33 * Acts::UnitConstants::mm; - // z of last layers to avoid iterations - std::pair zOutermostLayers{-2700 * Acts::UnitConstants::mm, - 2700 * Acts::UnitConstants::mm}; + /// Seeding parameters used to define the region of interest for middle + /// space-point - // radial range for middle SP - // variable range based on SP radius + /// Radial range for middle space-point + /// The range can be defined manually with (rMinMiddle, rMaxMiddle). If + /// useVariableMiddleSPRange is set to false and the vector rRangeMiddleSP is + /// empty, we use (rMinMiddle, rMaxMiddle) to cut the middle space-points + float rMinMiddle = 60.f * Acts::UnitConstants::mm; + float rMaxMiddle = 120.f * Acts::UnitConstants::mm; + /// If useVariableMiddleSPRange is set to false, the vector rRangeMiddleSP can + /// be used to define a fixed r range for each z bin: {{rMin, rMax}, ...} bool useVariableMiddleSPRange = true; + /// Range defined in vector for each z bin + std::vector> rRangeMiddleSP; + /// If useVariableMiddleSPRange is true, the radial range will be calculated + /// based on the maximum and minimum r values of the space-points in the event + /// and a deltaR (deltaRMiddleMinSPRange, deltaRMiddleMaxSPRange) float deltaRMiddleMinSPRange = 10. * Acts::UnitConstants::mm; float deltaRMiddleMaxSPRange = 10. * Acts::UnitConstants::mm; - // range defined in vector for each z region - std::vector> rRangeMiddleSP; - // range defined by rMinMiddle and rMaxMiddle - float rMinMiddle = 60.f * Acts::UnitConstants::mm; - float rMaxMiddle = 120.f * Acts::UnitConstants::mm; + /// Vector containing minimum and maximum z boundaries for cutting middle + /// space-points + std::pair zOutermostLayers{-2700 * Acts::UnitConstants::mm, + 2700 * Acts::UnitConstants::mm}; + + /// Seeding parameters used to define the cuts on space-point doublets + + /// Minimum radial distance between middle-outer doublet components + float deltaRMinTopSP = std::numeric_limits::quiet_NaN(); + /// Maximum radial distance between middle-outer doublet components + float deltaRMaxTopSP = std::numeric_limits::quiet_NaN(); + /// Minimum radial distance between inner-middle doublet components + float deltaRMinBottomSP = std::numeric_limits::quiet_NaN(); + /// Maximum radial distance between inner-middle doublet components + float deltaRMaxBottomSP = std::numeric_limits::quiet_NaN(); + + /// Shrink the phi range of middle space-point (analogous to phi bin size in + /// grid from default seeding + number of phi bins used in search) float deltaPhiMax = 0.085; - // cut to the maximum value of delta z between SPs + /// Maximum value of z-distance between space-points in doublet float deltaZMax = std::numeric_limits::infinity() * Acts::UnitConstants::mm; - // enable cut on the compatibility between interaction point and SPs + /// Maximum allowed cotTheta between two space-points in doublet, used to + /// check if forward angle is within bounds + float cotThetaMax = 7.40627; // equivalent to 2.7 eta (pseudorapidity) + + /// Limiting location of collision region in z-axis used to check if doublet + /// origin is within reasonable bounds + float collisionRegionMin = -150 * Acts::UnitConstants::mm; + float collisionRegionMax = +150 * Acts::UnitConstants::mm; + + /// Enable cut on the compatibility between interaction point and doublet, + /// this is an useful approximation to speed up the seeding bool interactionPointCut = false; - // seed confirmation + /// Seeding parameters used to define the cuts on space-point triplets + + /// Minimum transverse momentum (pT) used to check the r-z slope compatibility + /// of triplets with maximum multiple scattering effect (produced by the + /// minimum allowed pT particle) + a certain uncertainty term. Check the + /// documentation for more information + /// https://acts.readthedocs.io/en/latest/core/reconstruction/pattern_recognition/seeding.html + float minPt = 400. * Acts::UnitConstants::MeV; + /// Number of sigmas of scattering angle to be considered in the minimum pT + /// scattering term + float sigmaScattering = 5; + /// Term that accounts for the thickness of scattering medium in radiation + /// lengths in the Lynch & Dahl correction to the Highland equation default is + /// 5% + /// TODO: necessary to make amount of material dependent on detector region? + float radLengthPerSeed = 0.05; + /// Maximum transverse momentum for scattering calculation + float maxPtScattering = 10 * Acts::UnitConstants::GeV; + /// Maximum value of impact parameter estimation of the seed candidates + float impactMax = 20. * Acts::UnitConstants::mm; + /// Parameter which can loosen the tolerance of the track seed to form a + /// helix. This is useful for e.g. misaligned seeding. + float helixCutTolerance = 1.; + + /// Seeding parameters used for quality seed confirmation + + /// Enable quality seed confirmation, this is different than default seeding + /// confirmation because it can also be defined for different (r, z) regions + /// of the detector (e.g. forward or central region) by SeedConfirmationRange. + /// Seeds are classified as "high-quality" seeds and normal quality seeds. + /// Normal quality seeds are only selected if no other "high-quality" seeds + /// has been found for that inner-middle doublet. bool seedConfirmation = false; - // parameters for central seed confirmation + /// Contains parameters for central seed confirmation SeedConfirmationRangeConfig centralSeedConfirmationRange; - // parameters for forward seed confirmation + /// Contains parameters for forward seed confirmation SeedConfirmationRangeConfig forwardSeedConfirmationRange; + /// Maximum number (minus one) of accepted seeds per middle space-point + unsigned int maxSeedsPerSpM = 5; - // average radiation lengths of material on the length of a seed. used for - // scattering. - // default is 5% - // TODO: necessary to make amount of material dependent on detector region? - float radLengthPerSeed = 0.05; - - // Parameter which can loosen the tolerance of the track seed to form a - // helix. This is useful for e.g. misaligned seeding. - float helixCutTolerance = 1.; + /// Other parameters - // derived values, set on SeedFinder construction + /// derived values, set on SeedFinder construction float highland = 0; float maxScatteringAngle2 = 0; @@ -123,8 +154,8 @@ struct SeedFinderOrthogonalConfig { "calculateDerivedQuantities"); } SeedFinderOrthogonalConfig config = *this; - // calculation of scattering using the highland formula - // convert pT to p once theta angle is known + /// calculation of scattering using the highland formula + /// convert pT to p once theta angle is known config.highland = 13.6 * std::sqrt(radLengthPerSeed) * (1 + 0.038 * std::log(radLengthPerSeed)); config.maxScatteringAngle2 = std::pow(config.highland / config.minPt, 2); diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index cab69814517..66025969fa8 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -29,6 +29,7 @@ #include "Acts/Propagator/ConstrainedStep.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Propagator/StandardAborters.hpp" +#include "Acts/Propagator/detail/LoopProtection.hpp" #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp" #include "Acts/Surfaces/BoundaryCheck.hpp" #include "Acts/TrackFinding/CombinatorialKalmanFilterError.hpp" @@ -392,6 +393,13 @@ class CombinatorialKalmanFilter { ACTS_VERBOSE("CombinatorialKalmanFilter step"); + // Initialize path limit reached aborter + if (result.pathLimitReached.internalLimit == + std::numeric_limits::max()) { + detail::setupLoopProtection(state, stepper, result.pathLimitReached, + true, logger()); + } + if (!result.filtered && filterTargetReached(state, stepper, navigator, logger())) { navigator.navigationBreak(state.navigation, true); @@ -488,12 +496,19 @@ class CombinatorialKalmanFilter { } } - if (endOfWorldReached(state, stepper, navigator, logger()) || - result.pathLimitReached(state, stepper, navigator, logger())) { - navigator.targetReached(state.navigation, false); - if (result.activeTips.empty()) { - // we are already done - } else if (result.activeTips.size() == 1) { + bool isEndOfWorldReached = + endOfWorldReached(state, stepper, navigator, logger()); + bool isPathLimitReached = + result.pathLimitReached(state, stepper, navigator, logger()); + if (isEndOfWorldReached || isPathLimitReached) { + if (isEndOfWorldReached) { + ACTS_VERBOSE("End of world reached"); + } + if (isPathLimitReached) { + ACTS_VERBOSE("Path limit reached"); + } + + if (result.activeTips.size() <= 1) { // this was the last track - we are done ACTS_VERBOSE("Kalman filtering finds " << result.lastTrackIndices.size() << " tracks"); @@ -553,7 +568,7 @@ class CombinatorialKalmanFilter { // track parameters for found track indexed with iSmoothed bool isTargetReached = smoothingTargetReached(state, stepper, navigator, logger()); - bool isPathLimitReached = + isPathLimitReached = result.pathLimitReached(state, stepper, navigator, logger()); if (result.smoothed && (isTargetReached || isPathLimitReached)) { ACTS_VERBOSE( @@ -1357,6 +1372,20 @@ class CombinatorialKalmanFilter { } }; + /// Void path limit reached aborter to replace the default since the path + /// limit is handled in the CKF actor internally. + struct StubPathLimitReached { + double internalLimit{}; + + template + bool operator()(propagator_state_t& /*unused*/, const stepper_t& /*unused*/, + const navigator_t& /*unused*/, + const Logger& /*unused*/) const { + return false; + } + }; + public: /// Combinatorial Kalman Filter implementation, calls the Kalman filter /// and smoother @@ -1443,7 +1472,8 @@ class CombinatorialKalmanFilter { r.stateBuffer = stateBuffer; r.stateBuffer->clear(); - auto result = m_propagator.template propagate( + auto result = m_propagator.template propagate< + start_parameters_t, decltype(propOptions), PathLimitReached>( initialParameters, propOptions, false, std::move(inputResult)); if (!result.ok()) { diff --git a/Core/include/Acts/TrackFinding/TrackSelector.hpp b/Core/include/Acts/TrackFinding/TrackSelector.hpp index 34afb8db322..d73176b8796 100644 --- a/Core/include/Acts/TrackFinding/TrackSelector.hpp +++ b/Core/include/Acts/TrackFinding/TrackSelector.hpp @@ -46,6 +46,10 @@ class TrackSelector { double ptMax = inf; std::size_t minMeasurements = 0; + std::size_t maxHoles = std::numeric_limits::max(); + std::size_t maxOutliers = std::numeric_limits::max(); + std::size_t maxSharedHits = std::numeric_limits::max(); + double maxChi2 = inf; // Helper factory functions to produce a populated config object more // conveniently @@ -269,6 +273,10 @@ inline std::ostream& operator<<(std::ostream& os, print("eta", cuts.etaMin, cuts.etaMax); print("absEta", cuts.absEtaMin, cuts.absEtaMax); print("pt", cuts.ptMin, cuts.ptMax); + print("nHoles", 0, cuts.maxHoles); + print("nOutliers", 0, cuts.maxOutliers); + print("nSharedHits", 0, cuts.maxSharedHits); + print("chi2", 0.0, cuts.maxChi2); os << " - " << cuts.minMeasurements << " <= nMeasurements\n"; return os; @@ -342,6 +350,7 @@ void TrackSelector::selectTracks(const input_tracks_t& inputTracks, template bool TrackSelector::isValidTrack(const track_proxy_t& track) const { auto checkMin = [](auto x, auto min) { return min <= x; }; + auto checkMax = [](auto x, auto max) { return x <= max; }; auto within = [](double x, double min, double max) { return (min <= x) && (x < max); }; @@ -382,7 +391,11 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const { within(track.loc0(), cuts.loc0Min, cuts.loc0Max) && within(track.loc1(), cuts.loc1Min, cuts.loc1Max) && within(track.time(), cuts.timeMin, cuts.timeMax) && - checkMin(track.nMeasurements(), cuts.minMeasurements); + checkMin(track.nMeasurements(), cuts.minMeasurements) && + checkMax(track.nHoles(), cuts.maxHoles) && + checkMax(track.nOutliers(), cuts.maxOutliers) && + checkMax(track.nSharedHits(), cuts.maxSharedHits) && + checkMax(track.chi2(), cuts.maxChi2); } inline TrackSelector::TrackSelector( diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index 46f3d668bc5..6bb62a1d9b4 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -814,6 +814,10 @@ class Gx2Fitter { // TODO write test for calculateTrackQuantities calculateTrackQuantities(track); + // Set the chi2sum for the track summary manually, since we don't calculate + // it for each state + track.chi2() = chi2sum; + // Return the converted Track return track; } diff --git a/Core/include/Acts/Utilities/GridFwd.hpp b/Core/include/Acts/Utilities/GridFwd.hpp deleted file mode 100644 index 7cb2b5daa77..00000000000 --- a/Core/include/Acts/Utilities/GridFwd.hpp +++ /dev/null @@ -1,14 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019 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/. - -#pragma once - -namespace Acts { -template -class Grid; -} // namespace Acts diff --git a/Core/include/Acts/Utilities/Intersection.hpp b/Core/include/Acts/Utilities/Intersection.hpp index a4ef4990d52..2f82b531532 100644 --- a/Core/include/Acts/Utilities/Intersection.hpp +++ b/Core/include/Acts/Utilities/Intersection.hpp @@ -117,24 +117,9 @@ using MultiIntersection3D = boost::container::static_vector; -/// @brief class extensions to return also the object and a representation -template +template class ObjectIntersection { public: - /// Object intersection - symmetric setup - /// - /// @param intersection is the intersection - /// @param object is the object to be instersected - /// @param index is the intersection index - template ::value, int> = 0> - constexpr ObjectIntersection(const Intersection3D& intersection, - const object_t* object, std::uint8_t index = 0) - : m_intersection(intersection), - m_object(object), - m_representation(object), - m_index(index) {} - /// Object intersection /// /// @param intersection is the intersection @@ -142,13 +127,8 @@ class ObjectIntersection { /// @param representation is the object representation /// @param index is the intersection index constexpr ObjectIntersection(const Intersection3D& intersection, - const object_t* object, - const representation_t* representation, - std::uint8_t index = 0) - : m_intersection(intersection), - m_object(object), - m_representation(representation), - m_index(index) {} + const object_t* object, std::uint8_t index = 0) + : m_intersection(intersection), m_object(object), m_index(index) {} /// Returns whether the intersection was successful or not constexpr explicit operator bool() const { @@ -173,10 +153,6 @@ class ObjectIntersection { constexpr const object_t* object() const { return m_object; } - constexpr const representation_t* representation() const { - return m_representation; - } - constexpr std::uint8_t index() const { return m_index; } constexpr static ObjectIntersection invalid() { return ObjectIntersection(); } @@ -198,33 +174,18 @@ class ObjectIntersection { Intersection3D m_intersection = Intersection3D::invalid(); /// The object that was (tried to be) intersected const object_t* m_object = nullptr; - /// The representation of this object - const representation_t* m_representation = nullptr; /// The intersection index std::uint8_t m_index = 0; constexpr ObjectIntersection() = default; }; -/// @brief class extensions to return also the object and a representation -template +template class ObjectMultiIntersection { public: - using SplitIntersections = boost::container::static_vector< - ObjectIntersection, - s_maximumNumberOfIntersections>; - - /// Object intersection - symmetric setup - /// - /// @param intersections are the intersections - /// @param object is the object to be instersected - template ::value, int> = 0> - constexpr ObjectMultiIntersection(const MultiIntersection3D& intersections, - const object_t* object) - : m_intersections(intersections), - m_object(object), - m_representation(object) {} + using SplitIntersections = + boost::container::static_vector, + s_maximumNumberOfIntersections>; /// Object intersection /// @@ -232,25 +193,17 @@ class ObjectMultiIntersection { /// @param object is the object to be instersected /// @param representation is the object representation constexpr ObjectMultiIntersection(const MultiIntersection3D& intersections, - const object_t* object, - const representation_t* representation) - : m_intersections(intersections), - m_object(object), - m_representation(representation) {} - - constexpr ObjectIntersection operator[]( - std::uint8_t index) const { - return {m_intersections[index], m_object, m_representation, index}; + const object_t* object) + : m_intersections(intersections), m_object(object) {} + + constexpr ObjectIntersection operator[](std::uint8_t index) const { + return {m_intersections[index], m_object, index}; } constexpr std::size_t size() const { return m_intersections.size(); } constexpr const object_t* object() const { return m_object; } - constexpr const representation_t* representation() const { - return m_representation; - } - constexpr SplitIntersections split() const { SplitIntersections result; for (std::size_t i = 0; i < size(); ++i) { @@ -259,11 +212,11 @@ class ObjectMultiIntersection { return result; } - constexpr ObjectIntersection closest() const { + constexpr ObjectIntersection closest() const { auto splitIntersections = split(); - return *std::min_element( - splitIntersections.begin(), splitIntersections.end(), - ObjectIntersection::closestOrder); + return *std::min_element(splitIntersections.begin(), + splitIntersections.end(), + ObjectIntersection::closestOrder); } private: @@ -271,8 +224,6 @@ class ObjectMultiIntersection { MultiIntersection3D m_intersections; /// The object that was (tried to be) intersected const object_t* m_object = nullptr; - /// The representation of this object - const representation_t* m_representation = nullptr; }; namespace detail { diff --git a/Core/include/Acts/Utilities/detail/MPL/get_position.hpp b/Core/include/Acts/Utilities/detail/MPL/get_position.hpp deleted file mode 100644 index a3f9e8d2497..00000000000 --- a/Core/include/Acts/Utilities/detail/MPL/get_position.hpp +++ /dev/null @@ -1,42 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2016-2018 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/. - -#pragma once -namespace Acts { -/// @cond detail -namespace detail { -/** - * @brief get position of integral constant in template parameter pack - * - * @tparam T integral type of the values to be investigated - * @tparam target target value whose position in the template parameter pack - * should be determined - * @tparam values template parameter pack containing the list of values - * - * @return `get_position::value` yields the position of - * `target` inside `values`. - * If `target` is not in the list of `values`, a compile-time error is - * generated. - */ -template -struct get_position; - -/// @cond -template -struct get_position { - enum { value = 0 }; -}; - -template -struct get_position { - enum { value = get_position::value + 1 }; -}; -/// @endcond -} // namespace detail -/// @endcond -} // namespace Acts diff --git a/Core/include/Acts/Utilities/detail/MPL/is_contained.hpp b/Core/include/Acts/Utilities/detail/MPL/is_contained.hpp deleted file mode 100644 index 4b2ad62ed87..00000000000 --- a/Core/include/Acts/Utilities/detail/MPL/is_contained.hpp +++ /dev/null @@ -1,50 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2016-2018 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/. - -#pragma once -namespace Acts { -/// @cond detail -namespace detail { -/** - * @brief check whether given integral constant is contained in a template - * parameter pack - * - * @tparam T integral type of the values to be checked - * @tparam target target value to be looked for - * @tparam values template parameter pack containing the list of values - * - * @return `is_contained::value` is @c true if `target` is - * among `values`, otherwise @c false - */ -template -struct is_contained; - -/// @cond -template -struct is_contained { - enum { value = true }; -}; - -template -struct is_contained { - enum { value = true }; -}; - -template -struct is_contained { - enum { value = is_contained::value }; -}; - -template -struct is_contained { - enum { value = false }; -}; -/// @endcond -} // namespace detail -/// @endcond -} // namespace Acts diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index 2c2a3969180..89bb09ba2a2 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -307,9 +307,16 @@ Acts::Result Acts:: trkAtVtx.linearizedState = *result; trkAtVtx.isLinearized = true; } - // Update the vertex with the new track - KalmanVertexUpdater::updateVertexWithTrack(*vtx, - trkAtVtx); + // Update the vertex with the new track. The second template argument + // corresponds to the number of fitted vertex dimensions (i.e., 3 if we + // only fit spatial coordinates and 4 if we also fit time). + if (m_cfg.useTime) { + KalmanVertexUpdater::updateVertexWithTrack( + *vtx, trkAtVtx); + } else { + KalmanVertexUpdater::updateVertexWithTrack( + *vtx, trkAtVtx); + } } else { ACTS_VERBOSE("Track weight too low. Skip track."); } @@ -369,7 +376,15 @@ void Acts::AdaptiveMultiVertexFitter< for (const auto trk : state.vtxInfoMap[vtx].trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); if (trkAtVtx.trackWeight > m_cfg.minWeight) { - KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + // Update the new track under the assumption that it originates at the + // vertex. The second template argument corresponds to the number of + // fitted vertex dimensions (i.e., 3 if we only fit spatial coordinates + // and 4 if we also fit time). + if (m_cfg.useTime) { + KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + } else { + KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + } } } } diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp index 72b5f9b668d..e90de63943d 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp @@ -18,29 +18,44 @@ namespace Acts { namespace KalmanVertexTrackUpdater { /// KalmanVertexTrackUpdater -/// +/// Based on +/// Ref. (1): +/// R. Frühwirth et al. +/// Vertex reconstruction and track bundling at the lep collider using robust +/// algorithms +/// Computer Physics Comm.: 96 (1996) 189 +/// Chapter 2.1 + /// @brief Refits a single track with the knowledge of /// the vertex it has originated from /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param track Track to update /// @param vtx Vertex `track` belongs to -template +template void update(TrackAtVertex& track, const Vertex& vtx); namespace detail { -/// @brief reates a new covariance matrix for the -/// refitted track parameters +/// @brief Calculates a covariance matrix for the refitted track parameters +/// +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). /// -/// @param sMat Track ovariance in momentum space -/// @param newTrkCov New track covariance matrixs +/// @param wMat W_k matrix from Ref. (1) +/// @param crossCovVP Cross-covariance matrix between vertex position and track +/// momentum /// @param vtxCov Vertex covariance matrix -/// @param newTrkParams New track parameter -inline BoundMatrix createFullTrackCovariance(const SquareMatrix3& wMat, - const SquareMatrix3& crossCovVP, - const SquareMatrix3& vtxCov, - const BoundVector& newTrkParams); +/// @param newTrkParams Refitted track parameters +template +inline BoundMatrix calculateTrackCovariance( + const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, + const ActsSquareMatrix& vtxCov, + const BoundVector& newTrkParams); } // Namespace detail diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 8ff1cb323c6..5a183704a75 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -11,9 +11,20 @@ #include "Acts/Utilities/detail/periodic.hpp" #include "Acts/Vertexing/VertexingError.hpp" -template +template void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const Vertex& vtx) { + if constexpr (nDimVertex != 3 && nDimVertex != 4) { + throw std::invalid_argument( + "The vertex dimension must either be 3 (when fitting the spatial " + "coordinates) or 4 (when fitting the spatial coordinates + time)."); + } + + using VertexVector = ActsVector; + using VertexMatrix = ActsSquareMatrix; + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // Check if linearized state exists if (!track.isLinearized) { throw std::invalid_argument("TrackAtVertex object must be linearized."); @@ -21,30 +32,35 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, // Extract vertex position and covariance // \tilde{x_n} - const Vector3& vtxPos = vtx.position(); + const VertexVector vtxPos = vtx.fullPosition().template head(); // C_n - const SquareMatrix3& vtxCov = vtx.covariance(); + const VertexMatrix vtxCov = + vtx.fullCovariance().template block(0, 0); // Get the linearized track const LinearizedTrack& linTrack = track.linearizedState; - // Retrieve variables from the track linearization. The comments indicate the - // corresponding symbol used in the Ref. (1). + // corresponding symbol used in Ref. (1). // A_k - const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would use time - // G_k - const ActsSquareMatrix<5> trkParamWeight = - linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); // c_k - const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time. + // G_k + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); // Cache object filled with zeros - KalmanVertexUpdater::Cache cache; + KalmanVertexUpdater::Cache cache; // Calculate the update of the vertex position when the track is removed. This // might be unintuitive, but it is needed to compute a symmetric chi2. @@ -67,30 +83,34 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, newTrkParams(BoundIndices::eBoundQOverP) = newTrkMomentum(2); // qOverP // E_k^n - const SquareMatrix3 crossCovVP = + const ActsMatrix crossCovVP = -vtxCov * posJac.transpose() * trkParamWeight * momJac * cache.wMat; // Difference in positions. cache.newVertexPos corresponds to \tilde{x_k^{n*}} in Ref. (1). - Vector3 posDiff = vtxPos - cache.newVertexPos; + VertexVector posDiff = + vtxPos - cache.newVertexPos.template head(); // r_k^n - ActsVector<5> paramDiff = + ParameterVector paramDiff = trkParams - (constTerm + posJac * vtxPos + momJac * newTrkMomentum); // New chi2 to be set later - double chi2 = posDiff.dot(cache.newVertexWeight * posDiff) + - paramDiff.dot(trkParamWeight * paramDiff); + double chi2 = + posDiff.dot( + cache.newVertexWeight.template block(0, 0) * + posDiff) + + paramDiff.dot(trkParamWeight * paramDiff); - Acts::BoundMatrix fullPerTrackCov = detail::createFullTrackCovariance( + Acts::BoundMatrix newTrackCov = detail::calculateTrackCovariance( cache.wMat, crossCovVP, vtxCov, newTrkParams); // Create new refitted parameters std::shared_ptr perigeeSurface = - Surface::makeShared(vtxPos); + Surface::makeShared(vtxPos.template head<3>()); - BoundTrackParameters refittedPerigee = BoundTrackParameters( - perigeeSurface, newTrkParams, std::move(fullPerTrackCov), - track.fittedParams.particleHypothesis()); + BoundTrackParameters refittedPerigee = + BoundTrackParameters(perigeeSurface, newTrkParams, std::move(newTrackCov), + track.fittedParams.particleHypothesis()); // Set new properties track.fittedParams = refittedPerigee; @@ -100,42 +120,66 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, return; } +template inline Acts::BoundMatrix -Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( - const SquareMatrix3& wMat, const SquareMatrix3& crossCovVP, - const SquareMatrix3& vtxCov, const BoundVector& newTrkParams) { +Acts::KalmanVertexTrackUpdater::detail::calculateTrackCovariance( + const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, + const ActsSquareMatrix& vtxCov, + const BoundVector& newTrkParams) { // D_k^n ActsSquareMatrix<3> momCov = wMat + crossCovVP.transpose() * vtxCov.inverse() * crossCovVP; - // Full (x,y,z,phi, theta, q/p) covariance matrix - // To be made 7d again after switching to (x,y,z,phi, theta, q/p, t) - ActsSquareMatrix<6> fullTrkCov(ActsSquareMatrix<6>::Zero()); - - fullTrkCov.block<3, 3>(0, 0) = vtxCov; - fullTrkCov.block<3, 3>(0, 3) = crossCovVP; - fullTrkCov.block<3, 3>(3, 0) = crossCovVP.transpose(); - fullTrkCov.block<3, 3>(3, 3) = momCov; + // Full x, y, z, phi, theta, q/p, and, optionally, t covariance matrix. Note + // that we call this set of parameters "free" in the following even though + // that is not quite correct (free parameters correspond to + // x, y, z, t, px, py, pz) + constexpr unsigned int nFreeParams = nDimVertex + 3; + ActsSquareMatrix freeTrkCov( + ActsSquareMatrix::Zero()); + + freeTrkCov.template block<3, 3>(0, 0) = vtxCov.template block<3, 3>(0, 0); + freeTrkCov.template block<3, 3>(0, 3) = crossCovVP.template block<3, 3>(0, 0); + freeTrkCov.template block<3, 3>(3, 0) = + crossCovVP.template block<3, 3>(0, 0).transpose(); + freeTrkCov.template block<3, 3>(3, 3) = momCov; + + // Fill time (cross-)covariances + if constexpr (nFreeParams == 7) { + freeTrkCov.template block<3, 1>(0, 6) = vtxCov.template block<3, 1>(0, 3); + freeTrkCov.template block<3, 1>(3, 6) = + crossCovVP.template block<1, 3>(3, 0).transpose(); + freeTrkCov.template block<1, 3>(6, 0) = vtxCov.template block<1, 3>(3, 0); + freeTrkCov.template block<1, 3>(6, 3) = + crossCovVP.template block<1, 3>(3, 0); + freeTrkCov(6, 6) = vtxCov(3, 3); + } - // Combined track jacobian - ActsMatrix<5, 6> trkJac(ActsMatrix<5, 6>::Zero()); + // Jacobian relating "free" and bound track parameters + constexpr unsigned int nBoundParams = nDimVertex + 2; + ActsMatrix freeToBoundJac( + ActsMatrix::Zero()); + // TODO: Jacobian is not quite correct // First row - trkJac(0, 0) = -std::sin(newTrkParams[2]); - trkJac(0, 1) = std::cos(newTrkParams[2]); + freeToBoundJac(0, 0) = -std::sin(newTrkParams[2]); + freeToBoundJac(0, 1) = std::cos(newTrkParams[2]); double tanTheta = std::tan(newTrkParams[3]); // Second row - trkJac(1, 0) = -trkJac(0, 1) / tanTheta; - trkJac(1, 1) = trkJac(0, 0) / tanTheta; + freeToBoundJac(1, 0) = -freeToBoundJac(0, 1) / tanTheta; + freeToBoundJac(1, 1) = freeToBoundJac(0, 0) / tanTheta; - trkJac.block<4, 4>(1, 2) = ActsMatrix<4, 4>::Identity(); + // Dimension of the part of the Jacobian that is an identity matrix + constexpr unsigned int nDimIdentity = nFreeParams - 2; + freeToBoundJac.template block(1, 2) = + ActsMatrix::Identity(); // Full perigee track covariance - BoundMatrix fullPerTrackCov(BoundMatrix::Identity()); - fullPerTrackCov.block<5, 5>(0, 0) = - (trkJac * (fullTrkCov * trkJac.transpose())); + BoundMatrix boundTrackCov(BoundMatrix::Identity()); + boundTrackCov.block(0, 0) = + (freeToBoundJac * (freeTrkCov * freeToBoundJac.transpose())); - return fullPerTrackCov; + return boundTrackCov; } diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp index 29d4fe7920d..c7a2b587145 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp @@ -35,15 +35,20 @@ namespace KalmanVertexUpdater { /// Section 5.3.5 /// Cache object, the comments indicate the names of the variables in Ref. (1) +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +template struct Cache { + using VertexVector = ActsVector; + using VertexMatrix = ActsSquareMatrix; // \tilde{x_k} - Vector3 newVertexPos = Vector3::Zero(); + VertexVector newVertexPos = VertexVector::Zero(); // C_k - SquareMatrix3 newVertexCov = SquareMatrix3::Zero(); + VertexMatrix newVertexCov = VertexMatrix::Zero(); // C_k^-1 - SquareMatrix3 newVertexWeight = SquareMatrix3::Zero(); + VertexMatrix newVertexWeight = VertexMatrix::Zero(); // C_{k-1}^-1 - SquareMatrix3 oldVertexWeight = SquareMatrix3::Zero(); + VertexMatrix oldVertexWeight = VertexMatrix::Zero(); // W_k SquareMatrix3 wMat = SquareMatrix3::Zero(); }; @@ -53,15 +58,24 @@ struct Cache { /// However, it does not add the track to the TrackAtVertex list. This to be /// done manually after calling the method. /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param vtx Vertex to be updated /// @param trk Track to be used for updating the vertex -template +template void updateVertexWithTrack(Vertex& vtx, TrackAtVertex& trk); /// @brief Calculates updated vertex position and covariance as well as the /// updated track momentum when adding/removing linTrack. Saves the result in /// cache. +/// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param vtx Vertex /// @param linTrack Linearized track to be added or removed /// @param trackWeight Track weight @@ -69,44 +83,58 @@ void updateVertexWithTrack(Vertex& vtx, /// @note Tracks are removed during the smoothing procedure to compute /// the chi2 of the track wrt the updated vertex position /// @param[out] cache A cache to store the results of this function -template +template void calculateUpdate(const Acts::Vertex& vtx, const Acts::LinearizedTrack& linTrack, - const double trackWeight, const int sign, Cache& cache); + const double trackWeight, const int sign, + Cache& cache); namespace detail { /// @brief Calculates the update of the vertex position chi2 after /// adding/removing the track /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param oldVtx Vertex before the track was added/removed /// @param cache Cache containing updated vertex position /// /// @return Chi2 -template +template double vertexPositionChi2Update(const Vertex& oldVtx, - const Cache& cache); + const Cache& cache); /// @brief Calculates chi2 of refitted track parameters /// w.r.t. updated vertex /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param linTrack Linearized version of track /// @param cache Cache containing some quantities needed in /// this function /// /// @return Chi2 -template -double trackParametersChi2(const LinearizedTrack& linTrack, const Cache& cache); +template +double trackParametersChi2(const LinearizedTrack& linTrack, + const Cache& cache); /// @brief Adds or removes (depending on `sign`) tracks from vertex /// and updates the vertex /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param vtx Vertex to be updated /// @param trk Track to be added to/removed from vtx /// @param sign +1 (add track) or -1 (remove track) /// @note Tracks are removed during the smoothing procedure to compute /// the chi2 of the track wrt the updated vertex position -template +template void update(Vertex& vtx, TrackAtVertex& trk, int sign); } // Namespace detail diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index 8138eb5e7a8..b1886315fd3 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -10,19 +10,25 @@ #include -template +template void Acts::KalmanVertexUpdater::updateVertexWithTrack( Vertex& vtx, TrackAtVertex& trk) { - detail::update(vtx, trk, 1); + detail::update(vtx, trk, 1); } -template +template void Acts::KalmanVertexUpdater::detail::update( Vertex& vtx, TrackAtVertex& trk, int sign) { + if constexpr (nDimVertex != 3 && nDimVertex != 4) { + throw std::invalid_argument( + "The vertex dimension must either be 3 (when fitting the spatial " + "coordinates) or 4 (when fitting the spatial coordinates + time)."); + } + double trackWeight = trk.trackWeight; // Set up cache where entire content is set to 0 - Cache cache; + Cache cache; // Calculate update and save result in cache calculateUpdate(vtx, trk.linearizedState, trackWeight, sign, cache); @@ -47,8 +53,13 @@ void Acts::KalmanVertexUpdater::detail::update( ndf += sign * trackWeight * 2.; // Updating the vertex - vtx.setPosition(cache.newVertexPos); - vtx.setCovariance(cache.newVertexCov); + if constexpr (nDimVertex == 3) { + vtx.setPosition(cache.newVertexPos); + vtx.setCovariance(cache.newVertexCov); + } else if constexpr (nDimVertex == 4) { + vtx.setFullPosition(cache.newVertexPos); + vtx.setFullCovariance(cache.newVertexCov); + } vtx.setFitQuality(chi2, ndf); if (sign == 1) { @@ -65,22 +76,29 @@ void Acts::KalmanVertexUpdater::detail::update( } } -template +template void Acts::KalmanVertexUpdater::calculateUpdate( const Acts::Vertex& vtx, const Acts::LinearizedTrack& linTrack, const double trackWeight, - const int sign, Cache& cache) { + const int sign, Cache& cache) { + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // Retrieve variables from the track linearization. The comments indicate the // corresponding symbol used in Ref. (1). // A_k - const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); // c_k - const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would use time + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time. // G_k // Note that, when removing a track, G_k -> - G_k, see Ref. (1). // Further note that, as we use the weighted formalism, the track weight @@ -88,21 +106,25 @@ void Acts::KalmanVertexUpdater::calculateUpdate( // with the track weight from the AMVF formalism. Here, we choose to // consider these two multiplicative factors directly in the updates of // newVertexWeight and newVertexPos. - const ActsSquareMatrix<5> trkParamWeight = - linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); // Retrieve current position of the vertex and its current weight matrix - const Vector3& oldVtxPos = vtx.position(); + const ActsVector oldVtxPos = + vtx.fullPosition().template head(); // C_{k-1}^-1 - cache.oldVertexWeight = (vtx.covariance()).inverse(); + cache.oldVertexWeight = + (vtx.fullCovariance().template block(0, 0)) + .inverse(); // W_k cache.wMat = (momJac.transpose() * (trkParamWeight * momJac)).inverse(); // G_k^B = G_k - G_k*B_k*W_k*B_k^(T)*G_k - ActsSquareMatrix<5> gBMat = - trkParamWeight - trkParamWeight * momJac * cache.wMat * - momJac.transpose() * trkParamWeight; + ParameterMatrix gBMat = trkParamWeight - trkParamWeight * momJac * + cache.wMat * momJac.transpose() * + trkParamWeight; // C_k^-1 cache.newVertexWeight = cache.oldVertexWeight + sign * trackWeight * @@ -118,33 +140,42 @@ void Acts::KalmanVertexUpdater::calculateUpdate( (trkParams - constTerm)); } -template +template double Acts::KalmanVertexUpdater::detail::vertexPositionChi2Update( - const Vertex& oldVtx, const Cache& cache) { - Vector3 posDiff = cache.newVertexPos - oldVtx.position(); + const Vertex& oldVtx, const Cache& cache) { + ActsVector posDiff = + cache.newVertexPos - oldVtx.fullPosition().template head(); // Calculate and return corresponding chi2 return posDiff.transpose() * (cache.oldVertexWeight * posDiff); } -template +template double Acts::KalmanVertexUpdater::detail::trackParametersChi2( - const LinearizedTrack& linTrack, const Cache& cache) { + const LinearizedTrack& linTrack, const Cache& cache) { + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // A_k - const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); // c_k - const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would use time + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time. // G_k - const ActsSquareMatrix<5> trkParamWeight = - linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); // A_k * \tilde{x_k} - const ActsVector<5> posJacVtxPos = posJac * cache.newVertexPos; + const ParameterVector posJacVtxPos = posJac * cache.newVertexPos; // \tilde{q_k} Vector3 newTrkMom = cache.wMat * momJac.transpose() * trkParamWeight * @@ -161,11 +192,11 @@ double Acts::KalmanVertexUpdater::detail::trackParametersChi2( */ // \tilde{p_k} - ActsVector<5> linearizedTrackParameters = + ParameterVector linearizedTrackParameters = constTerm + posJacVtxPos + momJac * newTrkMom; // r_k - ActsVector<5> paramDiff = trkParams - linearizedTrackParameters; + ParameterVector paramDiff = trkParams - linearizedTrackParameters; // Return chi2 return paramDiff.transpose() * (trkParamWeight * paramDiff); diff --git a/Core/src/Detector/Portal.cpp b/Core/src/Detector/Portal.cpp index cda979faf01..2e1123d03e6 100644 --- a/Core/src/Detector/Portal.cpp +++ b/Core/src/Detector/Portal.cpp @@ -18,93 +18,103 @@ #include #include -namespace Acts { -namespace Experimental { -class DetectorVolume; -} // namespace Experimental -} // namespace Acts +namespace Acts::Experimental { -Acts::Experimental::Portal::Portal(std::shared_ptr surface) +Portal::Portal(std::shared_ptr surface) : m_surface(std::move(surface)) { throw_assert(m_surface, "Portal surface is nullptr"); } -std::shared_ptr -Acts::Experimental::Portal::makeShared( - std::shared_ptr surface) { - return std::shared_ptr(new Portal(std::move(surface))); -} - -const Acts::RegularSurface& Acts::Experimental::Portal::surface() const { +const Acts::RegularSurface& Portal::surface() const { return *m_surface.get(); } -Acts::RegularSurface& Acts::Experimental::Portal::surface() { +Acts::RegularSurface& Portal::surface() { return *m_surface.get(); } -const Acts::Experimental::Portal::DetectorVolumeUpdaters& -Acts::Experimental::Portal::detectorVolumeUpdaters() const { +const Portal::DetectorVolumeUpdaters& Portal::detectorVolumeUpdaters() const { return m_volumeUpdaters; } -Acts::Experimental::Portal::AttachedDetectorVolumes& -Acts::Experimental::Portal::attachedDetectorVolumes() { +Portal::AttachedDetectorVolumes& Portal::attachedDetectorVolumes() { return m_attachedVolumes; } -std::shared_ptr -Acts::Experimental::Portal::getSharedPtr() { - return shared_from_this(); +void Portal::assignGeometryId(const GeometryIdentifier& geometryId) { + m_surface->assignGeometryId(geometryId); } -std::shared_ptr -Acts::Experimental::Portal::getSharedPtr() const { - return shared_from_this(); -} +std::shared_ptr Portal::fuse(std::shared_ptr& aPortal, + std::shared_ptr& bPortal) { + if (aPortal == bPortal) { + return aPortal; + } -void Acts::Experimental::Portal::assignGeometryId( - const GeometryIdentifier& geometryId) { - m_surface->assignGeometryId(geometryId); -} + auto bothConnected = [](const auto& p) { + return p.m_volumeUpdaters[0u].connected() && + p.m_volumeUpdaters[1u].connected(); + }; -void Acts::Experimental::Portal::fuse(std::shared_ptr& other) { - Direction bDir = Direction::Backward; + auto noneConnected = [](const auto& p) { + return !p.m_volumeUpdaters[0u].connected() && + !p.m_volumeUpdaters[1u].connected(); + }; - // Determine this directioon - Direction tDir = (!m_volumeUpdaters[bDir.index()].connected()) - ? Direction::Forward - : Direction::Backward; + if (bothConnected(*aPortal) || bothConnected(*bPortal)) { + throw std::invalid_argument( + "Portal: trying to fuse two portals where at least one has links on " + "both sides."); + } - if (!m_volumeUpdaters[tDir.index()].connected()) { + if (noneConnected(*aPortal) || noneConnected(*bPortal)) { throw std::invalid_argument( - "Portal: trying to fuse portal (keep) with no links."); + "Portal: trying to fuse two portals where at least on has no links."); } + + // We checked they're not both empty, so one of them must be connected + Direction aDir = (aPortal->m_volumeUpdaters[0].connected()) + ? Direction::fromIndex(0) + : Direction::fromIndex(1); + Direction bDir = aDir.invert(); + // And now check other direction - Direction oDir = tDir.invert(); - if (!other->m_volumeUpdaters[oDir.index()].connected()) { + if (!bPortal->m_volumeUpdaters[bDir.index()].connected()) { throw std::runtime_error( - "Portal: trying to fuse portal (waste) with no links."); + "Portal: trying to fuse portal (discard) with no links."); } - if (m_surface->surfaceMaterial() != nullptr && - other->surface().surfaceMaterial() != nullptr) { + const auto& aSurface = aPortal->surface(); + const auto& bSurface = bPortal->surface(); + + // @TODO: There's no safety against fusing portals with different surfaces + std::shared_ptr fused = std::make_shared(aPortal->m_surface); + + if (aSurface.surfaceMaterial() != nullptr && + bSurface.surfaceMaterial() != nullptr) { throw std::runtime_error( "Portal: both surfaces have surface material, fusing will lead to " "information loss."); - } else if (other->surface().surfaceMaterial() != nullptr) { - m_surface->assignSurfaceMaterial( - other->surface().surfaceMaterialSharedPtr()); + } else if (aSurface.surfaceMaterial() != nullptr) { + fused->m_surface = aPortal->m_surface; + } else if (bSurface.surfaceMaterial() != nullptr) { + fused->m_surface = bPortal->m_surface; } - auto odx = oDir.index(); - m_volumeUpdaters[odx] = std::move(other->m_volumeUpdaters[odx]); - m_attachedVolumes[odx] = other->m_attachedVolumes[odx]; - // And finally overwrite the original portal - other = getSharedPtr(); + fused->m_volumeUpdaters[aDir.index()] = + std::move(aPortal->m_volumeUpdaters[aDir.index()]); + fused->m_attachedVolumes[aDir.index()] = + std::move(aPortal->m_attachedVolumes[aDir.index()]); + + fused->m_volumeUpdaters[bDir.index()] = + std::move(bPortal->m_volumeUpdaters[bDir.index()]); + fused->m_attachedVolumes[bDir.index()] = + std::move(bPortal->m_attachedVolumes[bDir.index()]); + + return fused; } -void Acts::Experimental::Portal::assignDetectorVolumeUpdater( +void Portal::assignDetectorVolumeUpdater( Direction dir, DetectorVolumeUpdater dVolumeUpdater, std::vector> attachedVolumes) { auto idx = dir.index(); @@ -112,7 +122,7 @@ void Acts::Experimental::Portal::assignDetectorVolumeUpdater( m_attachedVolumes[idx] = std::move(attachedVolumes); } -void Acts::Experimental::Portal::assignDetectorVolumeUpdater( +void Portal::assignDetectorVolumeUpdater( DetectorVolumeUpdater dVolumeUpdater, std::vector> attachedVolumes) { // Check and throw exceptions @@ -127,8 +137,8 @@ void Acts::Experimental::Portal::assignDetectorVolumeUpdater( m_attachedVolumes[idx] = std::move(attachedVolumes); } -void Acts::Experimental::Portal::updateDetectorVolume( - const GeometryContext& gctx, NavigationState& nState) const { +void Portal::updateDetectorVolume(const GeometryContext& gctx, + NavigationState& nState) const { const auto& position = nState.position; const auto& direction = nState.direction; const Vector3 normal = surface().normal(gctx, position); @@ -140,3 +150,5 @@ void Acts::Experimental::Portal::updateDetectorVolume( nState.currentVolume = nullptr; } } + +} // namespace Acts::Experimental diff --git a/Core/src/Detector/PortalGenerators.cpp b/Core/src/Detector/PortalGenerators.cpp index d950054ca90..f846c710475 100644 --- a/Core/src/Detector/PortalGenerators.cpp +++ b/Core/src/Detector/PortalGenerators.cpp @@ -33,7 +33,7 @@ Acts::Experimental::generatePortals( std::vector> portals; for (auto [i, oSurface] : enumerate(orientedSurfaces)) { // Create a portal from the surface - auto portal = Portal::makeShared(oSurface.first); + auto portal = std::make_shared(oSurface.first); // Create a shared link instance & delegate auto singleLinkImpl = std::make_unique(dVolume.get()); diff --git a/Core/src/Detector/detail/CylindricalDetectorHelper.cpp b/Core/src/Detector/detail/CylindricalDetectorHelper.cpp index c211dd7be9d..ce7973b3889 100644 --- a/Core/src/Detector/detail/CylindricalDetectorHelper.cpp +++ b/Core/src/Detector/detail/CylindricalDetectorHelper.cpp @@ -99,7 +99,7 @@ Acts::Experimental::PortalReplacement createDiscReplacement( const auto& stitchBoundaries = (stitchValue == Acts::binR) ? rBoundaries : phiBoundaries; return Acts::Experimental::PortalReplacement( - Acts::Experimental::Portal::makeShared(surface), index, dir, + std::make_shared(surface), index, dir, stitchBoundaries, stitchValue); } @@ -135,7 +135,7 @@ Acts::Experimental::PortalReplacement createCylinderReplacement( const auto& stitchBoundaries = (stitchValue == Acts::binZ) ? zBoundaries : phiBoundaries; return Acts::Experimental::PortalReplacement( - Acts::Experimental::Portal::makeShared(surface), index, dir, + std::make_shared(surface), index, dir, stitchBoundaries, stitchValue); } @@ -200,8 +200,8 @@ Acts::Experimental::PortalReplacement createSectorReplacement( transform, std::move(bounds)); // A make a portal and indicate the new link direction Acts::Experimental::PortalReplacement pRep = { - Acts::Experimental::Portal::makeShared(surface), index, dir, boundaries, - binning}; + std::make_shared(surface), index, dir, + boundaries, binning}; return pRep; } @@ -413,7 +413,7 @@ Acts::Experimental::detail::CylindricalDetectorHelper::connectInR( refValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector]; ActsScalar avgPhi = refValues[CylinderVolumeBounds::BoundValues::eAveragePhi]; - // Fuse the cylinders - portals can be reused for this operation + // Fuse the cylinders for (unsigned int iv = 1; iv < volumes.size(); ++iv) { refValues = volumes[iv]->volumeBounds().values(); // Keep on collecting the outside maximum r for the overall r boundaries @@ -423,15 +423,12 @@ Acts::Experimental::detail::CylindricalDetectorHelper::connectInR( ACTS_VERBOSE("Connect volume '" << volumes[iv - 1]->name() << "' to " << volumes[iv]->name() << "'."); - // When fusing volumes at a cylinder boundary, we *keep* one - // portal and transfer the portal link information from the other - // - // In this case the outer cylinder portal of the inner volume is kept, - // the inner cylinder of the outer portal goes to waste - auto& keepCylinder = volumes[iv - 1]->portalPtrs()[2u]; - auto& wasteCylinder = volumes[iv]->portalPtrs()[3u]; - keepCylinder->fuse(wasteCylinder); - volumes[iv]->updatePortal(keepCylinder, 3u); + // Fusing cylinders from inner and outer volume + auto innerCylinder = volumes[iv - 1]->portalPtrs()[2u]; + auto outerCylinder = volumes[iv]->portalPtrs()[3u]; + auto fusedCylinder = Portal::fuse(innerCylinder, outerCylinder); + volumes[iv - 1]->updatePortal(fusedCylinder, 2u); + volumes[iv]->updatePortal(fusedCylinder, 3u); } } @@ -583,30 +580,26 @@ Acts::Experimental::detail::CylindricalDetectorHelper::connectInZ( if (connectZ) { ACTS_VERBOSE("Connect volume '" << volumes[iv - 1]->name() << "' to " << volumes[iv]->name() << "'."); - // When fusing, one portal survives (keep) and gets the - // portal linking from the waste transferred - // - // In this case we keep the disc at positive z of the volume - // at lower relative z, and trash the disc at negative z of the - // following volume - auto& keepDisc = volumes[iv - 1]->portalPtrs()[1u]; - auto& wasteDisc = volumes[iv]->portalPtrs()[0u]; + // Fusing the discs: positive at lower z, negative at hgiher z + auto& pDisc = volumes[iv - 1]->portalPtrs()[1u]; + auto& nDisc = volumes[iv]->portalPtrs()[0u]; // Throw an exception if the discs are not at the same position - Vector3 keepPosition = keepDisc->surface().center(gctx); - Vector3 wastePosition = wasteDisc->surface().center(gctx); - if (!keepPosition.isApprox(wastePosition)) { + Vector3 pPosition = pDisc->surface().center(gctx); + Vector3 nPosition = nDisc->surface().center(gctx); + if (!pPosition.isApprox(nPosition)) { std::string message = "CylindricalDetectorHelper: '"; message += volumes[iv - 1]->name(); message += "' does not attach to '"; message += volumes[iv]->name(); message += "'\n"; message += " - along z with values "; - message += Acts::toString(keepPosition); - message += " / " + Acts::toString(wastePosition); + message += Acts::toString(pPosition); + message += " / " + Acts::toString(nPosition); throw std::runtime_error(message.c_str()); } - keepDisc->fuse(wasteDisc); - volumes[iv]->updatePortal(keepDisc, 0u); + auto fusedDisc = Portal::fuse(pDisc, nDisc); + volumes[iv - 1]->updatePortal(fusedDisc, 1u); + volumes[iv]->updatePortal(fusedDisc, 0u); } } @@ -758,16 +751,17 @@ Acts::Experimental::detail::CylindricalDetectorHelper::connectInPhi( phiBoundaries.push_back( refValues[CylinderVolumeBounds::BoundValues::eAveragePhi] + refValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector]); - // Fuse the sectors - portals can be reused for this operation + // Fuse the sectors for (unsigned int iv = 1; iv < volumes.size(); ++iv) { ACTS_VERBOSE("Connect volume '" << volumes[iv - 1]->name() << "' to " << volumes[iv]->name() << "'."); - // Fuse and swap - auto& keepSector = volumes[iv - 1]->portalPtrs()[iSecOffset + 1u]; - auto& wasteSector = volumes[iv]->portalPtrs()[iSecOffset]; - keepSector->fuse(wasteSector); - volumes[iv]->updatePortal(keepSector, iSecOffset); + // Fuse sector surfaces r handed at lower index, l handed at higher index + auto& rSector = volumes[iv - 1]->portalPtrs()[iSecOffset + 1u]; + auto& lSector = volumes[iv]->portalPtrs()[iSecOffset]; + auto fusedSector = Portal::fuse(rSector, lSector); + volumes[iv - 1]->updatePortal(fusedSector, iSecOffset + 1u); + volumes[iv]->updatePortal(fusedSector, iSecOffset); // The current values auto curValues = volumes[iv]->volumeBounds().values(); // Bail out if they do not match @@ -875,22 +869,25 @@ Acts::Experimental::detail::CylindricalDetectorHelper::wrapInZR( dShell[2u] = volumes[1u]->portalPtrs()[2u]; // Fuse outer cover of first with inner cylinder of wrapping volume - auto& keepCover = volumes[0u]->portalPtrs()[2u]; - auto& wasteCover = volumes[1u]->portalPtrs()[3u]; - keepCover->fuse(wasteCover); - volumes[1u]->updatePortal(keepCover, 3u); + auto& outerCover = volumes[0u]->portalPtrs()[2u]; + auto& innerCover = volumes[1u]->portalPtrs()[3u]; + auto fusedCover = Portal::fuse(outerCover, innerCover); + volumes[0u]->updatePortal(fusedCover, 2u); + volumes[1u]->updatePortal(fusedCover, 3u); // Stitch sides - negative - auto& keepDiscN = volumes[1u]->portalPtrs()[4u]; - auto& wasteDiscN = volumes[0u]->portalPtrs()[0u]; - keepDiscN->fuse(wasteDiscN); - volumes[0u]->updatePortal(keepDiscN, 0u); + auto& firstDiscN = volumes[1u]->portalPtrs()[4u]; + auto& secondDiscN = volumes[0u]->portalPtrs()[0u]; + auto fusedDiscN = Portal::fuse(firstDiscN, secondDiscN); + volumes[1u]->updatePortal(fusedDiscN, 4u); + volumes[0u]->updatePortal(fusedDiscN, 0u); // Stich sides - positive - auto& keepDiscP = volumes[0u]->portalPtrs()[1u]; - auto& wasteDiscP = volumes[1u]->portalPtrs()[5u]; - keepDiscP->fuse(wasteDiscP); - volumes[1u]->updatePortal(keepDiscP, 5u); + auto& firstDiscP = volumes[0u]->portalPtrs()[1u]; + auto& secondDiscP = volumes[1u]->portalPtrs()[5u]; + auto fusedDiscP = Portal::fuse(firstDiscP, secondDiscP); + volumes[0u]->updatePortal(fusedDiscP, 1u); + volumes[1u]->updatePortal(fusedDiscP, 5u); // If needed, insert new cylinder if (volumes[0u]->portalPtrs().size() == 4u && @@ -960,13 +957,27 @@ Acts::Experimental::detail::CylindricalDetectorHelper::connectInR( "not be connected in R"); } - // Fuse and swap - std::shared_ptr keepCylinder = containers[ic - 1].find(2u)->second; - std::shared_ptr wasteCylinder = containers[ic].find(3u)->second; - keepCylinder->fuse(wasteCylinder); - for (auto& av : wasteCylinder->attachedDetectorVolumes()[1u]) { - av->updatePortal(keepCylinder, 3u); - } + // Fuse containers, and update the attached volumes + std::shared_ptr innerCylinder = containers[ic - 1].find(2u)->second; + // Direction is explicitly addressed with a direction index + auto innerAttachedVolumes = + innerCylinder + ->attachedDetectorVolumes()[Direction(Direction::Backward).index()]; + std::shared_ptr outerCylinder = containers[ic].find(3u)->second; + auto outerAttachedVolume = + outerCylinder + ->attachedDetectorVolumes()[Direction(Direction::Forward).index()]; + auto fusedCylinder = Portal::fuse(innerCylinder, outerCylinder); + + // Update the attached volumes with the new portal + std::for_each(innerAttachedVolumes.begin(), innerAttachedVolumes.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedCylinder, 2u); + }); + std::for_each(outerAttachedVolume.begin(), outerAttachedVolume.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedCylinder, 3u); + }); } // Proto container refurbishment @@ -1017,13 +1028,26 @@ Acts::Experimental::detail::CylindricalDetectorHelper::connectInZ( "CylindricalDetectorHelper: proto container has no positive disc, " "can not be connected in Z"); } - std::shared_ptr keepDisc = formerContainer.find(1u)->second; - std::shared_ptr wasteDisc = currentContainer.find(0u)->second; - keepDisc->fuse(wasteDisc); - for (auto& av : wasteDisc->attachedDetectorVolumes()[1u]) { - ACTS_VERBOSE("Update portal of detector volume '" << av->name() << "'."); - av->updatePortal(keepDisc, 0u); - } + // Container attachment positive Disc of lower, negative Disc at higher + std::shared_ptr pDisc = formerContainer.find(1u)->second; + auto pAttachedVolumes = + pDisc + ->attachedDetectorVolumes()[Direction(Direction::Backward).index()]; + + std::shared_ptr nDisc = currentContainer.find(0u)->second; + auto nAttachedVolumes = + nDisc->attachedDetectorVolumes()[Direction(Direction::Forward).index()]; + + auto fusedDisc = Portal::fuse(pDisc, nDisc); + + std::for_each(pAttachedVolumes.begin(), pAttachedVolumes.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedDisc, 1u); + }); + std::for_each(nAttachedVolumes.begin(), nAttachedVolumes.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedDisc, 0u); + }); } // Proto container refurbishment @@ -1119,22 +1143,51 @@ Acts::Experimental::detail::CylindricalDetectorHelper::wrapInZR( dShell[2u] = wrappingVolume->portalPtrs()[2u]; // Fuse outer cover of first with inner cylinder of wrapping volume - auto& keepCover = innerContainer[2u]; - auto& wasteCover = wrappingVolume->portalPtrs()[3u]; - keepCover->fuse(wasteCover); - wrappingVolume->updatePortal(keepCover, 3u); + auto& innerCover = innerContainer[2u]; + auto innerAttachedVolumes = + innerCover + ->attachedDetectorVolumes()[Direction(Direction::Backward).index()]; + auto& innerTube = wrappingVolume->portalPtrs()[3u]; + auto fusedCover = Portal::fuse(innerCover, innerTube); + + std::for_each(innerAttachedVolumes.begin(), innerAttachedVolumes.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedCover, 2u); + }); + wrappingVolume->updatePortal(fusedCover, 3u); // Stitch sides - negative - auto& keepDiscN = innerContainer[0u]; - auto& wasteDiscN = wrappingVolume->portalPtrs()[4u]; - keepDiscN->fuse(wasteDiscN); - wrappingVolume->updatePortal(keepDiscN, 4u); + // positive disc of lower , negative disc of higher + auto& firstDiscN = innerContainer[0u]; + + auto firstNAttachedVolumes = + firstDiscN + ->attachedDetectorVolumes()[Direction(Direction::Forward).index()]; + + auto& secondDiscN = wrappingVolume->portalPtrs()[4u]; + auto fusedDiscN = Portal::fuse(firstDiscN, secondDiscN); + + std::for_each(firstNAttachedVolumes.begin(), firstNAttachedVolumes.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedDiscN, 0u); + }); + wrappingVolume->updatePortal(fusedDiscN, 4u); // Stich sides - positive - auto& keepDiscP = innerContainer[1u]; - auto& wasteDiscP = wrappingVolume->portalPtrs()[5u]; - keepDiscP->fuse(wasteDiscP); - wrappingVolume->updatePortal(keepDiscP, 5u); + auto& firstDiscP = innerContainer[1u]; + auto firstPAttachedVolumes = + firstDiscP + ->attachedDetectorVolumes()[Direction(Direction::Backward).index()]; + + auto& secondDiscP = wrappingVolume->portalPtrs()[5u]; + auto fusedDiscP = Portal::fuse(firstDiscP, secondDiscP); + + std::for_each(firstPAttachedVolumes.begin(), firstPAttachedVolumes.end(), + [&](std::shared_ptr& av) { + av->updatePortal(fusedDiscP, 1u); + }); + + wrappingVolume->updatePortal(fusedDiscP, 5u); // If inner stitching is necessary if (innerContainer.size() == 4u && diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp index af9f111446e..bc282f551b2 100644 --- a/Core/src/Geometry/TrackingVolume.cpp +++ b/Core/src/Geometry/TrackingVolume.cpp @@ -504,9 +504,7 @@ Acts::TrackingVolume::compatibleBoundaries( ACTS_VERBOSE("- intersection path length " << std::abs(sIntersection.pathLength()) << " < overstep limit " << std::abs(oLimit)); - return BoundaryIntersection(sIntersection.intersection(), bSurface, - sIntersection.object(), - sIntersection.index()); + return BoundaryIntersection(sIntersection, bSurface); } ACTS_VERBOSE("Can't force intersection: "); ACTS_VERBOSE("- intersection path length " @@ -520,14 +518,12 @@ Acts::TrackingVolume::compatibleBoundaries( decltype(logger)>( sIntersection.intersection(), pLimit, oLimit, s_onSurfaceTolerance, logger)) { - return BoundaryIntersection(sIntersection.intersection(), bSurface, - sIntersection.object(), - sIntersection.index()); + return BoundaryIntersection(sIntersection, bSurface); } } ACTS_VERBOSE("No intersection accepted"); - return BoundaryIntersection::invalid(); + return BoundaryIntersection(SurfaceIntersection::invalid(), bSurface); }; /// Helper function to process boundary surfaces @@ -548,7 +544,7 @@ Acts::TrackingVolume::compatibleBoundaries( options.boundaryCheck); // Intersect and continue auto bIntersection = checkIntersection(bCandidate, bsIter.get()); - if (bIntersection) { + if (bIntersection.first) { ACTS_VERBOSE(" - Proceed with surface"); bIntersections.push_back(bIntersection); } else { @@ -587,8 +583,8 @@ Acts::TrackingVolume::compatibleBoundaries( }; std::sort(bIntersections.begin(), bIntersections.end(), - [&](const auto& a, const auto& b) { - return comparator(a.pathLength(), b.pathLength()); + [&](const BoundaryIntersection& a, const BoundaryIntersection& b) { + return comparator(a.first.pathLength(), b.first.pathLength()); }); return bIntersections; } @@ -623,9 +619,7 @@ Acts::TrackingVolume::compatibleLayers( if (atIntersection && (atIntersection.object() != options.targetSurface) && withinLimit) { // create a layer intersection - lIntersections.push_back(LayerIntersection( - atIntersection.intersection(), tLayer, atIntersection.object(), - atIntersection.index())); + lIntersections.push_back(LayerIntersection(atIntersection, tLayer)); } } // move to next one or break because you reached the end layer @@ -634,7 +628,9 @@ Acts::TrackingVolume::compatibleLayers( : tLayer->nextLayer(gctx, position, direction); } std::sort(lIntersections.begin(), lIntersections.end(), - LayerIntersection::forwardOrder); + [](const LayerIntersection& a, const LayerIntersection& b) { + return SurfaceIntersection::forwardOrder(a.first, b.first); + }); } // and return return lIntersections; diff --git a/Core/src/TrackFinding/FasTrackConnector.cpp b/Core/src/TrackFinding/FasTrackConnector.cpp index 4bde548f85a..0d4ad1bacd6 100644 --- a/Core/src/TrackFinding/FasTrackConnector.cpp +++ b/Core/src/TrackFinding/FasTrackConnector.cpp @@ -193,5 +193,4 @@ FasTrackConnector::~FasTrackConnector() { } } } - } // namespace Acts diff --git a/Core/src/TrackFinding/RoiDescriptor.cpp b/Core/src/TrackFinding/RoiDescriptor.cpp index 4ae7288a7af..c76df027082 100644 --- a/Core/src/TrackFinding/RoiDescriptor.cpp +++ b/Core/src/TrackFinding/RoiDescriptor.cpp @@ -22,10 +22,21 @@ Acts::RoiDescriptor::RoiDescriptor(double eta, double etaMinus, double etaPlus, m_zed(zed), m_phiMinus(phiMinus), m_phiPlus(phiPlus), - m_etaMinus(etaMinus), + m_etaMinus(etaMinus), //-4.5 m_etaPlus(etaPlus), m_zedMinus(zedMinus), - m_zedPlus(zedPlus) {} + m_zedPlus(zedPlus) { + // catch in the athena roi code + // if ( std::isnan(m_etaPlus) ) throw std::invalid_argument( "RoiDescriptor: + // etaPlus nan" ); if ( std::isnan(m_etaMinus) ) throw std::invalid_argument( + // "RoiDescriptor: etaMinus nan" ); + + m_drdzMinus = std::tan(2 * std::atan(std::exp(-m_etaMinus))); //-0.02 + m_drdzPlus = std::tan(2 * std::atan(std::exp(-m_etaPlus))); // 0.02 + + m_dzdrMinus = 1 / m_drdzMinus; //-45 + m_dzdrPlus = 1 / m_drdzPlus; // 45 +} Acts::RoiDescriptor::~RoiDescriptor() = default; diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp index a1022d61566..c67999e573b 100644 --- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp @@ -14,6 +14,8 @@ #include "Acts/Seeding/SeedFilterConfig.hpp" #include "Acts/Seeding/SeedFinderFTF.hpp" #include "Acts/Seeding/SeedFinderFTFConfig.hpp" +#include "ActsExamples/EventData/Cluster.hpp" + // in core #include "Acts/Geometry/TrackingGeometry.hpp" #include "Acts/Seeding/SeedFinderConfig.hpp" @@ -50,6 +52,8 @@ class SeedingFTFAlgorithm final : public IAlgorithm { std::map, std::pair> ACTS_FTF_Map; bool fill_module_csv = false; + + std::string inputClusters; }; // constructor: @@ -87,6 +91,8 @@ class SeedingFTFAlgorithm final : public IAlgorithm { ReadDataHandle m_inputSourceLinks{ this, "InputSourceLinks"}; + + ReadDataHandle m_inputClusters{this, "InputClusters"}; }; } // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp index 4766127ad73..a092f75f4e6 100644 --- a/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp @@ -64,6 +64,8 @@ ActsExamples::SeedingFTFAlgorithm::SeedingFTFAlgorithm( m_inputSourceLinks.initialize(m_cfg.inputSourceLinks); + m_inputClusters.initialize(m_cfg.inputClusters); + m_cfg.seedFinderConfig.seedFilter = std::make_unique>( Acts::SeedFilter(m_cfg.seedFilterConfig)); @@ -92,6 +94,17 @@ ActsExamples::ProcessCode ActsExamples::SeedingFTFAlgorithm::execute( std::vector> FTF_spacePoints = Make_FTF_spacePoints(ctx, m_cfg.ACTS_FTF_Map); + // cluster width + // const ClusterContainer* clusters = &m_inputClusters(ctx) ; + + // for (const auto& sp : FTF_spacePoints){ + // const auto& sl = sp.SP->sourceLinks().front().get() ; + // const auto& cluster = clusters->at(sl.index()) ; + // std::cout << "testing 0: " << cluster.sizeLoc0 << " 1: " << + // cluster.sizeLoc1 << std::endl ; + + // } + for (auto sp : FTF_spacePoints) { ACTS_DEBUG("FTF space points: " << " FTF_id: " << sp.FTF_ID << " z: " << sp.SP->z() @@ -119,14 +132,15 @@ ActsExamples::ProcessCode ActsExamples::SeedingFTFAlgorithm::execute( finder.loadSpacePoints(FTF_spacePoints); + // trigFTF file : Acts::RoiDescriptor internalRoi(0, -4.5, 4.5, 0, -M_PI, M_PI, 0, -150.0, 150.0); + // ROI file: + // Acts::RoiDescriptor internalRoi(0, -5, 5, 0, -M_PI, M_PI, 0, -225.0, + // 225.0); - finder.createSeeds(internalRoi, *mGNNgeo); - - // still to develop - SimSeedContainer seeds = finder.createSeeds_old( - m_cfg.seedFinderOptions, FTF_spacePoints, create_coordinates); + // new version returns seeds + SimSeedContainer seeds = finder.createSeeds(internalRoi, *mGNNgeo); m_outputSeeds(ctx, std::move(seeds)); @@ -286,7 +300,8 @@ ActsExamples::SeedingFTFAlgorithm::LayerNumbering() const { std::make_pair(ACTS_joint_id, 0); // here the key needs to be pair of(vol*100+lay, 0) auto Find = m_cfg.ACTS_FTF_Map.find(key); - int FTF_id = Find->second.first; // new map, item is pair want first + int FTF_id = 0; // initialise first to avoid FLTUND later + FTF_id = Find->second.first; // new map, item is pair want first if (Find == m_cfg.ACTS_FTF_Map .end()) { // if end then make new key of (vol*100+lay, modid) diff --git a/Examples/Io/Json/include/ActsExamples/Io/Json/JsonSpacePointWriter.hpp b/Examples/Io/Json/include/ActsExamples/Io/Json/JsonSpacePointWriter.hpp deleted file mode 100644 index df6fff26baa..00000000000 --- a/Examples/Io/Json/include/ActsExamples/Io/Json/JsonSpacePointWriter.hpp +++ /dev/null @@ -1,111 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2017 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/. - -/// @file -/// @date 2016-05-23 Initial version -/// @date 2017-08-07 Rewrite with new interfaces - -#pragma once - -#include "ActsExamples/EventData/GeometryContainers.hpp" -#include "ActsExamples/Framework/WriterT.hpp" -#include "ActsExamples/Utilities/Paths.hpp" - -#include - -namespace ActsExamples { - -/// Write out a space point collection in JSON format. -/// -/// This writes one file per event into the configured output directory. By -/// default it writes to the current working directory. Files are named -/// using the following schema -/// -/// event000000001-spacepoints.json -/// event000000002-spacepoints.json -/// -template -class JsonSpacePointWriter : public WriterT> { - public: - struct Config { - std::string collection; ///< which collection to write - std::string outputDir; ///< where to place output files - std::size_t outputPrecision = 6; ///< floating point precision - }; - - JsonSpacePointWriter(const Config& cfg, - Acts::Logging::Level level = Acts::Logging::INFO); - - protected: - ActsExamples::ProcessCode writeT( - const ActsExamples::AlgorithmContext& context, - const GeometryIdMultimap& spacePoints) override; - - private: - // since class itself is templated, base class template must be fixed - using Base = WriterT>; - - Config m_cfg; -}; - -} // namespace ActsExamples - -template -ActsExamples::JsonSpacePointWriter::JsonSpacePointWriter( - const ActsExamples::JsonSpacePointWriter::Config& cfg, - Acts::Logging::Level level) - : Base(cfg.collection, "JsonSpacePointWriter", level), m_cfg(cfg) { - if (m_cfg.collection.empty()) { - throw std::invalid_argument("Missing input collection"); - } -} - -template -ActsExamples::ProcessCode ActsExamples::JsonSpacePointWriter::writeT( - const ActsExamples::AlgorithmContext& context, - const GeometryIdMultimap& spacePoints) { - // open per-event file - std::string path = perEventFilepath(m_cfg.outputDir, "spacepoints.json", - context.eventNumber); - std::ofstream os(path, std::ofstream::out | std::ofstream::trunc); - if (!os) { // NOLINT(readability-implicit-bool-conversion) - throw std::ios_base::failure("Could not open '" + path + "' to write"); - } - - os << std::setprecision(m_cfg.outputPrecision); - os << "{\n"; - - bool firstVolume = true; - for (auto& volumeData : spacePoints) { - geo_id_value volumeID = volumeData.first; - - if (!firstVolume) - os << ",\n"; - os << " \"SpacePoints_" << volumeID << "\" : [\n"; - - bool firstPoint = true; - for (auto& layerData : volumeData.second) { - for (auto& moduleData : layerData.second) { - for (auto& data : moduleData.second) { - // set the comma correctly - if (!firstPoint) - os << ",\n"; - // write the space point - os << " [" << data.x() << ", " << data.y() << ", " << data.z() - << "]"; - firstPoint = false; - } - } - } - os << "]"; - firstVolume = false; - } - os << "\n}\n"; - - return ProcessCode::SUCCESS; -} diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjSpacePointWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjSpacePointWriter.hpp deleted file mode 100644 index ef98f7d3eed..00000000000 --- a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjSpacePointWriter.hpp +++ /dev/null @@ -1,94 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2017-2018 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/. - -#pragma once - -#include "ActsExamples/EventData/GeometryContainers.hpp" -#include "ActsExamples/Framework/WriterT.hpp" -#include "ActsExamples/Utilities/Paths.hpp" - -#include - -namespace ActsExamples { - -/// Write out a space point collection in OBJ format. -/// -/// This writes one file per event into the configured output directory. By -/// default it writes to the current working directory. Files are named -/// using the following schema -/// -/// event000000001-spacepoints.obj -/// event000000002-spacepoints.obj -/// -/// One write call per thread and hence thread safe. -template -class ObjSpacePointWriter : public WriterT> { - public: - struct Config { - std::string collection; ///< which collection to write - std::string outputDir; ///< where to place output files - double outputScalor = 1.0; ///< scale output values - std::size_t outputPrecision = 6; ///< floating point precision - }; - - ObjSpacePointWriter(const Config& cfg, - Acts::Logging::Level level = Acts::Logging::INFO); - - protected: - ProcessCode writeT(const AlgorithmContext& context, - const GeometryIdMultimap& spacePoints); - - private: - // since class iitself is templated, base class template must be fixed - using Base = WriterT>; - - Config m_cfg; -}; - -} // namespace ActsExamples - -template -inline ActsExamples::ObjSpacePointWriter::ObjSpacePointWriter( - const ObjSpacePointWriter::Config& cfg, Acts::Logging::Level level) - : Base(cfg.collection, "ObjSpacePointWriter", level), m_cfg(cfg) { - if (m_cfg.collection.empty()) { - throw std::invalid_argument("Missing input collection"); - } -} - -template -inline ActsExamples::ProcessCode ActsExamples::ObjSpacePointWriter::writeT( - const ActsExamples::AlgorithmContext& context, - const ActsExamples::GeometryIdMultimap& spacePoints) { - // open per-event file - std::string path = ActsExamples::perEventFilepath( - m_cfg.outputDir, "spacepoints.obj", context.eventNumber); - std::ofstream os(path, std::ofstream::out | std::ofstream::trunc); - if (!os) { - throw std::ios_base::failure("Could not open '" + path + "' to write"); - } - - os << std::setprecision(m_cfg.outputPrecision); - // count the vertex - std::size_t vertex = 0; - // loop and fill the space point data - for (auto& volumeData : spacePoints) { - for (auto& layerData : volumeData.second) { - for (auto& moduleData : layerData.second) { - for (auto& data : moduleData.second) { - // write the space point - os << "v " << m_cfg.outputScalor * data.x() << " " - << m_cfg.outputScalor * data.y() << " " - << m_cfg.outputScalor * data.z() << '\n'; - os << "p " << ++vertex << '\n'; - } - } - } - } - return ProcessCode::SUCCESS; -} diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp index 99c53c97635..39ed0a986b9 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp @@ -154,6 +154,8 @@ ActsExamples::VertexPerformanceWriter::VertexPerformanceWriter( m_outputTree->Branch("trk_pullQOverP", &m_pullQOverP); m_outputTree->Branch("trk_pullQOverPFitted", &m_pullQOverPFitted); + m_outputTree->Branch("trk_weight", &m_trkWeight); + m_outputTree->Branch("nTracksTruthVtx", &m_nTracksOnTruthVertex); m_outputTree->Branch("nTracksRecoVtx", &m_nTracksOnRecoVertex); @@ -489,6 +491,8 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT( auto& innerPullThetaFitted = m_pullThetaFitted.emplace_back(); auto& innerPullQOverPFitted = m_pullQOverPFitted.emplace_back(); + auto& innerTrkWeight = m_trkWeight.emplace_back(); + for (std::size_t j = 0; j < associatedTruthParticles.size(); ++j) { const auto& particle = associatedTruthParticles[j]; int priVtxId = particle.particleId().vertexPrimary(); @@ -687,6 +691,8 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT( innerPullQOverPFitted.push_back( pull(diffMomFitted[2], momCovFitted(2, 2), "q/p")); + innerTrkWeight.push_back(trk.trackWeight); + const auto& recoUnitDirFitted = paramsAtVtxFitted->direction(); double overlapFitted = trueUnitDir.dot(recoUnitDirFitted); innerMomOverlapFitted.push_back(overlapFitted); @@ -752,6 +758,8 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT( m_pullQOverP.clear(); m_pullQOverPFitted.clear(); + m_trkWeight.clear(); + m_nTracksOnTruthVertex.clear(); m_nTracksOnRecoVertex.clear(); diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp index b8b639fd71a..c1386b809b4 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp @@ -184,6 +184,10 @@ class VertexPerformanceWriter final std::vector> m_pullQOverP; std::vector> m_pullQOverPFitted; + // Track weights from vertex fit, will be set to 1 if we do unweighted vertex + // fitting + std::vector> m_trkWeight; + // Number of tracks associated with truth/reconstructed vertex std::vector m_nTracksOnTruthVertex; std::vector m_nTracksOnRecoVertex; diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index f8a5788d0a0..661f63ea46c 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -885,6 +885,7 @@ def addFTFSeeding( inputSourceLinks="sourcelinks", trackingGeometry=trackingGeometry, fill_module_csv=False, + inputClusters="clusters", ) sequence.addAlgorithm(seedingAlg) diff --git a/Examples/Python/src/DetectorInspectorsStub.cpp b/Examples/Python/src/DetectorInspectorsStub.cpp deleted file mode 100644 index 59990c1ccb7..00000000000 --- a/Examples/Python/src/DetectorInspectorsStub.cpp +++ /dev/null @@ -1,17 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2021 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/. - -namespace Acts { -namespace Python { -struct Context; -} // namespace Python -} // namespace Acts - -namespace Acts::Python { -void addDetectorInspectors(Context& /*ctx*/) {} -} // namespace Acts::Python diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index 43968c37fa2..3881a6b7e3e 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -276,7 +276,7 @@ void addTrackFinding(Context& ctx) { ActsExamples::SeedingFTFAlgorithm, mex, "SeedingFTFAlgorithm", inputSpacePoints, outputSeeds, seedFilterConfig, seedFinderConfig, seedFinderOptions, layerMappingFile, geometrySelection, inputSourceLinks, - trackingGeometry, ACTS_FTF_Map, fill_module_csv); + trackingGeometry, ACTS_FTF_Map, fill_module_csv, inputClusters); ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::HoughTransformSeeder, mex, "HoughTransformSeeder", diff --git a/Examples/Run/Geant4/barrel_chambers00.png b/Examples/Run/Geant4/barrel_chambers00.png deleted file mode 100644 index f694d061537..00000000000 Binary files a/Examples/Run/Geant4/barrel_chambers00.png and /dev/null differ diff --git a/Examples/Scripts/Python/full_chain_itk_FTF.py b/Examples/Scripts/Python/full_chain_itk_FTF.py index 83f98a666a6..b4de671f5c7 100755 --- a/Examples/Scripts/Python/full_chain_itk_FTF.py +++ b/Examples/Scripts/Python/full_chain_itk_FTF.py @@ -99,4 +99,17 @@ outputDirRoot=outputDir, ) +addCKFTracks( + s, + trackingGeometry, + field, + TrackSelectorConfig( + pt=(1.0 * u.GeV if ttbar_pu200 else 0.0, None), + absEta=(None, 4.0), + nMeasurementsMin=6, + ), + outputDirRoot=outputDir, +) + + s.run() diff --git a/Plugins/DD4hep/src/DD4hepMaterialConverter.cpp b/Plugins/DD4hep/src/DD4hepMaterialConverter.cpp deleted file mode 100644 index b872084f353..00000000000 --- a/Plugins/DD4hep/src/DD4hepMaterialConverter.cpp +++ /dev/null @@ -1,131 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "Acts/Geometry/ApproachDescriptor.hpp" -#include "Acts/Geometry/Layer.hpp" -#include "Acts/Material/ProtoSurfaceMaterial.hpp" -#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp" -#include "Acts/Plugins/DD4hep/DD4hepConversionHelpers.hpp" -#include "Acts/Surfaces/Surface.hpp" -#include "Acts/Utilities/BinUtility.hpp" - -#include -#include -#include -#include -#include - -#include -#include - -std::shared_ptr Acts::createProtoMaterial( - const dd4hep::rec::VariantParameters& params, const std::string& valueTag, - const std::vector >& - binning, - const Logger& logger) { - using namespace std::string_literals; - - // Create the bin utility - Acts::BinUtility bu; - // Loop over the bins - for (auto& bin : binning) { - // finding the iterator position to determine the binning value - auto bit = std::find(Acts::binningValueNames().begin(), - Acts::binningValueNames().end(), bin.first); - std::size_t indx = std::distance(Acts::binningValueNames().begin(), bit); - Acts::BinningValue bval = Acts::BinningValue(indx); - Acts::BinningOption bopt = bin.second; - double min = 0.; - double max = 0.; - if (bopt == Acts::closed) { - min = -M_PI; - max = M_PI; - } - int bins = params.get(valueTag + "_"s + bin.first); - ACTS_VERBOSE(" - material binning for " << bin.first << " on " << valueTag - << ": " << bins); - if (bins >= 1) { - bu += Acts::BinUtility(bins, min, max, bopt, bval); - } - } - return std::make_shared(bu); -} - -void Acts::addLayerProtoMaterial( - const dd4hep::rec::VariantParameters& params, Layer& layer, - const std::vector >& - binning, - const Logger& logger) { - ACTS_VERBOSE("addLayerProtoMaterial"); - // Start with the representing surface - std::vector materialOptions = {"layer_material_representing"}; - std::vector materialSurfaces = { - &(layer.surfaceRepresentation())}; - // Now fill (optionally) with the approach surfaces - auto aDescriptor = layer.approachDescriptor(); - if (aDescriptor != nullptr and aDescriptor->containedSurfaces().size() >= 2) { - // Add the inner and outer approach surface - const std::vector& aSurfaces = - aDescriptor->containedSurfaces(); - materialOptions.push_back("layer_material_inner"); - materialSurfaces.push_back(aSurfaces[0]); - materialOptions.push_back("layer_material_outer"); - materialSurfaces.push_back(aSurfaces[1]); - } - - // Now loop over it and create the ProtoMaterial - for (unsigned int is = 0; is < materialOptions.size(); ++is) { - // if (actsExtension.hasValue(materialOptions[is])) { - ACTS_VERBOSE(" - checking material for: " << materialOptions[is]); - if (params.contains(materialOptions[is])) { - ACTS_VERBOSE(" - have material"); - // Create the material and assign it - auto psMaterial = - createProtoMaterial(params, materialOptions[is], binning, logger); - // const_cast (ugly - to be changed after internal geometry stored - // non-const) - Surface* surface = const_cast(materialSurfaces[is]); - surface->assignSurfaceMaterial(psMaterial); - } - } -} - -void Acts::addCylinderLayerProtoMaterial(dd4hep::DetElement detElement, - Layer& cylinderLayer, - const Logger& logger) { - ACTS_VERBOSE( - "Translating DD4hep material into Acts material for CylinderLayer : " - << detElement.name()); - if (hasParams(detElement)) { - ACTS_VERBOSE(" params: " << getParams(detElement)); - } else { - ACTS_VERBOSE(" NO params"); - } - if (getParamOr("layer_material", detElement, false)) { - addLayerProtoMaterial(getParams(detElement), cylinderLayer, - {{"binPhi", Acts::closed}, {"binZ", Acts::open}}, - logger); - } -} - -void Acts::addDiscLayerProtoMaterial(dd4hep::DetElement detElement, - Layer& discLayer, const Logger& logger) { - ACTS_VERBOSE("Translating DD4hep material into Acts material for DiscLayer : " - << detElement.name()); - - if (hasParams(detElement)) { - ACTS_VERBOSE(" params: " << getParams(detElement)); - } else { - ACTS_VERBOSE(" NO params"); - } - if (getParamOr("layer_material", detElement, false)) { - addLayerProtoMaterial(getParams(detElement), discLayer, - {{"binPhi", Acts::closed}, {"binR", Acts::open}}, - logger); - } -} diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXTiming.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXTiming.hpp deleted file mode 100644 index 8bbd9839098..00000000000 --- a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXTiming.hpp +++ /dev/null @@ -1,87 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2022 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/. - -#pragma once - -#include "Acts/Utilities/Logger.hpp" - -#include -#include -#include -#include -#include - -namespace Acts { - -/// @struct ExaTrkXTime -/// -/// @brief Collection of timing information of the Exa.TrkX algorithm -struct ExaTrkXTime { - float embedding = 0.0; - float building = 0.0; - float filtering = 0.0; - float gnn = 0.0; - float labeling = 0.0; - float total = 0.0; - - void summary(const Logger& logger) const { - ACTS_VERBOSE("1) embedding: " << embedding); - ACTS_VERBOSE("2) building: " << building); - ACTS_VERBOSE("3) filtering: " << filtering); - ACTS_VERBOSE("4) gnn: " << gnn); - ACTS_VERBOSE("5) labeling: " << labeling); - ACTS_VERBOSE("6) total: " << total); - } -}; - -/// @class ExaTrkXTimer -/// -/// A timer to allow easy timing -class ExaTrkXTimer { - public: - ExaTrkXTimer(bool disabled = false) : m_disabled(disabled) {} - - void start() { - if (!m_disabled) { - m_start = std::chrono::high_resolution_clock::now(); - m_running = true; - } - } - void stop() { - if (!m_disabled) { - m_end = std::chrono::high_resolution_clock::now(); - m_running = false; - } - } - double stopAndGetElapsedTime() { - stop(); - return elapsedSeconds(); - } - double elapsed() { - if (!m_disabled) { - std::chrono::time_point end; - if (m_running) { - end = std::chrono::high_resolution_clock::now(); - } else { - end = m_end; - } - return std::chrono::duration(end - m_start).count(); - } else { - return 0.0; - } - } - double elapsedSeconds() { return elapsed() / 1000.0; } - - private: - std::chrono::time_point m_start; - std::chrono::time_point m_end; - bool m_running = false; - bool m_disabled; -}; - -} // namespace Acts diff --git a/Plugins/Json/include/Acts/Plugins/Json/DetectorVolumeFinderJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/DetectorVolumeFinderJsonConverter.hpp index 1bc51a82696..1988e6ba536 100644 --- a/Plugins/Json/include/Acts/Plugins/Json/DetectorVolumeFinderJsonConverter.hpp +++ b/Plugins/Json/include/Acts/Plugins/Json/DetectorVolumeFinderJsonConverter.hpp @@ -43,7 +43,7 @@ void convert(nlohmann::json& jIndexedVolumes, auto castedDelegate = dynamic_cast(instance); if (castedDelegate != nullptr) { jIndexedVolumes = IndexedGridJsonHelper::convertImpl( - *castedDelegate, detray); + *castedDelegate, detray, false); if (detray) { jIndexedVolumes["volume_link"] = 1; nlohmann::json jAccLink; diff --git a/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp index 257d5d6f6b9..c9a50857e14 100644 --- a/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp +++ b/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp @@ -84,10 +84,14 @@ nlohmann::json toJson(const grid_type& grid) { /// /// @tparam grid_type the type of the grid /// @param grid the grid object +/// @param swapAxis - this is needed for detray +/// +/// @note detray has a different offset for the +/// local indices, it starts at 0 /// /// @return a json object to represent the grid template -nlohmann::json toJsonDetray(const grid_type& grid) { +nlohmann::json toJsonDetray(const grid_type& grid, bool swapAxis = false) { nlohmann::json jGrid; auto axes = grid.axes(); @@ -106,7 +110,7 @@ nlohmann::json toJsonDetray(const grid_type& grid) { if constexpr (grid_type::DIM == 1u) { for (unsigned int ib0 = 1u; ib0 <= axes[0u]->getNBins(); ++ib0) { typename grid_type::index_t lbin; - lbin[0u] = ib0; + lbin[0u] = ib0 - 1u; nlohmann::json jBin; jBin["loc_index"] = lbin; jBin["content"] = grid.atLocalBins(lbin); @@ -116,11 +120,13 @@ nlohmann::json toJsonDetray(const grid_type& grid) { // 2D connections if constexpr (grid_type::DIM == 2u) { - for (unsigned int ib0 = 1u; ib0 <= axes[0u]->getNBins(); ++ib0) { - for (unsigned int ib1 = 1u; ib1 <= axes[1u]->getNBins(); ++ib1) { + unsigned int iaxis0 = swapAxis ? 1u : 0u; + unsigned int iaxis1 = swapAxis ? 0u : 1u; + for (unsigned int ib0 = 1u; ib0 <= axes[iaxis0]->getNBins(); ++ib0) { + for (unsigned int ib1 = 1u; ib1 <= axes[iaxis1]->getNBins(); ++ib1) { typename grid_type::index_t lbin; - lbin[0u] = ib0; - lbin[1u] = ib1; + lbin[0u] = ib0 - 1u; + lbin[1u] = ib1 - 1u; nlohmann::json jBin; jBin["loc_index"] = lbin; jBin["content"] = grid.atLocalBins(lbin); diff --git a/Plugins/Json/include/Acts/Plugins/Json/IndexedGridJsonHelper.hpp b/Plugins/Json/include/Acts/Plugins/Json/IndexedGridJsonHelper.hpp index 5a396341e6f..6664b40e56f 100644 --- a/Plugins/Json/include/Acts/Plugins/Json/IndexedGridJsonHelper.hpp +++ b/Plugins/Json/include/Acts/Plugins/Json/IndexedGridJsonHelper.hpp @@ -22,10 +22,18 @@ using namespace Experimental::detail::GridAxisGenerators; namespace IndexedGridJsonHelper { /// @brief The actual conversion method +/// +/// @param indexGrid is the index grid to be written +/// @param detray is a flag indicating detray writout +/// @param checkSwap is a flag indicating if the axes should be swapped template -nlohmann::json convertImpl(const index_grid& indexGrid, bool detray = false) { +nlohmann::json convertImpl(const index_grid& indexGrid, bool detray = false, + bool checkSwap = false) { nlohmann::json jIndexedGrid; + // Axis swapping + bool swapAxes = checkSwap; + // Fill the casts nlohmann::json jCasts; // 1D casts @@ -36,12 +44,14 @@ nlohmann::json convertImpl(const index_grid& indexGrid, bool detray = false) { if constexpr (index_grid::grid_type::DIM == 2u) { jCasts.push_back(indexGrid.casts[0u]); jCasts.push_back(indexGrid.casts[1u]); + swapAxes = checkSwap && + (indexGrid.casts[0u] == binZ && indexGrid.casts[1u] == binPhi); } jIndexedGrid["casts"] = jCasts; jIndexedGrid["transform"] = Transform3JsonConverter::toJson(indexGrid.transform); if (detray) { - jIndexedGrid = GridJsonConverter::toJsonDetray(indexGrid.grid); + jIndexedGrid = GridJsonConverter::toJsonDetray(indexGrid.grid, swapAxes); } else { jIndexedGrid["grid"] = GridJsonConverter::toJson(indexGrid.grid); } diff --git a/Plugins/Json/src/DetectorJsonConverter.cpp b/Plugins/Json/src/DetectorJsonConverter.cpp index 816723a1b0e..54e96b1d06c 100644 --- a/Plugins/Json/src/DetectorJsonConverter.cpp +++ b/Plugins/Json/src/DetectorJsonConverter.cpp @@ -149,6 +149,23 @@ nlohmann::json Acts::DetectorJsonConverter::toJsonDetray( if (jSurfacesDelegate.is_null()) { continue; } + // Patch axes for cylindrical grid surfaces, axes are swapped + // at this point + auto jAccLink = jSurfacesDelegate["acc_link"]; + std::size_t accLinkType = jAccLink["type"]; + if (accLinkType == 4u) { + // Radial value to transfer phi to rphi + std::vector bValues = volume->volumeBounds().values(); + ActsScalar rRef = 0.5 * (bValues[1] + bValues[0]); + // Get the axes + auto& jAxes = jSurfacesDelegate["axes"]; + // r*phi axis is the first one + std::vector jAxesEdges = jAxes[0u]["edges"]; + std::for_each(jAxesEdges.begin(), jAxesEdges.end(), + [rRef](ActsScalar& phi) { phi *= rRef; }); + // Write back the patches axis edges + jSurfacesDelegate["axes"][0u]["edges"] = jAxesEdges; + } // Colplete the grid json for detray usage jSurfacesDelegate["volume_link"] = iv; // jSurfacesDelegate["acc_link"] = diff --git a/Plugins/Json/src/PortalJsonConverter.cpp b/Plugins/Json/src/PortalJsonConverter.cpp index 7051ae8e7e1..ae06f2ef26f 100644 --- a/Plugins/Json/src/PortalJsonConverter.cpp +++ b/Plugins/Json/src/PortalJsonConverter.cpp @@ -22,6 +22,7 @@ #include "Acts/Surfaces/RegularSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Enumerate.hpp" +#include "Acts/Utilities/VectorHelpers.hpp" #include #include @@ -134,11 +135,11 @@ std::vector Acts::PortalJsonConverter::toJsonDetray( const auto& cast = multiLink1D->indexedUpdater.casts[0u]; const auto& transform = multiLink1D->indexedUpdater.transform; const auto& volumes = multiLink1D->indexedUpdater.extractor.dVolumes; - if (!transform.isApprox(Transform3::Identity())) { - std::runtime_error( - "PortalJsonConverter: transformed boundary link implementation not " - "(yet) supported"); - } + + // Apply the correction from local to global boundaries + ActsScalar gCorr = VectorHelpers::cast(transform.translation(), cast); + std::for_each(boundaries.begin(), boundaries.end(), + [&gCorr](ActsScalar& b) { b -= gCorr; }); // Get the volume indices auto surfaceType = surfaceAdjusted->type(); @@ -305,7 +306,7 @@ std::shared_ptr Acts::PortalJsonConverter::fromJson( throw std::runtime_error( "PortalJsonConverter: surface is not a regular surface."); } - auto portal = Experimental::Portal::makeShared(regSurface); + auto portal = std::make_shared(regSurface); std::array normalDirs = {Direction::Backward, Direction::Forward}; diff --git a/Tests/IntegrationTests/PropagationTestsAtlasField.cpp b/Tests/IntegrationTests/PropagationTestsAtlasField.cpp deleted file mode 100644 index dbc1c532612..00000000000 --- a/Tests/IntegrationTests/PropagationTestsAtlasField.cpp +++ /dev/null @@ -1,125 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2017-2018 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 -#include - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/TrackParameters.hpp" -#include "Acts/MagneticField/BFieldMapUtils.hpp" -#include "Acts/MagneticField/InterpolatedBFieldMap.hpp" -#include "Acts/MagneticField/SharedBFieldMap.hpp" -#include "Acts/MagneticField/concept/AnyFieldLookup.hpp" -#include "Acts/Propagator/AtlasStepper.hpp" -#include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/Propagator.hpp" -#include "Acts/Surfaces/CylinderSurface.hpp" -#include "Acts/Surfaces/DiscSurface.hpp" -#include "Acts/Surfaces/PlaneSurface.hpp" -#include "Acts/Surfaces/StrawSurface.hpp" -#include "Acts/Utilities/Grid.hpp" -#include "Acts/Utilities/detail/Axis.hpp" - -#include "PropagationTestHelper.hpp" - -namespace Acts { - -/// Get the ATLAS field from : -/// http://acts.web.cern.ch/ACTS/data/AtlasField/AtlasField.tar.gz -/// to run this - -namespace IntegrationTest { - -// Create a mapper from a text file -InterpolatedBFieldMap::FieldMapper<3, 3> readFieldXYZ( - std::function binsXYZ, - std::array nBinsXYZ)> - localToGlobalBin, - std::string fieldMapFile = "Field.txt", double lengthUnit = 1., - double BFieldUnit = 1., std::size_t nPoints = 100000, - bool firstOctant = false) { - /// [1] Read in field map file - // Grid position points in x, y and z - std::vector xPos; - std::vector yPos; - std::vector zPos; - // components of magnetic field on grid points - std::vector bField; - // reserve estimated size - xPos.reserve(nPoints); - yPos.reserve(nPoints); - zPos.reserve(nPoints); - bField.reserve(nPoints); - // [1] Read in file and fill values - std::ifstream map_file(fieldMapFile.c_str(), std::ios::in); - std::string line; - double x = 0., y = 0., z = 0.; - double bx = 0., by = 0., bz = 0.; - while (std::getline(map_file, line)) { - if (line.empty() || line[0] == '%' || line[0] == '#' || - line.find_first_not_of(' ') == std::string::npos) - continue; - - std::istringstream tmp(line); - tmp >> x >> y >> z >> bx >> by >> bz; - xPos.push_back(x); - yPos.push_back(y); - zPos.push_back(z); - bField.push_back(Acts::Vector3(bx, by, bz)); - } - map_file.close(); - - return fieldMapperXYZ(localToGlobalBin, xPos, yPos, zPos, bField, lengthUnit, - BFieldUnit, firstOctant); -} - -// create a bfiel map from a mapper -std::shared_ptr atlasBField( - std::string fieldMapFile = "Field.txt") { - // Declare the mapper - Concepts ::AnyFieldLookup<> mapper; - double lengthUnit = UnitConstants::mm; - double BFieldUnit = UnitConstants::T; - // read the field x,y,z from a text file - mapper = readFieldXYZ( - [](std::array binsXYZ, - std::array nBinsXYZ) { - return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) + - binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2)); - }, - fieldMapFile, lengthUnit, BFieldUnit); - // create the config - InterpolatedBFieldMap::Config config; - config.scale = 1.; - config.mapper = std::move(mapper); - // make the interpolated field - return std::make_shared(std::move(config)); -} - -double Bz = 2_T; - -using BFieldType = InterpolatedBFieldMap; -using EigenStepperType = EigenStepper<>; -using AtlasStepperType = AtlasStepper; -using EigenPropagatorType = Propagator; -using AtlasPropagatorType = Propagator; - -auto bField = atlasBField("Field.txt"); - -EigenStepperType estepper(bField); -EigenPropagatorType epropagator(std::move(estepper)); -AtlasStepperType astepper(bField); -AtlasPropagatorType apropagator(std::move(astepper)); - -// The actual test - needs to be included to avoid -// template inside template definition through boost -#include "PropagationTestBase.hpp" - -} // namespace IntegrationTest -} // namespace Acts diff --git a/Tests/UnitTests/CMakeLists.txt b/Tests/UnitTests/CMakeLists.txt index 64c4e21f1bd..77de8ffdf6c 100644 --- a/Tests/UnitTests/CMakeLists.txt +++ b/Tests/UnitTests/CMakeLists.txt @@ -3,6 +3,9 @@ # the common libraries which are linked to every unittest can be # extended by setting the `unittest_extra_libraries` variables before # calling the macro. + +add_custom_target(unit_tests) + macro(add_unittest _name _source) # automatically prefix the target name set(_target "ActsUnitTest${_name}") @@ -25,6 +28,7 @@ macro(add_unittest _name _source) ${unittest_extra_libraries}) # register as unittest executable add_test(NAME ${_name} COMMAND ${_target}) + add_dependencies(unit_tests ${_target}) endmacro() # This function adds a non compile test. To this end it diff --git a/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp b/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp index 6c7b0ad91ce..a0c875ec6e7 100644 --- a/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp +++ b/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp @@ -131,7 +131,8 @@ BOOST_AUTO_TEST_CASE(UpdatePortal) { auto cylinderSurface = Acts::Surface::makeShared(nominal, 10., 100.); - auto cylinderPortal = Acts::Experimental::Portal::makeShared(cylinderSurface); + auto cylinderPortal = + std::make_shared(cylinderSurface); fullCylinderVolume->updatePortal(cylinderPortal, 2u); diff --git a/Tests/UnitTests/Core/Detector/PortalTests.cpp b/Tests/UnitTests/Core/Detector/PortalTests.cpp index 265ca1a6421..c2e89b7dacf 100644 --- a/Tests/UnitTests/Core/Detector/PortalTests.cpp +++ b/Tests/UnitTests/Core/Detector/PortalTests.cpp @@ -53,18 +53,6 @@ class LinkToVolumeImpl : public INavigationDelegate { } // namespace Experimental } // namespace Acts -/// Unpack to shared - simply to test the getSharedPtr mechanism -/// -/// @tparam referenced_type is the type of the referenced object -/// -/// @param rt is the referenced object -/// -/// @returns a shared pointer -template -std::shared_ptr unpackToShared(referenced_type& rt) { - return rt.getSharedPtr(); -} - using namespace Acts::Experimental; // A test context @@ -90,7 +78,7 @@ BOOST_AUTO_TEST_CASE(PortalTest) { Acts::Surface::makeShared(dTransform, rectangle); // Create a portal out of it - auto portalA = Portal::makeShared(surface); + auto portalA = std::make_shared(surface); BOOST_CHECK_EQUAL(&(portalA->surface()), surface.get()); @@ -98,9 +86,6 @@ BOOST_AUTO_TEST_CASE(PortalTest) { BOOST_CHECK_EQUAL(portalA->surface().geometryId(), Acts::GeometryIdentifier{5}); - BOOST_CHECK_EQUAL(portalA, unpackToShared(*portalA)); - BOOST_CHECK_EQUAL(portalA, unpackToShared(*portalA)); - // Create a links to volumes auto linkToAImpl = std::make_unique(volumeA); DetectorVolumeUpdater linkToA; @@ -124,7 +109,7 @@ BOOST_AUTO_TEST_CASE(PortalTest) { portalA->updateDetectorVolume(tContext, nState); BOOST_CHECK_EQUAL(nState.currentVolume, nullptr); - auto portalB = Portal::makeShared(surface); + auto portalB = std::make_shared(surface); DetectorVolumeUpdater linkToB; auto linkToBImpl = std::make_unique(volumeB); linkToB.connect<&LinkToVolumeImpl::link>(std::move(linkToBImpl)); @@ -139,8 +124,16 @@ BOOST_AUTO_TEST_CASE(PortalTest) { portalB->updateDetectorVolume(tContext, nState); BOOST_CHECK_EQUAL(nState.currentVolume, volumeB.get()); + Acts::GeometryContext gctx; + BOOST_CHECK_EQUAL(portalA->surface().center(gctx), + portalB->surface().center(gctx)); + + // Fuse with itself, nothing happens + BOOST_CHECK_EQUAL(portalA, Portal::fuse(portalA, portalA)); + // Now fuse the portals together, both links valid - portalA->fuse(portalB); + portalA = Portal::fuse(portalA, portalB); + nState.direction = Acts::Vector3(0., 0., 1.); portalA->updateDetectorVolume(tContext, nState); BOOST_CHECK_EQUAL(nState.currentVolume, volumeA.get()); @@ -148,26 +141,28 @@ BOOST_AUTO_TEST_CASE(PortalTest) { portalA->updateDetectorVolume(tContext, nState); BOOST_CHECK_EQUAL(nState.currentVolume, volumeB.get()); - // Portal A is now identical to portal B - BOOST_CHECK_EQUAL(portalA, portalB); + // Portal A retains identical position to B + BOOST_CHECK_EQUAL(portalA->surface().center(gctx), + portalB->surface().center(gctx)); // An invalid fusing setup auto linkToAIImpl = std::make_unique(volumeA); auto linkToBIImpl = std::make_unique(volumeB); - auto portalAI = Portal::makeShared(surface); + auto portalAI = std::make_shared(surface); DetectorVolumeUpdater linkToAI; linkToAI.connect<&LinkToVolumeImpl::link>(std::move(linkToAIImpl)); portalAI->assignDetectorVolumeUpdater(Acts::Direction::Positive, std::move(linkToAI), {volumeA}); - auto portalBI = Portal::makeShared(surface); + auto portalBI = std::make_shared(surface); DetectorVolumeUpdater linkToBI; linkToBI.connect<&LinkToVolumeImpl::link>(std::move(linkToBIImpl)); portalBI->assignDetectorVolumeUpdater(Acts::Direction::Positive, std::move(linkToBI), {volumeB}); - BOOST_CHECK_THROW(portalAI->fuse(portalBI), std::runtime_error); + BOOST_CHECK_THROW(Portal::fuse(portalAI, portalBI), std::runtime_error); + BOOST_CHECK_THROW(Portal::fuse(portalBI, portalAI), std::runtime_error); } BOOST_AUTO_TEST_CASE(PortalMaterialTest) { @@ -197,7 +192,7 @@ BOOST_AUTO_TEST_CASE(PortalMaterialTest) { auto surfaceA = Acts::Surface::makeShared( Acts::Transform3::Identity(), rectangle); surfaceA->assignSurfaceMaterial(materialA); - auto portalA = Acts::Experimental::Portal::makeShared(surfaceA); + auto portalA = std::make_shared(surfaceA); DetectorVolumeUpdater linkToA; auto linkToAImpl = std::make_unique(volumeA); @@ -207,7 +202,7 @@ BOOST_AUTO_TEST_CASE(PortalMaterialTest) { auto surfaceB = Acts::Surface::makeShared( Acts::Transform3::Identity(), rectangle); - auto portalB = Acts::Experimental::Portal::makeShared(surfaceB); + auto portalB = std::make_shared(surfaceB); DetectorVolumeUpdater linkToB; auto linkToBImpl = std::make_unique(volumeB); linkToB.connect<&LinkToVolumeImpl::link>(std::move(linkToBImpl)); @@ -215,13 +210,12 @@ BOOST_AUTO_TEST_CASE(PortalMaterialTest) { std::move(linkToB), {volumeB}); // Portal A fuses with B - // - has material and keeps it, portal B becomes portal A - portalA->fuse(portalB); + // - has material and keeps it + portalA = Portal::fuse(portalA, portalB); BOOST_CHECK_EQUAL(portalA->surface().surfaceMaterial(), materialA.get()); - BOOST_CHECK_EQUAL(portalA, portalB); // Remake portal B - portalB = Acts::Experimental::Portal::makeShared(surfaceB); + portalB = std::make_shared(surfaceB); DetectorVolumeUpdater linkToB2; auto linkToB2Impl = std::make_unique(volumeB); linkToB2.connect<&LinkToVolumeImpl::link>(std::move(linkToB2Impl)); @@ -229,14 +223,20 @@ BOOST_AUTO_TEST_CASE(PortalMaterialTest) { std::move(linkToB2), {volumeB}); // Portal B fuses with A - // - A has material and keeps it, portal B gets it from A, A becomes B + // - A has material, portal B gets it from A BOOST_REQUIRE_NE(portalA, portalB); - portalB->fuse(portalA); + + // This fails because A has accumulated volumes on both sides through fusing + BOOST_CHECK_THROW(Portal::fuse(portalB, portalA), std::invalid_argument); + // Remove Negative volume on A + portalA->assignDetectorVolumeUpdater(Acts::Direction::Negative, + DetectorVolumeUpdater{}, {}); + + portalB = Portal::fuse(portalB, portalA); BOOST_CHECK_EQUAL(portalB->surface().surfaceMaterial(), materialA.get()); - BOOST_CHECK_EQUAL(portalB, portalA); // Remake portal A and B, this time both with material - portalA = Acts::Experimental::Portal::makeShared(surfaceA); + portalA = std::make_shared(surfaceA); DetectorVolumeUpdater linkToA2; auto linkToA2Impl = std::make_unique(volumeA); linkToA2.connect<&LinkToVolumeImpl::link>(std::move(linkToA2Impl)); @@ -244,7 +244,7 @@ BOOST_AUTO_TEST_CASE(PortalMaterialTest) { std::move(linkToA2), {volumeA}); surfaceB->assignSurfaceMaterial(materialB); - portalB = Acts::Experimental::Portal::makeShared(surfaceB); + portalB = std::make_shared(surfaceB); DetectorVolumeUpdater linkToB3; auto linkToB3Impl = std::make_unique(volumeB); linkToB3.connect<&LinkToVolumeImpl::link>(std::move(linkToB3Impl)); @@ -252,9 +252,9 @@ BOOST_AUTO_TEST_CASE(PortalMaterialTest) { std::move(linkToB3), {volumeB}); // Portal A fuses with B - both have material, throw exception - BOOST_CHECK_THROW(portalA->fuse(portalB), std::runtime_error); + BOOST_CHECK_THROW(Portal::fuse(portalA, portalB), std::runtime_error); // Same in reverse - BOOST_CHECK_THROW(portalB->fuse(portalA), std::runtime_error); + BOOST_CHECK_THROW(Portal::fuse(portalB, portalA), std::runtime_error); } BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Core/Navigation/NavigationStateTests.cpp b/Tests/UnitTests/Core/Navigation/NavigationStateTests.cpp index fd5846ed9ed..0fce099f478 100644 --- a/Tests/UnitTests/Core/Navigation/NavigationStateTests.cpp +++ b/Tests/UnitTests/Core/Navigation/NavigationStateTests.cpp @@ -49,8 +49,8 @@ BOOST_AUTO_TEST_CASE(NavigationState) { Acts::Surface::makeShared(dTransform, rectangle); // Create a few fake portals out of it - auto portalA = Acts::Experimental::Portal::makeShared(pSurfaceA); - auto portalB = Acts::Experimental::Portal::makeShared(pSurfaceB); + auto portalA = std::make_shared(pSurfaceA); + auto portalB = std::make_shared(pSurfaceB); std::vector surfaces = {surfaceA.get(), surfaceB.get(), surfaceC.get()}; diff --git a/Tests/UnitTests/Core/Navigation/NavigationStateUpdatersTests.cpp b/Tests/UnitTests/Core/Navigation/NavigationStateUpdatersTests.cpp index 9a55a0364c8..f97244117e6 100644 --- a/Tests/UnitTests/Core/Navigation/NavigationStateUpdatersTests.cpp +++ b/Tests/UnitTests/Core/Navigation/NavigationStateUpdatersTests.cpp @@ -176,8 +176,8 @@ auto surfaceC = Acts::Surface::makeShared(); auto pSurfaceA = Acts::Surface::makeShared(); auto pSurfaceB = Acts::Surface::makeShared(); -auto portalA = Acts::Experimental::Portal::makeShared(pSurfaceA); -auto portalB = Acts::Experimental::Portal::makeShared(pSurfaceB); +auto portalA = std::make_shared(pSurfaceA); +auto portalB = std::make_shared(pSurfaceB); BOOST_AUTO_TEST_SUITE(Experimental) diff --git a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp index 96d9010efd3..9d507d5dbc5 100644 --- a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp +++ b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp @@ -508,7 +508,7 @@ BOOST_AUTO_TEST_CASE(Navigator_target_methods) { // The index should points to the begin BOOST_CHECK_EQUAL(state.navigation.navLayerIndex, 0); // Cache the beam pipe radius - double beamPipeR = perp(state.navigation.navLayer().position()); + double beamPipeR = perp(state.navigation.navLayer().first.position()); // step size has been updated CHECK_CLOSE_ABS(state.stepping.stepSize.value(), beamPipeR, s_onSurfaceTolerance); diff --git a/Tests/UnitTests/Core/Seeding/ATLASBottomBinFinder.hpp b/Tests/UnitTests/Core/Seeding/ATLASBottomBinFinder.hpp deleted file mode 100644 index a430a83801a..00000000000 --- a/Tests/UnitTests/Core/Seeding/ATLASBottomBinFinder.hpp +++ /dev/null @@ -1,38 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 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/. - -#pragma once - -#include "Acts/Seeding/IBinFinder.hpp" - -#include - -namespace Acts { - -// DEBUG: THIS REQUIRES THE BINS TO BE SET TO phi:41 z:11 - -/// @class ATLASBinFinder -/// The ATLASBinFinder is an implementation of the ATLAS bincombination behavior -/// satisfying the IBinFinder interface. Assumes the grid has 11 bins filled by -/// the same logic as ATLAS bins. -template -class ATLASBottomBinFinder : public IBinFinder { - public: - /// destructor - ~ATLASBottomBinFinder() = default; - - /// Return all bins that could contain space points that can be used with the - /// space points in the bin with the provided indices to create seeds. - /// @param phiBin phi index of bin with middle space points - /// @param zBin z index of bin with middle space points - /// @param binnedSP phi-z grid containing all bins - std::set findBins(std::size_t phiBin, std::size_t zBin, - const SpacePointGrid* binnedSP); -}; -} // namespace Acts -#include "ATLASBottomBinFinder.ipp" diff --git a/Tests/UnitTests/Core/Seeding/ATLASBottomBinFinder.ipp b/Tests/UnitTests/Core/Seeding/ATLASBottomBinFinder.ipp deleted file mode 100644 index 7b6a1d9c2aa..00000000000 --- a/Tests/UnitTests/Core/Seeding/ATLASBottomBinFinder.ipp +++ /dev/null @@ -1,46 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 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/. - -// DEBUG: THIS REQUIRES THE BINS TO BE SET TO phi:41 z:11 - -template -std::set Acts::ATLASBottomBinFinder::findBins( - std::size_t phiBin, std::size_t zBin, - const Acts::SpacePointGrid* binnedSP) { - std::set neighbourBins = - binnedSP->neighborHoodIndices({phiBin, zBin}, 1); - if (zBin == 6) { - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin - 1, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin + 1, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin, zBin - 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin - 1, zBin - 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin + 1, zBin - 1})); - return neighbourBins; - } - if (zBin > 6) { - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin - 1, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin + 1, zBin + 1})); - } else { - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin, zBin - 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin - 1, zBin - 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin + 1, zBin - 1})); - } - if (zBin == 4) { - neighbourBins.insert(binnedSP->getGlobalBinIndex({phiBin, zBin + 2})); - neighbourBins.insert(binnedSP->getGlobalBinIndex({phiBin - 1, zBin + 2})); - neighbourBins.insert(binnedSP->getGlobalBinIndex({phiBin + 1, zBin + 2})); - } - if (zBin == 8) { - neighbourBins.insert(binnedSP->getGlobalBinIndex({phiBin, zBin - 2})); - neighbourBins.insert(binnedSP->getGlobalBinIndex({phiBin - 1, zBin - 2})); - neighbourBins.insert(binnedSP->getGlobalBinIndex({phiBin + 1, zBin - 2})); - } - return neighbourBins; -} diff --git a/Tests/UnitTests/Core/Seeding/ATLASTopBinFinder.hpp b/Tests/UnitTests/Core/Seeding/ATLASTopBinFinder.hpp deleted file mode 100644 index d04a0ec5547..00000000000 --- a/Tests/UnitTests/Core/Seeding/ATLASTopBinFinder.hpp +++ /dev/null @@ -1,39 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 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/. - -#pragma once - -#include "Acts/Seeding/IBinFinder.hpp" - -#include - -namespace Acts { - -// DEBUG: THIS REQUIRES THE BINS TO BE SET TO phi:41 z:11 - -/// @class ATLASBinFinder -/// The ATLASBinFinder is an implementation of the ATLAS bincombination behavior -/// satisfying the IBinFinder interface. Assumes the grid has 11 bins filled by -/// the same logic as ATLAS bins. -template -class ATLASTopBinFinder : public IBinFinder { - public: - /// Virtual destructor - ~ATLASTopBinFinder() = default; - - /// Return all bins that could contain space points that can be used with the - /// space points in the bin with the provided indices to create seeds. - /// @param phiBin phi index of bin with middle space points - /// @param zBin z index of bin with middle space points - /// @param binnedSP phi-z grid containing all bins - std::set findBins(std::size_t phiBin, std::size_t zBin, - const SpacePointGrid* binnedSP); -}; -} // namespace Acts - -#include "ATLASTopBinFinder.ipp" diff --git a/Tests/UnitTests/Core/Seeding/ATLASTopBinFinder.ipp b/Tests/UnitTests/Core/Seeding/ATLASTopBinFinder.ipp deleted file mode 100644 index a3fdb6b47c8..00000000000 --- a/Tests/UnitTests/Core/Seeding/ATLASTopBinFinder.ipp +++ /dev/null @@ -1,31 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 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/. - -// DEBUG: THIS REQUIRES THE BINS TO BE SET TO phi:41 z:11 -template -std::set Acts::ATLASTopBinFinder::findBins( - std::size_t phiBin, std::size_t zBin, - const Acts::SpacePointGrid* binnedSP) { - std::set neighbourBins = - binnedSP->neighborHoodIndices({phiBin, zBin}, 1); - if (zBin == 6) { - return neighbourBins; - } - // <10 not necessary because grid doesn't return bins that don't exist (only - // works for 11 zbins) - if (zBin > 6) { - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin, zBin - 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin - 1, zBin - 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin + 1, zBin - 1})); - } else { - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin - 1, zBin + 1})); - neighbourBins.erase(binnedSP->getGlobalBinIndex({phiBin + 1, zBin + 1})); - } - return neighbourBins; -} diff --git a/Tests/UnitTests/Core/Seeding/TestCovarianceTool.hpp b/Tests/UnitTests/Core/Seeding/TestCovarianceTool.hpp deleted file mode 100644 index 33181404cb4..00000000000 --- a/Tests/UnitTests/Core/Seeding/TestCovarianceTool.hpp +++ /dev/null @@ -1,40 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 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/. - -#pragma once - -namespace Acts { -class CovarianceTool { - public: - /// Virtual destructor - ~CovarianceTool() = default; - - /// ICovarianceTool interface method - /// returns squared errors in z and r for the passed SpacePoint - /// @param sp is the SpacePoint from which the covariance values will be - /// retrieved - /// @param zAlign is the alignment uncertainty in z. - /// it is going to be squared and added to varianceZ. - /// @param rAlign is the alignment uncertainty in r. - /// it is going to be squared and added to varianceR. - /// @param sigma is multiplied with the combined alignment and covariance - /// errors - template - Acts::Vector2 getCovariances(const SpacePoint* sp, float zAlign = 0, - float rAlign = 0, float sigma = 1); -}; -template -inline Acts::Vector2 CovarianceTool::getCovariances(const SpacePoint* sp, - float zAlign, float rAlign, - float sigma) { - Acts::Vector2 cov; - cov[0] = ((*sp).varianceR + rAlign * rAlign) * sigma; - cov[1] = ((*sp).varianceZ + zAlign * zAlign) * sigma; - return cov; -} -} // namespace Acts diff --git a/Tests/UnitTests/Core/Surfaces/CMakeLists.txt b/Tests/UnitTests/Core/Surfaces/CMakeLists.txt index 653a42fbb84..65c8418c0bb 100644 --- a/Tests/UnitTests/Core/Surfaces/CMakeLists.txt +++ b/Tests/UnitTests/Core/Surfaces/CMakeLists.txt @@ -26,3 +26,4 @@ add_unittest(Surface SurfaceTests.cpp) add_unittest(TrapezoidBounds TrapezoidBoundsTests.cpp) add_unittest(VerticesHelper VerticesHelperTests.cpp) add_unittest(AlignmentHelper AlignmentHelperTests.cpp) +add_unittest(PolyhedronSurfacesTests PolyhedronSurfacesTests.cpp) diff --git a/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp b/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp index 7ea9b382d03..1b8a2709421 100644 --- a/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp +++ b/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2017-2018 CERN for the benefit of the Acts project +// Copyright (C) 2017-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -10,42 +10,26 @@ #include #include -// Helper -#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" - -// The class to test +#include "Acts/Definitions/Units.hpp" #include "Acts/Geometry/Extent.hpp" #include "Acts/Geometry/Polyhedron.hpp" - -// Cone surface +#include "Acts/Surfaces/AnnulusBounds.hpp" #include "Acts/Surfaces/ConeBounds.hpp" #include "Acts/Surfaces/ConeSurface.hpp" - -// Cylinder surface +#include "Acts/Surfaces/ConvexPolygonBounds.hpp" #include "Acts/Surfaces/CylinderBounds.hpp" #include "Acts/Surfaces/CylinderSurface.hpp" - -// Disc Surface -#include "Acts/Surfaces/AnnulusBounds.hpp" +#include "Acts/Surfaces/DiamondBounds.hpp" #include "Acts/Surfaces/DiscSurface.hpp" #include "Acts/Surfaces/DiscTrapezoidBounds.hpp" -#include "Acts/Surfaces/RadialBounds.hpp" - -// Plane Surface -#include "Acts/Surfaces/ConvexPolygonBounds.hpp" -#include "Acts/Surfaces/DiamondBounds.hpp" #include "Acts/Surfaces/EllipseBounds.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" #include "Acts/Surfaces/RectangleBounds.hpp" #include "Acts/Surfaces/TrapezoidBounds.hpp" +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" +#include "Acts/Utilities/Logger.hpp" -// Straw Surface -#include "Acts/Definitions/Units.hpp" -#include "Acts/Surfaces/LineBounds.hpp" -#include "Acts/Surfaces/StrawSurface.hpp" -#include "Acts/Visualization/ObjVisualization3D.hpp" - -#include #include #include #include @@ -54,475 +38,553 @@ namespace Acts { using namespace UnitLiterals; +Acts::Logging::Level logLevel = Acts::Logging::VERBOSE; + using IdentifiedPolyhedron = std::tuple; namespace Test { // Create a test context -GeometryContext tgContext = GeometryContext(); +const GeometryContext tgContext = GeometryContext(); -std::vector> testModes = { - {"", false, 72}, {"Triangulate", true, 72}, {"Extremas", false, 1}}; +const std::vector> testModes = { + {"Triangulate", 72}, {"Extrema", 1}}; -auto transform = std::make_shared(Transform3::Identity()); +const Transform3 transform = Transform3::Identity(); +const double epsAbs = 1e-12; -BOOST_AUTO_TEST_SUITE(Surfaces) +BOOST_AUTO_TEST_SUITE(PolyhedronSurfaces) /// Unit tests for Cone Surfaces BOOST_AUTO_TEST_CASE(ConeSurfacePolyhedrons) { - double hzpmin = 10_mm; - double hzpos = 35_mm; - double hzneg = -20_mm; - double alpha = 0.234; - double phiSector = 0.358; + ACTS_LOCAL_LOGGER( + Acts::getDefaultLogger("PolyhedronSurfacesTests", logLevel)); + ACTS_INFO("Test: ConeSurfacePolyhedrons"); + + const double hzPos = 35_mm; + const double hzNeg = -20_mm; + const double alpha = 0.234; + + const double rMax = hzPos * std::tan(alpha); for (const auto& mode : testModes) { - unsigned int segments = std::get(mode); - std::string modename = std::get(mode); - bool modetrg = std::get(mode); + ACTS_INFO("\tMode: " << std::get(mode)); + const unsigned int segments = std::get(mode); /// The full cone on one side - auto cone = std::make_shared(alpha, 0_mm, hzpos); - auto oneCone = Surface::makeShared(transform, cone); - auto oneConePh = oneCone->polyhedronRepresentation(tgContext, segments); - std::size_t expectedFaces = segments < 4 ? 4 : segments; - BOOST_CHECK_EQUAL(oneConePh.faces.size(), expectedFaces); - BOOST_CHECK_EQUAL(oneConePh.vertices.size(), expectedFaces + 1); - // Check the extent in space - double r = hzpos * std::tan(alpha); - auto extent = oneConePh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0_mm, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0_mm, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hzpos, 1e-6); + { + auto cone = std::make_shared(alpha, 0_mm, hzPos); + auto oneCone = Surface::makeShared(transform, cone); + auto oneConePh = oneCone->polyhedronRepresentation(tgContext, segments); + + const auto extent = oneConePh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0_mm, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0_mm, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), hzPos, epsAbs); + + const unsigned int expectedFaces = segments < 4 ? 4 : segments; + BOOST_CHECK_EQUAL(oneConePh.faces.size(), expectedFaces); + BOOST_CHECK_EQUAL(oneConePh.vertices.size(), expectedFaces + 1); + } /// The full cone on one side - auto conePiece = std::make_shared(alpha, hzpmin, hzpos); - auto oneConePiece = Surface::makeShared(transform, conePiece); - auto oneConePiecePh = - oneConePiece->polyhedronRepresentation(tgContext, segments); - expectedFaces = segments < 4 ? 4 : segments; - // Check the extent in space - double rmin = hzpmin * std::tan(alpha); - extent = oneConePiecePh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() rmin, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() hzpmin, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hzpos, 1e-6); - - // The full cone on both sides - auto coneBoth = std::make_shared(alpha, hzneg, hzpos); - auto twoCones = Surface::makeShared(transform, coneBoth); - auto twoConesPh = twoCones->polyhedronRepresentation(tgContext, segments); - expectedFaces = segments < 4 ? 8 : 2 * segments; - BOOST_CHECK_EQUAL(twoConesPh.faces.size(), expectedFaces); - BOOST_CHECK_EQUAL(twoConesPh.vertices.size(), expectedFaces + 1); - extent = twoConesPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0_mm, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() hzneg, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hzpos, 1e-6); - - // A centered sectoral cone on both sides - auto sectoralBoth = - std::make_shared(alpha, hzneg, hzpos, phiSector, 0.); - auto sectoralCones = - Surface::makeShared(transform, sectoralBoth); - auto sectoralConesPh = - sectoralCones->polyhedronRepresentation(tgContext, segments); - extent = sectoralConesPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0_mm, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() hzneg, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hzpos, 1e-6); + { + const double hzpMin = 10_mm; + const double rMin = hzpMin * std::tan(alpha); + + auto conePiece = std::make_shared(alpha, hzpMin, hzPos); + auto oneConePiece = + Surface::makeShared(transform, conePiece); + auto oneConePiecePh = + oneConePiece->polyhedronRepresentation(tgContext, segments); + + const auto extent = oneConePiecePh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), rMin, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), hzpMin, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), hzPos, epsAbs); + } + + /// The full cone on both sides + { + auto coneBoth = std::make_shared(alpha, hzNeg, hzPos); + auto twoCones = Surface::makeShared(transform, coneBoth); + auto twoConesPh = twoCones->polyhedronRepresentation(tgContext, segments); + + const auto extent = twoConesPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0_mm, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), hzNeg, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), hzPos, epsAbs); + + const unsigned int expectedFaces = segments < 4 ? 8 : 2 * segments; + BOOST_CHECK_EQUAL(twoConesPh.faces.size(), expectedFaces); + BOOST_CHECK_EQUAL(twoConesPh.vertices.size(), expectedFaces + 1); + } + + /// A centered sectoral cone on both sides + { + const double phiSector = 0.358; + + auto sectoralBoth = + std::make_shared(alpha, hzNeg, hzPos, phiSector, 0.); + auto sectoralCones = + Surface::makeShared(transform, sectoralBoth); + auto sectoralConesPh = + sectoralCones->polyhedronRepresentation(tgContext, segments); + + const auto extent = sectoralConesPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), 0, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rMax, epsAbs); + // CHECK_CLOSE_ABS(extent.range(binY).min(), ???, epsAbs); + // CHECK_CLOSE_ABS(extent.range(binY).max(), ???, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0_mm, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), rMax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), hzNeg, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), hzPos, epsAbs); + } } } /// Unit tests for Cylinder Surfaces BOOST_AUTO_TEST_CASE(CylinderSurfacePolyhedrons) { - double r = 25_mm; - double hZ = 35_mm; + ACTS_LOCAL_LOGGER( + Acts::getDefaultLogger("PolyhedronSurfacesTests", logLevel)); + ACTS_INFO("Test: CylinderSurfacePolyhedrons"); - double phiSector = 0.458; - double averagePhi = -1.345; + const double r = 25_mm; + const double hZ = 35_mm; for (const auto& mode : testModes) { - unsigned int segments = std::get(mode); - std::string modename = std::get(mode); - bool modetrg = std::get(mode); - - std::size_t expectedFaces = segments < 4 ? 4 : segments; - std::size_t expectedVertices = segments < 4 ? 8 : 2 * segments; + ACTS_INFO("\tMode: " << std::get(mode)); + const unsigned int segments = std::get(mode); /// The full cone on one side - auto cylinder = std::make_shared(r, hZ); - auto fullCylinder = - Surface::makeShared(transform, cylinder); - auto fullCylinderPh = - fullCylinder->polyhedronRepresentation(tgContext, segments); - - BOOST_CHECK_EQUAL(fullCylinderPh.faces.size(), expectedFaces); - BOOST_CHECK_EQUAL(fullCylinderPh.vertices.size(), expectedVertices); - // Check the extent in space - auto extent = fullCylinderPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() -hZ, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hZ, 1e-6); + { + auto cylinder = std::make_shared(r, hZ); + auto fullCylinder = + Surface::makeShared(transform, cylinder); + auto fullCylinderPh = + fullCylinder->polyhedronRepresentation(tgContext, segments); + + const auto extent = fullCylinderPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), -hZ, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), hZ, epsAbs); + + const unsigned int expectedFaces = segments < 4 ? 4 : segments; + const unsigned int expectedVertices = segments < 4 ? 8 : 2 * segments; + BOOST_CHECK_EQUAL(fullCylinderPh.faces.size(), expectedFaces); + BOOST_CHECK_EQUAL(fullCylinderPh.vertices.size(), expectedVertices); + } /// The full cone on one side - auto sectorCentered = std::make_shared(r, phiSector, hZ); - auto centerSectoredCylinder = - Surface::makeShared(transform, sectorCentered); - auto centerSectoredCylinderPh = - centerSectoredCylinder->polyhedronRepresentation(tgContext, segments); - - // Check the extent in space - extent = centerSectoredCylinderPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() r * std::cos(phiSector), 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -r * std::sin(phiSector), 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() r * std::sin(phiSector), 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() -hZ, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hZ, 1e-6); - - /// The full cone on one side - auto sectorShifted = - std::make_shared(r, averagePhi, phiSector, hZ); - auto shiftedSectoredCylinder = - Surface::makeShared(transform, sectorShifted); - auto shiftedSectoredCylinderPh = - shiftedSectoredCylinder->polyhedronRepresentation(tgContext, segments); - - // Check the extent in space - extent = shiftedSectoredCylinderPh.extent(); - CHECK_CLOSE_ABS((extent.range(binR).min() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() r, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() -hZ, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() hZ, 1e-6); + { + const double phiSector = 0.458; + + auto sectorCentered = std::make_shared(r, hZ, phiSector); + auto centerSectoredCylinder = + Surface::makeShared(transform, sectorCentered); + auto centerSectoredCylinderPh = + centerSectoredCylinder->polyhedronRepresentation(tgContext, segments); + + const auto extent = centerSectoredCylinderPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), r * std::cos(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -r * std::sin(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), r * std::sin(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), r, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), -hZ, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), hZ, epsAbs); + } } } /// Unit tests for Disc Surfaces BOOST_AUTO_TEST_CASE(DiscSurfacePolyhedrons) { - double innerR = 10_mm; - double outerR = 25_mm; - - double phiSector = 0.345; - double averagePhi = -1.0; + ACTS_LOCAL_LOGGER( + Acts::getDefaultLogger("PolyhedronSurfacesTests", logLevel)); + ACTS_INFO("Test: DiscSurfacePolyhedrons"); - double cphi = std::cos(phiSector); - double sphi = std::sin(phiSector); - - std::pair lineA = { - Vector3(0., 0., 0.), Vector3(outerR * cphi, outerR * sphi, 0.)}; - std::pair lineB = { - Vector3(0., 0., 0.), Vector3(outerR * cphi, -outerR * sphi, 0.)}; - - double minPhi = averagePhi - phiSector; - double maxPhi = averagePhi + phiSector; - lineA = {Vector3(0., 0., 0.), - Vector3(outerR * std::cos(minPhi), outerR * std::sin(minPhi), 0.)}; - lineB = {Vector3(0., 0., 0.), - Vector3(outerR * std::cos(maxPhi), outerR * std::sin(maxPhi), 0.)}; + const double innerR = 10_mm; + const double outerR = 25_mm; + const double phiSector = 0.345; for (const auto& mode : testModes) { - unsigned int segments = std::get(mode); - std::string modename = std::get(mode); - bool modetrg = std::get(mode); - - // Full disc - auto disc = std::make_shared(0_mm, outerR); - auto fullDisc = Surface::makeShared(transform, disc); - auto fullDiscPh = fullDisc->polyhedronRepresentation(tgContext, segments); - - unsigned int expectedVertices = segments > 4 ? segments : 4; - unsigned int expectedFaces = 1; - - BOOST_CHECK_EQUAL(fullDiscPh.faces.size(), expectedFaces); - BOOST_CHECK_EQUAL(fullDiscPh.vertices.size(), expectedVertices); - - auto extent = fullDiscPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - // Ring disc - auto radial = std::make_shared(innerR, outerR); - auto radialDisc = Surface::makeShared(transform, radial); - auto radialPh = radialDisc->polyhedronRepresentation(tgContext, segments); - extent = radialPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() innerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - // Sectoral disc - around 0. - auto sector = std::make_shared(0., outerR, phiSector); - auto sectorDisc = Surface::makeShared(transform, sector); - auto sectorPh = sectorDisc->polyhedronRepresentation(tgContext, segments); - extent = sectorPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -outerR * std::sin(phiSector), - 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() outerR * std::sin(phiSector), - 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - // Sectoral ring - around 0. - auto sectorRing = std::make_shared(innerR, outerR, phiSector); - auto sectorRingDisc = - Surface::makeShared(transform, sectorRing); - auto sectorRingDiscPh = - sectorRingDisc->polyhedronRepresentation(tgContext, segments); - extent = sectorRingDiscPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() innerR * std::cos(phiSector), - 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -outerR * std::sin(phiSector), - 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() outerR * std::sin(phiSector), - 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() innerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - // Sectoral disc - shifted - auto sectorRingShifted = - std::make_shared(innerR, outerR, averagePhi, phiSector); - auto sectorRingDiscShifted = - Surface::makeShared(transform, sectorRingShifted); - auto sectorRingDiscShiftedPh = - sectorRingDiscShifted->polyhedronRepresentation(tgContext, segments); - extent = sectorRingDiscShiftedPh.extent(); - CHECK_CLOSE_ABS((extent.range(binR).min() innerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - // Trapezoid for a disc - double halfXmin = 10_mm; - double halfXmax = 20_mm; - auto trapezoidDisc = std::make_shared( - halfXmin, halfXmax, innerR, outerR, 0.); - auto trapezoidDiscSf = - Surface::makeShared(transform, trapezoidDisc); - auto trapezoidDiscSfPh = - trapezoidDiscSf->polyhedronRepresentation(tgContext, segments); - extent = trapezoidDiscSfPh.extent(); - CHECK_CLOSE_ABS((extent.range(binR).min() innerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - auto trapezoidDiscShifted = std::make_shared( - halfXmin, halfXmax, innerR, outerR, averagePhi); - auto trapezoidDiscShiftedSf = - Surface::makeShared(transform, trapezoidDiscShifted); - auto trapezoidDiscShiftedSfPh = - trapezoidDiscShiftedSf->polyhedronRepresentation(tgContext, segments); - extent = trapezoidDiscShiftedSfPh.extent(); - CHECK_CLOSE_ABS((extent.range(binR).min() innerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() outerR, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - double minRadius = 7.; - double maxRadius = 12.; - double minPhiA = 0.75; - double maxPhiA = 1.4; - - Vector2 offset(-2., 2.); - - auto annulus = std::make_shared(minRadius, maxRadius, - minPhiA, maxPhiA, offset); - auto annulusDisc = Surface::makeShared(transform, annulus); - auto annulusDiscPh = - annulusDisc->polyhedronRepresentation(tgContext, segments); + ACTS_INFO("\tMode: " << std::get(mode)); + const unsigned int segments = std::get(mode); + + /// Full disc + { + auto disc = std::make_shared(0_mm, outerR); + auto fullDisc = Surface::makeShared(transform, disc); + auto fullDiscPh = fullDisc->polyhedronRepresentation(tgContext, segments); + + const auto extent = fullDiscPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + + const unsigned int expectedFaces = 1; + const unsigned int expectedVertices = segments > 4 ? segments + 1 : 4 + 1; + BOOST_CHECK_EQUAL(fullDiscPh.faces.size(), expectedFaces); + BOOST_CHECK_EQUAL(fullDiscPh.vertices.size(), expectedVertices); + } + + /// Ring disc + { + auto radial = std::make_shared(innerR, outerR); + auto radialDisc = Surface::makeShared(transform, radial); + auto radialPh = radialDisc->polyhedronRepresentation(tgContext, segments); + + const auto extent = radialPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), innerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } + + /// Sectoral disc - around 0. + { + auto sector = std::make_shared(0., outerR, phiSector); + auto sectorDisc = Surface::makeShared(transform, sector); + auto sectorPh = sectorDisc->polyhedronRepresentation(tgContext, segments); + + const auto extent = sectorPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -outerR * std::sin(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), outerR * std::sin(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } + + /// Sectoral ring - around 0. + { + auto sectorRing = + std::make_shared(innerR, outerR, phiSector); + auto sectorRingDisc = + Surface::makeShared(transform, sectorRing); + auto sectorRingDiscPh = + sectorRingDisc->polyhedronRepresentation(tgContext, segments); + + const auto extent = sectorRingDiscPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), innerR * std::cos(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -outerR * std::sin(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), outerR * std::sin(phiSector), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), innerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), outerR, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } + + /// Trapezoid for a disc + { + const double halfXmin = 10_mm; + const double halfXmax = 20_mm; + + auto trapezoidDisc = std::make_shared( + halfXmin, halfXmax, innerR, outerR, 0.); + auto trapezoidDiscSf = + Surface::makeShared(transform, trapezoidDisc); + auto trapezoidDiscSfPh = + trapezoidDiscSf->polyhedronRepresentation(tgContext, segments); + const auto extent = trapezoidDiscSfPh.extent(); + + CHECK_CLOSE_ABS(extent.range(binX).min(), -std::abs(outerR - innerR) / 2., + epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), std::abs(outerR - innerR) / 2., + epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -halfXmax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), halfXmax, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), + std::hypot(std::abs(outerR - innerR) / 2., halfXmax), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } + + /// AnnulusBounds for a disc + { + const double minRadius = 7.; + const double maxRadius = 12.; + const double minPhiA = 0.75; + const double maxPhiA = 1.4; + const Vector2 offset(0., 0.); + + auto annulus = std::make_shared(minRadius, maxRadius, + minPhiA, maxPhiA, offset); + auto annulusDisc = Surface::makeShared(transform, annulus); + auto annulusDiscPh = + annulusDisc->polyhedronRepresentation(tgContext, segments); + + const auto extent = annulusDiscPh.extent(); + // CHECK_CLOSE_ABS(extent.range(binX).min(), ???, epsAbs); + // CHECK_CLOSE_ABS(extent.range(binX).max(), ???, epsAbs); + // CHECK_CLOSE_ABS(extent.range(binY).min(), ???, epsAbs); + // CHECK_CLOSE_ABS(extent.range(binY).max(), ???, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), minRadius, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), maxRadius, epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } } } /// Unit tests for Plane Surfaces BOOST_AUTO_TEST_CASE(PlaneSurfacePolyhedrons) { - double rhX = 10_mm; - double rhY = 25_mm; - double shiftY = 50_mm; - auto rectangular = std::make_shared(rhX, rhY); - - // Special test for shifted plane to check rMin/rMax - Vector3 shift(0., shiftY, 0.); - auto shiftedTransform = std::make_shared(Transform3::Identity()); - shiftedTransform->pretranslate(shift); - auto shiftedPlane = - Surface::makeShared(shiftedTransform, rectangular); - auto shiftedPh = shiftedPlane->polyhedronRepresentation(tgContext, 1); - auto shiftedExtent = shiftedPh.extent(); - // Let's check the extent - CHECK_CLOSE_ABS(shifted(extent.range(binX).min() -rhX, 1e-6); - CHECK_CLOSE_ABS(shifted(extent.range(binX).max() rhX, 1e-6); - CHECK_CLOSE_ABS(shifted(extent.range(binY).min() -rhY + shiftY, 1e-6); - CHECK_CLOSE_ABS(shifted(extent.range(binY).max() rhY + shiftY, 1e-6); + ACTS_LOCAL_LOGGER( + Acts::getDefaultLogger("PolyhedronSurfacesTests", logLevel)); + ACTS_INFO("Test: PlaneSurfacePolyhedrons"); for (const auto& mode : testModes) { - unsigned int segments = std::get(mode); - std::string modename = std::get(mode); - bool modetrg = std::get(mode); + ACTS_INFO("\tMode: " << std::get(mode)); + const unsigned int segments = std::get(mode); /// Rectangular Plane - auto rectangularPlane = - Surface::makeShared(transform, rectangular); - auto rectangularPh = - rectangularPlane->polyhedronRepresentation(tgContext, segments); - auto extent = rectangularPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -rhX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() rhX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -rhY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() rhY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() - std::hypot(rhX, rhY), 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - BOOST_CHECK_EQUAL(rectangularPh.vertices.size() == 4); - BOOST_CHECK_EQUAL(rectangularPh.faces.size() == 1); - std::vector expectedRect = {0, 1, 2, 3}; - BOOST_CHECK_EQUAL(rectangularPh.faces[0] == expectedRect); + { + const double rhX = 10_mm; + const double rhY = 25_mm; + + auto rectangular = std::make_shared(rhX, rhY); + auto rectangularPlane = + Surface::makeShared(transform, rectangular); + auto rectangularPh = + rectangularPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = rectangularPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rhX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rhX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rhY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rhY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), std::hypot(rhX, rhY), epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + + BOOST_CHECK_EQUAL(rectangularPh.vertices.size(), 4); + BOOST_CHECK_EQUAL(rectangularPh.faces.size(), 1); + + const std::vector expectedRect = {0, 1, 2, 3}; + BOOST_CHECK(rectangularPh.faces[0] == expectedRect); + } /// Trapezoidal Plane - double thX1 = 10_mm; - double thX2 = 25_mm; - double thY = 35_mm; - - auto trapezoid = std::make_shared(thX1, thX2, thY); - auto trapezoidalPlane = - Surface::makeShared(transform, trapezoid); - auto trapezoidalPh = - trapezoidalPlane->polyhedronRepresentation(tgContext, segments); - extent = trapezoidalPh.extent(); - - double thX = std::max(thX1, thX2); - CHECK_CLOSE_ABS((extent.range(binX).min() -thX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() thX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -thY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() thY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() - std::hypot(thX, thY), 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - BOOST_CHECK_EQUAL(trapezoidalPh.vertices.size() == 4); - BOOST_CHECK_EQUAL(trapezoidalPh.faces.size() == 1); - std::vector expectedTra = {0, 1, 2, 3}; - BOOST_CHECK_EQUAL(trapezoidalPh.faces[0] == expectedTra); - - /// Ring-like ellispoidal Plane - double rMaxX = 30_mm; - double rMaxY = 40_mm; - auto ellipse = std::make_shared(0., 0., rMaxX, rMaxY); - auto ellipsoidPlane = Surface::makeShared(transform, ellipse); - auto ellispoidPh = - ellipsoidPlane->polyhedronRepresentation(tgContext, segments); - extent = ellispoidPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -rMaxX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() rMaxX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -rMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() rMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() rMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - double rMinX = 10_mm; - double rMinY = 20_mm; - auto ellipseRing = - std::make_shared(rMinX, rMaxX, rMinY, rMaxY); - auto ellipsoidRingPlane = - Surface::makeShared(transform, ellipseRing); - auto ellispoidRingPh = - ellipsoidRingPlane->polyhedronRepresentation(tgContext, segments); - - extent = ellispoidPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -rMaxX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() rMaxX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -rMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() rMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() rMinX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() rMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); - - /// ConvextPolygonBounds test - std::vector vtxs = { - Vector2(-40_mm, -10_mm), Vector2(-10_mm, -30_mm), - Vector2(30_mm, -20_mm), Vector2(10_mm, 20_mm), - Vector2(-20_mm, 50_mm), Vector2(-30_mm, 30_mm)}; - - auto sextagon = std::make_shared>(vtxs); - auto sextagonPlane = Surface::makeShared(transform, sextagon); - auto sextagonPlanePh = - sextagonPlane->polyhedronRepresentation(tgContext, segments); + { + const double thX1 = 10_mm; + const double thX2 = 25_mm; + const double thY = 35_mm; + + auto trapezoid = std::make_shared(thX1, thX2, thY); + auto trapezoidalPlane = + Surface::makeShared(transform, trapezoid); + auto trapezoidalPh = + trapezoidalPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = trapezoidalPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -std::max(thX1, thX2), epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), std::max(thX1, thX2), epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -thY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), thY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), + std::hypot(std::max(thX1, thX2), thY), epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + + BOOST_CHECK_EQUAL(trapezoidalPh.vertices.size(), 4); + BOOST_CHECK_EQUAL(trapezoidalPh.faces.size(), 1); + + const std::vector expectedTra = {0, 1, 2, 3}; + BOOST_CHECK(trapezoidalPh.faces[0] == expectedTra); + } + + /// Ring-like ellipsoidal plane + { + const double rMinX = 0_mm; + const double rMinY = 0_mm; + const double rMaxX = 30_mm; + const double rMaxY = 40_mm; + auto ellipse = + std::make_shared(rMinX, rMinY, rMaxX, rMaxY); + auto ellipsoidPlane = + Surface::makeShared(transform, ellipse); + auto ellipsoidPh = + ellipsoidPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = ellipsoidPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rMaxX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rMaxX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rMaxY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rMaxY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), std::min(rMinX, rMinY), epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), std::max(rMaxX, rMaxY), epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } + + { + const double rMinX = 10_mm; + const double rMinY = 20_mm; + const double rMaxX = 30_mm; + const double rMaxY = 40_mm; + auto ellipseRing = + std::make_shared(rMinX, rMinY, rMaxX, rMaxY); + auto ellipsoidRingPlane = + Surface::makeShared(transform, ellipseRing); + auto ellipsoidRingPh = + ellipsoidRingPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = ellipsoidRingPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rMaxX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rMaxX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rMaxY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rMaxY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), std::min(rMinX, rMinY), epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), std::max(rMaxX, rMaxY), epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } + + /// ConvexPolygonBounds test + { + std::vector vtxs = { + Vector2(-40_mm, -10_mm), Vector2(-10_mm, -30_mm), + Vector2(30_mm, -20_mm), Vector2(10_mm, 20_mm), + Vector2(-20_mm, 50_mm), Vector2(-30_mm, 30_mm)}; + + auto hexagon = std::make_shared>(vtxs); + auto hexagonPlane = Surface::makeShared(transform, hexagon); + auto hexagonPlanePh = + hexagonPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = hexagonPlanePh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -40, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), 30, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -30, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), 50, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), std::sqrt(2900), epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + } /// Diamond shaped plane - double hMinX = 10_mm; - double hMedX = 20_mm; - double hMaxX = 15_mm; - double hMinY = 40_mm; - double hMaxY = 50_mm; - auto diamond = - std::make_shared(hMinX, hMedX, hMaxX, hMinY, hMaxY); - auto diamondPlane = Surface::makeShared(transform, diamond); - auto diamondPh = - diamondPlane->polyhedronRepresentation(tgContext, segments); - BOOST_CHECK_EQUAL(diamondPh.vertices.size() == 6); - BOOST_CHECK_EQUAL(diamondPh.faces.size() == 1); - extent = diamondPh.extent(); - CHECK_CLOSE_ABS((extent.range(binX).min() -hMedX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binX).max() hMedX, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).min() -hMinY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binY).max() hMaxY, 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binR).max() - std::hypot(hMaxX, hMaxY), 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).min() 0., 1e-6); - CHECK_CLOSE_ABS((extent.range(binZ).max() 0., 1e-6); + { + const double hMinX = 10_mm; + const double hMedX = 20_mm; + const double hMaxX = 15_mm; + const double hMinY = 40_mm; + const double hMaxY = 50_mm; + + auto diamond = + std::make_shared(hMinX, hMedX, hMaxX, hMinY, hMaxY); + auto diamondPlane = Surface::makeShared(transform, diamond); + auto diamondPh = + diamondPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = diamondPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -hMedX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), hMedX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -hMinY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), hMaxY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), std::hypot(hMaxX, hMaxY), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + + BOOST_CHECK_EQUAL(diamondPh.vertices.size(), 6); + BOOST_CHECK_EQUAL(diamondPh.faces.size(), 1); + } } } -BOOST_AUTO_TEST_SUITE_END() +/// Unit tests shifted plane +BOOST_AUTO_TEST_CASE(ShiftedSurfacePolyhedrons) { + ACTS_LOCAL_LOGGER( + Acts::getDefaultLogger("PolyhedronSurfacesTests", logLevel)); + ACTS_INFO("Test: ShiftedSurfacePolyhedrons"); -} // namespace Test + const double shiftY = 50_mm; + Vector3 shift(0., shiftY, 0.); + Transform3 shiftedTransform = Transform3::Identity(); + shiftedTransform.pretranslate(shift); + for (const auto& mode : testModes) { + ACTS_INFO("\tMode: " << std::get(mode)); + const unsigned int segments = std::get(mode); + + /// Rectangular Plane + { + const double rhX = 10_mm; + const double rhY = 25_mm; + + auto rectangular = std::make_shared(rhX, rhY); + auto rectangularPlane = + Surface::makeShared(shiftedTransform, rectangular); + auto rectangularPh = + rectangularPlane->polyhedronRepresentation(tgContext, segments); + + const auto extent = rectangularPh.extent(); + CHECK_CLOSE_ABS(extent.range(binX).min(), -rhX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binX).max(), rhX, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).min(), -rhY + shiftY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binY).max(), rhY + shiftY, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).min(), 25, epsAbs); + CHECK_CLOSE_ABS(extent.range(binR).max(), std::hypot(rhX, rhY + shiftY), + epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).min(), 0., epsAbs); + CHECK_CLOSE_ABS(extent.range(binZ).max(), 0., epsAbs); + + BOOST_CHECK_EQUAL(rectangularPh.vertices.size(), 4); + BOOST_CHECK_EQUAL(rectangularPh.faces.size(), 1); + + const std::vector expectedRect = {0, 1, 2, 3}; + BOOST_CHECK(rectangularPh.faces[0] == expectedRect); + } + } +} +BOOST_AUTO_TEST_SUITE_END() +} // namespace Test } // namespace Acts diff --git a/Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp b/Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp index bdb851bf6e0..320044b03a4 100644 --- a/Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp +++ b/Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp @@ -23,7 +23,11 @@ struct MockTrack { double m_loc0; double m_loc1; double m_time; - double m_nMeasurements; + std::size_t m_nMeasurements; + std::size_t m_nHoles; + std::size_t m_nOutliers; + std::size_t m_nSharedHits; + float m_chi2; bool hasReferenceSurface() const { return true; } double theta() const { return m_theta; } @@ -32,7 +36,11 @@ struct MockTrack { double loc0() const { return m_loc0; } double loc1() const { return m_loc1; } double time() const { return m_time; } - double nMeasurements() const { return m_nMeasurements; } + std::size_t nMeasurements() const { return m_nMeasurements; } + std::size_t nHoles() const { return m_nHoles; } + std::size_t nOutliers() const { return m_nOutliers; } + std::size_t nSharedHits() const { return m_nSharedHits; } + float chi2() const { return m_chi2; } }; double thetaFromEta(double eta) { @@ -55,7 +63,11 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) { baseTrack.m_loc0 = 0.5; baseTrack.m_loc1 = 0.5; baseTrack.m_time = 0.5; - baseTrack.m_nMeasurements = 0.5; + baseTrack.m_nMeasurements = 1; + baseTrack.m_nHoles = 0; + baseTrack.m_nOutliers = 0; + baseTrack.m_nSharedHits = 0; + baseTrack.m_chi2 = 0.0; { TrackSelector selector{cfgBase}; @@ -257,6 +269,62 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) { track.m_nMeasurements = {0}; BOOST_CHECK(!selector.isValidTrack(track)); } + + { + BOOST_TEST_INFO_SCOPE("nHoles max"); + auto cfg = cfgBase; + cfg.cutSets.at(0).maxHoles = {3}; + TrackSelector selector{cfg}; + MockTrack track = baseTrack; + track.m_nHoles = {2}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_nHoles = {3}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_nHoles = {4}; + BOOST_CHECK(!selector.isValidTrack(track)); + } + + { + BOOST_TEST_INFO_SCOPE("nOutliers max"); + auto cfg = cfgBase; + cfg.cutSets.at(0).maxOutliers = {3}; + TrackSelector selector{cfg}; + MockTrack track = baseTrack; + track.m_nOutliers = {2}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_nOutliers = {3}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_nOutliers = {4}; + BOOST_CHECK(!selector.isValidTrack(track)); + } + + { + BOOST_TEST_INFO_SCOPE("nSharedHits max"); + auto cfg = cfgBase; + cfg.cutSets.at(0).maxSharedHits = {3}; + TrackSelector selector{cfg}; + MockTrack track = baseTrack; + track.m_nSharedHits = {2}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_nSharedHits = {3}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_nSharedHits = {4}; + BOOST_CHECK(!selector.isValidTrack(track)); + } + + { + BOOST_TEST_INFO_SCOPE("nSharedHits max"); + auto cfg = cfgBase; + cfg.cutSets.at(0).maxChi2 = {3}; + TrackSelector selector{cfg}; + MockTrack track = baseTrack; + track.m_chi2 = {2}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_chi2 = {3}; + BOOST_CHECK(selector.isValidTrack(track)); + track.m_chi2 = {4}; + BOOST_CHECK(!selector.isValidTrack(track)); + } } BOOST_AUTO_TEST_CASE(TestSingleBinEtaCutByBinEdge) { @@ -291,7 +359,11 @@ BOOST_AUTO_TEST_CASE(TestMultiBinCuts) { baseTrack.m_loc0 = 0.5; baseTrack.m_loc1 = 0.5; baseTrack.m_time = 0.5; - baseTrack.m_nMeasurements = 0.5; + baseTrack.m_nMeasurements = 1; + baseTrack.m_nHoles = 0; + baseTrack.m_nOutliers = 0; + baseTrack.m_nSharedHits = 0; + baseTrack.m_chi2 = 0.0; using Config = TrackSelector::Config; diff --git a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp index 08a4fbe9e14..7fc71e3b32d 100644 --- a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp @@ -58,9 +58,9 @@ Acts::CurvilinearTrackParameters makeParameters( stddev[Acts::eBoundPhi] = 2_degree; stddev[Acts::eBoundTheta] = 2_degree; stddev[Acts::eBoundQOverP] = 1 / 100_GeV; - Acts::BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal(); + const Acts::BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal(); // define a track in the transverse plane along x - Acts::Vector4 mPos4(x, y, z, w); + const Acts::Vector4 mPos4(x, y, z, w); return Acts::CurvilinearTrackParameters(mPos4, phi, theta, q / p, cov, Acts::ParticleHypothesis::pion()); } @@ -74,33 +74,43 @@ static std::vector prepareSourceLinks( return result; } +/// @brief Create a simple telescope detector +/// +/// @param geoCtx +/// @param nSurfaces Number of surfaces std::shared_ptr makeToyDetector( - const GeometryContext& tgContext, const std::size_t nSurfaces = 5) { + const Acts::GeometryContext& geoCtx, const std::size_t nSurfaces = 5) { if (nSurfaces < 1) { throw std::invalid_argument("At least 1 surfaces needs to be created."); } + + // Define the dimensions of the square surfaces + const double halfSizeSurface = 1_m; + + // Rotation of the surfaces + const double rotationAngle = M_PI * 0.5; + const Vector3 xPos(cos(rotationAngle), 0., sin(rotationAngle)); + const Vector3 yPos(0., 1., 0.); + const Vector3 zPos(-sin(rotationAngle), 0., cos(rotationAngle)); + // Construct builder CuboidVolumeBuilder cvb; // Create configurations for surfaces std::vector surfaceConfig; - for (unsigned int i = 1; i <= nSurfaces; i++) { + for (unsigned int surfPos = 1; surfPos <= nSurfaces; surfPos++) { // Position of the surfaces CuboidVolumeBuilder::SurfaceConfig cfg; - cfg.position = {i * UnitConstants::m, 0, 0.}; + cfg.position = {surfPos * UnitConstants::m, 0., 0.}; // Rotation of the surfaces - double rotationAngle = M_PI * 0.5; - Vector3 xPos(cos(rotationAngle), 0., sin(rotationAngle)); - Vector3 yPos(0., 1., 0.); - Vector3 zPos(-sin(rotationAngle), 0., cos(rotationAngle)); cfg.rotation.col(0) = xPos; cfg.rotation.col(1) = yPos; cfg.rotation.col(2) = zPos; - /// Shape of the surface - // Boundaries of the surfaces - cfg.rBounds = - std::make_shared(RectangleBounds(1_m, 1_m)); + + // Boundaries of the surfaces (shape) + cfg.rBounds = std::make_shared( + RectangleBounds(halfSizeSurface, halfSizeSurface)); // Material of the surfaces MaterialSlab matProp(makeBeryllium(), 0.5_mm); @@ -132,14 +142,16 @@ std::shared_ptr makeToyDetector( // Inner Volume - Build volume configuration CuboidVolumeBuilder::VolumeConfig volumeConfig; - volumeConfig.length = {(nSurfaces + 1) * 1_m, 2_m, 2_m}; + volumeConfig.length = {(nSurfaces + 1) * 1_m, 2 * halfSizeSurface, + 2 * halfSizeSurface}; volumeConfig.position = {volumeConfig.length.x() / 2, 0., 0.}; volumeConfig.layerCfg = layerConfig; volumeConfig.name = "Test volume"; // Outer volume - Build TrackingGeometry configuration CuboidVolumeBuilder::Config config; - config.length = {(nSurfaces + 1) * 1_m, 2_m, 2_m}; + config.length = {(nSurfaces + 1) * 1_m, 2 * halfSizeSurface, + 2 * halfSizeSurface}; config.position = {volumeConfig.length.x() / 2, 0., 0.}; config.volumeCfg = {volumeConfig}; @@ -155,7 +167,7 @@ std::shared_ptr makeToyDetector( TrackingGeometryBuilder tgb(tgbCfg); std::unique_ptr detector = - tgb.trackingGeometry(tgContext); + tgb.trackingGeometry(geoCtx); return detector; } @@ -165,43 +177,53 @@ struct Detector { }; BOOST_AUTO_TEST_SUITE(Gx2fTest) +ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)) + +// Context objects +const Acts::GeometryContext geoCtx; +const Acts::MagneticFieldContext magCtx; +const Acts::CalibrationContext calCtx; + +// Measurement resolutions +const MeasurementResolution resPixel = {MeasurementType::eLoc01, + {25_um, 50_um}}; +const MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {25_um}}; +const MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {50_um}}; +const MeasurementResolutionMap resMapAllPixel = { + {Acts::GeometryIdentifier().setVolume(0), resPixel}}; // This test checks if the call to the fitter works and no errors occur in the // framework, without fitting and updating any parameters BOOST_AUTO_TEST_CASE(NoFit) { - ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); ACTS_INFO("*** Test: NoFit -- Start"); - // Context objects - Acts::GeometryContext geoCtx; - Acts::MagneticFieldContext magCtx; - Acts::CalibrationContext calCtx; std::default_random_engine rng(42); - Detector detector; + ACTS_DEBUG("Create the detector"); const std::size_t nSurfaces = 5; + Detector detector; detector.geometry = makeToyDetector(geoCtx, nSurfaces); - auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, - 80_degree, 1_GeV, 1_e); - - MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; - MeasurementResolutionMap resolutions = { - {Acts::GeometryIdentifier().setVolume(0), resPixel}}; + ACTS_DEBUG("Set the start parameters for measurement creation and fit"); + const auto parametersMeasurements = makeParameters(); + const auto startParametersFit = makeParameters( + 7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e); - // propagator + ACTS_DEBUG("Create the measurements"); using SimPropagator = Acts::Propagator; - SimPropagator simPropagator = makeStraightPropagator(detector.geometry); - auto measurements = createMeasurements( - simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng); - auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + const SimPropagator simPropagator = makeStraightPropagator(detector.geometry); + const auto measurements = + createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements, + resMapAllPixel, rng); + const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size()); + ACTS_DEBUG("Set up the fitter"); // Reuse the SimPropagator, since we will not actually use it using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter fitter(simPropagator, gx2fLogger->clone()); + const Gx2Fitter fitter(simPropagator, gx2fLogger->clone()); const Surface* rSurface = ¶metersMeasurements.referenceSurface(); @@ -219,19 +241,31 @@ BOOST_AUTO_TEST_CASE(NoFit) { Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; - // Fit the track - auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), - startParametersFit, gx2fOptions, tracks); + ACTS_DEBUG("Fit the track"); + ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); + ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); + const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), + startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); - auto& track = *res; + const auto& track = *res; BOOST_CHECK_EQUAL(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid); BOOST_CHECK(track.hasReferenceSurface()); - BOOST_CHECK_EQUAL(track.nMeasurements(), 0u); + + // Track quantities + BOOST_CHECK_EQUAL(track.chi2(), 0.); + BOOST_CHECK_EQUAL(track.nDoF(), 0u); BOOST_CHECK_EQUAL(track.nHoles(), 0u); + BOOST_CHECK_EQUAL(track.nMeasurements(), 0u); + BOOST_CHECK_EQUAL(track.nSharedHits(), 0u); + BOOST_CHECK_EQUAL(track.nOutliers(), 0u); + + // Parameters BOOST_CHECK_EQUAL(track.parameters(), startParametersFit.parameters()); BOOST_CHECK_EQUAL(track.covariance(), BoundMatrix::Identity()); + + // Convergence BOOST_CHECK_EQUAL( (track.template component< std::size_t, @@ -242,47 +276,33 @@ BOOST_AUTO_TEST_CASE(NoFit) { } BOOST_AUTO_TEST_CASE(Fit5Iterations) { - ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); ACTS_INFO("*** Test: Fit5Iterations -- Start"); - // Create a test context - GeometryContext tgContext = GeometryContext(); + std::default_random_engine rng(42); - Detector detector; + ACTS_DEBUG("Create the detector"); const std::size_t nSurfaces = 5; - detector.geometry = makeToyDetector(tgContext, nSurfaces); - - ACTS_DEBUG("Go to propagator"); - - auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, - 80_degree, 1_GeV, 1_e); - - // Context objects - Acts::GeometryContext geoCtx; - Acts::MagneticFieldContext magCtx; - // Acts::CalibrationContext calCtx; - std::default_random_engine rng(42); + Detector detector; + detector.geometry = makeToyDetector(geoCtx, nSurfaces); - MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; - MeasurementResolutionMap resolutions = { - {Acts::GeometryIdentifier().setVolume(0), resPixel}}; + ACTS_DEBUG("Set the start parameters for measurement creation and fit"); + const auto parametersMeasurements = makeParameters(); + const auto startParametersFit = makeParameters( + 7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e); - // simulation propagator + ACTS_DEBUG("Create the measurements"); using SimPropagator = Acts::Propagator; - SimPropagator simPropagator = makeStraightPropagator(detector.geometry); - auto measurements = createMeasurements( - simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng); - auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + const SimPropagator simPropagator = makeStraightPropagator(detector.geometry); + const auto measurements = + createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements, + resMapAllPixel, rng); + const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size()); BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces); - ACTS_DEBUG("Start fitting"); - ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); - ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); - + ACTS_DEBUG("Set up the fitter"); const Surface* rSurface = ¶metersMeasurements.referenceSurface(); using RecoStepper = EigenStepper<>; @@ -292,7 +312,7 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { using RecoPropagator = decltype(recoPropagator); using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); + const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); Experimental::Gx2FitterExtensions extensions; extensions.calibrator @@ -301,27 +321,35 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { extensions.surfaceAccessor .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); - MagneticFieldContext mfContext; - CalibrationContext calContext; - const Experimental::Gx2FitterOptions gx2fOptions( - tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 5, true, 0); + geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface, + false, false, FreeToBoundCorrection(false), 5, true, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; - // Fit the track - auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), - startParametersFit, gx2fOptions, tracks); + ACTS_DEBUG("Fit the track"); + ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); + ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); + const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), + startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); - auto& track = *res; + const auto& track = *res; + BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1); BOOST_CHECK(track.hasReferenceSurface()); - BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + + // Track quantities + CHECK_CLOSE_ABS(track.chi2(), 8., 2.); + BOOST_CHECK_EQUAL(track.nDoF(), 10u); BOOST_CHECK_EQUAL(track.nHoles(), 0u); + BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + BOOST_CHECK_EQUAL(track.nSharedHits(), 0u); + BOOST_CHECK_EQUAL(track.nOutliers(), 0u); + + // Parameters // We need quite coarse checks here, since on different builds // the created measurements differ in the randomness BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0); @@ -331,6 +359,8 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1); BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0); + + // Convergence BOOST_CHECK_EQUAL( (track.template component< std::size_t, @@ -341,32 +371,22 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { } BOOST_AUTO_TEST_CASE(MixedDetector) { - ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); ACTS_INFO("*** Test: MixedDetector -- Start"); - // Create a test context - GeometryContext tgContext = GeometryContext(); + std::default_random_engine rng(42); - Detector detector; + ACTS_DEBUG("Create the detector"); const std::size_t nSurfaces = 7; - detector.geometry = makeToyDetector(tgContext, nSurfaces); - - ACTS_DEBUG("Go to propagator"); - - auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, - 80_degree, 1_GeV, 1_e); + Detector detector; + detector.geometry = makeToyDetector(geoCtx, nSurfaces); - // Context objects - Acts::GeometryContext geoCtx; - Acts::MagneticFieldContext magCtx; - // Acts::CalibrationContext calCtx; - std::default_random_engine rng(42); + ACTS_DEBUG("Set the start parameters for measurement creation and fit"); + const auto parametersMeasurements = makeParameters(); + const auto startParametersFit = makeParameters( + 7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e); - MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; - MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {25_um}}; - MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {50_um}}; - MeasurementResolutionMap resolutions = { + ACTS_DEBUG("Create the measurements"); + const MeasurementResolutionMap resMap = { {Acts::GeometryIdentifier().setVolume(2).setLayer(2), resPixel}, {Acts::GeometryIdentifier().setVolume(2).setLayer(4), resStrip0}, {Acts::GeometryIdentifier().setVolume(2).setLayer(6), resStrip1}, @@ -376,21 +396,17 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { {Acts::GeometryIdentifier().setVolume(2).setLayer(14), resPixel}, }; - // simulation propagator using SimPropagator = Acts::Propagator; - SimPropagator simPropagator = makeStraightPropagator(detector.geometry); - auto measurements = createMeasurements( - simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng); - auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + const SimPropagator simPropagator = makeStraightPropagator(detector.geometry); + const auto measurements = createMeasurements( + simPropagator, geoCtx, magCtx, parametersMeasurements, resMap, rng); + const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size()); BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces); - ACTS_DEBUG("Start fitting"); - ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); - ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); - + ACTS_DEBUG("Set up the fitter"); const Surface* rSurface = ¶metersMeasurements.referenceSurface(); using RecoStepper = EigenStepper<>; @@ -400,7 +416,7 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { using RecoPropagator = decltype(recoPropagator); using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); + const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); Experimental::Gx2FitterExtensions extensions; extensions.calibrator @@ -409,27 +425,34 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { extensions.surfaceAccessor .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); - MagneticFieldContext mfContext; - CalibrationContext calContext; - const Experimental::Gx2FitterOptions gx2fOptions( - tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 5, true, 0); + geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface, + false, false, FreeToBoundCorrection(false), 5, true, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; - // Fit the track - auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), - startParametersFit, gx2fOptions, tracks); + ACTS_DEBUG("Fit the track"); + ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); + ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); + const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), + startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); - auto& track = *res; + const auto& track = *res; BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1); BOOST_CHECK(track.hasReferenceSurface()); - BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + + // Track quantities + CHECK_CLOSE_ABS(track.chi2(), 8.5, 4.); + BOOST_CHECK_EQUAL(track.nDoF(), 10u); BOOST_CHECK_EQUAL(track.nHoles(), 0u); + BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + BOOST_CHECK_EQUAL(track.nSharedHits(), 0u); + BOOST_CHECK_EQUAL(track.nOutliers(), 0u); + + // Parameters // We need quite coarse checks here, since on different builds // the created measurements differ in the randomness BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0); @@ -439,6 +462,8 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1); BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); BOOST_CHECK_CLOSE(track.covariance().determinant(), 2e-28, 1e0); + + // Convergence BOOST_CHECK_EQUAL( (track.template component< std::size_t, @@ -450,55 +475,42 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { // This test checks if we can fit QOverP, when a magnetic field is introduced BOOST_AUTO_TEST_CASE(FitWithBfield) { - ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); ACTS_INFO("*** Test: FitWithBfield -- Start"); - // Create a test context - GeometryContext tgContext = GeometryContext(); + std::default_random_engine rng(42); - Detector detector; + ACTS_DEBUG("Create the detector"); const std::size_t nSurfaces = 5; - detector.geometry = makeToyDetector(tgContext, nSurfaces); - - ACTS_DEBUG("Go to propagator"); - - auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, - 80_degree, 1_GeV, 1_e); - - // Context objects - Acts::GeometryContext geoCtx; - Acts::MagneticFieldContext magCtx; - // Acts::CalibrationContext calCtx; - std::default_random_engine rng(42); + Detector detector; + detector.geometry = makeToyDetector(geoCtx, nSurfaces); - MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; - MeasurementResolutionMap resolutions = { - {Acts::GeometryIdentifier().setVolume(0), resPixel}}; + ACTS_DEBUG("Set the start parameters for measurement creation and fit"); + const auto parametersMeasurements = makeParameters(); + const auto startParametersFit = makeParameters( + 7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e); + ACTS_DEBUG("Create the measurements"); using SimStepper = EigenStepper<>; const auto simPropagator = makeConstantFieldPropagator(detector.geometry, 0.3_T); - auto measurements = createMeasurements( - simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng); + const auto measurements = + createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements, + resMapAllPixel, rng); - auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size()); BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces); - ACTS_DEBUG("Start fitting"); - ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); - ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); - + ACTS_DEBUG("Set up the fitter"); const Surface* rSurface = ¶metersMeasurements.referenceSurface(); // Reuse the SimPropagator, since it already uses the EigenStepper<> using SimPropagator = decltype(simPropagator); using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter fitter(simPropagator, gx2fLogger->clone()); + const Gx2Fitter fitter(simPropagator, gx2fLogger->clone()); Experimental::Gx2FitterExtensions extensions; extensions.calibrator @@ -507,27 +519,34 @@ BOOST_AUTO_TEST_CASE(FitWithBfield) { extensions.surfaceAccessor .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); - MagneticFieldContext mfContext; - CalibrationContext calContext; - const Experimental::Gx2FitterOptions gx2fOptions( - tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 5, false, 0); + geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface, + false, false, FreeToBoundCorrection(false), 5, false, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; - // Fit the track - auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), - startParametersFit, gx2fOptions, tracks); + ACTS_DEBUG("Fit the track"); + ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); + ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); + const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), + startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); - auto& track = *res; + const auto& track = *res; BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1); BOOST_CHECK(track.hasReferenceSurface()); - BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + + // Track quantities + CHECK_CLOSE_ABS(track.chi2(), 7.5, 1.5); + BOOST_CHECK_EQUAL(track.nDoF(), 10u); BOOST_CHECK_EQUAL(track.nHoles(), 0u); + BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + BOOST_CHECK_EQUAL(track.nSharedHits(), 0u); + BOOST_CHECK_EQUAL(track.nOutliers(), 0u); + + // Parameters // We need quite coarse checks here, since on different builds // the created measurements differ in the randomness // TODO investigate further the reference values for eBoundPhi and det(cov) @@ -538,6 +557,8 @@ BOOST_AUTO_TEST_CASE(FitWithBfield) { BOOST_CHECK_CLOSE(track.parameters()[eBoundQOverP], 0.5, 2e-1); BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); BOOST_CHECK_CLOSE(track.covariance().determinant(), 8e-35, 4e0); + + // Convergence BOOST_CHECK_EQUAL( (track.template component< std::size_t, @@ -548,47 +569,34 @@ BOOST_AUTO_TEST_CASE(FitWithBfield) { } BOOST_AUTO_TEST_CASE(relChi2changeCutOff) { - ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); ACTS_INFO("*** Test: relChi2changeCutOff -- Start"); - // Create a test context - GeometryContext tgContext = GeometryContext(); + std::default_random_engine rng(42); - Detector detector; + ACTS_DEBUG("Create the detector"); const std::size_t nSurfaces = 5; - detector.geometry = makeToyDetector(tgContext, nSurfaces); - - ACTS_DEBUG("Go to propagator"); - - auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, - 80_degree, 1_GeV, 1_e); - - // Context objects - Acts::GeometryContext geoCtx; - Acts::MagneticFieldContext magCtx; - // Acts::CalibrationContext calCtx; - std::default_random_engine rng(42); + Detector detector; + detector.geometry = makeToyDetector(geoCtx, nSurfaces); - MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; - MeasurementResolutionMap resolutions = { - {Acts::GeometryIdentifier().setVolume(0), resPixel}}; + ACTS_DEBUG("Set the start parameters for measurement creation and fit"); + const auto parametersMeasurements = makeParameters(); + const auto startParametersFit = makeParameters( + 7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e); + ACTS_DEBUG("Create the measurements"); // simulation propagator using SimPropagator = Acts::Propagator; - SimPropagator simPropagator = makeStraightPropagator(detector.geometry); - auto measurements = createMeasurements( - simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng); - auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + const SimPropagator simPropagator = makeStraightPropagator(detector.geometry); + const auto measurements = + createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements, + resMapAllPixel, rng); + const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size()); BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces); - ACTS_DEBUG("Start fitting"); - ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); - ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); - + ACTS_DEBUG("Set up the fitter"); const Surface* rSurface = ¶metersMeasurements.referenceSurface(); using RecoStepper = EigenStepper<>; @@ -598,7 +606,7 @@ BOOST_AUTO_TEST_CASE(relChi2changeCutOff) { using RecoPropagator = decltype(recoPropagator); using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); + const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); Experimental::Gx2FitterExtensions extensions; extensions.calibrator @@ -607,27 +615,34 @@ BOOST_AUTO_TEST_CASE(relChi2changeCutOff) { extensions.surfaceAccessor .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); - MagneticFieldContext mfContext; - CalibrationContext calContext; - const Experimental::Gx2FitterOptions gx2fOptions( - tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 500, true, 1e-5); + geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface, + false, false, FreeToBoundCorrection(false), 500, true, 1e-5); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; - // Fit the track - auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), - startParametersFit, gx2fOptions, tracks); + ACTS_DEBUG("Fit the track"); + ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); + ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); + const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), + startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); - auto& track = *res; + const auto& track = *res; BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1); BOOST_CHECK(track.hasReferenceSurface()); - BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + + // Track quantities + CHECK_CLOSE_ABS(track.chi2(), 8., 2.); + BOOST_CHECK_EQUAL(track.nDoF(), 10u); BOOST_CHECK_EQUAL(track.nHoles(), 0u); + BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + BOOST_CHECK_EQUAL(track.nSharedHits(), 0u); + BOOST_CHECK_EQUAL(track.nOutliers(), 0u); + + // Parameters // We need quite coarse checks here, since on different builds // the created measurements differ in the randomness BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0); @@ -637,6 +652,8 @@ BOOST_AUTO_TEST_CASE(relChi2changeCutOff) { BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1); BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0); + + // Convergence BOOST_CHECK_LT( (track.template component< std::size_t, diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp index a6a01858ba2..7b23de2f213 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp @@ -82,12 +82,18 @@ std::uniform_real_distribution phiDist(-M_PI, M_PI); std::uniform_real_distribution thetaDist(1.0, M_PI - 1.0); // Track charge helper distribution std::uniform_real_distribution qDist(-1, 1); +// Distribution of track time (relative to vertex time). Values are unrealistic +// and only used for testing purposes. +std::uniform_real_distribution relTDist(-4_ps, 4_ps); // Track IP resolution distribution std::uniform_real_distribution resIPDist(0., 100._um); // Track angular distribution std::uniform_real_distribution resAngDist(0., 0.1); // Track q/p resolution distribution std::uniform_real_distribution resQoPDist(-0.1, 0.1); +// Track time resolution distribution. Values are unrealistic and only used for +// testing purposes. +std::uniform_real_distribution resTDist(0_ps, 8_ps); /// @brief Unit test for AdaptiveMultiVertexFitter /// @@ -314,6 +320,146 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test) { << vtxList.at(2).fullPosition()); } +/// @brief Unit test for fitting a 4D vertex position +/// +BOOST_AUTO_TEST_CASE(time_fitting) { + // Set up RNG + int mySeed = 31415; + std::mt19937 gen(mySeed); + + // Set up constant B-Field + auto bField = std::make_shared(Vector3{0.0, 0.0, 1_T}); + + // Set up EigenStepper + EigenStepper<> stepper(bField); + + // Set up propagator with void navigator + auto propagator = std::make_shared(stepper); + + VertexingOptions vertexingOptions(geoContext, + magFieldContext); + + // IP 3D Estimator + using IPEstimator = ImpactPointEstimator; + + IPEstimator::Config ip3dEstCfg(bField, propagator); + IPEstimator ip3dEst(ip3dEstCfg); + + AdaptiveMultiVertexFitter::Config fitterCfg( + ip3dEst); + + // Linearizer for BoundTrackParameters type test + Linearizer::Config ltConfig(bField, propagator); + Linearizer linearizer(ltConfig); + + // Test smoothing + fitterCfg.doSmoothing = true; + // Do time fit + fitterCfg.useTime = true; + + AdaptiveMultiVertexFitter fitter( + std::move(fitterCfg)); + + // Vertex position + double trueVtxTime = 40.0_ps; + Vector3 trueVtxPos(-0.15_mm, -0.1_mm, -1.5_mm); + + // Seed position of the vertex + Vector4 vtxSeedPos(0.0_mm, 0.0_mm, -1.4_mm, 0.0_ps); + + Vertex vtx(vtxSeedPos); + // Set initial covariance matrix to a large value + SquareMatrix4 initialCovariance(SquareMatrix4::Identity() * 1e+8); + vtx.setFullCovariance(initialCovariance); + + // Vector of all tracks + std::vector trks; + + unsigned int nTracks = 4; + for (unsigned int _ = 0; _ < nTracks; _++) { + // Construct positive or negative charge randomly + double q = qDist(gen) < 0 ? -1. : 1.; + + // Track resolution + double resD0 = resIPDist(gen); + double resZ0 = resIPDist(gen); + double resPh = resAngDist(gen); + double resTh = resAngDist(gen); + double resQp = resQoPDist(gen); + double resT = resTDist(gen); + + // Random diagonal covariance matrix + Covariance covMat; + + // clang-format off + covMat << + resD0 * resD0, 0., 0., 0., 0., 0., + 0., resZ0 * resZ0, 0., 0., 0., 0., + 0., 0., resPh * resPh, 0., 0., 0., + 0., 0., 0., resTh * resTh, 0., 0., + 0., 0., 0., 0., resQp * resQp, 0., + 0., 0., 0., 0., 0., resT * resT; + // clang-format on + + // Random track parameters + BoundTrackParameters::ParametersVector paramVec; + paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen), + q / pTDist(gen), trueVtxTime + relTDist(gen); + + std::shared_ptr perigeeSurface = + Surface::makeShared(trueVtxPos); + + trks.emplace_back(perigeeSurface, paramVec, std::move(covMat), + ParticleHypothesis::pion()); + } + + std::vector trksPtr; + for (const auto& trk : trks) { + trksPtr.push_back(&trk); + } + + // Prepare fitter state + AdaptiveMultiVertexFitter::State state( + *bField, magFieldContext); + + for (const auto& trk : trks) { + ACTS_DEBUG("Track parameters:\n" << trk); + // Index of current vertex + state.vtxInfoMap[&vtx].trackLinks.push_back(&trk); + state.tracksAtVerticesMap.insert( + std::make_pair(std::make_pair(&trk, &vtx), + TrackAtVertex(1., trk, &trk))); + } + + state.addVertexToMultiMap(vtx); + + auto res = fitter.addVtxToFit(state, vtx, linearizer, vertexingOptions); + + BOOST_CHECK(res.ok()); + + ACTS_DEBUG("Truth vertex position: " << trueVtxPos.transpose()); + ACTS_DEBUG("Fitted vertex position: " << vtx.position().transpose()); + + ACTS_DEBUG("Truth vertex time: " << trueVtxTime); + ACTS_DEBUG("Fitted vertex time: " << vtx.time()); + + // Check that true vertex position and time approximately recovered + CHECK_CLOSE_ABS(trueVtxPos, vtx.position(), 60_um); + CHECK_CLOSE_ABS(trueVtxTime, vtx.time(), 2_ps); + + const SquareMatrix4& vtxCov = vtx.fullCovariance(); + + ACTS_DEBUG("Vertex 4D covariance after the fit:\n" << vtxCov); + + // Check that variances of the vertex position/time are positive + for (std::size_t i = 0; i <= 3; i++) { + BOOST_CHECK_GT(vtxCov(i, i), 0.); + } + + // Check that the covariance matrix is approximately symmetric + CHECK_CLOSE_ABS(vtxCov - vtxCov.transpose(), SquareMatrix4::Zero(), 1e-5); +} + /// @brief Unit test for AdaptiveMultiVertexFitter /// based on Athena unit test, i.e. same setting and /// test values are used here diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp index 871899777ac..83d5fcb8935 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp @@ -178,7 +178,7 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) { Vertex vtx(vtxPos); // Update trkAtVertex with assumption of originating from vtx - KalmanVertexTrackUpdater::update(trkAtVtx, vtx); + KalmanVertexTrackUpdater::update(trkAtVtx, vtx); // The old distance double oldDistance = diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp index 0360756c319..6d2ddea556d 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp @@ -172,8 +172,8 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) { vtx.setFullCovariance(SquareMatrix4::Identity() * 0.01); // Update trkAtVertex with assumption of originating from vtx - KalmanVertexUpdater::updateVertexWithTrack(vtx, - trkAtVtx); + KalmanVertexUpdater::updateVertexWithTrack( + vtx, trkAtVtx); if (debug) { std::cout << "Old vertex position: " << vtxPos << std::endl; diff --git a/Tests/UnitTests/Plugins/Json/PortalJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/PortalJsonConverterTests.cpp index 04f592ea4e6..ae59c4bcc97 100644 --- a/Tests/UnitTests/Plugins/Json/PortalJsonConverterTests.cpp +++ b/Tests/UnitTests/Plugins/Json/PortalJsonConverterTests.cpp @@ -38,7 +38,8 @@ BOOST_AUTO_TEST_CASE(PortalUnconnected) { auto surface = Acts::Surface::makeShared( Acts::Vector3(0., 0., 0.), Acts::Vector3(0., 1., 0.)); - auto portal = Acts::Experimental::Portal::makeShared(std::move(surface)); + auto portal = + std::make_shared(std::move(surface)); BOOST_CHECK_NE(portal, nullptr); @@ -70,7 +71,8 @@ BOOST_AUTO_TEST_CASE(PortalSingleConnected) { auto surface = Acts::Surface::makeShared( Acts::Vector3(0., 0., 0.), Acts::Vector3(0., 1., 0.)); - auto portal = Acts::Experimental::Portal::makeShared(std::move(surface)); + auto portal = + std::make_shared(std::move(surface)); BOOST_CHECK_NE(portal, nullptr); // Attaching the portals Acts::Experimental::detail::PortalHelper::attachDetectorVolumeUpdater( @@ -115,7 +117,8 @@ BOOST_AUTO_TEST_CASE(PortalMultiConnected) { auto surface = Acts::Surface::makeShared( Acts::Vector3(0., 0., 0.), Acts::Vector3(0., 1., 0.)); - auto portal = Acts::Experimental::Portal::makeShared(std::move(surface)); + auto portal = + std::make_shared(std::move(surface)); BOOST_CHECK_NE(portal, nullptr); // Attaching the portals diff --git a/docs/core/figures/material/MaterialAssignment.jpeg b/docs/core/figures/material/MaterialAssignment.jpg similarity index 100% rename from docs/core/figures/material/MaterialAssignment.jpeg rename to docs/core/figures/material/MaterialAssignment.jpg diff --git a/docs/core/figures/material/MaterialAveraging.jpeg b/docs/core/figures/material/MaterialAveraging.jpg similarity index 100% rename from docs/core/figures/material/MaterialAveraging.jpeg rename to docs/core/figures/material/MaterialAveraging.jpg diff --git a/docs/core/geometry/figures/DiscLayerAB.png b/docs/core/geometry/figures/DiscLayerAB.png deleted file mode 100644 index 54676e2a735..00000000000 Binary files a/docs/core/geometry/figures/DiscLayerAB.png and /dev/null differ diff --git a/docs/core/material.md b/docs/core/material.md index cc6f76b6441..05efff08748 100644 --- a/docs/core/material.md +++ b/docs/core/material.md @@ -16,5 +16,5 @@ Link to instructions how to run the material mapping in the examples Improve material mapping technical description ::: -![](figures/material/MaterialAssignment.jpeg) -![](figures/material/MaterialAveraging.jpeg) +![](figures/material/MaterialAssignment.jpg) +![](figures/material/MaterialAveraging.jpg) diff --git a/docs/examples/howto/figures/materialMapping/ActsMaterialMappingAutoTuning.png b/docs/examples/howto/figures/materialMapping/ActsMaterialMappingAutoTuning.png deleted file mode 100644 index 9bf06a09dff..00000000000 Binary files a/docs/examples/howto/figures/materialMapping/ActsMaterialMappingAutoTuning.png and /dev/null differ diff --git a/docs/figures/doxygen/RotatedTrapezoidBounds.gif b/docs/figures/doxygen/RotatedTrapezoidBounds.gif deleted file mode 100644 index c9c6c0fbe4d..00000000000 Binary files a/docs/figures/doxygen/RotatedTrapezoidBounds.gif and /dev/null differ diff --git a/docs/figures/performance/CKF/duplicationrate_vs_eta_ttbar_pu200.png b/docs/figures/performance/CKF/duplicationrate_vs_eta_ttbar_pu200.png deleted file mode 100644 index 219281a0df9..00000000000 Binary files a/docs/figures/performance/CKF/duplicationrate_vs_eta_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/CKF/fakerate_vs_eta_ttbar_pu200.png b/docs/figures/performance/CKF/fakerate_vs_eta_ttbar_pu200.png deleted file mode 100644 index 421b3d222ff..00000000000 Binary files a/docs/figures/performance/CKF/fakerate_vs_eta_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/CKF/trackeff_vs_eta_ttbar_pu200.png b/docs/figures/performance/CKF/trackeff_vs_eta_ttbar_pu200.png deleted file mode 100644 index 86077f64e71..00000000000 Binary files a/docs/figures/performance/CKF/trackeff_vs_eta_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/fitter/nHoles_vs_eta_ttbar_pu200.png b/docs/figures/performance/fitter/nHoles_vs_eta_ttbar_pu200.png deleted file mode 100644 index 4c38da8824a..00000000000 Binary files a/docs/figures/performance/fitter/nHoles_vs_eta_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/fitter/nMeasurements_vs_eta_ttbar_pu200.png b/docs/figures/performance/fitter/nMeasurements_vs_eta_ttbar_pu200.png deleted file mode 100644 index cfa5af10d4b..00000000000 Binary files a/docs/figures/performance/fitter/nMeasurements_vs_eta_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/fitter/pull_perigee_parameters_ttbar_pu200.png b/docs/figures/performance/fitter/pull_perigee_parameters_ttbar_pu200.png deleted file mode 100644 index 4aecbd6cc2b..00000000000 Binary files a/docs/figures/performance/fitter/pull_perigee_parameters_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/fitter/trackeff_vs_eta_ttbar_pu200.png b/docs/figures/performance/fitter/trackeff_vs_eta_ttbar_pu200.png deleted file mode 100644 index 232cb1306c9..00000000000 Binary files a/docs/figures/performance/fitter/trackeff_vs_eta_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/fitter/trackeff_vs_pT_ttbar_pu200.png b/docs/figures/performance/fitter/trackeff_vs_pT_ttbar_pu200.png deleted file mode 100644 index 760a636fd01..00000000000 Binary files a/docs/figures/performance/fitter/trackeff_vs_pT_ttbar_pu200.png and /dev/null differ diff --git a/docs/figures/performance/seeding/seeding_eff_vs_eta.png b/docs/figures/performance/seeding/seeding_eff_vs_eta.png deleted file mode 100644 index 304813f03bc..00000000000 Binary files a/docs/figures/performance/seeding/seeding_eff_vs_eta.png and /dev/null differ diff --git a/docs/figures/performance/seeding/seeding_eff_vs_pt.png b/docs/figures/performance/seeding/seeding_eff_vs_pt.png deleted file mode 100644 index bfc8e29882b..00000000000 Binary files a/docs/figures/performance/seeding/seeding_eff_vs_pt.png and /dev/null differ diff --git a/docs/figures/tracking/layer_geo.svg b/docs/figures/tracking/layer_geo.svg deleted file mode 100644 index ccd1d344816..00000000000 --- a/docs/figures/tracking/layer_geo.svg +++ /dev/null @@ -1 +0,0 @@ -xyzbarrelpositive endcapnegative endcap \ No newline at end of file diff --git a/docs/figures/tracking/navigation.svg b/docs/figures/tracking/navigation.svg deleted file mode 100644 index c47e633544a..00000000000 --- a/docs/figures/tracking/navigation.svg +++ /dev/null @@ -1 +0,0 @@ -pADpDHpCDABCDEFGHxyz \ No newline at end of file diff --git a/docs/figures/tracking/strip_ghosts.svg b/docs/figures/tracking/strip_ghosts.svg deleted file mode 100644 index 2239ddea9b6..00000000000 --- a/docs/figures/tracking/strip_ghosts.svg +++ /dev/null @@ -1 +0,0 @@ -(a)particleAparticleB(b) \ No newline at end of file diff --git a/docs/plugins/figures/sycl/count_duplets.png b/docs/plugins/figures/sycl/count_duplets.png deleted file mode 100644 index 41b21818b24..00000000000 Binary files a/docs/plugins/figures/sycl/count_duplets.png and /dev/null differ diff --git a/docs/plugins/figures/sycl/duplet_search_matrix.png b/docs/plugins/figures/sycl/duplet_search_matrix.png deleted file mode 100644 index aaf5cc8e500..00000000000 Binary files a/docs/plugins/figures/sycl/duplet_search_matrix.png and /dev/null differ diff --git a/docs/plugins/figures/sycl/flat_matrix.png b/docs/plugins/figures/sycl/flat_matrix.png deleted file mode 100644 index 643776b4945..00000000000 Binary files a/docs/plugins/figures/sycl/flat_matrix.png and /dev/null differ diff --git a/docs/plugins/figures/sycl/middle_duplet_indices.png b/docs/plugins/figures/sycl/middle_duplet_indices.png deleted file mode 100644 index 20a7ace1116..00000000000 Binary files a/docs/plugins/figures/sycl/middle_duplet_indices.png and /dev/null differ diff --git a/docs/plugins/figures/sycl/prefix_sum_array.png b/docs/plugins/figures/sycl/prefix_sum_array.png deleted file mode 100644 index b6e01718b4a..00000000000 Binary files a/docs/plugins/figures/sycl/prefix_sum_array.png and /dev/null differ diff --git a/docs/plugins/figures/tgeo/TGeoArb8_PlaneSurface.png b/docs/plugins/figures/tgeo/TGeoArb8_PlaneSurface.png deleted file mode 100644 index 316543d961d..00000000000 Binary files a/docs/plugins/figures/tgeo/TGeoArb8_PlaneSurface.png and /dev/null differ diff --git a/docs/whitepaper-template.png b/docs/whitepaper-template.png deleted file mode 100644 index e8e4c8acf86..00000000000 Binary files a/docs/whitepaper-template.png and /dev/null differ