diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 9a17048cdd5..c2ece0d20f2 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,41 +1,7 @@ -PLEASE FOLLOW THE CHECKLIST BELOW WHEN CREATING A NEW PULL REQUEST. THE -CHECKLIST IS FOR YOUR INFORMATION AND MUST BE REMOVED BEFORE SUBMITTING THE PULL -REQUEST. +PLEASE DESCRIBE YOUR CHANGES. +THIS MESSAGE ENDS UP AS THE COMMIT MESSAGE. +DO NOT USE @-MENTIONS HERE! -## Checklist +--- END COMMIT MESSAGE --- -- [ ] Does the PR title follow the `: title` scheme? - - The prefix must be one of: - - - `fix`: for a bugfix - - `feat`: for a new feature - - `refactor`: for an improvement of an existing feature - - `perf`, `test`: for performance- or test-related changes - - `docs`: for documentation-related changes - - `build`, `ci`, `chore`: as appropriated for infrastructure changes - -- [ ] Does this modify the public API as defined in `docs/versioning.rst`? - - - [ ] Does the PR title contain a `!` to indicate a breaking change? - - [ ] Is there section starting with `BREAKING CHANGE:` in the PR body - that explains the breaking change? - -- [ ] Is the PR ready to be merged? - - - [ ] If not: is it marked as a draft PR? - -- [ ] Does this PR close an existing issue? - - - [ ] Is the issue correctly linked so it will be automatically closed - upon successful merge (See closing keywords link in the sidebar)? - -- The CI will initially report a missing milestone. One of the maintainers will - handle assigning a milestone for book-keeping. - -- An automated workflow will assign labels based on changed files, and whether - or not reference files were changed. These do not have to be set manually. - -- If you push updates, and you know they will be superseded later on, consider adding - `[skip ci]` in the commit message. This will instruct the CI system not to run any - jobs on this commit. +Any further description goes here, @-mentions are ok here! diff --git a/.kodiak.toml b/.kodiak.toml index b5be47df1ee..651905f7ce0 100644 --- a/.kodiak.toml +++ b/.kodiak.toml @@ -2,7 +2,7 @@ version = 1 merge.message.title = "pull_request_title" merge.message.body = "pull_request_body" -merge.message.cut_body_before = "--- BEGIN COMMIT MESSAGE ---" +merge.message.cut_body_after = "--- END COMMIT MESSAGE ---" merge.message.cut_body_and_text = true merge.method = "squash" merge.message.include_coauthors = true diff --git a/CI/physmon/workflows/physmon_simulation.py b/CI/physmon/workflows/physmon_simulation.py index 746c68c39aa..37874663187 100755 --- a/CI/physmon/workflows/physmon_simulation.py +++ b/CI/physmon/workflows/physmon_simulation.py @@ -64,6 +64,14 @@ ) ) + s.addWriter( + acts.examples.RootParticleWriter( + level=acts.logging.INFO, + inputParticles="particles_input", + filePath=tp / "particles.root", + ) + ) + addFatras( s, setup.trackingGeometry, @@ -100,6 +108,7 @@ s.run() for file, name in [ + (tp / "particles.root", "particles_gun.root"), (tp / "fatras" / "particles_simulation.root", "particles_fatras.root"), (tp / "geant4" / "particles_simulation.root", "particles_geant4.root"), ]: @@ -135,8 +144,8 @@ s.run() for file, name in [ - (tp / "pythia8_particles.root", "particles_ttbar.root"), - (tp / "pythia8_vertices.root", "vertices_ttbar.root"), + (tp / "particles.root", "particles_ttbar.root"), + (tp / "vertices.root", "vertices_ttbar.root"), ]: assert file.exists(), "file not found" shutil.copy(file, setup.outdir / name) diff --git a/CMakeLists.txt b/CMakeLists.txt index f81002cab83..f5209e0d210 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -73,7 +73,6 @@ option(ACTS_BUILD_EXAMPLES_HEPMC3 "Build HepMC3-based code in the examples" OFF) option(ACTS_BUILD_EXAMPLES_HASHING "Build Hashing-based code in the examples" OFF) option(ACTS_BUILD_EXAMPLES_PYTHIA8 "Build Pythia8-based code in the examples" OFF) option(ACTS_BUILD_EXAMPLES_PYTHON_BINDINGS "Build python bindings for the examples" OFF) -option(ACTS_USE_EXAMPLES_TBB "Use Threading Building Blocks library in the examples" ON) option(ACTS_BUILD_ANALYSIS_APPS "Build Analysis applications in the examples" OFF) # test related options option(ACTS_BUILD_BENCHMARKS "Build benchmarks" OFF) diff --git a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp index f624cd88aee..38f297764b1 100644 --- a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Definitions/Tolerance.hpp" +#include "Acts/EventData/TrackParameterHelpers.hpp" #include "Acts/EventData/TransformationHelpers.hpp" #include "Acts/EventData/detail/PrintParameters.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -18,7 +19,6 @@ #include #include #include -#include namespace Acts { @@ -61,7 +61,10 @@ class GenericBoundTrackParameters { m_cov(std::move(cov)), m_surface(std::move(surface)), m_particleHypothesis(std::move(particleHypothesis)) { - assert(m_surface); + // TODO set `validateAngleRange` to `true` after fixing caller code + assert(isBoundVectorValid(m_params, false) && + "Invalid bound parameters vector"); + assert(m_surface != nullptr && "Reference surface must not be null"); normalizePhiTheta(); } diff --git a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp index 214aa2e1551..7a2b4b46522 100644 --- a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/TrackParameterHelpers.hpp" #include "Acts/EventData/TrackParametersConcept.hpp" #include "Acts/EventData/TransformationHelpers.hpp" #include "Acts/EventData/detail/PrintParameters.hpp" @@ -18,10 +19,8 @@ #include "Acts/Utilities/UnitVectors.hpp" #include "Acts/Utilities/VectorHelpers.hpp" -#include #include #include -#include namespace Acts { @@ -55,7 +54,9 @@ class GenericFreeTrackParameters { ParticleHypothesis particleHypothesis) : m_params(params), m_cov(std::move(cov)), - m_particleHypothesis(std::move(particleHypothesis)) {} + m_particleHypothesis(std::move(particleHypothesis)) { + assert(isFreeVectorValid(m_params) && "Invalid free parameters vector"); + } /// Construct from four-position, direction, absolute momentum, and charge. /// @@ -78,6 +79,8 @@ class GenericFreeTrackParameters { m_params[eFreeDir1] = dir[eMom1]; m_params[eFreeDir2] = dir[eMom2]; m_params[eFreeQOverP] = qOverP; + + assert(isFreeVectorValid(m_params) && "Invalid free parameters vector"); } /// Construct from four-position, angles, absolute momentum, and charge. @@ -103,6 +106,8 @@ class GenericFreeTrackParameters { m_params[eFreeDir1] = dir[eMom1]; m_params[eFreeDir2] = dir[eMom2]; m_params[eFreeQOverP] = qOverP; + + assert(isFreeVectorValid(m_params) && "Invalid free parameters vector"); } /// Converts a free track parameter with a different hypothesis. diff --git a/Core/include/Acts/EventData/Seed.hpp b/Core/include/Acts/EventData/Seed.hpp index db93b319dcb..f3c6707b2e1 100644 --- a/Core/include/Acts/EventData/Seed.hpp +++ b/Core/include/Acts/EventData/Seed.hpp @@ -14,8 +14,9 @@ namespace Acts { template - requires(N >= 3ul) class Seed { + static_assert(N >= 3ul); + public: using value_type = external_spacepoint_t; static constexpr std::size_t DIM = N; diff --git a/Core/include/Acts/EventData/Seed.ipp b/Core/include/Acts/EventData/Seed.ipp index 74f2c0e4169..e977a3d67d7 100644 --- a/Core/include/Acts/EventData/Seed.ipp +++ b/Core/include/Acts/EventData/Seed.ipp @@ -9,7 +9,6 @@ namespace Acts { template - requires(N >= 3ul) template requires(sizeof...(args_t) == N) && (std::same_as && ...) @@ -17,32 +16,27 @@ Seed::Seed(const args_t&... points) : m_spacepoints({&points...}) {} template - requires(N >= 3ul) void Seed::setVertexZ(float vertex) { m_vertexZ = vertex; } template - requires(N >= 3ul) void Seed::setQuality(float seedQuality) { m_seedQuality = seedQuality; } template - requires(N >= 3ul) const std::array& Seed::sp() const { return m_spacepoints; } template - requires(N >= 3ul) float Seed::z() const { return m_vertexZ; } template - requires(N >= 3ul) float Seed::seedQuality() const { return m_seedQuality; } diff --git a/Core/include/Acts/EventData/TrackParameterHelpers.hpp b/Core/include/Acts/EventData/TrackParameterHelpers.hpp index c9089ff9ab9..df40c934824 100644 --- a/Core/include/Acts/EventData/TrackParameterHelpers.hpp +++ b/Core/include/Acts/EventData/TrackParameterHelpers.hpp @@ -11,8 +11,38 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Utilities/detail/periodic.hpp" +#include + namespace Acts { +/// Check if a bound vector is valid. This checks the following: +/// - All values are finite +/// - (optionally) The phi value is in the range [-pi, pi) +/// - (optionally) The theta value is in the range [0, pi] +/// +/// @param v The bound vector to check +/// @param validateAngleRange If true, the phi and theta values are range checked +/// @param epsilon The epsilon to use for the checks +/// @param maxAbsEta The maximum allowed eta value +/// +/// @return True if the bound vector is valid +bool isBoundVectorValid( + const BoundVector& v, bool validateAngleRange, double epsilon = 1e-6, + double maxAbsEta = std::numeric_limits::infinity()); + +/// Check if a free vector is valid. This checks the following: +/// - All values are finite +/// - Direction is normalized +/// +/// @param v The free vector to check +/// @param epsilon The epsilon to use for the checks +/// @param maxAbsEta The maximum allowed eta value +/// +/// @return True if the free vector is valid +bool isFreeVectorValid( + const FreeVector& v, double epsilon = 1e-6, + double maxAbsEta = std::numeric_limits::infinity()); + /// Normalize the bound parameter angles /// /// @param boundParams The bound parameters to normalize diff --git a/Core/include/Acts/EventData/TrackParametersConcept.hpp b/Core/include/Acts/EventData/TrackParametersConcept.hpp index b39255e5ef7..5e4982545a4 100644 --- a/Core/include/Acts/EventData/TrackParametersConcept.hpp +++ b/Core/include/Acts/EventData/TrackParametersConcept.hpp @@ -12,9 +12,7 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Geometry/GeometryContext.hpp" -#include #include -#include namespace Acts { @@ -68,4 +66,5 @@ concept BoundTrackParametersConcept = { p.position(c) } -> std::same_as; }; }; + } // namespace Acts diff --git a/Core/include/Acts/EventData/detail/GenerateParameters.hpp b/Core/include/Acts/EventData/detail/GenerateParameters.hpp index 023314e2527..64e7bc3405b 100644 --- a/Core/include/Acts/EventData/detail/GenerateParameters.hpp +++ b/Core/include/Acts/EventData/detail/GenerateParameters.hpp @@ -143,6 +143,9 @@ struct GenerateQoverPOptions { /// Indicate if the momentum referse to transverse momentum bool pTransverse = true; + /// Indicate if the momentum should be uniformly distributed in log space. + bool pLogUniform = false; + /// Charge of the parameters. double charge = 1; @@ -157,6 +160,19 @@ inline double generateQoverP(generator_t& rng, using UniformIndex = std::uniform_int_distribution; using UniformReal = std::uniform_real_distribution; + auto drawP = [&options](generator_t& rng_, double theta_) -> double { + const double pTransverseScaling = + options.pTransverse ? 1. / std::sin(theta_) : 1.; + + if (options.pLogUniform) { + UniformReal pLogDist(std::log(options.pMin), std::log(options.pMax)); + return std::exp(pLogDist(rng_)) * pTransverseScaling; + } + + UniformReal pDist(options.pMin, options.pMax); + return pDist(rng_) * pTransverseScaling; + }; + // choose between particle/anti-particle if requested // the upper limit of the distribution is inclusive UniformIndex particleTypeChoice(0u, options.randomizeCharge ? 1u : 0u); @@ -165,14 +181,12 @@ inline double generateQoverP(generator_t& rng, options.charge, -options.charge, }; - UniformReal pDist(options.pMin, options.pMax); // draw parameters const std::uint8_t type = particleTypeChoice(rng); const double q = qChoices[type]; - const double p = - pDist(rng) * (options.pTransverse ? 1. / std::sin(theta) : 1.); + const double p = drawP(rng, theta); const double qOverP = (q != 0) ? q / p : 1 / p; return qOverP; @@ -195,6 +209,10 @@ struct GenerateBoundParametersOptions { GenerateQoverPOptions qOverP; }; +inline BoundVector someBoundParametersA() { + return {0.0, 0.0, 0.0, std::numbers::pi / 2, 1.0, 0.0}; +} + template inline BoundVector generateBoundParameters( generator_t& rng, const GenerateBoundParametersOptions& options) { @@ -242,6 +260,10 @@ struct GenerateFreeParametersOptions { GenerateQoverPOptions qOverP; }; +inline FreeVector someFreeParametersA() { + return {0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0}; +} + template inline FreeVector generateFreeParameters( generator_t& rng, const GenerateFreeParametersOptions& options) { diff --git a/Core/include/Acts/EventData/detail/TestTrackState.hpp b/Core/include/Acts/EventData/detail/TestTrackState.hpp index c1471566c5b..514ae073c2f 100644 --- a/Core/include/Acts/EventData/detail/TestTrackState.hpp +++ b/Core/include/Acts/EventData/detail/TestTrackState.hpp @@ -38,11 +38,11 @@ struct TestTrackState { : surface( CurvilinearSurface(Vector3::Zero(), Vector3::UnitZ()).surface()), // set bogus parameters first since they are not default-constructible - predicted(surface, BoundVector::Zero(), std::nullopt, + predicted(surface, someBoundParametersA(), std::nullopt, ParticleHypothesis::pion()), - filtered(surface, BoundVector::Zero(), std::nullopt, + filtered(surface, someBoundParametersA(), std::nullopt, ParticleHypothesis::pion()), - smoothed(surface, BoundVector::Zero(), std::nullopt, + smoothed(surface, someBoundParametersA(), std::nullopt, ParticleHypothesis::pion()), jacobian(BoundMatrix::Identity()), chi2(std::chi_squared_distribution(measdim)(rng)), diff --git a/Core/include/Acts/Geometry/NavigationPolicyFactory.hpp b/Core/include/Acts/Geometry/NavigationPolicyFactory.hpp index 4d9542af826..ed9acce9aa2 100644 --- a/Core/include/Acts/Geometry/NavigationPolicyFactory.hpp +++ b/Core/include/Acts/Geometry/NavigationPolicyFactory.hpp @@ -54,17 +54,12 @@ class NavigationPolicyFactory { namespace detail { template -concept NavigationPolicyIsolatedFactoryConcept = requires(F f) { - { - f(std::declval(), - std::declval(), std::declval(), - std::declval()...) - } -> std::derived_from; - - requires NavigationPolicyConcept(), - std::declval(), std::declval(), - std::declval()...))>; +concept NavigationPolicyIsolatedFactoryConcept = requires( + F f, const GeometryContext& gctx, const TrackingVolume& volume, + const Logger& logger, Args&&... args) { + { f(gctx, volume, logger, args...) } -> std::derived_from; + + requires NavigationPolicyConcept; requires(std::is_copy_constructible_v && ...); }; @@ -205,7 +200,7 @@ class NavigationPolicyFactoryImpl : public NavigationPolicyFactory { template friend class NavigationPolicyFactoryImpl; - NavigationPolicyFactoryImpl(std::tuple&& factories) + explicit NavigationPolicyFactoryImpl(std::tuple&& factories) : m_factories(std::move(factories)) {} std::tuple m_factories; diff --git a/Core/include/Acts/Geometry/TrackingVolume.hpp b/Core/include/Acts/Geometry/TrackingVolume.hpp index f47f8b440fc..664faa39d46 100644 --- a/Core/include/Acts/Geometry/TrackingVolume.hpp +++ b/Core/include/Acts/Geometry/TrackingVolume.hpp @@ -117,8 +117,8 @@ class TrackingVolume : public Volume { ~TrackingVolume() override; TrackingVolume(const TrackingVolume&) = delete; TrackingVolume& operator=(const TrackingVolume&) = delete; - TrackingVolume(TrackingVolume&&); - TrackingVolume& operator=(TrackingVolume&&); + TrackingVolume(TrackingVolume&&) noexcept; + TrackingVolume& operator=(TrackingVolume&&) noexcept; /// Constructor for a container Volume /// - vacuum filled volume either as a for other tracking volumes diff --git a/Core/include/Acts/Navigation/INavigationPolicy.hpp b/Core/include/Acts/Navigation/INavigationPolicy.hpp index 6794dbdea09..de98f8d26c3 100644 --- a/Core/include/Acts/Navigation/INavigationPolicy.hpp +++ b/Core/include/Acts/Navigation/INavigationPolicy.hpp @@ -42,9 +42,11 @@ class INavigationPolicy { public: /// Noop update function that is suitable as a default for default navigation /// delegates. - static void noopInitializeCandidates(const NavigationArguments& /*unused*/, - AppendOnlyNavigationStream& /*unused*/, - const Logger& /*unused*/) {} + static void noopInitializeCandidates( + const NavigationArguments& /*unused*/, + const AppendOnlyNavigationStream& /*unused*/, const Logger& /*unused*/) { + // This is a noop + } /// Virtual destructor so policies can be held through this base class. virtual ~INavigationPolicy() = default; diff --git a/Core/include/Acts/Navigation/MultiNavigationPolicy.hpp b/Core/include/Acts/Navigation/MultiNavigationPolicy.hpp index c5c8a397aea..e66872fc39e 100644 --- a/Core/include/Acts/Navigation/MultiNavigationPolicy.hpp +++ b/Core/include/Acts/Navigation/MultiNavigationPolicy.hpp @@ -30,7 +30,7 @@ class MultiNavigationPolicy final : public MultiNavigationPolicyBase { public: /// Constructor from a set of child policies. /// @param policies The child policies - MultiNavigationPolicy(Policies&&... policies) + explicit MultiNavigationPolicy(Policies&&... policies) : m_policies{std::move(policies)...} {} /// Implementation of the connection to a navigation delegate. diff --git a/Core/include/Acts/Navigation/SurfaceArrayNavigationPolicy.hpp b/Core/include/Acts/Navigation/SurfaceArrayNavigationPolicy.hpp index 9143cf0087d..731122caff7 100644 --- a/Core/include/Acts/Navigation/SurfaceArrayNavigationPolicy.hpp +++ b/Core/include/Acts/Navigation/SurfaceArrayNavigationPolicy.hpp @@ -60,13 +60,14 @@ class SurfaceArrayNavigationPolicy : public INavigationPolicy { friend std::ostream& operator<<(std::ostream& os, const LayerType& layerType) { switch (layerType) { - case LayerType::Cylinder: + using enum LayerType; + case Cylinder: os << "Cylinder"; break; - case LayerType::Disc: + case Disc: os << "Disc"; break; - case LayerType::Plane: + case Plane: os << "Plane"; break; } diff --git a/Core/include/Acts/Seeding/GbtsTrackingFilter.hpp b/Core/include/Acts/Seeding/GbtsTrackingFilter.hpp index 86640eb946b..c2111482947 100644 --- a/Core/include/Acts/Seeding/GbtsTrackingFilter.hpp +++ b/Core/include/Acts/Seeding/GbtsTrackingFilter.hpp @@ -10,6 +10,7 @@ // TODO: update to C++17 style #include "Acts/Seeding/GbtsDataStorage.hpp" //includes geo which has trigindetsilayer, may move this to trigbase +#include "Acts/Utilities/Logger.hpp" #include #include @@ -112,8 +113,11 @@ template class GbtsTrackingFilter { public: GbtsTrackingFilter(const std::vector& g, - std::vector>& sb) - : m_geo(g), m_segStore(sb) {} + std::vector>& sb, + std::unique_ptr logger = + Acts::getDefaultLogger("Filter", + Acts::Logging::Level::INFO)) + : m_geo(g), m_segStore(sb), m_logger(std::move(logger)) {} void followTrack(Acts::GbtsEdge* pS, GbtsEdgeState& output) { @@ -233,11 +237,15 @@ class GbtsTrackingFilter { const float add_hit = 14.0; if (ts.m_Cx[2][2] < 0.0 || ts.m_Cx[1][1] < 0.0 || ts.m_Cx[0][0] < 0.0) { - std::cout << "Negative cov_x" << std::endl; + ACTS_WARNING("Negative covariance detected in X components: " + << "cov[2][2]=" << ts.m_Cx[2][2] << " cov[1][1]=" + << ts.m_Cx[1][1] << " cov[0][0]=" << ts.m_Cx[0][0]); } if (ts.m_Cy[1][1] < 0.0 || ts.m_Cy[0][0] < 0.0) { - std::cout << "Negative cov_y" << std::endl; + ACTS_WARNING("Negative covariance detected in Y components: " + << "cov[1][1]=" << ts.m_Cy[1][1] + << " cov[0][0]=" << ts.m_Cy[0][0]); } // add ms. @@ -380,4 +388,7 @@ class GbtsTrackingFilter { GbtsEdgeState m_stateStore[MAX_EDGE_STATE]; int m_globalStateCounter{0}; + + const Acts::Logger& logger() const { return *m_logger; } + std::unique_ptr m_logger{nullptr}; }; diff --git a/Core/include/Acts/Seeding/SeedFilter.hpp b/Core/include/Acts/Seeding/SeedFilter.hpp index 3c8720c42f0..d6e1479a9cd 100644 --- a/Core/include/Acts/Seeding/SeedFilter.hpp +++ b/Core/include/Acts/Seeding/SeedFilter.hpp @@ -13,6 +13,7 @@ #include "Acts/Seeding/CandidatesForMiddleSp.hpp" #include "Acts/Seeding/IExperimentCuts.hpp" #include "Acts/Seeding/SeedFilterConfig.hpp" +#include "Acts/Utilities/Logger.hpp" #include #include @@ -37,8 +38,15 @@ struct SeedFilterState { template class SeedFilter final { public: - SeedFilter(SeedFilterConfig config, + SeedFilter(const SeedFilterConfig& config, IExperimentCuts* expCuts = nullptr); + SeedFilter(const SeedFilterConfig& config, + std::unique_ptr logger, + IExperimentCuts* expCuts = nullptr); + SeedFilter(const SeedFilter&) = delete; + SeedFilter& operator=(const SeedFilter&) = delete; + SeedFilter(SeedFilter&&) noexcept = default; + SeedFilter& operator=(SeedFilter&&) noexcept = default; SeedFilter() = delete; ~SeedFilter() = default; @@ -95,7 +103,11 @@ class SeedFilter final { } private: + const Logger& logger() const { return *m_logger; } + const SeedFilterConfig m_cfg; + std::unique_ptr m_logger = + Acts::getDefaultLogger("Filter", Logging::Level::INFO); const IExperimentCuts* m_experimentCuts; }; } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFilter.ipp b/Core/include/Acts/Seeding/SeedFilter.ipp index 5e54789b714..66498916a7e 100644 --- a/Core/include/Acts/Seeding/SeedFilter.ipp +++ b/Core/include/Acts/Seeding/SeedFilter.ipp @@ -16,7 +16,7 @@ namespace Acts { // constructor template SeedFilter::SeedFilter( - SeedFilterConfig config, + const SeedFilterConfig& config, IExperimentCuts* expCuts /* = 0*/) : m_cfg(config), m_experimentCuts(expCuts) { if (!config.isInInternalUnits) { @@ -24,6 +24,18 @@ SeedFilter::SeedFilter( "SeedFilterConfig not in ACTS internal units in SeedFilter"); } } + +template +SeedFilter::SeedFilter( + const SeedFilterConfig& config, std::unique_ptr logger, + IExperimentCuts* expCuts /* = 0*/) + : m_cfg(config), m_logger(std::move(logger)), m_experimentCuts(expCuts) { + if (!config.isInInternalUnits) { + throw std::runtime_error( + "SeedFilterConfig not in ACTS internal units in SeedFilter"); + } +} + // function to filter seeds based on all seeds with same bottom- and // middle-spacepoint. // return vector must contain weight of each seed @@ -295,10 +307,14 @@ void SeedFilter::filterSeeds_1SpFixed( seed.setVertexZ(zOrigin); seed.setQuality(bestSeedQuality); + ACTS_VERBOSE("Adding seed: [b=" << bottom->index() << ", m=" + << medium->index() << ", t=" << top->index() + << "], quality=" << bestSeedQuality + << ", vertexZ=" << zOrigin); Acts::detail::pushBackOrInsertAtEnd(outputCollection, std::move(seed)); - ++numTotalSeeds; } + ACTS_VERBOSE("Identified " << numTotalSeeds << " seeds"); } } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinder.hpp b/Core/include/Acts/Seeding/SeedFinder.hpp index cfe9908574c..35a101015a5 100644 --- a/Core/include/Acts/Seeding/SeedFinder.hpp +++ b/Core/include/Acts/Seeding/SeedFinder.hpp @@ -18,6 +18,7 @@ #include "Acts/Seeding/SeedFinderUtils.hpp" #include "Acts/Seeding/SpacePointGrid.hpp" #include "Acts/Seeding/detail/UtilityFunctions.hpp" +#include "Acts/Utilities/Logger.hpp" #include #include @@ -37,11 +38,11 @@ concept GridBinCollection = std::ranges::random_access_range && std::same_as; -template -concept CollectionStoresSeedsTo = requires(Coll coll, external_t sp) { - Acts::detail::pushBackOrInsertAtEnd(coll, - Acts::Seed(sp, sp, sp)); -}; +template +concept CollectionStoresSeedsTo = + requires(collection_t coll, Acts::Seed seed) { + Acts::detail::pushBackOrInsertAtEnd(coll, seed); + }; enum class SpacePointCandidateType : short { eBottom, eTop }; @@ -87,7 +88,14 @@ class SeedFinder { /// The only constructor. Requires a config object. /// @param config the configuration for the SeedFinder - SeedFinder(const Acts::SeedFinderConfig& config); + /// @param logger the ACTS logger + SeedFinder(const Acts::SeedFinderConfig& config, + std::unique_ptr logger = + getDefaultLogger("Finder", Logging::Level::INFO)); + SeedFinder(SeedFinder&&) noexcept = + default; + SeedFinder& operator=(SeedFinder&&) noexcept = default; ~SeedFinder() = default; /** @name Disallow default instantiation, copy, assignment */ //@{ @@ -95,7 +103,7 @@ class SeedFinder { SeedFinder(const SeedFinder&) = delete; SeedFinder& operator=( - const SeedFinder&) = default; + const SeedFinder&) = delete; //@} /// Create all seeds from the space points in the three iterators. @@ -172,7 +180,10 @@ class SeedFinder { SeedingState& state) const; private: + const Logger& logger() const { return *m_logger; } + Acts::SeedFinderConfig m_config; + std::unique_ptr m_logger{nullptr}; }; } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinder.ipp b/Core/include/Acts/Seeding/SeedFinder.ipp index 45b95ec05d5..e6ccb9f589f 100644 --- a/Core/include/Acts/Seeding/SeedFinder.ipp +++ b/Core/include/Acts/Seeding/SeedFinder.ipp @@ -15,8 +15,9 @@ namespace Acts { template SeedFinder::SeedFinder( - const Acts::SeedFinderConfig& config) - : m_config(config) { + const Acts::SeedFinderConfig& config, + std::unique_ptr logger) + : m_config(config), m_logger(std::move(logger)) { if (!config.isInInternalUnits) { throw std::runtime_error( "SeedFinderConfig not in ACTS internal units in SeedFinder"); @@ -110,6 +111,12 @@ void SeedFinder::createSeedsForGroup( // same z-bin auto [minRadiusRangeForMiddle, maxRadiusRangeForMiddle] = retrieveRadiusRangeForMiddle(*middleSPs.front(), rMiddleSPRange); + ACTS_VERBOSE("Current global bin: " << middleSPsIdx << ", z value of " + << middleSPs.front()->z()); + ACTS_VERBOSE("Validity range (radius) for the middle space point is [" + << minRadiusRangeForMiddle << ", " << maxRadiusRangeForMiddle + << "]"); + for (const external_spacepoint_t* spM : middleSPs) { const float rM = spM->radius(); @@ -136,6 +143,7 @@ void SeedFinder::createSeedsForGroup( // no top SP found -> try next spM if (state.compatTopSP.empty()) { + ACTS_VERBOSE("No compatible Tops, moving to next middle candidate"); continue; } @@ -157,6 +165,10 @@ void SeedFinder::createSeedsForGroup( seedFilterState.rMaxSeedConf = seedConfRange.rMaxSeedConf; // continue if number of top SPs is smaller than minimum if (state.compatTopSP.size() < seedFilterState.nTopSeedConf) { + ACTS_VERBOSE( + "Number of top SPs is " + << state.compatTopSP.size() + << " and is smaller than minimum, moving to next middle candidate"); continue; } } @@ -170,9 +182,14 @@ void SeedFinder::createSeedsForGroup( // no bottom SP found -> try next spM if (state.compatBottomSP.empty()) { + ACTS_VERBOSE("No compatible Bottoms, moving to next middle candidate"); continue; } + ACTS_VERBOSE("Candidates: " << state.compatBottomSP.size() + << " bottoms and " << state.compatTopSP.size() + << " tops for middle candidate indexed " + << spM->index()); // filter candidates if (m_config.useDetailedDoubleMeasurementInfo) { filterCandidates( diff --git a/Core/include/Acts/Seeding/SeedFinderGbts.hpp b/Core/include/Acts/Seeding/SeedFinderGbts.hpp index aa8646c09d7..a7ca8f55268 100644 --- a/Core/include/Acts/Seeding/SeedFinderGbts.hpp +++ b/Core/include/Acts/Seeding/SeedFinderGbts.hpp @@ -15,7 +15,7 @@ #include "Acts/Seeding/SeedFinderConfig.hpp" #include "Acts/Seeding/SeedFinderGbtsConfig.hpp" #include "Acts/TrackFinding/RoiDescriptor.hpp" -#include "Acts/Utilities/KDTree.hpp" +#include "Acts/Utilities/Logger.hpp" #include #include @@ -46,14 +46,15 @@ class SeedFinderGbts { static constexpr std::size_t NDims = 3; using seed_t = Seed; - // using internal_sp_t = InternalSpacePoint; - // using tree_t = KDTree; // constructors SeedFinderGbts(const SeedFinderGbtsConfig &config, - const GbtsGeometry &gbtsgeo); + const GbtsGeometry &gbtsgeo, + std::unique_ptr logger = + Acts::getDefaultLogger("Finder", + Acts::Logging::Level::INFO)); - ~SeedFinderGbts(); //!!! is it dangerous not to use default? got def in ipp + ~SeedFinderGbts() = default; SeedFinderGbts() = default; SeedFinderGbts(const SeedFinderGbts &) = delete; SeedFinderGbts &operator=( @@ -85,10 +86,14 @@ class SeedFinderGbts { const Acts::GbtsGeometry &gbtsgeo); // needs to be member of class so can accessed by all member functions - GbtsDataStorage *m_storage; + std::unique_ptr> m_storage{nullptr}; // for create seeds: std::vector> m_triplets; + + const Acts::Logger &logger() const { return *m_logger; } + std::unique_ptr m_logger = + Acts::getDefaultLogger("Finder", Acts::Logging::Level::INFO); }; } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinderGbts.ipp b/Core/include/Acts/Seeding/SeedFinderGbts.ipp index 91f5f3f64e4..ab557d86977 100644 --- a/Core/include/Acts/Seeding/SeedFinderGbts.ipp +++ b/Core/include/Acts/Seeding/SeedFinderGbts.ipp @@ -34,22 +34,18 @@ namespace Acts { template SeedFinderGbts::SeedFinderGbts( const SeedFinderGbtsConfig& config, - const GbtsGeometry& gbtsGeo) - : m_config(config) { - m_storage = new GbtsDataStorage(gbtsGeo); -} - -template -SeedFinderGbts::~SeedFinderGbts() { - delete m_storage; - - m_storage = nullptr; -} + const GbtsGeometry& gbtsGeo, + std::unique_ptr logger) + : m_config(config), + m_storage( + std::make_unique>(gbtsGeo)), + m_logger(std::move(logger)) {} // define loadspace points function template void SeedFinderGbts::loadSpacePoints( const std::vector>& gbtsSPvect) { + ACTS_VERBOSE("Loading space points"); for (const auto& gbtssp : gbtsSPvect) { bool is_Pixel = gbtssp.isPixel(); if (!is_Pixel) { @@ -70,6 +66,7 @@ void SeedFinderGbts::runGbts_TrackFinder( std::vector>& vTracks, const Acts::RoiDescriptor& roi, const Acts::GbtsGeometry& gbtsGeo) { + ACTS_VERBOSE("Running GBTS Track Finder"); const float min_z0 = roi.zedMinus(); const float max_z0 = roi.zedPlus(); const float cut_zMinU = min_z0 + m_config.maxOuterRadius * roi.dzdrMinus(); @@ -159,7 +156,7 @@ void SeedFinderGbts::runGbts_TrackFinder( .m_phiSliceWidth; // the default sliding window along phi if (m_config.m_useEtaBinning) { - deltaPhi = 0.001f + m_maxCurv * std::fabs(rb2 - rb1); + deltaPhi = 0.001f + m_maxCurv * std::abs(rb2 - rb1); } unsigned int first_it = 0; @@ -219,7 +216,7 @@ void SeedFinderGbts::runGbts_TrackFinder( float dz = z2 - z1; float tau = dz / dr; - float ftau = std::fabs(tau); + float ftau = std::abs(tau); if (ftau > 36.0) { continue; } @@ -288,17 +285,18 @@ void SeedFinderGbts::runGbts_TrackFinder( float tau2 = edgeStorage.at(n2_in_idx).m_p[0]; float tau_ratio = tau2 * uat_1 - 1.0f; - if (std::fabs(tau_ratio) > - m_config.cut_tau_ratio_max) { // bad - // match + // bad match + if (std::abs(tau_ratio) > m_config.cut_tau_ratio_max) { continue; } - isGood = true; // good match found + + // good match found + isGood = true; break; } } if (!isGood) { - continue; // no moatch found, skip creating [n1 <- n2] edge + continue; // no match found, skip creating [n1 <- n2] edge } float curv = D * std::sqrt(L2); // signed curvature @@ -327,6 +325,7 @@ void SeedFinderGbts::runGbts_TrackFinder( m_storage->getConnectingNodes(vNodes); if (vNodes.empty()) { + ACTS_VERBOSE("No nodes"); return; } @@ -505,8 +504,9 @@ void SeedFinderGbts::runGbts_TrackFinder( // backtracking - GbtsTrackingFilter tFilter(m_config.m_layerGeometry, - edgeStorage); + GbtsTrackingFilter tFilter( + m_config.m_layerGeometry, edgeStorage, + logger().cloneWithSuffix("GbtsFilter")); for (auto pS : vSeeds) { if (pS->m_level == -1) { @@ -650,6 +650,7 @@ void SeedFinderGbts::createSeeds( const Acts::RoiDescriptor& roi, const Acts::GbtsGeometry& gbtsGeo, output_container_t& out_cont) { + ACTS_VERBOSE("Creating seeds"); std::vector> vTracks; // make empty vector diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp b/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp index 2221bff6626..594e35b44cc 100644 --- a/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp @@ -13,6 +13,7 @@ #include "Acts/Seeding/SeedFinderConfig.hpp" #include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp" #include "Acts/Utilities/KDTree.hpp" +#include "Acts/Utilities/Logger.hpp" #include #include @@ -53,20 +54,26 @@ class SeedFinderOrthogonal { * @brief Construct a new orthogonal seed finder. * * @param config The configuration parameters for this seed finder. + * @param logger The ACTS logger. */ SeedFinderOrthogonal( - const Acts::SeedFinderOrthogonalConfig &config); - + const Acts::SeedFinderOrthogonalConfig &config, + std::unique_ptr logger = + Acts::getDefaultLogger("Finder", Logging::Level::INFO)); /** * @brief Destroy the orthogonal seed finder object. */ ~SeedFinderOrthogonal() = default; - SeedFinderOrthogonal() = default; + SeedFinderOrthogonal() = delete; SeedFinderOrthogonal(const SeedFinderOrthogonal &) = delete; SeedFinderOrthogonal &operator=( - const SeedFinderOrthogonal &) = default; + const SeedFinderOrthogonal &) = delete; + SeedFinderOrthogonal( + SeedFinderOrthogonal &&) noexcept = default; + SeedFinderOrthogonal &operator=( + SeedFinderOrthogonal &&) noexcept = default; /** * @brief Perform seed finding, appending seeds to a container. @@ -238,6 +245,16 @@ class SeedFinderOrthogonal { * @brief The configuration for the seeding algorithm. */ Acts::SeedFinderOrthogonalConfig m_config; + + /** + * @brief Get the logger. + */ + const Logger &logger() const { return *m_logger; } + + /** + * @brief The logger + */ + std::unique_ptr m_logger{nullptr}; }; } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp b/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp index edd8e0f3074..a725d68e92d 100644 --- a/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp @@ -226,7 +226,7 @@ bool SeedFinderOrthogonal::validTuple( * Cut: Ensure that the forward angle (z / r) lies within reasonable bounds, * which is to say the absolute value must be smaller than the max cot(θ). */ - if (std::fabs(cotTheta) > m_config.cotThetaMax) { + if (std::abs(cotTheta) > m_config.cotThetaMax) { return false; } @@ -244,8 +244,9 @@ bool SeedFinderOrthogonal::validTuple( template SeedFinderOrthogonal::SeedFinderOrthogonal( - const SeedFinderOrthogonalConfig &config) - : m_config(config) { + const SeedFinderOrthogonalConfig &config, + std::unique_ptr logger) + : m_config(config), m_logger(std::move(logger)) { if (!config.isInInternalUnits) { throw std::runtime_error( "SeedFinderOrthogonalConfig not in ACTS internal units in " @@ -687,6 +688,7 @@ auto SeedFinderOrthogonal::createTree( const std::vector &spacePoints) const -> tree_t { std::vector points; + points.reserve(spacePoints.size()); /* * For every input point, we create a coordinate-pointer pair, which we then @@ -703,6 +705,8 @@ auto SeedFinderOrthogonal::createTree( points.emplace_back(point, sp); } + ACTS_VERBOSE("Created k-d tree populated with " << points.size() + << " space points"); return tree_t(std::move(points)); } @@ -711,6 +715,7 @@ template void SeedFinderOrthogonal::createSeeds( const Acts::SeedFinderOptions &options, const input_container_t &spacePoints, output_container_t &out_cont) const { + ACTS_VERBOSE("Creating seeds with Orthogonal strategy"); if (!options.isInInternalUnits) { throw std::runtime_error( "SeedFinderOptions not in ACTS internal units in " @@ -734,6 +739,7 @@ void SeedFinderOrthogonal::createSeeds( * take each external spacepoint, allocate a corresponding internal space * point, and save it in a vector. */ + ACTS_VERBOSE("Running on " << spacePoints.size() << " input space points"); Acts::Extent rRangeSPExtent; std::vector internal_sps; internal_sps.reserve(spacePoints.size()); @@ -746,6 +752,8 @@ void SeedFinderOrthogonal::createSeeds( rRangeSPExtent.extend({p.x(), p.y(), p.z()}); internal_sps.push_back(&p); } + ACTS_VERBOSE(rRangeSPExtent); + // variable middle SP radial region of interest const Acts::Range1D rMiddleSPRange( std::floor(rRangeSPExtent.min(Acts::BinningValue::binR) / 2) * 2 + diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp index 4410f1c125b..a2a214fa39b 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp @@ -11,6 +11,7 @@ #include "Acts/Seeding/BinnedGroup.hpp" #include "Acts/Seeding/SeedFinderConfig.hpp" #include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/Logger.hpp" #include #include @@ -158,7 +159,8 @@ class CylindricalSpacePointGridCreator { template static Acts::CylindricalSpacePointGrid createGrid( const Acts::CylindricalSpacePointGridConfig& _config, - const Acts::CylindricalSpacePointGridOptions& _options); + const Acts::CylindricalSpacePointGridOptions& _options, + const Acts::Logger& logger = Acts::getDummyLogger()); template @@ -167,7 +169,8 @@ class CylindricalSpacePointGridCreator { const Acts::SeedFinderOptions& options, Acts::CylindricalSpacePointGrid& grid, external_spacepoint_iterator_t spBegin, - external_spacepoint_iterator_t spEnd); + external_spacepoint_iterator_t spEnd, + const Acts::Logger& logger = Acts::getDummyLogger()); template requires std::ranges::range && @@ -177,7 +180,8 @@ class CylindricalSpacePointGridCreator { const Acts::SeedFinderConfig& config, const Acts::SeedFinderOptions& options, Acts::CylindricalSpacePointGrid& grid, - const external_collection_t& collection); + const external_collection_t& collection, + const Acts::Logger& logger = Acts::getDummyLogger()); }; } // namespace Acts diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp index 523d470eae5..b99cdbdfaf9 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp @@ -12,7 +12,8 @@ template Acts::CylindricalSpacePointGrid Acts::CylindricalSpacePointGridCreator::createGrid( const Acts::CylindricalSpacePointGridConfig& config, - const Acts::CylindricalSpacePointGridOptions& options) { + const Acts::CylindricalSpacePointGridOptions& options, + const Acts::Logger& logger) { if (!config.isInInternalUnits) { throw std::runtime_error( "CylindricalSpacePointGridConfig not in ACTS internal units in " @@ -30,6 +31,9 @@ Acts::CylindricalSpacePointGridCreator::createGrid( // for no magnetic field, create 100 phi-bins if (options.bFieldInZ == 0) { phiBins = 100; + ACTS_VERBOSE( + "B-Field is 0 (z-coordinate), setting the number of bins in phi to " + << phiBins); } else { // calculate circle intersections of helix and max detector radius float minHelixRadius = @@ -39,7 +43,7 @@ Acts::CylindricalSpacePointGridCreator::createGrid( // = pT[MeV] / (300 *Bz[kT]) // sanity check: if yOuter takes the square root of a negative number - if (minHelixRadius < config.rMax / 2) { + if (minHelixRadius < config.rMax * 0.5) { throw std::domain_error( "The value of minHelixRadius cannot be smaller than rMax / 2. Please " "check the configuration of bFieldInZ and minPt"); @@ -141,6 +145,12 @@ Acts::CylindricalSpacePointGridCreator::createGrid( Axis zAxis(std::move(zValues)); Axis rAxis(std::move(rValues)); + + ACTS_VERBOSE("Defining Grid:"); + ACTS_VERBOSE("- Phi Axis: " << phiAxis); + ACTS_VERBOSE("- Z axis : " << zAxis); + ACTS_VERBOSE("- R axis : " << rAxis); + return Acts::CylindricalSpacePointGrid( std::make_tuple(std::move(phiAxis), std::move(zAxis), std::move(rAxis))); } @@ -152,7 +162,7 @@ void Acts::CylindricalSpacePointGridCreator::fillGrid( const Acts::SeedFinderOptions& options, Acts::CylindricalSpacePointGrid& grid, external_spacepoint_iterator_t spBegin, - external_spacepoint_iterator_t spEnd) { + external_spacepoint_iterator_t spEnd, const Acts::Logger& logger) { if (!config.isInInternalUnits) { throw std::runtime_error( "SeedFinderConfig not in ACTS internal units in BinnedSPGroup"); @@ -180,9 +190,10 @@ void Acts::CylindricalSpacePointGridCreator::fillGrid( std::vector rBinsIndex; rBinsIndex.reserve(grid.size()); + ACTS_VERBOSE("Fetching " << std::distance(spBegin, spEnd) + << " space points to the grid"); std::size_t counter = 0ul; - for (external_spacepoint_iterator_t it = spBegin; it != spEnd; - it++, ++counter) { + for (external_spacepoint_iterator_t it = spBegin; it != spEnd; ++it) { const external_spacepoint_t& sp = *it; // remove SPs according to experiment specific cuts @@ -199,6 +210,7 @@ void Acts::CylindricalSpacePointGridCreator::fillGrid( std::size_t globIndex = grid.globalBinFromPosition(position); auto& rbin = grid.at(globIndex); rbin.push_back(&sp); + ++counter; // keep track of the bins we modify so that we can later sort the SPs in // those bins only @@ -213,6 +225,9 @@ void Acts::CylindricalSpacePointGridCreator::fillGrid( auto& rbin = grid.atPosition(binIndex); std::ranges::sort(rbin, {}, [](const auto& rb) { return rb->radius(); }); } + + ACTS_VERBOSE( + "Number of space points inserted (within grid range): " << counter); } template @@ -223,8 +238,8 @@ void Acts::CylindricalSpacePointGridCreator::fillGrid( const Acts::SeedFinderConfig& config, const Acts::SeedFinderOptions& options, Acts::CylindricalSpacePointGrid& grid, - const external_collection_t& collection) { + const external_collection_t& collection, const Acts::Logger& logger) { Acts::CylindricalSpacePointGridCreator::fillGrid( config, options, grid, std::ranges::begin(collection), - std::ranges::end(collection)); + std::ranges::end(collection), logger); } diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index 05e236e1e29..35acbc49feb 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -399,7 +399,7 @@ class CombinatorialKalmanFilter { measurementSelector(trackStateCandidates, isOutlier, logger); if (!selectorResult.ok()) { ACTS_ERROR("Selection of calibrated measurements failed: " - << selectorResult.error()); + << selectorResult.error().message()); resultTrackStateList = ResultTrackStateList::failure(selectorResult.error()); } else { diff --git a/Core/include/Acts/TrackFinding/TrackSelector.hpp b/Core/include/Acts/TrackFinding/TrackSelector.hpp index bca7fa8d2a0..459ab996857 100644 --- a/Core/include/Acts/TrackFinding/TrackSelector.hpp +++ b/Core/include/Acts/TrackFinding/TrackSelector.hpp @@ -427,8 +427,9 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const { const Config* cutsPtr{nullptr}; if (!m_isUnbinned) { - if (absEta() < m_cfg.absEtaEdges.front() || - _absEta >= m_cfg.absEtaEdges.back()) { + // return false if |eta| is outside its range, or nan. + if (!(absEta() >= m_cfg.absEtaEdges.front() && + _absEta < m_cfg.absEtaEdges.back())) { return false; } cutsPtr = &m_cfg.getCuts(_eta); diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index 92a6a4f76ba..968843fa479 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -18,6 +18,7 @@ #include "Acts/EventData/SourceLink.hpp" #include "Acts/EventData/TrackContainerFrontendConcept.hpp" #include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/TrackProxyConcept.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/EventData/VectorTrackContainer.hpp" #include "Acts/Geometry/GeometryContext.hpp" @@ -372,30 +373,30 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, ACTS_VERBOSE( "Contributions in addMeasurementToGx2fSums:\n" - << "kMeasDim: " << kMeasDim << "\n" - << "predicted" << predicted.transpose() << "\n" - << "measurement: " << measurement.transpose() << "\n" - << "covarianceMeasurement:\n" + << " kMeasDim: " << kMeasDim << "\n" + << " predicted: " << predicted.transpose() << "\n" + << " measurement: " << measurement.transpose() << "\n" + << " covarianceMeasurement:\n" << covarianceMeasurement << "\n" - << "projector:\n" + << " projector:\n" << projector.eval() << "\n" - << "projJacobian:\n" + << " projJacobian:\n" << projJacobian.eval() << "\n" - << "projPredicted: " << (projPredicted.transpose()).eval() << "\n" - << "residual: " << (residual.transpose()).eval() << "\n" - << "extendedJacobian:\n" + << " projPredicted: " << (projPredicted.transpose()).eval() << "\n" + << " residual: " << (residual.transpose()).eval() << "\n" + << " extendedJacobian:\n" << extendedJacobian << "\n" - << "aMatrixMeas:\n" + << " aMatrix contribution:\n" << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) .eval() << "\n" - << "bVectorMeas: " + << " bVector contribution: " << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval() << "\n" - << "chi2sumMeas: " + << " chi2sum contribution: " << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) << "\n" - << "safeInvCovMeasurement:\n" + << " safeInvCovMeasurement:\n" << (*safeInvCovMeasurement)); return; @@ -459,8 +460,8 @@ void addMaterialToGx2fSums( ACTS_VERBOSE( "Contributions in addMaterialToGx2fSums:\n" - << " invCov: " << invCov << "\n" - << " sinThetaLoc: " << sinThetaLoc << "\n" + << " invCov: " << invCov << "\n" + << " sinThetaLoc: " << sinThetaLoc << "\n" << " deltaPosition: " << deltaPosition << "\n" << " Phi:\n" << " scattering angle: " << scatteringAngles[eBoundPhi] << "\n" @@ -485,6 +486,103 @@ void addMaterialToGx2fSums( return; } +/// @brief Fill the GX2F system with data from a track +/// +/// This function processes a track proxy and updates the aMatrix, bVector, and +/// chi2 values for the GX2F fitting system. It considers material only if +/// multiple scattering is enabled. +/// +/// @tparam track_proxy_t The type of the track proxy +/// +/// @param track A mutable track proxy to operate on +/// @param aMatrixExtended The aMatrix, summing over second derivatives for the fitting system +/// @param bVectorExtended The bVector, summing over first derivatives for the fitting system +/// @param chi2sum Accumulated chi2 value of the system +/// @param countNdf The number of degrees of freedom counted so far +/// @param multipleScattering Flag to consider multiple scattering in the calculation +/// @param scatteringMap Map of geometry identifiers to scattering properties, containing all scattering angles and covariances +/// @param geoIdVector A vector to store geometry identifiers for tracking processed elements +/// @param logger A logger instance +template +void fillGx2fSystem( + const track_proxy_t track, Eigen::MatrixXd& aMatrixExtended, + Eigen::VectorXd& bVectorExtended, double& chi2sum, std::size_t& countNdf, + const bool multipleScattering, + const std::unordered_map& + scatteringMap, + std::vector& geoIdVector, const Logger& logger) { + std::vector jacobianFromStart; + jacobianFromStart.emplace_back(BoundMatrix::Identity()); + + for (const auto& trackState : track.trackStates()) { + // Get and store geoId for the current surface + const GeometryIdentifier geoId = trackState.referenceSurface().geometryId(); + ACTS_DEBUG("Start to investigate trackState on surface " << geoId); + const auto typeFlags = trackState.typeFlags(); + const bool stateHasMeasurement = + typeFlags.test(TrackStateFlag::MeasurementFlag); + const bool stateHasMaterial = typeFlags.test(TrackStateFlag::MaterialFlag); + + // First we figure out, if we would need to look into material + // surfaces at all. Later, we also check, if the material slab is + // valid, otherwise we modify this flag to ignore the material + // completely. + bool doMaterial = multipleScattering && stateHasMaterial; + if (doMaterial) { + const auto scatteringMapId = scatteringMap.find(geoId); + assert(scatteringMapId != scatteringMap.end() && + "No scattering angles found for material surface."); + doMaterial = doMaterial && scatteringMapId->second.materialIsValid(); + } + + // We only consider states with a measurement (and/or material) + if (!stateHasMeasurement && !doMaterial) { + ACTS_DEBUG(" Skip state."); + continue; + } + + // update all Jacobians from start + for (auto& jac : jacobianFromStart) { + jac = trackState.jacobian() * jac; + } + + // Handle measurement + if (stateHasMeasurement) { + ACTS_DEBUG(" Handle measurement."); + + const auto measDim = trackState.calibratedSize(); + + if (measDim < 1 || 6 < measDim) { + ACTS_ERROR("Can not process state with measurement with " + << measDim << " dimensions."); + throw std::domain_error( + "Found measurement with less than 1 or more than 6 dimension(s)."); + } + + countNdf += measDim; + + visit_measurement(measDim, [&](auto N) { + addMeasurementToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, + jacobianFromStart, trackState, logger); + }); + } + + // Handle material + if (doMaterial) { + ACTS_DEBUG(" Handle material"); + // Add for this material a new Jacobian, starting from this surface. + jacobianFromStart.emplace_back(BoundMatrix::Identity()); + + // Add the material contribution to the system + addMaterialToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, + geoIdVector.size(), scatteringMap, trackState, + logger); + + geoIdVector.emplace_back(geoId); + } + } +} + /// @brief Calculate and update the covariance of the fitted parameters /// /// This function calculates the covariance of the fitted parameters using @@ -727,13 +825,15 @@ class Gx2Fitter { if (doMaterial) { ACTS_DEBUG(" Update parameters with scattering angles."); const auto scatteringMapId = scatteringMap->find(geoId); - ACTS_VERBOSE(" scatteringAngles:\n" - << scatteringMapId->second.scatteringAngles() - << "\n boundParams before the update:\n" - << boundParams); + ACTS_VERBOSE( + " scatteringAngles: " + << scatteringMapId->second.scatteringAngles().transpose()); + ACTS_VERBOSE(" boundParams before the update: " + << boundParams.parameters().transpose()); boundParams.parameters() += scatteringMapId->second.scatteringAngles(); - ACTS_VERBOSE(" boundParams after the update:\n" << boundParams); + ACTS_VERBOSE(" boundParams after the update: " + << boundParams.parameters().transpose()); } // Fill the track state @@ -829,13 +929,15 @@ class Gx2Fitter { // multipleScattering and have material ACTS_DEBUG(" Update parameters with scattering angles."); const auto scatteringMapId = scatteringMap->find(geoId); - ACTS_VERBOSE(" scatteringAngles:\n" - << scatteringMapId->second.scatteringAngles() - << "\n boundParams before the update:\n" - << boundParams); + ACTS_VERBOSE( + " scatteringAngles: " + << scatteringMapId->second.scatteringAngles().transpose()); + ACTS_VERBOSE(" boundParams before the update: " + << boundParams.parameters().transpose()); boundParams.parameters() += scatteringMapId->second.scatteringAngles(); - ACTS_VERBOSE(" boundParams after the update:\n" << boundParams); + ACTS_VERBOSE(" boundParams after the update: " + << boundParams.parameters().transpose()); // Fill the track state trackStateProxy.smoothed() = boundParams.parameters(); @@ -998,7 +1100,7 @@ class Gx2Fitter { requires(!isDirectNavigator) { // Preprocess Measurements (SourceLinks -> map) - // To be able to find measurements later, we put them into a map + // To be able to find measurements later, we put them into a map. // We need to copy input SourceLinks anyway, so the map can own them. ACTS_VERBOSE("Preparing " << std::distance(it, end) << " input measurements"); @@ -1009,7 +1111,6 @@ class Gx2Fitter { auto geoId = gx2fOptions.extensions.surfaceAccessor(sl)->geometryId(); inputMeasurements.emplace(geoId, std::move(sl)); } - ACTS_VERBOSE("inputMeasurements.size() = " << inputMeasurements.size()); // Store, if we want to do multiple scattering. We still need to pass this // option to the Actor. @@ -1056,21 +1157,21 @@ class Gx2Fitter { // track parameters. BoundMatrix fullCovariancePredicted = BoundMatrix::Identity(); - ACTS_VERBOSE("params:\n" << params); + ACTS_VERBOSE("Initial parameters: " << params.parameters().transpose()); /// Actual Fitting ///////////////////////////////////////////////////////// ACTS_DEBUG("Start to iterate"); // Iterate the fit and improve result. Abort after n steps or after - // convergence - // nUpdate is initialized outside to save its state for the track + // convergence. + // nUpdate is initialized outside to save its state for the track. std::size_t nUpdate = 0; for (nUpdate = 0; nUpdate < gx2fOptions.nUpdateMax; nUpdate++) { ACTS_DEBUG("nUpdate = " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax); // update params params.parameters() += deltaParams; - ACTS_VERBOSE("updated params:\n" << params); + ACTS_VERBOSE("Updated parameters: " << params.parameters().transpose()); // set up propagator and co Acts::GeometryContext geoCtx = gx2fOptions.geoContext; @@ -1130,9 +1231,8 @@ class Gx2Fitter { tipIndex = gx2fResult.lastMeasurementIndex; // It could happen, that no measurements were found. Then the track would - // be empty and the following operations would be invalid. - // Usually, this only happens during the first iteration, due to bad - // initial parameters. + // be empty and the following operations would be invalid. Usually, this + // only happens during the first iteration, due to bad initial parameters. if (tipIndex == Acts::MultiTrajectoryTraits::kInvalid) { ACTS_INFO("Did not find any measurements in nUpdate " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax); @@ -1185,85 +1285,15 @@ class Gx2Fitter { Eigen::VectorXd bVectorExtended = Eigen::VectorXd::Zero(dimsExtendedParams); - std::vector jacobianFromStart; - jacobianFromStart.emplace_back(BoundMatrix::Identity()); - // This vector stores the IDs for each visited material in order. We use // it later for updating the scattering angles. We cannot use // scatteringMap directly, since we cannot guarantee, that we will visit // all stored material in each propagation. std::vector geoIdVector; - for (const auto& trackState : track.trackStates()) { - // Get and store geoId for the current surface - const GeometryIdentifier geoId = - trackState.referenceSurface().geometryId(); - ACTS_DEBUG("Start to investigate trackState on surface " << geoId); - const auto typeFlags = trackState.typeFlags(); - const bool stateHasMeasurement = - typeFlags.test(TrackStateFlag::MeasurementFlag); - const bool stateHasMaterial = - typeFlags.test(TrackStateFlag::MaterialFlag); - - // First we figure out, if we would need to look into material surfaces - // at all. Later, we also check, if the material slab is valid, - // otherwise we modify this flag to ignore the material completely. - bool doMaterial = multipleScattering && stateHasMaterial; - if (doMaterial) { - const auto scatteringMapId = scatteringMap.find(geoId); - assert(scatteringMapId != scatteringMap.end() && - "No scattering angles found for material surface."); - doMaterial = doMaterial && scatteringMapId->second.materialIsValid(); - } - - // We only consider states with a measurement (and/or material) - if (!stateHasMeasurement && !doMaterial) { - ACTS_DEBUG(" Skip state."); - continue; - } - - // update all Jacobians from start - for (auto& jac : jacobianFromStart) { - jac = trackState.jacobian() * jac; - } - - // Handle measurement - if (stateHasMeasurement) { - ACTS_DEBUG(" Handle measurement."); - - const auto measDim = trackState.calibratedSize(); - - if (measDim < 1 || 6 < measDim) { - ACTS_ERROR("Can not process state with measurement with " - << measDim << " dimensions."); - throw std::domain_error( - "Found measurement with less than 1 or more than 6 " - "dimension(s)."); - } - - countNdf += measDim; - - visit_measurement(measDim, [&](auto N) { - addMeasurementToGx2fSums(aMatrixExtended, bVectorExtended, - chi2sum, jacobianFromStart, trackState, - *m_addToSumLogger); - }); - } - - // Handle material - if (doMaterial) { - ACTS_DEBUG(" Handle material"); - // Add for this material a new Jacobian, starting from this surface. - jacobianFromStart.emplace_back(BoundMatrix::Identity()); - - // Add the material contribution to the system - addMaterialToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, - geoIdVector.size(), scatteringMap, trackState, - *m_addToSumLogger); - - geoIdVector.emplace_back(geoId); - } - } + fillGx2fSystem(track, aMatrixExtended, bVectorExtended, chi2sum, countNdf, + multipleScattering, scatteringMap, geoIdVector, + *m_addToSumLogger); // Get required number of degrees of freedom ndfSystem. // We have only 3 cases, because we always have l0, l1, phi, theta @@ -1279,12 +1309,11 @@ class Gx2Fitter { } // This check takes into account the evaluated dimensions of the - // measurements. To fit, we need at least NDF+1 measurements. However, - // we count n-dimensional measurements for n measurements, reducing the - // effective number of needed measurements. - // We might encounter the case, where we cannot use some (parts of a) - // measurements, maybe if we do not support that kind of measurement. This - // is also taken into account here. + // measurements. To fit, we need at least NDF+1 measurements. However, we + // count n-dimensional measurements for n measurements, reducing the + // effective number of needed measurements. We might encounter the case, + // where we cannot use some (parts of a) measurements, maybe if we do not + // support that kind of measurement. This is also taken into account here. // We skip the check during the first iteration, since we cannot guarantee // to hit all/enough measurement surfaces with the initial parameter // guess. @@ -1371,7 +1400,7 @@ class Gx2Fitter { oldChi2sum = chi2sum; } ACTS_DEBUG("Finished to iterate"); - ACTS_VERBOSE("final params:\n" << params); + ACTS_VERBOSE("Final parameters: " << params.parameters().transpose()); /// Finish Fitting ///////////////////////////////////////////////////////// ACTS_VERBOSE("Final scattering angles:"); @@ -1384,7 +1413,7 @@ class Gx2Fitter { << " )"); } - ACTS_VERBOSE("final covariance:\n" << fullCovariancePredicted); + ACTS_VERBOSE("Final covariance:\n" << fullCovariancePredicted); // Propagate again with the final covariance matrix. This is necessary to // obtain the propagated covariance for each state. @@ -1392,7 +1421,7 @@ class Gx2Fitter { // step, we will not ignore the boundary checks for measurement surfaces. We // want to create trackstates only on surfaces, that we actually hit. if (gx2fOptions.nUpdateMax > 0) { - ACTS_VERBOSE("final deltaParams:\n" << deltaParams); + ACTS_VERBOSE("Final delta parameters: " << deltaParams.transpose()); ACTS_VERBOSE("Propagate with the final covariance."); // update covariance params.covariance() = fullCovariancePredicted; diff --git a/Core/include/Acts/Utilities/DelegateChainBuilder.hpp b/Core/include/Acts/Utilities/DelegateChainBuilder.hpp index 6e4f2c89650..a24ab98ddb8 100644 --- a/Core/include/Acts/Utilities/DelegateChainBuilder.hpp +++ b/Core/include/Acts/Utilities/DelegateChainBuilder.hpp @@ -102,7 +102,7 @@ class DelegateChainBuilder, template static constexpr auto invoke(result_ptr result, const tuple_type* payloads, - callable_args... args) { + callable_args&&... args) { const auto& callable = findCallable(); if constexpr (!std::is_same_v, @@ -134,7 +134,7 @@ class DelegateChainBuilder, tuple_type m_payloads{}; - auto dispatch(callable_args... args) const { + auto dispatch(callable_args&&... args) const { if constexpr (std::is_same_v) { invoke(nullptr, &m_payloads, std::forward(args)...); } else { diff --git a/Core/src/EventData/CMakeLists.txt b/Core/src/EventData/CMakeLists.txt index 08a69711378..c425d533649 100644 --- a/Core/src/EventData/CMakeLists.txt +++ b/Core/src/EventData/CMakeLists.txt @@ -8,4 +8,5 @@ target_sources( TrackStatePropMask.cpp VectorMultiTrajectory.cpp VectorTrackContainer.cpp + TrackParameterHelpers.cpp ) diff --git a/Core/src/EventData/TrackParameterHelpers.cpp b/Core/src/EventData/TrackParameterHelpers.cpp new file mode 100644 index 00000000000..9e7ddb40abf --- /dev/null +++ b/Core/src/EventData/TrackParameterHelpers.cpp @@ -0,0 +1,64 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include "Acts/EventData/TrackParameterHelpers.hpp" + +#include "Acts/Utilities/AngleHelpers.hpp" +#include "Acts/Utilities/VectorHelpers.hpp" + +#include + +bool Acts::isBoundVectorValid(const BoundVector& v, bool validateAngleRange, + double epsilon, double maxAbsEta) { + constexpr auto pi = std::numbers::pi_v; + + bool loc0Valid = std::isfinite(v[eBoundLoc0]); + bool loc1Valid = std::isfinite(v[eBoundLoc1]); + bool phiValid = std::isfinite(v[eBoundPhi]); + bool thetaValid = std::isfinite(v[eBoundTheta]); + bool qOverPValid = std::isfinite(v[eBoundQOverP]); + bool timeValid = std::isfinite(v[eBoundTime]); + + if (validateAngleRange) { + phiValid = phiValid && (v[eBoundPhi] + epsilon >= -pi) && + (v[eBoundPhi] - epsilon < pi); + thetaValid = thetaValid && (v[eBoundTheta] + epsilon >= 0) && + (v[eBoundTheta] - epsilon <= pi); + } + + bool etaValid = true; + if (std::isfinite(maxAbsEta)) { + double eta = AngleHelpers::etaFromTheta(v[eBoundTheta]); + etaValid = std::isfinite(eta) && (std::abs(eta) - epsilon <= maxAbsEta); + } + + return loc0Valid && loc1Valid && phiValid && thetaValid && qOverPValid && + timeValid && etaValid; +} + +bool Acts::isFreeVectorValid(const FreeVector& v, double epsilon, + double maxAbsEta) { + bool pos0Valid = std::isfinite(v[eFreePos0]); + bool pos1Valid = std::isfinite(v[eFreePos1]); + bool pos2Valid = std::isfinite(v[eFreePos2]); + bool dir0Valid = std::isfinite(v[eFreeDir0]); + bool dir1Valid = std::isfinite(v[eFreeDir1]); + bool dir2Valid = std::isfinite(v[eFreeDir2]); + bool dirValid = (std::abs(v.segment<3>(eFreeDir0).norm() - 1) - epsilon <= 0); + bool qOverPValid = std::isfinite(v[eFreeQOverP]); + bool timeValid = std::isfinite(v[eFreeTime]); + + bool etaValid = true; + if (std::isfinite(maxAbsEta)) { + double eta = VectorHelpers::eta(v.segment<3>(eFreeDir0)); + etaValid = std::isfinite(eta) && (std::abs(eta) - epsilon <= maxAbsEta); + } + + return pos0Valid && pos1Valid && pos2Valid && dir0Valid && dir1Valid && + dir2Valid && dirValid && qOverPValid && timeValid && etaValid; +} diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp index 611904537f0..bbe9b60fc98 100644 --- a/Core/src/Geometry/TrackingVolume.cpp +++ b/Core/src/Geometry/TrackingVolume.cpp @@ -69,8 +69,8 @@ TrackingVolume::TrackingVolume(const Transform3& transform, {}, volumeName) {} TrackingVolume::~TrackingVolume() = default; -TrackingVolume::TrackingVolume(TrackingVolume&&) = default; -TrackingVolume& TrackingVolume::operator=(TrackingVolume&&) = default; +TrackingVolume::TrackingVolume(TrackingVolume&&) noexcept = default; +TrackingVolume& TrackingVolume::operator=(TrackingVolume&&) noexcept = default; const TrackingVolume* TrackingVolume::lowestTrackingVolume( const GeometryContext& gctx, const Vector3& position, diff --git a/Core/src/MagneticField/BFieldMapUtils.cpp b/Core/src/MagneticField/BFieldMapUtils.cpp index e7216aa3bf9..b5b92191589 100644 --- a/Core/src/MagneticField/BFieldMapUtils.cpp +++ b/Core/src/MagneticField/BFieldMapUtils.cpp @@ -58,8 +58,8 @@ Acts::fieldMapRZ( double zMax = zPos[nBinsZ - 1]; // calculate maxima (add one last bin, because bin value always corresponds to // left boundary) - double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1); - double stepR = std::fabs(rMax - rMin) / (nBinsR - 1); + double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); + double stepR = std::abs(rMax - rMin) / (nBinsR - 1); rMax += stepR; zMax += stepZ; if (firstQuadrant) { @@ -172,9 +172,9 @@ Acts::fieldMapXYZ( double zMax = zPos[nBinsZ - 1]; // calculate maxima (add one last bin, because bin value always corresponds to // left boundary) - double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1); - double stepY = std::fabs(yMax - yMin) / (nBinsY - 1); - double stepX = std::fabs(xMax - xMin) / (nBinsX - 1); + double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); + double stepY = std::abs(yMax - yMin) / (nBinsY - 1); + double stepX = std::abs(xMax - xMin) / (nBinsX - 1); xMax += stepX; yMax += stepY; zMax += stepZ; diff --git a/Core/src/Material/MaterialGridHelper.cpp b/Core/src/Material/MaterialGridHelper.cpp index a6d1a90416a..83c7946f77c 100644 --- a/Core/src/Material/MaterialGridHelper.cpp +++ b/Core/src/Material/MaterialGridHelper.cpp @@ -32,8 +32,8 @@ Acts::Grid2D Acts::createGrid(Acts::MaterialGridAxisData gridAxis1, // calculate maxima (add one last bin, because bin value always corresponds // to // left boundary) - double stepAxis1 = std::fabs(maxAxis1 - minAxis1) / (nBinsAxis1 - 1); - double stepAxis2 = std::fabs(maxAxis2 - minAxis2) / (nBinsAxis2 - 1); + double stepAxis1 = std::abs(maxAxis1 - minAxis1) / (nBinsAxis1 - 1); + double stepAxis2 = std::abs(maxAxis2 - minAxis2) / (nBinsAxis2 - 1); maxAxis1 += stepAxis1; maxAxis2 += stepAxis2; @@ -64,11 +64,11 @@ Acts::Grid3D Acts::createGrid(Acts::MaterialGridAxisData gridAxis1, // to // left boundary) double stepAxis1 = - std::fabs(maxAxis1 - minAxis1) / std::max(nBinsAxis1 - 1, std::size_t{1}); + std::abs(maxAxis1 - minAxis1) / std::max(nBinsAxis1 - 1, std::size_t{1}); double stepAxis2 = - std::fabs(maxAxis2 - minAxis2) / std::max(nBinsAxis2 - 1, std::size_t{1}); + std::abs(maxAxis2 - minAxis2) / std::max(nBinsAxis2 - 1, std::size_t{1}); double stepAxis3 = - std::fabs(maxAxis3 - minAxis3) / std::max(nBinsAxis3 - 1, std::size_t{1}); + std::abs(maxAxis3 - minAxis3) / std::max(nBinsAxis3 - 1, std::size_t{1}); maxAxis1 += stepAxis1; maxAxis2 += stepAxis2; maxAxis3 += stepAxis3; diff --git a/Core/src/Material/MaterialMapUtils.cpp b/Core/src/Material/MaterialMapUtils.cpp index a9b04275925..b51b1a3e223 100644 --- a/Core/src/Material/MaterialMapUtils.cpp +++ b/Core/src/Material/MaterialMapUtils.cpp @@ -64,8 +64,8 @@ auto Acts::materialMapperRZ( double zMax = *minMaxZ.second; // calculate maxima (add one last bin, because bin value always corresponds to // left boundary) - double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1); - double stepR = std::fabs(rMax - rMin) / (nBinsR - 1); + double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); + double stepR = std::abs(rMax - rMin) / (nBinsR - 1); rMax += stepR; zMax += stepZ; @@ -156,9 +156,9 @@ auto Acts::materialMapperXYZ( double zMax = *minMaxZ.second; // calculate maxima (add one last bin, because bin value always corresponds to // left boundary) - double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1); - double stepY = std::fabs(yMax - yMin) / (nBinsY - 1); - double stepX = std::fabs(xMax - xMin) / (nBinsX - 1); + double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); + double stepY = std::abs(yMax - yMin) / (nBinsY - 1); + double stepX = std::abs(xMax - xMin) / (nBinsX - 1); xMax += stepX; yMax += stepY; zMax += stepZ; diff --git a/Core/src/Navigation/SurfaceArrayNavigationPolicy.cpp b/Core/src/Navigation/SurfaceArrayNavigationPolicy.cpp index f40621fdd06..f60bbbf6839 100644 --- a/Core/src/Navigation/SurfaceArrayNavigationPolicy.cpp +++ b/Core/src/Navigation/SurfaceArrayNavigationPolicy.cpp @@ -45,8 +45,6 @@ SurfaceArrayNavigationPolicy::SurfaceArrayNavigationPolicy( auto [binsPhi, binsZ] = config.bins; m_surfaceArray = sac.surfaceArrayOnCylinder(gctx, std::move(surfaces), binsPhi, binsZ); - // m_surfaces = sac.createCylinderSurfaces(config.bins.first, - // config.bins.second); } else if (config.layerType == LayerType::Plane) { ACTS_ERROR("Plane layers are not yet supported"); throw std::invalid_argument("Plane layers are not yet supported"); diff --git a/Core/src/Surfaces/CylinderBounds.cpp b/Core/src/Surfaces/CylinderBounds.cpp index 27123d29ec0..64d3e9f547f 100644 --- a/Core/src/Surfaces/CylinderBounds.cpp +++ b/Core/src/Surfaces/CylinderBounds.cpp @@ -69,8 +69,8 @@ bool Acts::CylinderBounds::inside( double localx = lposition[0] > radius ? 2 * radius - lposition[0] : lposition[0]; Vector2 shiftedlposition = shifted(lposition); - if ((std::fabs(shiftedlposition[0]) <= halfPhi && - std::fabs(shiftedlposition[1]) <= halfLengthZ)) { + if ((std::abs(shiftedlposition[0]) <= halfPhi && + std::abs(shiftedlposition[1]) <= halfLengthZ)) { return true; } diff --git a/Core/src/Surfaces/CylinderSurface.cpp b/Core/src/Surfaces/CylinderSurface.cpp index a39dce4a4fa..4b3cd2db7be 100644 --- a/Core/src/Surfaces/CylinderSurface.cpp +++ b/Core/src/Surfaces/CylinderSurface.cpp @@ -180,7 +180,7 @@ double Acts::CylinderSurface::pathCorrection( const Acts::Vector3& direction) const { Vector3 normalT = normal(gctx, position); double cosAlpha = normalT.dot(direction); - return std::fabs(1. / cosAlpha); + return std::abs(1. / cosAlpha); } const Acts::CylinderBounds& Acts::CylinderSurface::bounds() const { diff --git a/Core/src/TrackFinding/MeasurementSelector.cpp b/Core/src/TrackFinding/MeasurementSelector.cpp index 130c8381536..f3de466f5df 100644 --- a/Core/src/TrackFinding/MeasurementSelector.cpp +++ b/Core/src/TrackFinding/MeasurementSelector.cpp @@ -29,6 +29,11 @@ MeasurementSelector::MeasurementSelector(const MeasurementSelectorCuts& cuts) : MeasurementSelector({{GeometryIdentifier(), cuts}}) {} MeasurementSelector::MeasurementSelector(const Config& config) { + if (config.empty()) { + throw std::invalid_argument( + "MeasurementSelector: Configuration must not be empty"); + } + std::vector tmp; tmp.reserve(config.size()); for (std::size_t i = 0; i < config.size(); ++i) { diff --git a/Core/src/Utilities/SpacePointUtility.cpp b/Core/src/Utilities/SpacePointUtility.cpp index 4b79fc03cc3..7018db37830 100644 --- a/Core/src/Utilities/SpacePointUtility.cpp +++ b/Core/src/Utilities/SpacePointUtility.cpp @@ -204,8 +204,8 @@ Result SpacePointUtility::calculateStripSPPosition( } // Check if m and n can be resolved in the interval (-1, 1) - if (fabs(spParams.m) <= spParams.limit && - fabs(spParams.n) <= spParams.limit) { + if (std::abs(spParams.m) <= spParams.limit && + std::abs(spParams.n) <= spParams.limit) { return Result::success(); } return Result::failure(m_error); @@ -226,7 +226,7 @@ Result SpacePointUtility::recoverSpacePoint( spParams.limit + stripLengthGapTolerance / spParams.mag_firstBtmToTop; // Check if m is just slightly outside - if (fabs(spParams.m) > spParams.limitExtended) { + if (std::abs(spParams.m) > spParams.limitExtended) { return Result::failure(m_error); } // Calculate n if not performed previously @@ -236,7 +236,7 @@ Result SpacePointUtility::recoverSpacePoint( spParams.secondBtmToTop.dot(spParams.firstBtmToTopXvtxToFirstMid2); } // Check if n is just slightly outside - if (fabs(spParams.n) > spParams.limitExtended) { + if (std::abs(spParams.n) > spParams.limitExtended) { return Result::failure(m_error); } /// The following code considers an overshoot of m and n in the same direction @@ -275,8 +275,8 @@ Result SpacePointUtility::recoverSpacePoint( spParams.n -= (biggerOvershoot / secOnFirstScale); // Check if this recovered the space point - if (fabs(spParams.m) < spParams.limit && - fabs(spParams.n) < spParams.limit) { + if (std::abs(spParams.m) < spParams.limit && + std::abs(spParams.n) < spParams.limit) { return Result::success(); } else { return Result::failure(m_error); @@ -294,8 +294,8 @@ Result SpacePointUtility::recoverSpacePoint( spParams.m += biggerOvershoot; spParams.n += (biggerOvershoot / secOnFirstScale); // Check if this recovered the space point - if (fabs(spParams.m) < spParams.limit && - fabs(spParams.n) < spParams.limit) { + if (std::abs(spParams.m) < spParams.limit && + std::abs(spParams.n) < spParams.limit) { return Result::success(); } } @@ -325,7 +325,7 @@ Result SpacePointUtility::calcPerpendicularProjection( double qr = (spParams.firstBtmToTop).dot(spParams.secondBtmToTop); double denom = spParams.firstBtmToTop.dot(spParams.firstBtmToTop) - qr * qr; // Check for numerical stability - if (fabs(denom) > 1e-6) { + if (std::abs(denom) > 1e-6) { // Return lambda0 return Result::success( (ac.dot(spParams.secondBtmToTop) * qr - diff --git a/Examples/Algorithms/Digitization/CMakeLists.txt b/Examples/Algorithms/Digitization/CMakeLists.txt index 14729f2784f..de98c7e74d4 100644 --- a/Examples/Algorithms/Digitization/CMakeLists.txt +++ b/Examples/Algorithms/Digitization/CMakeLists.txt @@ -3,6 +3,7 @@ add_library( SHARED src/DigitizationAlgorithm.cpp src/DigitizationConfig.cpp + src/DigitizationCoordinatesConverter.cpp src/MeasurementCreation.cpp src/DigitizationConfigurator.cpp src/ModuleClusters.cpp diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationCoordinatesConverter.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationCoordinatesConverter.hpp new file mode 100644 index 00000000000..21c4f619b8a --- /dev/null +++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationCoordinatesConverter.hpp @@ -0,0 +1,39 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "ActsExamples/Digitization/DigitizationConfig.hpp" + +namespace ActsExamples { + +/// Utility that converts hit coordinates. +class DigitizationCoordinatesConverter final { + public: + /// Construct the converter + /// + /// @param config is the configuration + explicit DigitizationCoordinatesConverter(DigitizationConfig config); + + /// Get const access to the config + const DigitizationConfig& config() const { return m_cfg; } + + /// Convert the hit coordinates to the local frame. + std::tuple globalToLocal(std::uint64_t moduleId, double x, + double y, double z) const; + + /// Convert the hit coordinates to the global frame. + std::tuple localToGlobal(std::uint64_t moduleId, + double x, double y) const; + + private: + /// Configuration + DigitizationConfig m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/Digitization/src/DigitizationCoordinatesConverter.cpp b/Examples/Algorithms/Digitization/src/DigitizationCoordinatesConverter.cpp new file mode 100644 index 00000000000..e83e7be995b --- /dev/null +++ b/Examples/Algorithms/Digitization/src/DigitizationCoordinatesConverter.cpp @@ -0,0 +1,63 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include "ActsExamples/Digitization/DigitizationCoordinatesConverter.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/GeometryContext.hpp" + +#include +#include +#include + +ActsExamples::DigitizationCoordinatesConverter:: + DigitizationCoordinatesConverter(DigitizationConfig config) + : m_cfg(std::move(config)) { + if (m_cfg.surfaceByIdentifier.empty()) { + throw std::invalid_argument("Missing Surface-GeometryID association map"); + } +} + +std::tuple +ActsExamples::DigitizationCoordinatesConverter::globalToLocal( + std::uint64_t moduleId, double x, double y, double z) const { + const Acts::GeometryIdentifier moduleGeoId = moduleId; + auto surfaceItr = m_cfg.surfaceByIdentifier.find(moduleGeoId); + if (surfaceItr == m_cfg.surfaceByIdentifier.end()) { + throw std::runtime_error("Surface not found for moduleGeoId"); + } + + // Transform the position into the local surface frame + const Acts::Surface* surfacePtr = surfaceItr->second; + const auto& invTransform = + surfacePtr->transform(Acts::GeometryContext()).inverse(); + + const Acts::Vector3 pos(x, y, z); + Acts::Vector2 pos2Local = (invTransform * pos).segment<2>(0); + + return {pos2Local.x(), pos2Local.y()}; +} + +std::tuple +ActsExamples::DigitizationCoordinatesConverter::localToGlobal( + std::uint64_t moduleId, double x, double y) const { + const Acts::GeometryIdentifier moduleGeoId = moduleId; + auto surfaceItr = m_cfg.surfaceByIdentifier.find(moduleGeoId); + if (surfaceItr == m_cfg.surfaceByIdentifier.end()) { + throw std::runtime_error("Surface not found for moduleGeoId"); + } + + // Transform the position into the global frame + const Acts::Surface* surfacePtr = surfaceItr->second; + const auto& transform = surfacePtr->transform(Acts::GeometryContext()); + + const Acts::Vector3 pos(x, y, 0); + Acts::Vector3 pos2Global = (transform * pos); + + return {pos2Global.x(), pos2Global.y(), pos2Global.z()}; +} diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp index dcc7ac573f3..240669274b0 100644 --- a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp +++ b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp @@ -9,15 +9,12 @@ #include "ActsExamples/Generators/ParametricParticleGenerator.hpp" #include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/ParticleData.hpp" #include "Acts/Utilities/AngleHelpers.hpp" #include "ActsFatras/EventData/Barcode.hpp" #include "ActsFatras/EventData/Particle.hpp" -#include #include -#include #include namespace ActsExamples { @@ -25,41 +22,71 @@ namespace ActsExamples { ParametricParticleGenerator::ParametricParticleGenerator(const Config& cfg) : m_cfg(cfg), m_charge(cfg.charge.value_or(Acts::findCharge(m_cfg.pdg).value_or(0))), - m_mass(cfg.mass.value_or(Acts::findMass(m_cfg.pdg).value_or(0))), - // since we want to draw the direction uniform on the unit sphere, we must - // draw from cos(theta) instead of theta. see e.g. - // https://mathworld.wolfram.com/SpherePointPicking.html - m_cosThetaMin(std::cos(m_cfg.thetaMin)), - // ensure upper bound is included. see e.g. - // https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution - m_cosThetaMax(std::nextafter(std::cos(m_cfg.thetaMax), - std::numeric_limits::max())), - // in case we force uniform eta generation - m_etaMin(Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMin)), - m_etaMax(Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMax)) {} - -std::pair -ParametricParticleGenerator::operator()(RandomEngine& rng) { - using UniformIndex = std::uniform_int_distribution; - using UniformReal = std::uniform_real_distribution; - - // choose between particle/anti-particle if requested - // the upper limit of the distribution is inclusive - UniformIndex particleTypeChoice(0u, m_cfg.randomizeCharge ? 1u : 0u); - // (anti-)particle choice is one random draw but defines two properties - const Acts::PdgParticle pdgChoices[] = { + m_mass(cfg.mass.value_or(Acts::findMass(m_cfg.pdg).value_or(0))) { + m_pdgChoices = { m_cfg.pdg, static_cast(-m_cfg.pdg), }; - const double qChoices[] = { + m_qChoices = { m_charge, -m_charge, }; - UniformReal phiDist(m_cfg.phiMin, m_cfg.phiMax); - UniformReal cosThetaDist(m_cosThetaMin, m_cosThetaMax); - UniformReal etaDist(m_etaMin, m_etaMax); - UniformReal pDist(m_cfg.pMin, m_cfg.pMax); + // choose between particle/anti-particle if requested + // the upper limit of the distribution is inclusive + m_particleTypeChoice = UniformIndex(0u, m_cfg.randomizeCharge ? 1u : 0u); + m_phiDist = UniformReal(m_cfg.phiMin, m_cfg.phiMax); + + if (m_cfg.etaUniform) { + double etaMin = Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMin); + double etaMax = Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMax); + + UniformReal etaDist(etaMin, etaMax); + + m_sinCosThetaDist = + [=](RandomEngine& rng) mutable -> std::pair { + const double eta = etaDist(rng); + const double theta = Acts::AngleHelpers::thetaFromEta(eta); + return {std::sin(theta), std::cos(theta)}; + }; + } else { + // since we want to draw the direction uniform on the unit sphere, we must + // draw from cos(theta) instead of theta. see e.g. + // https://mathworld.wolfram.com/SpherePointPicking.html + double cosThetaMin = std::cos(m_cfg.thetaMin); + // ensure upper bound is included. see e.g. + // https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution + double cosThetaMax = std::nextafter(std::cos(m_cfg.thetaMax), + std::numeric_limits::max()); + + UniformReal cosThetaDist(cosThetaMin, cosThetaMax); + + m_sinCosThetaDist = + [=](RandomEngine& rng) mutable -> std::pair { + const double cosTheta = cosThetaDist(rng); + return {std::sqrt(1 - cosTheta * cosTheta), cosTheta}; + }; + } + + if (m_cfg.pLogUniform) { + // distributes p or pt uniformly in log space + UniformReal dist(std::log(m_cfg.pMin), std::log(m_cfg.pMax)); + + m_somePDist = [=](RandomEngine& rng) mutable -> double { + return std::exp(dist(rng)); + }; + } else { + // distributes p or pt uniformly + UniformReal dist(m_cfg.pMin, m_cfg.pMax); + + m_somePDist = [=](RandomEngine& rng) mutable -> double { + return dist(rng); + }; + } +} + +std::pair +ParametricParticleGenerator::operator()(RandomEngine& rng) { SimVertexContainer::sequence_type vertices; SimParticleContainer::sequence_type particles; @@ -74,33 +101,22 @@ ParametricParticleGenerator::operator()(RandomEngine& rng) { primaryVertex.outgoing.insert(pid); // draw parameters - const unsigned int type = particleTypeChoice(rng); - const Acts::PdgParticle pdg = pdgChoices[type]; - const double q = qChoices[type]; - const double phi = phiDist(rng); - double p = pDist(rng); - - // we already have sin/cos theta. they can be used directly to - Acts::Vector3 dir; - double cosTheta = 0.; - double sinTheta = 0.; - if (!m_cfg.etaUniform) { - cosTheta = cosThetaDist(rng); - sinTheta = std::sqrt(1 - cosTheta * cosTheta); - } else { - const double eta = etaDist(rng); - const double theta = Acts::AngleHelpers::thetaFromEta(eta); - sinTheta = std::sin(theta); - cosTheta = std::cos(theta); - } - dir[Acts::eMom0] = sinTheta * std::cos(phi); - dir[Acts::eMom1] = sinTheta * std::sin(phi); - dir[Acts::eMom2] = cosTheta; + const unsigned int type = m_particleTypeChoice(rng); + const Acts::PdgParticle pdg = m_pdgChoices[type]; + const double q = m_qChoices[type]; + const double phi = m_phiDist(rng); + const double someP = m_somePDist(rng); + + const auto [sinTheta, cosTheta] = m_sinCosThetaDist(rng); + // we already have sin/cos theta. they can be used directly + const Acts::Vector3 dir = {sinTheta * std::cos(phi), + sinTheta * std::sin(phi), cosTheta}; + + const double p = someP * (m_cfg.pTransverse ? 1. / sinTheta : 1.); // construct the particle; ActsFatras::Particle particle(pid, pdg, q, m_mass); particle.setDirection(dir); - p *= m_cfg.pTransverse ? 1. / sinTheta : 1.; particle.setAbsoluteMomentum(p); // generated particle ids are already ordered and should end up at the end diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp index a7d9bfe3e08..0efc47a0c5c 100644 --- a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp +++ b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp @@ -15,10 +15,10 @@ #include "ActsExamples/Generators/EventGenerator.hpp" #include -#include #include #include #include +#include namespace ActsExamples { @@ -49,8 +49,10 @@ class ParametricParticleGenerator : public EventGenerator::ParticlesGenerator { /// Low, high (exclusive) for absolute/transverse momentum. double pMin = 1 * Acts::UnitConstants::GeV; double pMax = 10 * Acts::UnitConstants::GeV; - /// Indicate if the momentum referse to transverse momentum + /// Indicate if the momentum refers to transverse momentum. bool pTransverse = false; + /// Indicate if the momentum should be uniformly distributed in log space. + bool pLogUniform = false; /// (Absolute) PDG particle number to identify the particle type. Acts::PdgParticle pdg = Acts::PdgParticle::eMuon; /// Randomize the charge and flip the PDG particle number sign accordingly. @@ -71,14 +73,23 @@ class ParametricParticleGenerator : public EventGenerator::ParticlesGenerator { RandomEngine& rng) override; private: + using UniformIndex = std::uniform_int_distribution; + using UniformReal = std::uniform_real_distribution; + Config m_cfg; + // will be automatically set from PDG data tables - double m_charge; - double m_mass; - double m_cosThetaMin; - double m_cosThetaMax; - double m_etaMin; - double m_etaMax; + double m_charge{}; + double m_mass{}; + + // (anti-)particle choice is one random draw but defines two properties + std::array m_pdgChoices{}; + std::array m_qChoices{}; + + UniformIndex m_particleTypeChoice; + UniformReal m_phiDist; + std::function(RandomEngine& rng)> m_sinCosThetaDist; + std::function m_somePDist; }; } // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp index 527860ee556..963c179476b 100644 --- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp @@ -70,7 +70,7 @@ class SeedingOrthogonalAlgorithm final : public IAlgorithm { private: Config m_cfg; - Acts::SeedFinderOrthogonal m_finder; + std::unique_ptr> m_finder{nullptr}; std::vector>> m_inputSpacePoints{}; diff --git a/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp index 385a3db7ec1..48f1490946e 100644 --- a/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp @@ -88,20 +88,20 @@ ActsExamples::ProcessCode ActsExamples::GbtsSeedingAlgorithm::execute( MakeGbtsSpacePoints(ctx, m_cfg.ActsGbtsMap); for (auto sp : GbtsSpacePoints) { - ACTS_DEBUG("Gbts space points: " << " Gbts_id: " << sp.gbtsID << " z: " - << sp.SP->z() << " r: " << sp.SP->r() - << " ACTS volume: " - << sp.SP->sourceLinks() - .front() - .get() - .geometryId() - .volume() - << "\n"); + ACTS_DEBUG("Gbts space points: Gbts_id: " + << sp.gbtsID << " z: " << sp.SP->z() << " r: " << sp.SP->r() + << " ACTS volume: " + << sp.SP->sourceLinks() + .front() + .get() + .geometryId() + .volume()); } // this is now calling on a core algorithm - Acts::SeedFinderGbts finder(m_cfg.seedFinderConfig, - *m_gbtsGeo); + Acts::SeedFinderGbts finder( + m_cfg.seedFinderConfig, *m_gbtsGeo, + logger().cloneWithSuffix("GbtdFinder")); // output of function needed for seed diff --git a/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp index ed5b11e4df2..2cea4b7b505 100644 --- a/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp @@ -45,9 +45,8 @@ ActsExamples::SeedingAlgorithm::SeedingAlgorithm( // internal units m_cfg.seedFilterConfig = m_cfg.seedFilterConfig.toInternalUnits(); m_cfg.seedFinderConfig.seedFilter = - std::make_shared>( - m_cfg.seedFilterConfig); - + std::make_unique>( + m_cfg.seedFilterConfig, logger().cloneWithSuffix("SeedFilter")); m_cfg.seedFinderConfig = m_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities(); m_cfg.seedFinderOptions = @@ -195,13 +194,10 @@ ActsExamples::SeedingAlgorithm::SeedingAlgorithm( m_topBinFinder = std::make_unique>( m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsTop, 0); - m_cfg.seedFinderConfig.seedFilter = - std::make_unique>( - m_cfg.seedFilterConfig); m_seedFinder = Acts::SeedFinder>( - m_cfg.seedFinderConfig); + m_cfg.seedFinderConfig, logger().cloneWithSuffix("SeedFinder")); } ActsExamples::ProcessCode ActsExamples::SeedingAlgorithm::execute( @@ -244,10 +240,11 @@ ActsExamples::ProcessCode ActsExamples::SeedingAlgorithm::execute( Acts::CylindricalSpacePointGrid grid = Acts::CylindricalSpacePointGridCreator::createGrid( - m_cfg.gridConfig, m_cfg.gridOptions); + m_cfg.gridConfig, m_cfg.gridOptions, logger()); Acts::CylindricalSpacePointGridCreator::fillGrid( - m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid, spContainer); + m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid, spContainer, + logger()); // Compute radius Range // we rely on the fact the grid is storing the proxies diff --git a/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp index dcfd995f4e9..6a1673fe6f6 100644 --- a/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp @@ -67,9 +67,11 @@ ActsExamples::SeedingOrthogonalAlgorithm::SeedingOrthogonalAlgorithm( // construct seed filter m_cfg.seedFinderConfig.seedFilter = - std::make_unique>(m_cfg.seedFilterConfig); + std::make_unique>( + m_cfg.seedFilterConfig, logger().cloneWithSuffix("Filter")); - m_finder = Acts::SeedFinderOrthogonal(m_cfg.seedFinderConfig); + m_finder = std::make_unique>( + m_cfg.seedFinderConfig, logger().cloneWithSuffix("Finder")); } ActsExamples::ProcessCode ActsExamples::SeedingOrthogonalAlgorithm::execute( @@ -95,9 +97,8 @@ ActsExamples::ProcessCode ActsExamples::SeedingOrthogonalAlgorithm::execute( ACTS_INFO("About to process " << spContainer.size() << " space points ..."); - Acts::SeedFinderOrthogonal finder(m_cfg.seedFinderConfig); std::vector> seeds = - finder.createSeeds(m_cfg.seedFinderOptions, spContainer); + m_finder->createSeeds(m_cfg.seedFinderOptions, spContainer); ACTS_INFO("Created " << seeds.size() << " track seeds from " << spacePoints.size() << " space points"); diff --git a/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp index 6a5d2e0059b..6f835c5c967 100644 --- a/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp @@ -62,7 +62,7 @@ struct SimpleReverseFilteringLogic { bool doBackwardFiltering( Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const { - auto momentum = fabs(1 / trackState.filtered()[Acts::eBoundQOverP]); + auto momentum = std::abs(1 / trackState.filtered()[Acts::eBoundQOverP]); return (momentum <= momentumThreshold); } }; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedSelector.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedSelector.cpp deleted file mode 100644 index 49640e5079f..00000000000 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedSelector.cpp +++ /dev/null @@ -1,89 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#include "ActsExamples/TruthTracking/TruthSeedSelector.hpp" - -#include "Acts/Utilities/MultiIndex.hpp" -#include "Acts/Utilities/VectorHelpers.hpp" -#include "ActsExamples/EventData/Index.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/Utilities/Range.hpp" -#include "ActsFatras/EventData/Barcode.hpp" -#include "ActsFatras/EventData/Particle.hpp" - -#include -#include -#include - -namespace ActsExamples { -struct AlgorithmContext; -} // namespace ActsExamples - -using namespace ActsExamples; - -TruthSeedSelector::TruthSeedSelector(const Config& config, - Acts::Logging::Level level) - : IAlgorithm("TruthSeedSelector", level), m_cfg(config) { - if (m_cfg.inputParticles.empty()) { - throw std::invalid_argument("Missing input truth particles collection"); - } - if (m_cfg.inputMeasurementParticlesMap.empty()) { - throw std::invalid_argument("Missing input hit-particles map collection"); - } - if (m_cfg.outputParticles.empty()) { - throw std::invalid_argument("Missing output truth particles collection"); - } - - m_inputParticles.initialize(m_cfg.inputParticles); - m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap); - m_outputParticles.initialize(m_cfg.outputParticles); -} - -ProcessCode TruthSeedSelector::execute(const AlgorithmContext& ctx) const { - // prepare input collections - const auto& inputParticles = m_inputParticles(ctx); - const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx); - // compute particle_id -> {hit_id...} map from the - // hit_id -> {particle_id...} map on the fly. - const auto& particleHitsMap = invertIndexMultimap(hitParticlesMap); - - // prepare output collection - SimParticleContainer selectedParticles; - selectedParticles.reserve(inputParticles.size()); - - auto within = [](double x, double min, double max) { - return (min <= x) && (x < max); - }; - auto isValidparticle = [&](const auto& p) { - const auto eta = Acts::VectorHelpers::eta(p.direction()); - const auto phi = Acts::VectorHelpers::phi(p.direction()); - const auto rho = Acts::VectorHelpers::perp(p.position()); - // find the corresponding hits for this particle - const auto& hits = makeRange(particleHitsMap.equal_range(p.particleId())); - // number of recorded hits - std::size_t nHits = hits.size(); - return within(rho, 0., m_cfg.rhoMax) && - within(p.position().z(), m_cfg.zMin, m_cfg.zMax) && - within(std::abs(eta), m_cfg.absEtaMin, m_cfg.absEtaMax) && - within(eta, m_cfg.etaMin, m_cfg.etaMax) && - within(phi, m_cfg.phiMin, m_cfg.phiMax) && - within(p.transverseMomentum(), m_cfg.ptMin, m_cfg.ptMax) && - within(nHits, m_cfg.nHitsMin, m_cfg.nHitsMax) && - (m_cfg.keepNeutral || (p.charge() != 0)); - }; - - // create prototracks for all input particles - for (const auto& particle : inputParticles) { - if (isValidparticle(particle)) { - selectedParticles.insert(particle); - } - } - - m_outputParticles(ctx, std::move(selectedParticles)); - return ProcessCode::SUCCESS; -} diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedSelector.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedSelector.hpp deleted file mode 100644 index 1bb12042aa3..00000000000 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedSelector.hpp +++ /dev/null @@ -1,92 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/SimHit.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/Framework/DataHandle.hpp" -#include "ActsExamples/Framework/IAlgorithm.hpp" -#include "ActsExamples/Framework/ProcessCode.hpp" - -#include -#include -#include - -namespace ActsFatras { -class Barcode; -} // namespace ActsFatras - -namespace ActsExamples { -struct AlgorithmContext; - -/// Select truth particles to be used as 'seeds' of reconstruction algorithms, -/// e.g. track fitting and track finding. -/// -/// This pre-selection could help guarantee quality of the 'seeds', i.e. to -/// avoid empty proto track (no recorded hits for the particle). In addition, it -/// could help save unnecessary reconstruction time. For instance, when -/// investigating performance of CombinatorialKalmanFilter (CKF), we might be -/// interested in its performance for only truth particles with pT and number of -/// recorded hits (on sensitive detectors) safistying provided criteria (input -/// measurements of CKF are still recorded hits from all possible particles). -/// Then we could use particles only satisfying provided criteria as the 'seeds' -/// of CKF instead of handling all the truth particles. -// -class TruthSeedSelector final : public IAlgorithm { - public: - struct Config { - /// The input truth particles that should be used to create proto tracks. - std::string inputParticles; - /// The input hit-particles map collection. - std::string inputMeasurementParticlesMap; - /// The output proto tracks collection. - std::string outputParticles; - /// Maximum distance from the origin in the transverse plane - double rhoMin = 0.; - double rhoMax = std::numeric_limits::max(); - /// Minimum/Maximum absolute distance from the origin along z - double zMin = std::numeric_limits::lowest(); - double zMax = std::numeric_limits::max(); - // Truth particle kinematic cuts - double phiMin = std::numeric_limits::lowest(); - double phiMax = std::numeric_limits::max(); - double etaMin = std::numeric_limits::lowest(); - double etaMax = std::numeric_limits::max(); - double absEtaMin = std::numeric_limits::lowest(); - double absEtaMax = std::numeric_limits::max(); - double ptMin = 0.0; - double ptMax = std::numeric_limits::max(); - /// Keep neutral particles - bool keepNeutral = false; - /// Requirement on number of recorded hits - //@TODO: implement detector-specific requirements - std::size_t nHitsMin = 0; - std::size_t nHitsMax = std::numeric_limits::max(); - }; - - TruthSeedSelector(const Config& config, Acts::Logging::Level level); - - ProcessCode execute(const AlgorithmContext& ctx) const final; - - /// Get readonly access to the config parameters - const Config& config() const { return m_cfg; } - - private: - Config m_cfg; - - ReadDataHandle m_inputParticles{this, "InputParticles"}; - ReadDataHandle m_inputMeasurementParticlesMap{ - this, "InputMeasurementParticlesMap"}; - - WriteDataHandle m_outputParticles{this, - "OutputParticles"}; -}; - -} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/CMakeLists.txt b/Examples/Algorithms/TruthTracking/CMakeLists.txt index 7c1eff1ed51..4594cdb7480 100644 --- a/Examples/Algorithms/TruthTracking/CMakeLists.txt +++ b/Examples/Algorithms/TruthTracking/CMakeLists.txt @@ -6,7 +6,6 @@ add_library( ActsExamples/TruthTracking/TrackParameterSelector.cpp ActsExamples/TruthTracking/TrackModifier.cpp ActsExamples/TruthTracking/TrackTruthMatcher.cpp - ActsExamples/TruthTracking/TruthSeedSelector.cpp ActsExamples/TruthTracking/TruthTrackFinder.cpp ActsExamples/TruthTracking/TruthVertexFinder.cpp ActsExamples/TruthTracking/TruthSeedingAlgorithm.cpp diff --git a/Examples/Framework/CMakeLists.txt b/Examples/Framework/CMakeLists.txt index efe06302800..db45cffc925 100644 --- a/Examples/Framework/CMakeLists.txt +++ b/Examples/Framework/CMakeLists.txt @@ -46,30 +46,15 @@ target_compile_definitions( PRIVATE BOOST_FILESYSTEM_NO_DEPRECATED ) -if(ACTS_USE_EXAMPLES_TBB) - # newer DD4hep version require TBB and search internally for TBB in - # config-only mode. to avoid mismatches we explicitly search using - # config-only mode first to be sure that we find the same version. - find_package(TBB ${_acts_tbb_version} CONFIG) - if(NOT TBB_FOUND) - # no version check possible when using the find module - find_package(TBB ${_acts_tbb_version} MODULE) - endif() -else() - set(TBB_FOUND FALSE) -endif() -if(TBB_FOUND) - target_link_libraries(ActsExamplesFramework PUBLIC TBB::tbb) -else() - message( - STATUS - "disable TBB for Examples/Framework - only single-threaded running will be supported" - ) - target_compile_definitions( - ActsExamplesFramework - PUBLIC -DACTS_EXAMPLES_NO_TBB - ) +# newer DD4hep version require TBB and search internally for TBB in +# config-only mode. to avoid mismatches we explicitly search using +# config-only mode first to be sure that we find the same version. +find_package(TBB ${_acts_tbb_version} CONFIG) +if(NOT TBB_FOUND) + # no version check possible when using the find module + find_package(TBB ${_acts_tbb_version} MODULE REQUIRED) endif() +target_link_libraries(ActsExamplesFramework PUBLIC TBB::tbb) install( TARGETS ActsExamplesFramework diff --git a/Examples/Framework/include/ActsExamples/Utilities/tbbWrap.hpp b/Examples/Framework/include/ActsExamples/Utilities/tbbWrap.hpp index 84784c55d31..79969369221 100644 --- a/Examples/Framework/include/ActsExamples/Utilities/tbbWrap.hpp +++ b/Examples/Framework/include/ActsExamples/Utilities/tbbWrap.hpp @@ -8,20 +8,11 @@ #pragma once -// uncomment to remove all use of tbb library. -// #define ACTS_EXAMPLES_NO_TBB - -#ifdef ACTS_EXAMPLES_NO_TBB -#define ACTS_EXAMPLES_WITH_TBB(a) -#include -#else -#define ACTS_EXAMPLES_WITH_TBB(a) a #include #include #include #include -#endif /// Wrapper for most of the tbb functions that we use in Sequencer. /// @@ -30,34 +21,9 @@ /// tbb::blocked_range (which doesn't require any thread setup) is still taken /// from the tbb library. /// -/// However, if ACTS_EXAMPLES_NO_TBB is defined, then don't use tbb library at -/// all (requires nthreads=1 or -1). This allows the ACTS Examples to be built -/// without the tbb library (and reduces the dependency on ROOT). -/// In this case, we provide our own minimal implementation of -/// tbb::blocked_range. -/// /// Based on an idea from /// https://stackoverflow.com/questions/59736661/how-to-completely-switch-off-threading-in-tbb-code -#ifdef ACTS_EXAMPLES_NO_TBB -namespace ActsExamples::tbb { -namespace task_arena { -constexpr int automatic = -1; -} // namespace task_arena - -template -struct blocked_range { - blocked_range(Value begin_, Value end_) : my_end(end_), my_begin(begin_) {} - Value begin() const { return my_begin; } - Value end() const { return my_end; } - - private: - Value my_end; - Value my_begin; -}; -} // namespace ActsExamples::tbb -#endif - namespace ActsExamples::tbbWrap { /// enableTBB keeps a record of whether we are multi-threaded (nthreads!=1) or /// not. This is set once in task_arena and stored globally. @@ -67,17 +33,10 @@ namespace ActsExamples::tbbWrap { static bool enableTBB(int nthreads = -99) { static bool setting = false; if (nthreads != -99) { -#ifdef ACTS_EXAMPLES_NO_TBB - if (nthreads > 1) { - throw std::runtime_error( - "tbb is not available, so can't do multi-threading."); - } -#else bool newSetting = (nthreads != 1); if (!setting && newSetting) { setting = newSetting; } -#endif } return setting; } @@ -87,28 +46,20 @@ static bool enableTBB(int nthreads = -99) { /// That should be fine because the task_arena is initialised before spawning /// any threads. class task_arena { -#ifndef ACTS_EXAMPLES_NO_TBB std::optional tbb; -#endif public: - task_arena(int nthreads = tbb::task_arena::automatic, - unsigned ACTS_EXAMPLES_WITH_TBB(res) = 1) { + task_arena(int nthreads = tbb::task_arena::automatic, unsigned res = 1) { if (enableTBB(nthreads)) { -#ifndef ACTS_EXAMPLES_NO_TBB tbb.emplace(nthreads, res); -#endif } } template void execute(const F& f) { -#ifndef ACTS_EXAMPLES_NO_TBB if (tbb) { tbb->execute(f); - } else -#endif - { + } else { f(); } } @@ -119,12 +70,9 @@ class parallel_for { public: template parallel_for(const R& r, const F& f) { -#ifndef ACTS_EXAMPLES_NO_TBB if (enableTBB()) { tbb::parallel_for(r, f); - } else -#endif - { + } else { for (auto i = r.begin(); i != r.end(); ++i) { // use default grainsize=1 f(R(i, i + 1)); } @@ -134,39 +82,29 @@ class parallel_for { /// Small wrapper for tbb::queuing_mutex and tbb::queuing_mutex::scoped_lock. class queuing_mutex { -#ifndef ACTS_EXAMPLES_NO_TBB std::optional tbb; -#endif public: queuing_mutex() { -#ifndef ACTS_EXAMPLES_NO_TBB if (enableTBB()) { tbb.emplace(); } -#endif } class scoped_lock { -#ifndef ACTS_EXAMPLES_NO_TBB std::optional tbb; -#endif public: scoped_lock() { -#ifndef ACTS_EXAMPLES_NO_TBB if (enableTBB()) { tbb.emplace(); } -#endif } - explicit scoped_lock(queuing_mutex& ACTS_EXAMPLES_WITH_TBB(m)) { -#ifndef ACTS_EXAMPLES_NO_TBB + explicit scoped_lock(queuing_mutex& m) { if (enableTBB()) { tbb.emplace(*m.tbb); } -#endif } }; }; diff --git a/Examples/Framework/src/Framework/Sequencer.cpp b/Examples/Framework/src/Framework/Sequencer.cpp index 5c1a76b5bd4..eb9bf8fd1c0 100644 --- a/Examples/Framework/src/Framework/Sequencer.cpp +++ b/Examples/Framework/src/Framework/Sequencer.cpp @@ -42,15 +42,11 @@ #include #include -#include - -#ifndef ACTS_EXAMPLES_NO_TBB #include -#endif - #include #include #include +#include namespace ActsExamples { @@ -108,16 +104,12 @@ Sequencer::Sequencer(const Sequencer::Config& cfg) m_taskArena((m_cfg.numThreads < 0) ? tbb::task_arena::automatic : m_cfg.numThreads), m_logger(Acts::getDefaultLogger("Sequencer", m_cfg.logLevel)) { -#ifndef ACTS_EXAMPLES_NO_TBB if (m_cfg.numThreads == 1) { -#endif ACTS_INFO("Create Sequencer (single-threaded)"); -#ifndef ACTS_EXAMPLES_NO_TBB } else { ROOT::EnableThreadSafety(); ACTS_INFO("Create Sequencer with " << m_cfg.numThreads << " threads"); } -#endif const char* envvar = std::getenv("ACTS_SEQUENCER_DISABLE_FPEMON"); if (envvar != nullptr) { @@ -267,10 +259,13 @@ void Sequencer::addWhiteboardAlias(const std::string& aliasName, auto [it, success] = m_whiteboardObjectAliases.insert({objectName, aliasName}); if (!success) { - throw std::invalid_argument("Alias to '" + aliasName + "' -> '" + - objectName + "' already set"); + ACTS_INFO("Key '" << objectName << "' aliased to '" << aliasName + << "' already set"); + return; } + ACTS_INFO("Key '" << objectName << "' aliased to '" << aliasName << "'"); + if (auto oit = m_whiteBoardState.find(objectName); oit != m_whiteBoardState.end()) { m_whiteBoardState[aliasName] = oit->second; diff --git a/Examples/Io/Csv/src/CsvSeedWriter.cpp b/Examples/Io/Csv/src/CsvSeedWriter.cpp index 96dd04b07a4..fe9d881c24f 100644 --- a/Examples/Io/Csv/src/CsvSeedWriter.cpp +++ b/Examples/Io/Csv/src/CsvSeedWriter.cpp @@ -123,10 +123,10 @@ ActsExamples::ProcessCode ActsExamples::CsvSeedWriter::writeT( // Compute the distance between the truth and estimated directions float truthPhi = phi(truthUnitDir); float truthEta = std::atanh(std::cos(theta(truthUnitDir))); - float dEta = fabs(truthEta - seedEta); - float dPhi = fabs(truthPhi - seedPhi) < M_PI - ? fabs(truthPhi - seedPhi) - : fabs(truthPhi - seedPhi) - M_PI; + float dEta = std::abs(truthEta - seedEta); + float dPhi = std::abs(truthPhi - seedPhi) < M_PI + ? std::abs(truthPhi - seedPhi) + : std::abs(truthPhi - seedPhi) - M_PI; truthDistance = sqrt(dPhi * dPhi + dEta * dEta); // If the seed is truth matched, check if it is the closest one for the // contributing particle diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index fd02c6bddac..bc1b1d4df9d 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -1403,8 +1403,14 @@ def addCKFTracks( TrackSelector configuration. Each range is specified as a tuple of (min,max). Specify as a list(TrackSelectorConfig) for eta-dependent cuts, with binning specified by absEta[1]. Defaults of no cuts specified in Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackSelector.hpp - writeTrajectories : bool, True - write trackstates_ckf.root and tracksummary_ckf.root ntuples? These can be quite large. + writeTrackSummary : bool, True + write tracksummary_ckf.root ntuple? + writeTrackStates : bool, False + write trackstates_ckf.root ntuple? This can be quite large. + writePerformance : bool, True + write performance_fitting_ckf.root and performance_finding_ckf.root ntuples? + writeCovMat : bool, False + write covaraiance matrices to tracksummary_ckf.root ntuple? """ customLogLevel = acts.examples.defaultLogging(s, logLevel) @@ -1741,21 +1747,6 @@ def addExaTrkX( ) -> None: customLogLevel = acts.examples.defaultLogging(s, logLevel) - # Run the particle selection - # The pre-selection will select truth particles satisfying provided criteria - # from all particles read in by particle reader for further processing. It - # has no impact on the truth hits themselves - s.addAlgorithm( - acts.examples.TruthSeedSelector( - level=customLogLevel(), - ptMin=500 * u.MeV, - nHitsMin=9, - inputParticles="particles_initial", - inputMeasurementParticlesMap="measurement_particles_map", - outputParticles="particles_seed_selected", - ) - ) - # Create space points s.addAlgorithm( acts.examples.SpacePointMaker( @@ -1917,8 +1908,8 @@ def addAmbiguityResolution( tracks=alg.config.outputTracks, outputDirCsv=outputDirCsv, outputDirRoot=outputDirRoot, - writeSummary=writeTrackStates, - writeStates=writeTrackSummary, + writeSummary=writeTrackSummary, + writeStates=writeTrackStates, writeFitterPerformance=writePerformance, writeFinderPerformance=writePerformance, writeCovMat=writeCovMat, diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py index 8f6b6d0e603..b72d1f456df 100644 --- a/Examples/Python/python/acts/examples/simulation.py +++ b/Examples/Python/python/acts/examples/simulation.py @@ -18,8 +18,8 @@ # Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp MomentumConfig = namedtuple( "MomentumConfig", - ["min", "max", "transverse"], - defaults=[None, None, None], + ["min", "max", "transverse", "logUniform"], + defaults=[None, None, None, None], ) EtaConfig = namedtuple( "EtaConfig", ["min", "max", "uniform"], defaults=[None, None, None] @@ -80,8 +80,8 @@ def addParticleGun( the output folder for the Csv output, None triggers no output outputDirRoot : Path|str, path, None the output folder for the Root output, None triggers no output - momentumConfig : MomentumConfig(min, max, transverse) - momentum configuration: minimum momentum, maximum momentum, transverse + momentumConfig : MomentumConfig(min, max, transverse, logUniform) + momentum configuration: minimum momentum, maximum momentum, transverse, log-uniform etaConfig : EtaConfig(min, max, uniform) pseudorapidity configuration: eta min, eta max, uniform phiConfig : PhiConfig(min, max) @@ -118,6 +118,7 @@ def addParticleGun( **acts.examples.defaultKWArgs( p=(momentumConfig.min, momentumConfig.max), pTransverse=momentumConfig.transverse, + pLogUniform=momentumConfig.logUniform, eta=(etaConfig.min, etaConfig.max), phi=(phiConfig.min, phiConfig.max), etaUniform=etaConfig.uniform, @@ -335,7 +336,7 @@ def addPythia8( acts.examples.RootParticleWriter( level=customLogLevel(), inputParticles=evGen.config.outputParticles, - filePath=str(outputDirRoot / "pythia8_particles.root"), + filePath=str(outputDirRoot / "particles.root"), ) ) @@ -343,7 +344,7 @@ def addPythia8( acts.examples.RootVertexWriter( level=customLogLevel(), inputVertices=evGen.config.outputVertices, - filePath=str(outputDirRoot / "pythia8_vertices.root"), + filePath=str(outputDirRoot / "vertices.root"), ) ) diff --git a/Examples/Python/src/Digitization.cpp b/Examples/Python/src/Digitization.cpp index c392d6dcc8c..948bcc9593b 100644 --- a/Examples/Python/src/Digitization.cpp +++ b/Examples/Python/src/Digitization.cpp @@ -13,6 +13,7 @@ #include "ActsExamples/Digitization/DigitizationAlgorithm.hpp" #include "ActsExamples/Digitization/DigitizationConfig.hpp" #include "ActsExamples/Digitization/DigitizationConfigurator.hpp" +#include "ActsExamples/Digitization/DigitizationCoordinatesConverter.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" #include "ActsExamples/Io/Json/JsonDigitizationConfig.hpp" @@ -101,6 +102,19 @@ void addDigitization(Context& ctx) { ACTS_PYTHON_MEMBER(outputDigiComponents); ACTS_PYTHON_STRUCT_END(); } + + { + py::class_>( + mex, "DigitizationCoordinatesConverter") + .def(py::init(), py::arg("config")) + .def_property_readonly( + "config", &ActsExamples::DigitizationCoordinatesConverter::config) + .def("globalToLocal", + &ActsExamples::DigitizationCoordinatesConverter::globalToLocal) + .def("localToGlobal", + &ActsExamples::DigitizationCoordinatesConverter::localToGlobal); + } } } // namespace Acts::Python diff --git a/Examples/Python/src/Generators.cpp b/Examples/Python/src/Generators.cpp index 800c912aa89..454227ddbf6 100644 --- a/Examples/Python/src/Generators.cpp +++ b/Examples/Python/src/Generators.cpp @@ -176,6 +176,7 @@ void addGenerators(Context& ctx) { .def_readwrite("pMin", &Config::pMin) .def_readwrite("pMax", &Config::pMax) .def_readwrite("pTransverse", &Config::pTransverse) + .def_readwrite("pLogUniform", &Config::pLogUniform) .def_readwrite("pdg", &Config::pdg) .def_readwrite("randomizeCharge", &Config::randomizeCharge) .def_readwrite("numParticles", &Config::numParticles) diff --git a/Examples/Python/src/Navigation.cpp b/Examples/Python/src/Navigation.cpp index ac2f1a77ef6..733bbe1a151 100644 --- a/Examples/Python/src/Navigation.cpp +++ b/Examples/Python/src/Navigation.cpp @@ -42,7 +42,7 @@ struct AnyNavigationPolicyFactory : public Acts::NavigationPolicyFactory { template , typename... Policies> struct NavigationPolicyFactoryT : public AnyNavigationPolicyFactory { - NavigationPolicyFactoryT(Factory impl) + explicit NavigationPolicyFactoryT(Factory impl) requires(sizeof...(Policies) > 0) : m_impl(std::move(impl)) {} @@ -74,7 +74,7 @@ struct NavigationPolicyFactoryT : public AnyNavigationPolicyFactory { private: template std::unique_ptr add(Args&&... args) { - if constexpr (!((std::is_same_v || ...))) { + if constexpr (!(std::is_same_v || ...)) { auto impl = std::move(m_impl).template add(std::forward(args)...); return std::make_unique< diff --git a/Examples/Python/src/TruthTracking.cpp b/Examples/Python/src/TruthTracking.cpp index 2c43b3cfc4d..7278b3f821f 100644 --- a/Examples/Python/src/TruthTracking.cpp +++ b/Examples/Python/src/TruthTracking.cpp @@ -13,15 +13,11 @@ #include "ActsExamples/TruthTracking/TrackModifier.hpp" #include "ActsExamples/TruthTracking/TrackParameterSelector.hpp" #include "ActsExamples/TruthTracking/TrackTruthMatcher.hpp" -#include "ActsExamples/TruthTracking/TruthSeedSelector.hpp" #include "ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp" #include "ActsExamples/TruthTracking/TruthTrackFinder.hpp" #include "ActsExamples/TruthTracking/TruthVertexFinder.hpp" #include "ActsExamples/Utilities/HitSelector.hpp" -#include "ActsExamples/Utilities/Range.hpp" -#include -#include #include #include @@ -45,48 +41,6 @@ void addTruthTracking(Context& ctx) { ActsExamples::TruthTrackFinder, mex, "TruthTrackFinder", inputParticles, inputMeasurementParticlesMap, outputProtoTracks); - { - using Alg = ActsExamples::TruthSeedSelector; - using Config = Alg::Config; - - auto alg = py::class_>( - mex, "TruthSeedSelector") - .def(py::init(), - py::arg("config"), py::arg("level")) - .def_property_readonly("config", &Alg::config); - - auto c = py::class_(alg, "Config").def(py::init<>()); - - ACTS_PYTHON_STRUCT_BEGIN(c, Config); - ACTS_PYTHON_MEMBER(inputParticles); - ACTS_PYTHON_MEMBER(inputMeasurementParticlesMap); - ACTS_PYTHON_MEMBER(outputParticles); - ACTS_PYTHON_MEMBER(rhoMin); - ACTS_PYTHON_MEMBER(rhoMax); - ACTS_PYTHON_MEMBER(zMin); - ACTS_PYTHON_MEMBER(zMax); - ACTS_PYTHON_MEMBER(phiMin); - ACTS_PYTHON_MEMBER(phiMax); - ACTS_PYTHON_MEMBER(etaMin); - ACTS_PYTHON_MEMBER(etaMax); - ACTS_PYTHON_MEMBER(absEtaMin); - ACTS_PYTHON_MEMBER(absEtaMax); - ACTS_PYTHON_MEMBER(ptMin); - ACTS_PYTHON_MEMBER(ptMax); - ACTS_PYTHON_MEMBER(keepNeutral); - ACTS_PYTHON_MEMBER(nHitsMin); - ACTS_PYTHON_MEMBER(nHitsMax); - ACTS_PYTHON_STRUCT_END(); - - pythonRangeProperty(c, "rho", &Config::rhoMin, &Config::rhoMax); - pythonRangeProperty(c, "z", &Config::zMin, &Config::zMax); - pythonRangeProperty(c, "phi", &Config::phiMin, &Config::phiMax); - pythonRangeProperty(c, "eta", &Config::etaMin, &Config::etaMax); - pythonRangeProperty(c, "absEta", &Config::absEtaMin, &Config::absEtaMax); - pythonRangeProperty(c, "pt", &Config::ptMin, &Config::ptMax); - pythonRangeProperty(c, "nHits", &Config::nHitsMin, &Config::nHitsMax); - } - ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::ParticleSmearing, mex, "ParticleSmearing", inputParticles, outputTrackParameters, sigmaD0, sigmaD0PtA, sigmaD0PtB, sigmaZ0, diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 3ade64a68ed..2ecf8b57c94 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -1,4 +1,4 @@ -test_pythia8__pythia8_particles.root: 91c852f3e0e20bcd382c616a7b643985d092decd42bdd653deae67ed8652e8d8 +test_pythia8__particles.root: 91c852f3e0e20bcd382c616a7b643985d092decd42bdd653deae67ed8652e8d8 test_fatras__particles_simulation.root: bc970873fef0c2efd86ed5413623802353d2cd04abea72de14e8cdfc0e40076f test_fatras__hits.root: 6e4beb045fa1712c4d14c280ba33c3fa13e4aff9de88d55c3e32f62ad226f724 test_geant4__particles_simulation.root: 49926c71a9b54e13aa1cc7596d3302baf3c87d8e2c1d0267cb4523f6abdc0ac2 diff --git a/Examples/Python/tests/test_algorithms.py b/Examples/Python/tests/test_algorithms.py index ede2610f604..4cca7199291 100644 --- a/Examples/Python/tests/test_algorithms.py +++ b/Examples/Python/tests/test_algorithms.py @@ -11,7 +11,6 @@ EventGenerator, FatrasSimulation, MaterialMapping, - TruthSeedSelector, TruthTrackFinder, ParticleSelector, TruthVertexFinder, @@ -42,7 +41,6 @@ EventGenerator, FatrasSimulation, MaterialMapping, - TruthSeedSelector, TruthTrackFinder, ParticleSelector, TruthVertexFinder, diff --git a/Examples/Python/tests/test_detectors.py b/Examples/Python/tests/test_detectors.py index e99406c69d6..c4d2c6ea9ed 100644 --- a/Examples/Python/tests/test_detectors.py +++ b/Examples/Python/tests/test_detectors.py @@ -1,4 +1,5 @@ import pytest +from pathlib import Path from helpers import dd4hepEnabled @@ -157,3 +158,34 @@ def eq(self, other): v = Volume(**{key: (4, None)}) assert getattr(v, key) == Interval(4, None) + + +def test_coordinate_converter(trk_geo): + digiCfg = acts.examples.DigitizationConfig( + acts.examples.readDigiConfigFromJson( + str( + Path(__file__).parent.parent.parent.parent + / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json" + ) + ), + surfaceByIdentifier=trk_geo.geoIdSurfaceMap(), + ) + converter = acts.examples.DigitizationCoordinatesConverter(digiCfg) + + def test_surface(surface): + gctx = acts.GeometryContext() + geo_id = surface.geometryId().value() + geo_center = surface.center(gctx) + x, y, z = geo_center[0], geo_center[1], geo_center[2] + + # test if surface center can be reproduced + assert converter.globalToLocal(geo_id, x, y, z) == (0, 0) + assert converter.localToGlobal(geo_id, 0, 0) == (x, y, z) + + # test if we can get back to the same local coordinates + global_shifted = converter.localToGlobal(geo_id, 5, 5) + local_shifted = converter.globalToLocal(geo_id, *global_shifted) + assert abs(local_shifted[0] - 5) / 5 < 1e-6 + assert abs(local_shifted[1] - 5) / 5 < 1e-6 + + trk_geo.visitSurfaces(test_surface) diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py index 26f69f4aa82..69f06101fdd 100644 --- a/Examples/Python/tests/test_examples.py +++ b/Examples/Python/tests/test_examples.py @@ -87,14 +87,14 @@ def test_pythia8(tmp_path, seq, assert_root_hash): (tmp_path / "csv").mkdir() - assert not (tmp_path / "pythia8_particles.root").exists() + assert not (tmp_path / "particles.root").exists() assert len(list((tmp_path / "csv").iterdir())) == 0 events = seq.config.events runPythia8(str(tmp_path), outputRoot=True, outputCsv=True, s=seq).run() - fp = tmp_path / "pythia8_particles.root" + fp = tmp_path / "particles.root" assert fp.exists() assert fp.stat().st_size > 2**10 * 50 assert_entries(fp, "particles", events) diff --git a/Examples/Scripts/MaterialMapping/Mat_map.C b/Examples/Scripts/MaterialMapping/Mat_map.C index 9f58e605244..60c70fa2339 100644 --- a/Examples/Scripts/MaterialMapping/Mat_map.C +++ b/Examples/Scripts/MaterialMapping/Mat_map.C @@ -57,7 +57,7 @@ void Draw_ratio(TCanvas* c, TProfile* h1, TProfile* h2, TLegend* leg, std::strin h5->SetStats(0); // No statistics on lower plot h5->Divide(h1); - double maxi = min( max( fabs(h5->GetMinimum()-0.1*h5->GetMinimum()),h5->GetMaximum()+0.1*h5->GetMaximum() ), 10. ); + double maxi = min( max( std::abs(h5->GetMinimum()-0.1*h5->GetMinimum()),h5->GetMaximum()+0.1*h5->GetMaximum() ), 10. ); h5->SetMinimum( 0.5 ); // Define Y .. h5->SetMaximum( 1.1 ); // .. range @@ -145,7 +145,7 @@ void Mat_map(std::string Val = "", std::string geantino = "", std::string name = // 2D map for Validation input TCanvas *VM = new TCanvas("VM","Validation Map") ; - Val_file->Draw("mat_y:mat_z","fabs(mat_x)<1"); + Val_file->Draw("mat_y:mat_z","std::abs(mat_x)<1"); eta_0->Draw("Same"); eta_1p->Draw("Same"); @@ -206,7 +206,7 @@ void Mat_map(std::string Val = "", std::string geantino = "", std::string name = // 2D map for Geantino input TCanvas *GM = new TCanvas("GM","Geantino Map") ; - geantino_file->Draw("mat_y:mat_z","fabs(mat_x)<1"); + geantino_file->Draw("mat_y:mat_z","std::abs(mat_x)<1"); eta_0->Draw("Same"); eta_1p->Draw("Same"); diff --git a/Examples/Scripts/Optimization/ckf.py b/Examples/Scripts/Optimization/ckf.py index c1761141262..ba18291735c 100755 --- a/Examples/Scripts/Optimization/ckf.py +++ b/Examples/Scripts/Optimization/ckf.py @@ -263,7 +263,7 @@ def runCKFTracks( field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T)) - inputParticlePath = Path(Inputdir) / "pythia8_particles.root" + inputParticlePath = Path(Inputdir) / "particles.root" if not inputParticlePath.exists(): inputParticlePath = None diff --git a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp index dea8a0ae334..b1066919dc7 100644 --- a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp +++ b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp @@ -41,7 +41,7 @@ struct NuclearInteraction { /// The storage of the parameterisation detail::MultiParticleNuclearInteractionParametrisation multiParticleParameterisation; - /// The number of trials to match momenta and inveriant masses + /// The number of trials to match momenta and invariant masses //~ unsigned int nMatchingTrials = std::numeric_limits::max(); unsigned int nMatchingTrials = 100; unsigned int nMatchingTrialsTotal = 1000; @@ -56,7 +56,7 @@ struct NuclearInteraction { template std::pair generatePathLimits(generator_t& generator, const Particle& particle) const { - // Fast exit: No paramtrization provided + // Fast exit: No parameterisation provided if (multiParticleParameterisation.empty()) { return std::make_pair(std::numeric_limits::infinity(), std::numeric_limits::infinity()); @@ -416,7 +416,7 @@ Acts::ActsDynamicVector NuclearInteraction::sampleInvariantMasses( for (unsigned int i = 0; i < size; i++) { float variance = parametrisation.eigenvaluesInvariantMass[i]; std::normal_distribution dist{ - parametrisation.meanInvariantMass[i], sqrtf(variance)}; + parametrisation.meanInvariantMass[i], std::sqrt(variance)}; parameters[i] = dist(generator); } // Transform to multivariate normal distribution @@ -446,7 +446,7 @@ Acts::ActsDynamicVector NuclearInteraction::sampleMomenta( for (unsigned int i = 0; i < size; i++) { float variance = parametrisation.eigenvaluesMomentum[i]; std::normal_distribution dist{ - parametrisation.meanMomentum[i], sqrtf(variance)}; + parametrisation.meanMomentum[i], std::sqrt(variance)}; parameters[i] = dist(generator); } diff --git a/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp b/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp index ae80bb1a0c7..848de200db0 100644 --- a/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp +++ b/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp @@ -429,7 +429,8 @@ std::shared_ptr volumeBuilder_dd4hep( plbConfig.layerIdentification = subDetector.name(); plbConfig.centralLayerRadii = std::vector(1, 0.5 * (rMax + rMin)); plbConfig.centralLayerHalflengthZ = std::vector(1, halfZ); - plbConfig.centralLayerThickness = std::vector(1, fabs(rMax - rMin)); + plbConfig.centralLayerThickness = + std::vector(1, std::abs(rMax - rMin)); plbConfig.centralLayerMaterial = {plMaterial}; auto pcLayerBuilder = std::make_shared( plbConfig, logger.clone(std::string("D2A_PL:") + subDetector.name())); diff --git a/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp b/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp index 3b85391ab84..c7cffaaf5aa 100644 --- a/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp +++ b/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp @@ -208,8 +208,8 @@ const Acts::LayerVector Acts::DD4hepLayerBuilder::endcapLayers( // create the share disc bounds auto dBounds = std::make_shared( pl.min(Acts::BinningValue::binR), pl.max(Acts::BinningValue::binR)); - double thickness = std::fabs(pl.max(Acts::BinningValue::binZ) - - pl.min(Acts::BinningValue::binZ)); + double thickness = std::abs(pl.max(Acts::BinningValue::binZ) - + pl.min(Acts::BinningValue::binZ)); // Create the layer containing the sensitive surface endcapLayer = DiscLayer::create(transform, dBounds, std::move(sArray), thickness, nullptr, Acts::active); @@ -357,8 +357,8 @@ const Acts::LayerVector Acts::DD4hepLayerBuilder::centralLayers( double layerR = (pl.min(Acts::BinningValue::binR) + pl.max(Acts::BinningValue::binR)) * 0.5; - double thickness = std::fabs(pl.max(Acts::BinningValue::binR) - - pl.min(Acts::BinningValue::binR)); + double thickness = std::abs(pl.max(Acts::BinningValue::binR) - + pl.min(Acts::BinningValue::binR)); auto cBounds = std::make_shared(layerR, halfZ); // Create the layer containing the sensitive surface centralLayer = diff --git a/Plugins/GeoModel/src/detail/GeoPolygonConverter.cpp b/Plugins/GeoModel/src/detail/GeoPolygonConverter.cpp index 6b0642d25e8..11a0f983f1a 100644 --- a/Plugins/GeoModel/src/detail/GeoPolygonConverter.cpp +++ b/Plugins/GeoModel/src/detail/GeoPolygonConverter.cpp @@ -47,9 +47,9 @@ Acts::detail::GeoPolygonConverter::operator()( // sort based on the y-coordinate std::ranges::sort(vertices, {}, [](const auto& v) { return v[1]; }); if (nVertices == 4) { - double hlxnegy = fabs(vertices[0][0] - vertices[1][0]) / 2; - double hlxposy = fabs(vertices[2][0] - vertices[3][0]) / 2; - double hly = fabs(vertices[0][1] - vertices[3][1]) / 2; + double hlxnegy = std::abs(vertices[0][0] - vertices[1][0]) / 2; + double hlxposy = std::abs(vertices[2][0] - vertices[3][0]) / 2; + double hly = std::abs(vertices[0][1] - vertices[3][1]) / 2; std::vector halfLengths = {hlxnegy, hlxposy, hly}; // Create the surface @@ -78,10 +78,10 @@ Acts::detail::GeoPolygonConverter::operator()( // Return the detector element and surface return std::make_tuple(detectorElement, surface); } else if (nVertices == 6) { - double hlxnegy = fabs(vertices[0][0] - vertices[1][0]) / 2; - double hlxzeroy = fabs(vertices[2][0] - vertices[3][0]) / 2; - double hlxposy = fabs(vertices[4][0] - vertices[5][0]) / 2; - double hly = fabs(vertices[0][1] - vertices[4][1]) / 2; + double hlxnegy = std::abs(vertices[0][0] - vertices[1][0]) / 2; + double hlxzeroy = std::abs(vertices[2][0] - vertices[3][0]) / 2; + double hlxposy = std::abs(vertices[4][0] - vertices[5][0]) / 2; + double hly = std::abs(vertices[0][1] - vertices[4][1]) / 2; std::vector halfLengths = {hlxnegy, hlxzeroy, hlxposy, hly, hly}; diff --git a/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.hpp b/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.hpp index 12b3137682c..b7e060a1f80 100644 --- a/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.hpp +++ b/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.hpp @@ -318,7 +318,7 @@ inline SPForSeed* AtlasSeedFinder::newSpacePoint( if (m_checketa) { // filter SP outside of eta-range - float z = (fabs(r[2]) + m_zmax); + float z = std::abs(r[2]) + m_zmax; float x = r[0] * m_dzdrmin; float y = r[1] * m_dzdrmin; if ((z * z) < (x * x + y * y)) { diff --git a/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.ipp b/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.ipp index ce4f4d5525a..4b5b5c60875 100644 --- a/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.ipp +++ b/Plugins/Legacy/include/Acts/Seeding/AtlasSeedFinder.ipp @@ -205,7 +205,7 @@ void Acts::Legacy::AtlasSeedFinder::findNext() { /////////////////////////////////////////////////////////////////// template void Acts::Legacy::AtlasSeedFinder::buildFrameWork() { - m_ptmin = fabs(m_ptmin); + m_ptmin = std::abs(m_ptmin); if (m_ptmin < 100.) { m_ptmin = 100.; @@ -218,7 +218,7 @@ void Acts::Legacy::AtlasSeedFinder::buildFrameWork() { m_divermax = m_diversss; } - if (fabs(m_etamin) < .1) { + if (std::abs(m_etamin) < 0.1) { m_etamin = -m_etamax; } m_dzdrmax0 = 1. / tan(2. * atan(exp(-m_etamax))); @@ -227,7 +227,7 @@ void Acts::Legacy::AtlasSeedFinder::buildFrameWork() { // scattering factor. depends on error, forward direction and distance between // SP m_COF = 134 * .05 * 9.; - m_ipt = 1. / fabs(.9 * m_ptmin); + m_ipt = 1. / std::abs(0.9 * m_ptmin); m_ipt2 = m_ipt * m_ipt; m_K = 0.; @@ -646,7 +646,7 @@ void Acts::Legacy::AtlasSeedFinder::production3Sp( } // forward direction of SP duplet float Tz = (Z - (*r)->z()) / dR; - float aTz = fabs(Tz); + float aTz = std::abs(Tz); // why also exclude seeds with small pseudorapidity?? if (aTz < m_dzdrmin || aTz > m_dzdrmax) { continue; @@ -690,7 +690,7 @@ void Acts::Legacy::AtlasSeedFinder::production3Sp( } float Tz = ((*r)->z() - Z) / dR; - float aTz = fabs(Tz); + float aTz = std::abs(Tz); if (aTz < m_dzdrmin || aTz > m_dzdrmax) { continue; @@ -793,14 +793,14 @@ void Acts::Legacy::AtlasSeedFinder::production3Sp( continue; } - float Im = fabs((A - B * R) * R); + float Im = std::abs((A - B * R) * R); if (Im <= imax) { // Add penalty factor dependent on difference between cot(theta) to // the quality Im (previously Impact) float dr = 0; m_R[t] < m_R[b] ? dr = m_R[t] : dr = m_R[b]; - Im += fabs((Tzb - m_Tz[t]) / (dr * sTzb2)); + Im += std::abs((Tzb - m_Tz[t]) / (dr * sTzb2)); // B/sqrt(S2) = 1/helixradius m_CmSp.push_back(std::make_pair(B / sqrt(S2), m_SP[t])); m_SP[t]->setParam(Im); @@ -926,7 +926,7 @@ void Acts::Legacy::AtlasSeedFinder:: } // Compared seeds should have at least deltaRMin distance float Rj = (*j).second->radius(); - if (fabs(Rj - Ri) < m_drmin) { + if (std::abs(Rj - Ri) < m_drmin) { continue; } diff --git a/README.md b/README.md index 6cb20bfbe6c..7b75acee190 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ or *A Common Tracking Software* if you do not like recursive acronyms [![10.5281/zenodo.5141418](https://zenodo.org/badge/DOI/10.5281/zenodo.5141418.svg)](https://doi.org/10.5281/zenodo.5141418) [![Chat on Mattermost](https://badgen.net/badge/chat/on%20mattermost/cyan)](https://mattermost.web.cern.ch/acts/) -[![codecov](https://codecov.io/gh/acts-project/acts/graph/badge.svg)](https://codecov.io/gh/acts-project/acts) +[![Coverage](https://sonarcloud.io/api/project_badges/measure?project=acts-project_acts&metric=coverage)](https://sonarcloud.io/summary/new_code?id=acts-project_acts) [![Latest release](https://badgen.net/github/release/acts-project/acts)](https://github.com/acts-project/acts/releases) [![Status](https://badgen.net/github/checks/acts-project/acts/main)](https://github.com/acts-project/acts/actions) [![Metrics](https://badgen.net/badge/metric/tracker/purple)](https://acts-project.github.io/metrics/) diff --git a/Tests/UnitTests/Core/EventData/TrackParameterHelpersTests.cpp b/Tests/UnitTests/Core/EventData/TrackParameterHelpersTests.cpp index 0a6f3064594..bba77ff8ef2 100644 --- a/Tests/UnitTests/Core/EventData/TrackParameterHelpersTests.cpp +++ b/Tests/UnitTests/Core/EventData/TrackParameterHelpersTests.cpp @@ -14,6 +14,16 @@ BOOST_AUTO_TEST_SUITE(TrackParameterHelpers) +BOOST_AUTO_TEST_CASE(isBoundVectorValid) { + BOOST_CHECK(!Acts::isBoundVectorValid({1, 2, 3, 4, 5, 6}, true)); + BOOST_CHECK(Acts::isBoundVectorValid({1, 2, 1, 1, 5, 6}, true)); +} + +BOOST_AUTO_TEST_CASE(isFreeVectorValid) { + BOOST_CHECK(!Acts::isFreeVectorValid({1, 2, 3, 4, 5, 6, 7, 8})); + BOOST_CHECK(Acts::isFreeVectorValid({1, 2, 3, 4, 1, 0, 0, 8})); +} + BOOST_AUTO_TEST_CASE(normalizeBoundParameters) { CHECK_CLOSE_OR_SMALL(Acts::normalizeBoundParameters({1, 2, 3, 4, 5, 6}), Acts::BoundVector(1, 2, -0.141593, 2.28319, 5, 6), 1e-3, diff --git a/Tests/UnitTests/Core/Material/MaterialGridHelperTests.cpp b/Tests/UnitTests/Core/Material/MaterialGridHelperTests.cpp index 514da1ddd74..dde3a05edba 100644 --- a/Tests/UnitTests/Core/Material/MaterialGridHelperTests.cpp +++ b/Tests/UnitTests/Core/Material/MaterialGridHelperTests.cpp @@ -60,10 +60,8 @@ BOOST_AUTO_TEST_CASE(Square_Grid_test) { BOOST_CHECK_EQUAL(Grid.minPosition()[0], bd[0].min); BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min); - float max1 = - bd[0].max + std::fabs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); - float max2 = - bd[1].max + std::fabs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); + float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); + float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1); BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2); @@ -152,10 +150,8 @@ BOOST_AUTO_TEST_CASE(PhiZ_Grid_test) { BOOST_CHECK_EQUAL(Grid.minPosition()[0], bd[0].min); BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min); - float max1 = - bd[0].max + std::fabs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); - float max2 = - bd[1].max + std::fabs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); + float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); + float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1); BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2); @@ -244,12 +240,9 @@ BOOST_AUTO_TEST_CASE(Cubic_Grid_test) { BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min); BOOST_CHECK_EQUAL(Grid.minPosition()[2], bd[2].min); - float max1 = - bd[0].max + std::fabs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); - float max2 = - bd[1].max + std::fabs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); - float max3 = - bd[2].max + std::fabs(bd[2].max - bd[2].min) / (bd[2].bins() - 1); + float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); + float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); + float max3 = bd[2].max + std::abs(bd[2].max - bd[2].min) / (bd[2].bins() - 1); BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1); BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2); @@ -341,12 +334,9 @@ BOOST_AUTO_TEST_CASE(Cylindrical_Grid_test) { BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min); BOOST_CHECK_EQUAL(Grid.minPosition()[2], bd[2].min); - float max1 = - bd[0].max + std::fabs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); - float max2 = - bd[1].max + std::fabs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); - float max3 = - bd[2].max + std::fabs(bd[2].max - bd[2].min) / (bd[2].bins() - 1); + float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1); + float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1); + float max3 = bd[2].max + std::abs(bd[2].max - bd[2].min) / (bd[2].bins() - 1); BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1); BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2); diff --git a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp index b9842b55420..e7edefb1f03 100644 --- a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -124,9 +125,30 @@ auto makeDefaultBoundPars(bool cov = true, std::size_t n = 4, return c; }; + // note that we are using the default random device + std::mt19937 gen; + std::uniform_real_distribution<> locDis(-10.0, 10.0); + std::uniform_real_distribution<> phiDis(-M_PI, M_PI); + std::uniform_real_distribution<> thetaDis(0, M_PI); + std::uniform_real_distribution<> qOverPDis(-10.0, 10.0); + std::uniform_real_distribution<> timeDis(0.0, 100.0); + for (auto i = 0ul; i < n; ++i) { - cmps.push_back({1. / n, ext_pars ? *ext_pars : BoundVector::Random(), - cov ? Opt{make_random_sym_matrix()} : Opt{}}); + BoundVector params = BoundVector::Zero(); + + if (ext_pars) { + params = *ext_pars; + } else { + params[eBoundLoc0] = locDis(gen); + params[eBoundLoc1] = locDis(gen); + params[eBoundPhi] = phiDis(gen); + params[eBoundTheta] = thetaDis(gen); + params[eBoundQOverP] = qOverPDis(gen); + params[eBoundTime] = timeDis(gen); + } + + cmps.push_back( + {1. / n, params, cov ? Opt{make_random_sym_matrix()} : Opt{}}); } auto surface = Acts::CurvilinearSurface(Vector3::Zero(), Vector3{1., 0., 0.}) @@ -430,7 +452,8 @@ void test_multi_stepper_surface_status_update() { std::vector>> cmps(2, {0.5, BoundVector::Zero(), std::nullopt}); std::get(cmps[0])[eBoundTheta] = M_PI_2; - std::get(cmps[1])[eBoundTheta] = -M_PI_2; + std::get(cmps[1])[eBoundPhi] = M_PI; + std::get(cmps[1])[eBoundTheta] = M_PI_2; std::get(cmps[0])[eBoundQOverP] = 1.0; std::get(cmps[1])[eBoundQOverP] = 1.0; @@ -541,7 +564,8 @@ void test_component_bound_state() { std::vector>> cmps(2, {0.5, BoundVector::Zero(), std::nullopt}); std::get(cmps[0])[eBoundTheta] = M_PI_2; - std::get(cmps[1])[eBoundTheta] = -M_PI_2; + std::get(cmps[1])[eBoundPhi] = M_PI; + std::get(cmps[1])[eBoundTheta] = M_PI_2; std::get(cmps[0])[eBoundQOverP] = 1.0; std::get(cmps[1])[eBoundQOverP] = 1.0; @@ -703,18 +727,7 @@ void test_single_component_interface_function() { using MultiState = typename multi_stepper_t::State; using MultiStepper = multi_stepper_t; - std::vector>> - cmps; - for (int i = 0; i < 4; ++i) { - cmps.push_back({0.25, BoundVector::Random(), BoundSquareMatrix::Random()}); - } - - auto surface = - Acts::CurvilinearSurface(Vector3::Zero(), Vector3::Ones().normalized()) - .planeSurface(); - - MultiComponentBoundTrackParameters multi_pars(surface, cmps, - particleHypothesis); + MultiComponentBoundTrackParameters multi_pars = makeDefaultBoundPars(true, 4); MultiState multi_state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize); diff --git a/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp b/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp index 9d89c1c96c6..8368bdbc3f2 100644 --- a/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp +++ b/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp @@ -186,8 +186,8 @@ int main(int argc, char** argv) { Acts::SeedFilterConfig sfconf; Acts::ATLASCuts atlasCuts = Acts::ATLASCuts(); - config.seedFilter = std::make_unique>( - Acts::SeedFilter(sfconf, &atlasCuts)); + config.seedFilter = + std::make_unique>(sfconf, &atlasCuts); Acts::SeedFinder> a; // test creation of unconfigured finder a = Acts::SeedFinder>( diff --git a/Tests/UnitTests/Core/Surfaces/LineSurfaceTests.cpp b/Tests/UnitTests/Core/Surfaces/LineSurfaceTests.cpp index 0175578a423..15d930929af 100644 --- a/Tests/UnitTests/Core/Surfaces/LineSurfaceTests.cpp +++ b/Tests/UnitTests/Core/Surfaces/LineSurfaceTests.cpp @@ -14,6 +14,7 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/ParticleHypothesis.hpp" #include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/detail/GenerateParameters.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Material/HomogeneousSurfaceMaterial.hpp" #include "Acts/Propagator/Propagator.hpp" @@ -343,8 +344,9 @@ BOOST_AUTO_TEST_CASE(LineSurfaceIntersection) { .closest(); CHECK_CLOSE_ABS(intersection.pathLength(), pathLimit, eps); - BoundTrackParameters endParameters{surface, BoundVector::Zero(), std::nullopt, - ParticleHypothesis::pion()}; + BoundTrackParameters endParameters{surface, + detail::Test::someBoundParametersA(), + std::nullopt, ParticleHypothesis::pion()}; { PropagatorOptions options(tgContext, {}); options.direction = Acts::Direction::Forward; diff --git a/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp b/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp index 6692974ef63..269370d4fc4 100644 --- a/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp +++ b/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp @@ -77,7 +77,7 @@ struct TestReverseFilteringLogic { template bool operator()(typename traj_t::ConstTrackStateProxy state) const { // can't determine an outlier w/o a measurement or predicted parameters - auto momentum = fabs(1 / state.filtered()[Acts::eBoundQOverP]); + auto momentum = std::abs(1 / state.filtered()[Acts::eBoundQOverP]); std::cout << "momentum : " << momentum << std::endl; return (momentum <= momentumMax); } diff --git a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp index c9743e90d7f..279f2c1a739 100644 --- a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -371,6 +372,9 @@ BOOST_DATA_TEST_CASE(VertexCompatibility4D, IPs* vertices, d0, l0, vx0, vy0, BoundVector paramVecClose = BoundVector::Zero(); paramVecClose[eBoundLoc0] = d0; paramVecClose[eBoundLoc1] = l0; + paramVecClose[eBoundPhi] = 0; + paramVecClose[eBoundTheta] = std::numbers::pi / 2; + paramVecClose[eBoundQOverP] = 0; paramVecClose[eBoundTime] = vt0 + sgnClose * timeDiffClose; BoundVector paramVecFar = paramVecClose; diff --git a/Tests/UnitTests/Core/Vertexing/SingleSeedVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/SingleSeedVertexFinderTests.cpp index 95ba93c0efb..7be7e174e6a 100644 --- a/Tests/UnitTests/Core/Vertexing/SingleSeedVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/SingleSeedVertexFinderTests.cpp @@ -294,10 +294,10 @@ BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_full_planes_test) { double x1 = (D * dirY + sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) / dirR2; double y1 = - (-D * dirX + std::fabs(dirY) * std::sqrt(r * r * dirR2 - D * D)) / + (-D * dirX + std::abs(dirY) * std::sqrt(r * r * dirR2 - D * D)) / dirR2; // how many units from the vertex to the intersection - double zDist = std::fabs((x1 - posX) / dirX); + double zDist = std::abs((x1 - posX) / dirX); // position of the new spacepoint posX = x1; @@ -411,11 +411,11 @@ BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_full_rays_test) { double x1 = (D * dirY + part * sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) / dirR2; - double y1 = (-D * dirX + part * std::fabs(dirY) * - std::sqrt(r * r * dirR2 - D * D)) / + double y1 = (-D * dirX + + part * std::abs(dirY) * std::sqrt(r * r * dirR2 - D * D)) / dirR2; // how many units from the vertex to the intersection - double zDist = std::fabs((x1 - vtxX) / dirX); + double zDist = std::abs((x1 - vtxX) / dirX); // use the same amount of units for distance in Z inputSpacepoints.emplace_back(x1, y1, zDist * dirZ + vtxZ); } diff --git a/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp b/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp index f281c7c1095..033f2c9cb4a 100644 --- a/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp +++ b/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp @@ -101,9 +101,9 @@ GeoGeometry constructGeoModel() { {123, 50}, {-123, 50}, {-153, 0}}; geoDims.tube = {5, 6, 100}; geoDims.trapHls = { - fabs(geoDims.trapVerts[0][0] - geoDims.trapVerts[1][0]) / 2, - fabs(geoDims.trapVerts[2][0] - geoDims.trapVerts[3][0]) / 2, - fabs(geoDims.trapVerts[0][1] - geoDims.trapVerts[2][1]) / 2}; + std::abs(geoDims.trapVerts[0][0] - geoDims.trapVerts[1][0]) / 2, + std::abs(geoDims.trapVerts[2][0] - geoDims.trapVerts[3][0]) / 2, + std::abs(geoDims.trapVerts[0][1] - geoDims.trapVerts[2][1]) / 2}; // create shapes GeoIntrusivePtr boxXY( diff --git a/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp index 566464a1889..5c0882cb12c 100644 --- a/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp +++ b/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp @@ -8,12 +8,12 @@ #include +#include "Acts/EventData/ParticleHypothesis.hpp" +#include "Acts/EventData/TrackParameters.hpp" #include "Acts/Plugins/Json/TrackParametersJsonConverter.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" #include "Acts/Surfaces/RectangleBounds.hpp" -#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" -#include #include #include diff --git a/docs/getting_started.md b/docs/getting_started.md index f17ec5abcf1..aff8e38a3b0 100644 --- a/docs/getting_started.md +++ b/docs/getting_started.md @@ -301,7 +301,6 @@ components. | ACTS_BUILD_EXAMPLES_HASHING | Build Hashing-based code in the examples
type: `bool`, default: `OFF` | | ACTS_BUILD_EXAMPLES_PYTHIA8 | Build Pythia8-based code in the examples
type: `bool`, default: `OFF` | | ACTS_BUILD_EXAMPLES_PYTHON_BINDINGS | Build python bindings for the examples
type: `bool`, default: `OFF` | -| ACTS_USE_EXAMPLES_TBB | Use Threading Building Blocks library in
the examples
type: `bool`, default: `ON` | | ACTS_BUILD_ANALYSIS_APPS | Build Analysis applications in the
examples
type: `bool`, default: `OFF` | | ACTS_BUILD_BENCHMARKS | Build benchmarks
type: `bool`, default: `OFF` | | ACTS_BUILD_INTEGRATIONTESTS | Build integration tests
type: `bool`, default: `OFF` |