diff --git a/CI/physmon/fpe_masks.yml b/CI/physmon/fpe_masks.yml index 1c8d227d8b6..33fd9a6381c 100644 --- a/CI/physmon/fpe_masks.yml +++ b/CI/physmon/fpe_masks.yml @@ -1,8 +1,8 @@ "Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheHeitler.hpp:66": FLTUND: 1 -"Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp:177": - FLTUND: 1 "Core/include/Acts/TrackFitting/detail/GsfUtils.hpp:197": - FLTUND: 1 + FLTUND: 1 +"Core/include/Acts/TrackFitting/detail/GsfUtils.hpp:201": + FLTUND: 1 "Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp:120": FLTUND: 1 diff --git a/CI/physmon/reference/performance_ambi_ttbar.root b/CI/physmon/reference/performance_ambi_ttbar.root index 0ad0b43cd6f..6bfe9250999 100644 Binary files a/CI/physmon/reference/performance_ambi_ttbar.root and b/CI/physmon/reference/performance_ambi_ttbar.root differ diff --git a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root index a9a393e8afb..60f6d54088e 100644 Binary files a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root differ diff --git a/CI/physmon/reference/performance_amvf_ttbar_hist.root b/CI/physmon/reference/performance_amvf_ttbar_hist.root index 16504e57111..33e77040382 100644 Binary files a/CI/physmon/reference/performance_amvf_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_ttbar_hist.root differ diff --git a/CI/physmon/reference/performance_ckf_ttbar.root b/CI/physmon/reference/performance_ckf_ttbar.root index dd7dbcdffa7..8bce7567c63 100644 Binary files a/CI/physmon/reference/performance_ckf_ttbar.root and b/CI/physmon/reference/performance_ckf_ttbar.root differ diff --git a/CI/physmon/reference/performance_gsf.root b/CI/physmon/reference/performance_gsf.root index 5e60795b569..05f33a8f130 100644 Binary files a/CI/physmon/reference/performance_gsf.root and b/CI/physmon/reference/performance_gsf.root differ diff --git a/CI/physmon/reference/performance_seeding_ttbar.root b/CI/physmon/reference/performance_seeding_ttbar.root index 8aeb7ae7379..38a6722ece5 100644 Binary files a/CI/physmon/reference/performance_seeding_ttbar.root and b/CI/physmon/reference/performance_seeding_ttbar.root differ diff --git a/CI/physmon/reference/performance_truth_tracking.root b/CI/physmon/reference/performance_truth_tracking.root index 069dacb1b94..28ffe8c18d7 100644 Binary files a/CI/physmon/reference/performance_truth_tracking.root and b/CI/physmon/reference/performance_truth_tracking.root differ diff --git a/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root index 4e63e8dcdbd..54c2b202133 100644 Binary files a/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root and b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root differ diff --git a/CI/physmon/workflows/physmon_track_finding_ttbar.py b/CI/physmon/workflows/physmon_track_finding_ttbar.py index 7107a6a492c..bbc023bb7cb 100755 --- a/CI/physmon/workflows/physmon_track_finding_ttbar.py +++ b/CI/physmon/workflows/physmon_track_finding_ttbar.py @@ -43,6 +43,12 @@ acts.FpeType.FLTINV, 1, ), + acts.examples.Sequencer.FpeMask( + "Examples/Algorithms/Fatras/src/FatrasSimulation.cpp", + (235, 236), + acts.FpeType.FLTINV, + 1, + ), acts.examples.Sequencer.FpeMask( "Examples/Algorithms/Fatras/src/FatrasSimulation.cpp", (235, 236), diff --git a/CI/physmon/workflows/physmon_truth_tracking_gsf.py b/CI/physmon/workflows/physmon_truth_tracking_gsf.py index 518fc04af95..ca92f7b183a 100755 --- a/CI/physmon/workflows/physmon_truth_tracking_gsf.py +++ b/CI/physmon/workflows/physmon_truth_tracking_gsf.py @@ -12,7 +12,7 @@ with tempfile.TemporaryDirectory() as temp: s = acts.examples.Sequencer( - events=500, + events=10000, numThreads=-1, logLevel=acts.logging.INFO, fpeMasks=acts.examples.Sequencer.FpeMask.fromFile( @@ -23,8 +23,8 @@ tp = Path(temp) runTruthTrackingGsf( setup.trackingGeometry, - setup.digiConfig, setup.field, + setup.digiConfig, outputDir=tp, s=s, ) diff --git a/CI/physmon/workflows/physmon_truth_tracking_kalman.py b/CI/physmon/workflows/physmon_truth_tracking_kalman.py index 325ba8968ad..c01e1a0e182 100755 --- a/CI/physmon/workflows/physmon_truth_tracking_kalman.py +++ b/CI/physmon/workflows/physmon_truth_tracking_kalman.py @@ -24,7 +24,7 @@ runTruthTrackingKalman( setup.trackingGeometry, setup.field, - digiConfigFile=setup.digiConfig, + setup.digiConfig, outputDir=tp, s=s, ) diff --git a/Core/include/Acts/Detector/Blueprint.hpp b/Core/include/Acts/Detector/Blueprint.hpp index ead405ff8cc..a6d6fccbc50 100644 --- a/Core/include/Acts/Detector/Blueprint.hpp +++ b/Core/include/Acts/Detector/Blueprint.hpp @@ -20,7 +20,9 @@ namespace Acts { namespace Experimental { +class IGeometryIdGenerator; class IInternalStructureBuilder; +class IRootVolumeFinderBuilder; /// A Blueprint is an instruction tree that allows you to defina a tree sequence /// of volume building using the provided tools. @@ -87,6 +89,12 @@ struct Node final { std::vector> children = {}; /// Branch definition binning std::vector binning = {}; + + /// Builders and helper tools that can be attached + std::shared_ptr rootVolumeFinderBuilder = + nullptr; + /// Geometry id generator + std::shared_ptr geoIdGenerator = nullptr; /// Internal structure builder - for leaf nodes std::shared_ptr internalsBuilder = nullptr; @@ -144,6 +152,22 @@ struct Node final { ss << name << " -> " << name + "_int" << ";" << '\n'; } + + if (geoIdGenerator != nullptr) { + ss << name + "_geoID" + << " [shape=\"note\";style=\"filled\";fillcolor=\"azure\"];" << '\n'; + ss << name << " -> " << name + "_geoID" + << ";" << '\n'; + } + + if (rootVolumeFinderBuilder != nullptr) { + ss << name + "_roots" + << " [shape=\"Msquare\";style=\"filled\";fillcolor=\"darkkhaki\"];" + << '\n'; + ss << name << " -> " << name + "_roots" + << ";" << '\n'; + } + if (isRoot()) { ss << "}" << '\n'; } diff --git a/Core/include/Acts/Detector/CylindricalContainerBuilder.hpp b/Core/include/Acts/Detector/CylindricalContainerBuilder.hpp index 8d213b589f6..0d7b546dfef 100644 --- a/Core/include/Acts/Detector/CylindricalContainerBuilder.hpp +++ b/Core/include/Acts/Detector/CylindricalContainerBuilder.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/Detector/Blueprint.hpp" #include "Acts/Detector/DetectorComponents.hpp" #include "Acts/Detector/interface/IDetectorComponentBuilder.hpp" #include "Acts/Geometry/GeometryContext.hpp" @@ -30,6 +31,10 @@ class IGeometryIdGenerator; /// and allows for DetectorVolume attachment in z/r/phi, such as wrapping /// of bevelled cylinder objects in z/r /// +/// There exists an option to create this container builder (recursively) +/// from a blueprint tree, which attempts to fill in the gap volumes +/// accordingly. +/// /// @note the builder expects a fully consistent set of sub volume builders /// that will be executed in a chain /// @@ -45,7 +50,8 @@ class CylindricalContainerBuilder : public IDetectorComponentBuilder { /// Binning prescription of attachment std::vector binning = {}; /// The root volume finder - std::shared_ptr rootVolumeFinderBuilder = nullptr; + std::shared_ptr rootVolumeFinderBuilder = + nullptr; /// The geometry id generator std::shared_ptr geoIdGenerator = nullptr; /// An eventual reverse geometry id generation @@ -54,7 +60,7 @@ class CylindricalContainerBuilder : public IDetectorComponentBuilder { std::string auxiliary = ""; }; - /// Constructor with configuration arguments + /// Constructor with configuration struct /// /// @param cfg is the configuration struct /// @param logger logging instance for screen output @@ -63,6 +69,24 @@ class CylindricalContainerBuilder : public IDetectorComponentBuilder { std::unique_ptr logger = getDefaultLogger("CylindricalContainerBuilder", Logging::INFO)); + /// Constructor from blueprint and logging level + /// + /// It will create recursively the builders of sub volumes + /// + /// @param bpNode is the entry blue print node + /// @param logLevel is the logging output level for the builder tools + /// + /// @note no checking is being done on consistency of the blueprint, + /// it is assumed it has passed first through gap filling via the + /// blueprint helper. + /// + /// @note that the naming of the builders is taken from the bluprint nodes + /// + /// @return a cylindrical container builder representing this blueprint + CylindricalContainerBuilder( + const Acts::Experimental::Blueprint::Node& bpNode, + Acts::Logging::Level logLevel = Acts::Logging::INFO); + /// The final implementation of the cylindrical container builder /// /// @param gctx The geometry context for this call diff --git a/Core/include/Acts/Detector/DetectorVolumeBuilder.hpp b/Core/include/Acts/Detector/DetectorVolumeBuilder.hpp index 35fba51948e..a5ca7b5fb75 100644 --- a/Core/include/Acts/Detector/DetectorVolumeBuilder.hpp +++ b/Core/include/Acts/Detector/DetectorVolumeBuilder.hpp @@ -10,8 +10,6 @@ #include "Acts/Detector/DetectorComponents.hpp" #include "Acts/Detector/interface/IDetectorComponentBuilder.hpp" -#include "Acts/Detector/interface/IExternalStructureBuilder.hpp" -#include "Acts/Detector/interface/IInternalStructureBuilder.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Utilities/Logger.hpp" diff --git a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp index 833b7851717..e690fa8a4da 100644 --- a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/Definitions/Tolerance.hpp" #include "Acts/EventData/detail/PrintParameters.hpp" #include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -82,6 +83,7 @@ class GenericBoundTrackParameters { /// @param qOverP Charge over momentum /// @param cov Bound parameters covariance matrix /// @param particleHypothesis Particle hypothesis + /// @param tolerance Tolerance used for globalToLocal /// /// @note The returned result indicates whether the free parameters could /// successfully be converted to on-surface parameters. @@ -89,9 +91,11 @@ class GenericBoundTrackParameters { std::shared_ptr surface, const GeometryContext& geoCtx, const Vector4& pos4, const Vector3& dir, Scalar qOverP, std::optional cov, - ParticleHypothesis particleHypothesis) { + ParticleHypothesis particleHypothesis, + ActsScalar tolerance = s_onSurfaceTolerance) { Result bound = detail::transformFreeToBoundParameters( - pos4.segment<3>(ePos0), pos4[eTime], dir, qOverP, *surface, geoCtx); + pos4.segment<3>(ePos0), pos4[eTime], dir, qOverP, *surface, geoCtx, + tolerance); if (!bound.ok()) { return bound.error(); diff --git a/Core/include/Acts/EventData/TrackProxy.hpp b/Core/include/Acts/EventData/TrackProxy.hpp index 0664042128e..d9d8747a69f 100644 --- a/Core/include/Acts/EventData/TrackProxy.hpp +++ b/Core/include/Acts/EventData/TrackProxy.hpp @@ -673,19 +673,33 @@ class TrackProxy { /// "innermost" track state /// @note This is dangerous with branching track state sequences, as it will break them /// @note This also automatically forward-links the track! + /// @param invertJacobians Whether to invert the Jacobians of the track states template > - void reverseTrackStates() { + void reverseTrackStates(bool invertJacobians = false) { IndexType current = tipIndex(); IndexType next = kInvalid; IndexType prev = kInvalid; stemIndex() = tipIndex(); + // @TODO: Maybe refactor to not need this variable if invertJacobians == false + BoundMatrix nextJacobian; + while (current != kInvalid) { auto ts = m_container->trackStateContainer().getTrackState(current); prev = ts.previous(); ts.template component(hashString("next")) = prev; ts.previous() = next; + if (invertJacobians) { + if (next != kInvalid) { + BoundMatrix curJacobian = ts.jacobian(); + ts.jacobian() = nextJacobian.inverse(); + nextJacobian = curJacobian; + } else { + nextJacobian = ts.jacobian(); + ts.jacobian().setZero(); + } + } next = current; tipIndex() = current; current = prev; diff --git a/Core/include/Acts/Vertexing/AMVFInfo.hpp b/Core/include/Acts/Vertexing/AMVFInfo.hpp index 247fc8e4fd0..534c72cdbab 100644 --- a/Core/include/Acts/Vertexing/AMVFInfo.hpp +++ b/Core/include/Acts/Vertexing/AMVFInfo.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -21,31 +21,34 @@ template struct VertexInfo { VertexInfo() = default; - VertexInfo(const Acts::Vertex& vtx, const Acts::Vector4& pos) - : constraintVertex(vtx), + VertexInfo(const Acts::Vertex& constr, + const Acts::Vector4& pos) + : constraint(constr), linPoint(pos), oldPosition(pos), seedPosition(pos) {} - // The constraint vertex - Acts::Vertex constraintVertex; + // Vertex constraint + Acts::Vertex constraint; - // The linearization point + // Point where all associated tracks are linearized Acts::Vector4 linPoint{Acts::Vector4::Zero()}; - // Old position from last iteration + // Vertex position from the last iteration of the fit Acts::Vector4 oldPosition{Acts::Vector4::Zero()}; - // The seed position + // The seed position (i.e., the first estimate for the vertex position as + // obtained by the vertex seed finder) Acts::Vector4 seedPosition{Acts::Vector4::Zero()}; - // Needs relinearization bool + // If set to true, the associated tracks need to be relinearized at a more + // recent vertex position bool relinearize = true; - // Vector of all track currently held by vertex + // Vector of all tracks that are currently assigned to vertex std::vector trackLinks; - std::map ip3dParams; + std::map impactParams3D; }; } // namespace Acts diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp index 2d5c43c2576..c2e23bd9959 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -135,6 +135,9 @@ class AdaptiveMultiVertexFitter { // Do smoothing after multivertex fit bool doSmoothing{false}; + + // Use time information when calculating the vertex compatibility + bool useTime{false}; }; /// @brief Constructor used if InputTrack_t type == BoundTrackParameters @@ -167,12 +170,12 @@ class AdaptiveMultiVertexFitter { m_extractParameters(func), m_logger(std::move(logger)) {} - /// @brief The actual fit function, performs a simultaneous - /// fit of all vertices in `verticesToFit` by invoking `fitImpl` + /// @brief Performs a simultaneous fit of all vertices in `verticesToFit` + /// by invoking `fitImpl`. /// - /// @param state The state object + /// @param state Fitter state /// @param verticesToFit Vector containing all vertices to be fitted - /// @param linearizer The track linearizer + /// @param linearizer Track linearizer /// @param vertexingOptions Vertexing options /// /// @return Result object @@ -181,26 +184,18 @@ class AdaptiveMultiVertexFitter { const Linearizer_t& linearizer, const VertexingOptions& vertexingOptions) const; - /// @brief Adds new vertex to an existing multi-vertex fit - /// and fits everything together (by invoking the fit_impl method): - /// 1. The new vertex is added to the fit: all associated tracks get - /// initialized, i.e. ParamsAtIP3d are created (from ImpactPointEstimator) - /// to be later able to estimate in a fast way the compatibility of the tracks - /// to their respective vertices. - /// 2. All tracks belonging to the new vertex are scanned and all the vertices - /// which share tracks with the new vertex to be fit are also added to the - /// fit. - /// 3. The multivertex fit is performed with all involved vertices. - /// - /// This has the advantage that only vertices that are affected by adding the - /// new vertex are refitted. + /// @brief Adds a new vertex to an existing multi-vertex fit. + /// 1. The 3D impact parameters are calculated for all tracks associated + /// with newVertex. + /// 2. A list of all vertices that are connected with newVertex via shared + /// tracks is created. It also considers indirect connections (e.g., vertex A + /// which shares a track with vertex B which, in turn, shares a track with + /// newVertex). + /// 3. The multivertex fit is performed for all vertices on said list. /// - /// Note: newVertex has to be properly initialized (seed vertex, - /// constraint vertex, list of MAV) - /// - /// @param state The state object - /// @param newVertex New vertex to be added to fit - /// @param linearizer The track linearizer + /// @param state Fitter state + /// @param newVertex Vertex to be added to fit + /// @param linearizer Track linearizer /// @param vertexingOptions Vertexing options /// /// @return Result object @@ -224,11 +219,11 @@ class AdaptiveMultiVertexFitter { /// Private access to logging instance const Logger& logger() const { return *m_logger; } - /// @brief The actual fit function, performs a simultaneous - /// fit of all vertices in state.vertexCollection + /// @brief Performs a simultaneous fit of all vertices in + /// state.vertexCollection /// - /// @param state The state object - /// @param linearizer The track linearizer + /// @param state Fitter state + /// @param linearizer Track linearizer /// @param vertexingOptions Vertexing options /// /// @return Result object @@ -246,23 +241,22 @@ class AdaptiveMultiVertexFitter { Vertex* vtx, const std::vector*>& verticesVec) const; - /// @brief Prepares vertex object for the actual fit, i.e. - /// all TrackAtVertex objects at current vertex will obtain - /// `ip3dParams` from ImpactPointEstimator::estimate3DImpactParameters - /// in order to later faster estimate compatibilities of track - /// with different vertices + /// @brief 1) Calls ImpactPointEstimator::estimate3DImpactParameters + /// for all tracks that are associated with vtx (i.e., all elements + /// of the trackLinks vector in the VertexInfo of vtx). + /// 2) Saves the 3D impact parameters in the VertexInfo of vtx. /// - /// @param state The state to operate on - /// @param vtx The vertex object + /// @param state Vertex fitter state + /// @param vtx Vertex object /// @param vertexingOptions Vertexing options Result prepareVertexForFit( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const; - /// @brief Sets vertexCompatibility for all TrackAtVertex objects - /// at current vertex + /// @brief Sets the vertexCompatibility for all TrackAtVertex objects + /// at the current vertex /// - /// @param state The state object + /// @param state Fitter state /// @param currentVtx Current vertex /// @param vertexingOptions Vertexing options Result setAllVertexCompatibilities( @@ -272,36 +266,35 @@ class AdaptiveMultiVertexFitter { /// @brief Sets weights to the track according to Eq.(5.46) in Ref.(1) /// and updates the vertices by calling the VertexUpdater /// - /// @param state The state object + /// @param state Fitter state /// @param linearizer The track linearizer /// @param vertexingOptions Vertexing options Result setWeightsAndUpdate( State& state, const Linearizer_t& linearizer, const VertexingOptions& vertexingOptions) const; - /// @brief Collects all compatibility values of the track `trk` - /// at all vertices it is currently attached to and outputs - /// these values in a vector + /// @brief Collects the compatibility values of the track `trk` + /// wrt to all of its associated vertices /// - /// @param state The state object - /// @param trk The track + /// @param state Fitter state + /// @param trk Track /// /// @return Vector of compatibility values std::vector collectTrackToVertexCompatibilities( State& state, const InputTrack_t* trk) const; - /// @brief Determines if vertex position has shifted more than - /// m_cfg.maxRelativeShift in last iteration + /// @brief Determines if any vertex position has shifted more than + /// m_cfg.maxRelativeShift in the last iteration /// - /// @param state The state object + /// @param state Fitter state /// - /// @return False if shift was larger than maxRelativeShift + /// @return False if at least one shift was larger than maxRelativeShift bool checkSmallShift(State& state) const; /// @brief Updates tracks for current vertex with knowledge /// of current vertex position /// - /// @param state The state object + /// @param state Fitter state void doVertexSmoothing(State& state) const; }; diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index b137d421d9c..8f59ac534dc 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -37,12 +37,11 @@ Acts::AdaptiveMultiVertexFitter::fitImpl( // Reset annealing tool state.annealingState = AnnealingUtility::State(); - // Indicates how much the vertex positions have shifted - // in last fit iteration. Will be false if vertex position - // shift was too big. Needed if equilibrium is reached in - // annealing procedure but fitter has not fully converged - // yet and needs some more iterations until vertex position - // shifts between iterations are small (converged). + // Boolean indicating whether any of the vertices has moved more than + // m_cfg.maxRelativeShift during the last iteration. We will keep iterating + // until the equilibrium (i.e., the lowest temperature) is reached in + // the annealing procedure and isSmallShift is true (or until the maximum + // number of iterations is exceeded). bool isSmallShift = true; // Number of iterations counter @@ -52,48 +51,54 @@ Acts::AdaptiveMultiVertexFitter::fitImpl( while (nIter < m_cfg.maxIterations && (!state.annealingState.equilibriumReached || !isSmallShift)) { // Initial loop over all vertices in state.vertexCollection - for (auto currentVtx : state.vertexCollection) { - VertexInfo& currentVtxInfo = state.vtxInfoMap[currentVtx]; - currentVtxInfo.relinearize = false; + for (auto vtx : state.vertexCollection) { + VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; + vtxInfo.relinearize = false; // Store old position of vertex, i.e. seed position // in case of first iteration or position determined // in previous iteration afterwards - currentVtxInfo.oldPosition = currentVtx->fullPosition(); - - Vector4 dist = currentVtxInfo.oldPosition - currentVtxInfo.linPoint; - double perpDist = std::hypot(dist[0], dist[1]); - // Determine if relinearization is needed - if (perpDist > m_cfg.maxDistToLinPoint) { - // Relinearization needed, distance too big - currentVtxInfo.relinearize = true; - // Prepare for fit with new vertex position - prepareVertexForFit(state, currentVtx, vertexingOptions); + vtxInfo.oldPosition = vtx->fullPosition(); + + // Calculate the x-y-distance between the current vertex position + // and the linearization point of the tracks. If it is too large, + // we relinearize the tracks and recalculate their 3D impact + // parameters. + ActsVector<2> xyDiff = vtxInfo.oldPosition.template head<2>() - + vtxInfo.linPoint.template head<2>(); + if (xyDiff.norm() > m_cfg.maxDistToLinPoint) { + // Set flag for relinearization + vtxInfo.relinearize = true; + // Recalculate the track impact parameters at the current vertex + // position + prepareVertexForFit(state, vtx, vertexingOptions); } - // Determine if constraint vertex exist - if (state.vtxInfoMap[currentVtx].constraintVertex.fullCovariance() != + + // Check if we use the constraint during the vertex fit + if (state.vtxInfoMap[vtx].constraint.fullCovariance() != SquareMatrix4::Zero()) { - currentVtx->setFullPosition( - state.vtxInfoMap[currentVtx].constraintVertex.fullPosition()); - currentVtx->setFitQuality( - state.vtxInfoMap[currentVtx].constraintVertex.fitQuality()); - currentVtx->setFullCovariance( - state.vtxInfoMap[currentVtx].constraintVertex.fullCovariance()); - } else if (currentVtx->fullCovariance() == SquareMatrix4::Zero()) { + const Acts::Vertex& constraint = + state.vtxInfoMap[vtx].constraint; + vtx->setFullPosition(constraint.fullPosition()); + vtx->setFitQuality(constraint.fitQuality()); + vtx->setFullCovariance(constraint.fullCovariance()); + } else if (vtx->fullCovariance() == SquareMatrix4::Zero()) { return VertexingError::NoCovariance; } // Set vertexCompatibility for all TrackAtVertex objects - // at current vertex - setAllVertexCompatibilities(state, currentVtx, vertexingOptions); + // at the current vertex + setAllVertexCompatibilities(state, vtx, vertexingOptions); } // End loop over vertex collection - // Now after having estimated all compatibilities of all tracks at - // all vertices, run again over all vertices to set track weights - // and update the vertex + // Recalculate all track weights and update vertices setWeightsAndUpdate(state, linearizer, vertexingOptions); + + // Cool the system down, i.e., reduce the temperature parameter. At lower + // temperatures, outlying tracks are downweighted more. if (!state.annealingState.equilibriumReached) { m_cfg.annealingTool.anneal(state.annealingState); } + isSmallShift = checkSmallShift(state); ++nIter; } @@ -119,46 +124,47 @@ Acts::AdaptiveMultiVertexFitter::addVtxToFit( std::vector*> verticesToFit; - // Prepares vtx and tracks for fast estimation method of their - // compatibility with vertex + // Save the 3D impact parameters of all tracks associated with newVertex. auto res = prepareVertexForFit(state, &newVertex, vertexingOptions); if (!res.ok()) { return res.error(); } + // List of vertices added in last iteration std::vector*> lastIterAddedVertices = {&newVertex}; // List of vertices added in current iteration std::vector*> currentIterAddedVertices; - // Loop as long as new vertices are found that share tracks with - // previously added vertices + // Fill verticesToFit with vertices that are connected to newVertex (via + // tracks and/or other vertices). while (!lastIterAddedVertices.empty()) { - for (auto& lastVtxIter : lastIterAddedVertices) { - // Loop over all track at current lastVtxIter + for (auto& lastIterAddedVertex : lastIterAddedVertices) { + // Loop over all tracks at lastIterAddedVertex const std::vector& trks = - state.vtxInfoMap[lastVtxIter].trackLinks; + state.vtxInfoMap[lastIterAddedVertex].trackLinks; for (const auto& trk : trks) { - // Retrieve list of links to all vertices that currently use the current - // track - auto range = state.trackToVerticesMultiMap.equal_range(trk); - - // Loop over all attached vertices and add those to vertex fit - // which are not already in `verticesToFit` - for (auto vtxIter = range.first; vtxIter != range.second; ++vtxIter) { - auto newVtxIter = vtxIter->second; - if (!isAlreadyInList(newVtxIter, verticesToFit)) { - // Add newVtxIter to verticesToFit - verticesToFit.push_back(newVtxIter); - - // Add newVtxIter vertex to currentIterAddedVertices - // if vertex != lastVtxIter - if (newVtxIter != lastVtxIter) { - currentIterAddedVertices.push_back(newVtxIter); + // Range of vertices that are associated with trk. The range is + // represented via its bounds: begin refers to the first iterator of the + // range; end refers to the iterator after the last iterator of the + // range. + auto [begin, end] = state.trackToVerticesMultiMap.equal_range(trk); + + for (auto it = begin; it != end; ++it) { + // it->first corresponds to trk, it->second to one of its associated + // vertices + auto vtxToFit = it->second; + // Add vertex to the fit if it is not already included + if (!isAlreadyInList(vtxToFit, verticesToFit)) { + verticesToFit.push_back(vtxToFit); + + // Collect vertices that were added this iteration + if (vtxToFit != lastIterAddedVertex) { + currentIterAddedVertices.push_back(vtxToFit); } } - } // End for loop over linksToVertices - } - } // End loop over lastIterAddedVertices + } // End for loop over range of associated vertices + } // End loop over trackLinks + } // End loop over lastIterAddedVertices lastIterAddedVertices = currentIterAddedVertices; currentIterAddedVertices.clear(); @@ -177,11 +183,9 @@ Acts::AdaptiveMultiVertexFitter::addVtxToFit( template bool Acts::AdaptiveMultiVertexFitter:: - isAlreadyInList( - Vertex* vtx, - const std::vector*>& verticesVec) const { - return std::find(verticesVec.begin(), verticesVec.end(), vtx) != - verticesVec.end(); + isAlreadyInList(Vertex* vtx, + const std::vector*>& vertices) const { + return std::find(vertices.begin(), vertices.end(), vtx) != vertices.end(); } template @@ -189,21 +193,21 @@ Acts::Result Acts:: AdaptiveMultiVertexFitter::prepareVertexForFit( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const { - // The current vertex info object - auto& currentVtxInfo = state.vtxInfoMap[vtx]; - // The seed position - const Vector3& seedPos = currentVtxInfo.seedPosition.template head<3>(); + // Vertex info object + auto& vtxInfo = state.vtxInfoMap[vtx]; + // Vertex seed position + const Vector3& seedPos = vtxInfo.seedPosition.template head<3>(); - // Loop over all tracks at current vertex - for (const auto& trk : currentVtxInfo.trackLinks) { + // Loop over all tracks at the vertex + for (const auto& trk : vtxInfo.trackLinks) { auto res = m_cfg.ipEst.estimate3DImpactParameters( vertexingOptions.geoContext, vertexingOptions.magFieldContext, m_extractParameters(*trk), seedPos, state.ipState); if (!res.ok()) { return res.error(); } - // Set ip3dParams for current trackAtVertex - currentVtxInfo.ip3dParams.emplace(trk, res.value()); + // Save 3D impact parameters of the track + vtxInfo.impactParams3D.emplace(trk, res.value()); } return {}; } @@ -212,37 +216,42 @@ template Acts::Result Acts::AdaptiveMultiVertexFitter:: setAllVertexCompatibilities( - State& state, Vertex* currentVtx, + State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const { - VertexInfo& currentVtxInfo = state.vtxInfoMap[currentVtx]; + VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; - // Loop over tracks at current vertex and - // estimate compatibility with vertex - for (const auto& trk : currentVtxInfo.trackLinks) { - auto& trkAtVtx = - state.tracksAtVerticesMap.at(std::make_pair(trk, currentVtx)); + // Loop over all tracks that are associated with vtx and estimate their + // compatibility + for (const auto& trk : vtxInfo.trackLinks) { + auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); // Recover from cases where linearization point != 0 but // more tracks were added later on - if (currentVtxInfo.ip3dParams.find(trk) == - currentVtxInfo.ip3dParams.end()) { + if (vtxInfo.impactParams3D.find(trk) == vtxInfo.impactParams3D.end()) { auto res = m_cfg.ipEst.estimate3DImpactParameters( vertexingOptions.geoContext, vertexingOptions.magFieldContext, - m_extractParameters(*trk), - VectorHelpers::position(currentVtxInfo.linPoint), state.ipState); + m_extractParameters(*trk), VectorHelpers::position(vtxInfo.linPoint), + state.ipState); if (!res.ok()) { return res.error(); } - // Set ip3dParams for current trackAtVertex - currentVtxInfo.ip3dParams.emplace(trk, res.value()); + // Set impactParams3D for current trackAtVertex + vtxInfo.impactParams3D.emplace(trk, res.value()); } // Set compatibility with current vertex - auto compRes = m_cfg.ipEst.template getVertexCompatibility<3>( - vertexingOptions.geoContext, &(currentVtxInfo.ip3dParams.at(trk)), - VectorHelpers::position(currentVtxInfo.oldPosition)); - if (!compRes.ok()) { - return compRes.error(); + Acts::Result compatibilityResult(0.); + if (m_cfg.useTime) { + compatibilityResult = m_cfg.ipEst.template getVertexCompatibility<4>( + vertexingOptions.geoContext, &(vtxInfo.impactParams3D.at(trk)), + vtxInfo.oldPosition); + } else { + compatibilityResult = m_cfg.ipEst.template getVertexCompatibility<3>( + vertexingOptions.geoContext, &(vtxInfo.impactParams3D.at(trk)), + VectorHelpers::position(vtxInfo.oldPosition)); + } + if (!compatibilityResult.ok()) { + return compatibilityResult.error(); } - trkAtVtx.vertexCompatibility = *compRes; + trkAtVtx.vertexCompatibility = *compatibilityResult; } return {}; } @@ -253,36 +262,36 @@ Acts::Result Acts:: State& state, const linearizer_t& linearizer, const VertexingOptions& vertexingOptions) const { for (auto vtx : state.vertexCollection) { - VertexInfo& currentVtxInfo = state.vtxInfoMap[vtx]; + VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; + + if (vtxInfo.relinearize) { + vtxInfo.linPoint = vtxInfo.oldPosition; + } const std::shared_ptr vtxPerigeeSurface = Surface::makeShared( - VectorHelpers::position(state.vtxInfoMap[vtx].oldPosition)); + VectorHelpers::position(vtxInfo.linPoint)); - for (const auto& trk : currentVtxInfo.trackLinks) { + for (const auto& trk : vtxInfo.trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); // Set trackWeight for current track - double currentTrkWeight = m_cfg.annealingTool.getWeight( + trkAtVtx.trackWeight = m_cfg.annealingTool.getWeight( state.annealingState, trkAtVtx.vertexCompatibility, collectTrackToVertexCompatibilities(state, trk)); - trkAtVtx.trackWeight = currentTrkWeight; if (trkAtVtx.trackWeight > m_cfg.minWeight) { - // Check if linearization state exists or need to be relinearized - if (not trkAtVtx.isLinearized || state.vtxInfoMap[vtx].relinearize) { + // Check if track is already linearized and whether we need to + // relinearize + if (!trkAtVtx.isLinearized || vtxInfo.relinearize) { auto result = linearizer.linearizeTrack( - m_extractParameters(*trk), state.vtxInfoMap[vtx].oldPosition[3], + m_extractParameters(*trk), vtxInfo.linPoint[3], *vtxPerigeeSurface, vertexingOptions.geoContext, vertexingOptions.magFieldContext, state.linearizerState); if (!result.ok()) { return result.error(); } - if (trkAtVtx.isLinearized) { - state.vtxInfoMap[vtx].linPoint = state.vtxInfoMap[vtx].oldPosition; - } - trkAtVtx.linearizedState = *result; trkAtVtx.isLinearized = true; } @@ -304,13 +313,22 @@ std::vector Acts::AdaptiveMultiVertexFitter:: collectTrackToVertexCompatibilities(State& state, const input_track_t* trk) const { + // Compatibilities of trk wrt all of its associated vertices std::vector trkToVtxCompatibilities; - trkToVtxCompatibilities.reserve(state.vertexCollection.size()); - auto range = state.trackToVerticesMultiMap.equal_range(trk); - for (auto vtxIter = range.first; vtxIter != range.second; ++vtxIter) { + // Range of vertices that are associated with trk. The range is + // represented via its bounds: begin refers to the first iterator of the + // range; end refers to the iterator after the last iterator of the range. + auto [begin, end] = state.trackToVerticesMultiMap.equal_range(trk); + + // Allocate space in memory for the vector of compatibilities + trkToVtxCompatibilities.reserve(std::distance(begin, end)); + + for (auto it = begin; it != end; ++it) { + // it->first corresponds to trk, it->second to one of its associated + // vertices trkToVtxCompatibilities.push_back( - state.tracksAtVerticesMap.at(std::make_pair(trk, vtxIter->second)) + state.tracksAtVerticesMap.at(std::make_pair(trk, it->second)) .vertexCompatibility); } @@ -321,11 +339,10 @@ template bool Acts::AdaptiveMultiVertexFitter< input_track_t, linearizer_t>::checkSmallShift(State& state) const { for (auto vtx : state.vertexCollection) { - Vector3 diff = state.vtxInfoMap[vtx].oldPosition.template head<3>() - - vtx->fullPosition().template head<3>(); - SquareMatrix3 vtxWgt = - (vtx->fullCovariance().template block<3, 3>(0, 0)).inverse(); - double relativeShift = diff.dot(vtxWgt * diff); + Vector3 diff = + state.vtxInfoMap[vtx].oldPosition.template head<3>() - vtx->position(); + const SquareMatrix3& vtxCov = vtx->covariance(); + double relativeShift = diff.dot(vtxCov.inverse() * diff); if (relativeShift > m_cfg.maxRelativeShift) { return false; } diff --git a/Core/src/Detector/CylindricalContainerBuilder.cpp b/Core/src/Detector/CylindricalContainerBuilder.cpp index e4784790fd4..bbeb67a7d04 100644 --- a/Core/src/Detector/CylindricalContainerBuilder.cpp +++ b/Core/src/Detector/CylindricalContainerBuilder.cpp @@ -9,6 +9,8 @@ #include "Acts/Detector/CylindricalContainerBuilder.hpp" #include "Acts/Detector/DetectorComponents.hpp" +#include "Acts/Detector/DetectorVolumeBuilder.hpp" +#include "Acts/Detector/VolumeStructureBuilder.hpp" #include "Acts/Detector/detail/CylindricalDetectorHelper.hpp" #include "Acts/Detector/interface/IGeometryIdGenerator.hpp" #include "Acts/Detector/interface/IRootVolumeFinderBuilder.hpp" @@ -114,6 +116,50 @@ Acts::Experimental::CylindricalContainerBuilder::CylindricalContainerBuilder( } } +Acts::Experimental::CylindricalContainerBuilder::CylindricalContainerBuilder( + const Acts::Experimental::Blueprint::Node& bpNode, + Acts::Logging::Level logLevel) + : IDetectorComponentBuilder(), + m_logger(getDefaultLogger(bpNode.name + "_cont", logLevel)) { + if (bpNode.boundsType != VolumeBounds::BoundsType::eCylinder) { + throw std::invalid_argument( + "CylindricalContainerBuilder: boundary type must be cylinder - for " + "building from a blueprint node."); + } + + std::vector> builders; + for (const auto& child : bpNode.children) { + if (child->isLeaf()) { + // Volume structure + VolumeStructureBuilder::Config vsCfg; + vsCfg.transform = child->transform; + vsCfg.boundsType = child->boundsType; + vsCfg.boundValues = child->boundaryValues; + vsCfg.auxiliary = "*** acts auto-generated shape builder ***"; + auto vsBuilder = std::make_shared( + vsCfg, getDefaultLogger(child->name + "_shape", logLevel)); + // Detector volume builder + DetectorVolumeBuilder::Config dvCfg; + dvCfg.name = child->name; + dvCfg.externalsBuilder = vsBuilder; + dvCfg.internalsBuilder = child->internalsBuilder; + dvCfg.auxiliary = "*** acts auto-generated volume builder ***"; + // Add the builder + m_cfg.builders.push_back(std::make_shared( + dvCfg, getDefaultLogger(child->name, logLevel))); + } else { + // This evokes the recursive stepping down the tree + m_cfg.builders.push_back( + std::make_shared(*child, logLevel)); + } + } + + m_cfg.binning = bpNode.binning; + m_cfg.auxiliary = "*** acts auto-generated from proxy ***"; + m_cfg.geoIdGenerator = bpNode.geoIdGenerator; + m_cfg.rootVolumeFinderBuilder = bpNode.rootVolumeFinderBuilder; +} + Acts::Experimental::DetectorComponent Acts::Experimental::CylindricalContainerBuilder::construct( const GeometryContext& gctx) const { diff --git a/Core/src/Detector/DetectorVolumeBuilder.cpp b/Core/src/Detector/DetectorVolumeBuilder.cpp index fecc0326eda..8ffcbceed3f 100644 --- a/Core/src/Detector/DetectorVolumeBuilder.cpp +++ b/Core/src/Detector/DetectorVolumeBuilder.cpp @@ -42,7 +42,7 @@ Acts::Experimental::DetectorVolumeBuilder::construct( if (not m_cfg.auxiliary.empty()) { ACTS_DEBUG(m_cfg.auxiliary); } - ACTS_DEBUG("Building a volume with name " << m_cfg.name); + ACTS_DEBUG("Building a volume with name '" << m_cfg.name << "'."); // Get transform and bounds from the volume auto [transform, bounds, portalGenerator] = diff --git a/Core/src/Detector/VolumeStructureBuilder.cpp b/Core/src/Detector/VolumeStructureBuilder.cpp index 6a8eacec776..26e575fe7ee 100644 --- a/Core/src/Detector/VolumeStructureBuilder.cpp +++ b/Core/src/Detector/VolumeStructureBuilder.cpp @@ -17,6 +17,7 @@ #include "Acts/Geometry/TrapezoidVolumeBounds.hpp" #include "Acts/Utilities/Enumerate.hpp" #include "Acts/Utilities/Helpers.hpp" +#include "Acts/Utilities/StringHelpers.hpp" Acts::Experimental::VolumeStructureBuilder::VolumeStructureBuilder( const Acts::Experimental::VolumeStructureBuilder::Config& cfg, @@ -140,6 +141,10 @@ Acts::Experimental::VolumeStructureBuilder::construct( boundValues.push_back(M_PI); boundValues.push_back(0.); } + ACTS_VERBOSE(" - cylindrical shape with [iR, oR, hZ, sPhi, mPhi] = " + << boundValues[0] << ", " << boundValues[1] << ", " + << boundValues[2] << ", " << boundValues[3] << ", " + << boundValues[4]); auto bArray = to_array(boundValues); volumeBounds = std::make_unique(bArray); @@ -175,8 +180,14 @@ Acts::Experimental::VolumeStructureBuilder::construct( default: break; } + + Transform3 fTransform = m_cfg.transform * eTransform; + ACTS_VERBOSE(" - translation: " << Acts::toString(fTransform.translation())); + if (not fTransform.rotation().isApprox( + Acts::Transform3::Identity().rotation())) { + ACTS_VERBOSE(" - rotation: " << Acts::toString(fTransform.rotation())); + } // Return the transform, the volume bounds, and some default portal // generators - return {Transform3(m_cfg.transform * eTransform), std::move(volumeBounds), - defaultPortalGenerator()}; + return {fTransform, std::move(volumeBounds), defaultPortalGenerator()}; } diff --git a/Core/src/Detector/detail/BlueprintHelper.cpp b/Core/src/Detector/detail/BlueprintHelper.cpp index 7e240f0efae..1bde8053b3a 100644 --- a/Core/src/Detector/detail/BlueprintHelper.cpp +++ b/Core/src/Detector/detail/BlueprintHelper.cpp @@ -97,7 +97,8 @@ void Acts::Experimental::detail::BlueprintHelper::fillGaps( unsigned int igap = 0; for (auto& child : node.children) { auto [neg, pos] = cylEndpointsZ(*child); - if ((neg - negC).norm() > s_onSurfaceTolerance) { + ActsScalar gapSpan = (neg - negC).norm(); + if (gapSpan > s_onSurfaceTolerance) { // Fill a gap node auto gapName = node.name + "_gap_" + std::to_string(igap); auto gapTransform = Transform3::Identity(); @@ -105,7 +106,7 @@ void Acts::Experimental::detail::BlueprintHelper::fillGaps( gapTransform.translate(0.5 * (neg + negC)); auto gap = std::make_unique( gapName, gapTransform, VolumeBounds::eCylinder, - std::vector{cInnerR, cOuterR, negC.z() - neg.z()}); + std::vector{cInnerR, cOuterR, 0.5 * gapSpan}); gaps.push_back(std::move(gap)); ++igap; } @@ -113,7 +114,8 @@ void Acts::Experimental::detail::BlueprintHelper::fillGaps( negC = pos; } // Check if a last one needs to be filled - if ((negC - posC).norm() > s_onSurfaceTolerance) { + ActsScalar gapSpan = (negC - posC).norm(); + if (gapSpan > s_onSurfaceTolerance) { // Fill a gap node auto gapName = node.name + "_gap_" + std::to_string(igap); auto gapTransform = Transform3::Identity(); @@ -121,7 +123,7 @@ void Acts::Experimental::detail::BlueprintHelper::fillGaps( gapTransform.translate(0.5 * (negC + posC)); auto gap = std::make_unique( gapName, gapTransform, VolumeBounds::eCylinder, - std::vector{cInnerR, cOuterR, posC.z() - negC.z()}); + std::vector{cInnerR, cOuterR, 0.5 * gapSpan}); gaps.push_back(std::move(gap)); } diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp index d8a19c97641..0e9ff902d50 100644 --- a/Core/src/Geometry/TrackingVolume.cpp +++ b/Core/src/Geometry/TrackingVolume.cpp @@ -474,7 +474,7 @@ Acts::TrackingVolume::compatibleBoundaries( // The Limits: current, path & overstepping double pLimit = options.pathLimit; - double oLimit = options.overstepLimit; + double oLimit = 0; // Helper function to test intersection auto checkIntersection = diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp index de22b03914c..e96045bbb15 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp @@ -10,6 +10,8 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Utilities/VectorHelpers.hpp" +#include "ActsExamples/EventData/GeometryContainers.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" #include "ActsFatras/EventData/Particle.hpp" @@ -50,14 +52,38 @@ ActsExamples::ParticleSelector::ParticleSelector(const Config& config, ACTS_DEBUG("remove charged particles " << m_cfg.removeCharged); ACTS_DEBUG("remove neutral particles " << m_cfg.removeNeutral); ACTS_DEBUG("remove secondary particles " << m_cfg.removeSecondaries); + + // We only initialize this if we actually select on this + if (m_cfg.measurementsMin > 0 or + m_cfg.measurementsMax < std::numeric_limits::max()) { + m_inputMap.initialize(m_cfg.inputMeasurementParticlesMap); + ACTS_DEBUG("selection particle number of measurements [" + << m_cfg.measurementsMin << "," << m_cfg.measurementsMax << ")"); + } } ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( const AlgorithmContext& ctx) const { + using ParticlesMeasurmentMap = + boost::container::flat_multimap; + + // prepare input/ output types + const auto& inputParticles = m_inputParticles(ctx); + + // Make global particles measurement map if necessary + std::optional particlesMeasMap; + if (m_inputMap.isInitialized()) { + particlesMeasMap = invertIndexMultimap(m_inputMap(ctx)); + } + + std::size_t nInvalidCharge = 0; + std::size_t nInvalidMeasurementCount = 0; + // helper functions to select tracks - auto within = [](double x, double min, double max) { + auto within = [](auto x, auto min, auto max) { return (min <= x) and (x < max); }; + auto isValidParticle = [&](const ActsFatras::Particle& p) { const auto eta = Acts::VectorHelpers::eta(p.direction()); const auto phi = Acts::VectorHelpers::phi(p.direction()); @@ -67,7 +93,26 @@ ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( const bool validCharged = (p.charge() != 0) and not m_cfg.removeCharged; const bool validCharge = validNeutral or validCharged; const bool validSecondary = not m_cfg.removeSecondaries or !p.isSecondary(); - return validCharge and validSecondary and + + nInvalidCharge += static_cast(not validCharge); + + // default valid measurement count to true and only change if we have loaded + // the measurement particles map + bool validMeasurementCount = true; + if (particlesMeasMap) { + auto [b, e] = particlesMeasMap->equal_range(p.particleId()); + validMeasurementCount = + within(static_cast(std::distance(b, e)), + m_cfg.measurementsMin, m_cfg.measurementsMax); + + ACTS_VERBOSE("Found " << std::distance(b, e) << " measurements for " + << p.particleId()); + } + + nInvalidMeasurementCount += + static_cast(not validMeasurementCount); + + return validCharge and validSecondary and validMeasurementCount and within(p.transverseMomentum(), m_cfg.ptMin, m_cfg.ptMax) and within(std::abs(eta), m_cfg.absEtaMin, m_cfg.absEtaMax) and within(eta, m_cfg.etaMin, m_cfg.etaMax) and @@ -79,9 +124,6 @@ ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( within(p.mass(), m_cfg.mMin, m_cfg.mMax); }; - // prepare input/ output types - const auto& inputParticles = m_inputParticles(ctx); - SimParticleContainer outputParticles; outputParticles.reserve(inputParticles.size()); @@ -97,6 +139,9 @@ ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( ACTS_DEBUG("event " << ctx.eventNumber << " selected " << outputParticles.size() << " from " << inputParticles.size() << " particles"); + ACTS_DEBUG("filtered out because of charge: " << nInvalidCharge); + ACTS_DEBUG("filtered out because of measurement count: " + << nInvalidMeasurementCount); m_outputParticles(ctx, std::move(outputParticles)); return ProcessCode::SUCCESS; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp index 92ce7445dc6..d6c9fba3bc9 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp @@ -12,7 +12,10 @@ #pragma once +#include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/Index.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Framework/IAlgorithm.hpp" @@ -30,6 +33,8 @@ class ParticleSelector final : public IAlgorithm { struct Config { /// The input particles collection. std::string inputParticles; + /// Input measurement particles map (Optional) + std::string inputMeasurementParticlesMap; /// The output particles collection. std::string outputParticles; // Minimum/maximum distance from the origin in the transverse plane. @@ -54,6 +59,9 @@ class ParticleSelector final : public IAlgorithm { // Rest mass cuts double mMin = 0; double mMax = std::numeric_limits::infinity(); + /// Measurement number cuts + std::size_t measurementsMin = 0; + std::size_t measurementsMax = std::numeric_limits::max(); /// Remove charged particles. bool removeCharged = false; /// Remove neutral particles. @@ -74,6 +82,8 @@ class ParticleSelector final : public IAlgorithm { Config m_cfg; ReadDataHandle m_inputParticles{this, "InputParticles"}; + ReadDataHandle> m_inputMap{ + this, "InputMeasurementParticlesMap"}; WriteDataHandle m_outputParticles{this, "OutputParticles"}; diff --git a/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp new file mode 100644 index 00000000000..b9042452399 --- /dev/null +++ b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp @@ -0,0 +1,95 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" +#include "ActsFatras/EventData/Barcode.hpp" + +#include +#include + +namespace ActsExamples { + +/// Simple algorithm, that allows to extract a subset of a +/// measurment-particles-map. +/// This allows to conveniently work on subsets of the geometry. +/// +class MeasurementMapSelector final : public IAlgorithm { + using Map = IndexMultimap; + + public: + struct Config { + /// Input spacepoints collection. + std::string inputMeasurementParticleMap; + + /// Input source links + std::string inputSourceLinks; + + /// Output protoTracks collection. + std::string outputMeasurementParticleMap; + + /// What spacepoints to keep + std::vector geometrySelection; + }; + + /// Constructor of the track finding algorithm + /// + /// @param cfg is the config struct to configure the algorithm + /// @param level is the logging level + MeasurementMapSelector(Config cfg, Acts::Logging::Level lvl) + : IAlgorithm("MeasurementMapSelector", lvl), m_cfg(cfg) { + m_inputSourceLinks.initialize(m_cfg.inputSourceLinks); + m_inputMap.initialize(m_cfg.inputMeasurementParticleMap); + m_outputMap.initialize(m_cfg.outputMeasurementParticleMap); + } + + virtual ~MeasurementMapSelector() = default; + + /// Filter the measurements + /// + /// @param ctx is the algorithm context that holds event-wise information + /// @return a process code to steer the algorithm flow + ActsExamples::ProcessCode execute( + const ActsExamples::AlgorithmContext& ctx) const final { + const auto& inputSourceLinks = m_inputSourceLinks(ctx); + const auto& inputMap = m_inputMap(ctx); + + Map outputMap; + + for (const auto geoId : m_cfg.geometrySelection) { + auto range = selectLowestNonZeroGeometryObject(inputSourceLinks, geoId); + for (const auto& sl : range) { + const auto [begin, end] = inputMap.equal_range(sl.index()); + outputMap.insert(begin, end); + } + } + + m_outputMap(ctx, std::move(outputMap)); + + return ProcessCode::SUCCESS; + } + + const Config& config() const { return m_cfg; } + + private: + // configuration + Config m_cfg; + + ReadDataHandle m_inputSourceLinks{ + this, "InputSourceLinks"}; + ReadDataHandle m_inputMap{this, "InputMeasurementParticleMap"}; + WriteDataHandle m_outputMap{this, "OutputMeasurementParticleMap"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Io/Root/src/RootTrajectoryStatesWriter.cpp b/Examples/Io/Root/src/RootTrajectoryStatesWriter.cpp index 172b1de3a98..7e14671d540 100644 --- a/Examples/Io/Root/src/RootTrajectoryStatesWriter.cpp +++ b/Examples/Io/Root/src/RootTrajectoryStatesWriter.cpp @@ -545,12 +545,11 @@ ActsExamples::ProcessCode ActsExamples::RootTrajectoryStatesWriter::writeT( m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME); // track parameters error + // MARK: fpeMaskBegin(FLTINV, 1, #2348) m_err_eLOC0[ipar].push_back( - std::sqrt(covariance( // MARK: fpeMask(FLTINV, 1, #2348) - Acts::eBoundLoc0, Acts::eBoundLoc0))); + std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0))); m_err_eLOC1[ipar].push_back( - std::sqrt(covariance( // MARK: fpeMask(FLTINV, 1, #2348) - Acts::eBoundLoc1, Acts::eBoundLoc1))); + std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1))); m_err_ePHI[ipar].push_back( std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi))); m_err_eTHETA[ipar].push_back( @@ -559,6 +558,7 @@ ActsExamples::ProcessCode ActsExamples::RootTrajectoryStatesWriter::writeT( std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP))); m_err_eT[ipar].push_back( std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime))); + // MARK: fpeMaskEnd(FLTINV) // track parameters pull m_pull_eLOC0[ipar].push_back( diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 27e2f3731d5..2724d04933b 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -1008,6 +1008,7 @@ def addTruthTrackingGsf( calibrator=acts.examples.makePassThroughCalibrator(), ) s.addAlgorithm(gsfAlg) + s.addWhiteboardAlias("tracks", gsfAlg.config.outputTracks) trackConverter = acts.examples.TracksToTrajectories( level=customLogLevel(), diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py index 8faf65e6f5a..f0a9854337e 100644 --- a/Examples/Python/python/acts/examples/simulation.py +++ b/Examples/Python/python/acts/examples/simulation.py @@ -40,11 +40,12 @@ "absEta", # (min,max) "pt", # (min,max) "m", # (min,max) + "measurements", # (min,max) "removeCharged", # bool "removeNeutral", # bool "removeSecondaries", # bool ], - defaults=[(None, None)] * 8 + [None] * 3, + defaults=[(None, None)] * 9 + [None] * 3, ) @@ -333,6 +334,7 @@ def addParticleSelection( config: ParticleSelectorConfig, inputParticles: str, outputParticles: str, + inputMeasurementParticlesMap: str = "", logLevel: Optional[acts.logging.Level] = None, ) -> None: """ @@ -370,6 +372,8 @@ def addParticleSelection( ptMax=config.pt[1], mMin=config.m[0], mMax=config.m[1], + measurementsMin=config.measurements[0], + measurementsMax=config.measurements[1], removeCharged=config.removeCharged, removeNeutral=config.removeNeutral, removeSecondaries=config.removeSecondaries, @@ -377,6 +381,7 @@ def addParticleSelection( level=customLogLevel(), inputParticles=inputParticles, outputParticles=outputParticles, + inputMeasurementParticlesMap=inputMeasurementParticlesMap, ) ) diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index 66eab7d4157..897a21c98e1 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -26,6 +26,7 @@ #include "ActsExamples/TrackFinding/SpacePointMaker.hpp" #include "ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp" #include "ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp" +#include "ActsExamples/Utilities/MeasurementMapSelector.hpp" #include "ActsExamples/Utilities/PrototracksToSeeds.hpp" #include "ActsExamples/Utilities/SeedsToPrototracks.hpp" #include "ActsExamples/Utilities/TracksToTrajectories.hpp" @@ -377,6 +378,11 @@ void addTrackFinding(Context& ctx) { ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::PrototracksToSeeds, mex, "PrototracksToSeeds", inputProtoTracks, inputSpacePoints, outputSeeds, outputProtoTracks); + + ACTS_PYTHON_DECLARE_ALGORITHM( + ActsExamples::MeasurementMapSelector, mex, "MeasurementMapSelector", + inputMeasurementParticleMap, inputSourceLinks, + outputMeasurementParticleMap, geometrySelection); } } // namespace Acts::Python diff --git a/Examples/Python/src/TruthTracking.cpp b/Examples/Python/src/TruthTracking.cpp index af4f36a2b4c..5b37e7b3753 100644 --- a/Examples/Python/src/TruthTracking.cpp +++ b/Examples/Python/src/TruthTracking.cpp @@ -106,6 +106,7 @@ void addTruthTracking(Context& ctx) { 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); @@ -123,6 +124,8 @@ void addTruthTracking(Context& ctx) { ACTS_PYTHON_MEMBER(mMax); ACTS_PYTHON_MEMBER(ptMin); ACTS_PYTHON_MEMBER(ptMax); + ACTS_PYTHON_MEMBER(measurementsMin); + ACTS_PYTHON_MEMBER(measurementsMax); ACTS_PYTHON_MEMBER(removeCharged); ACTS_PYTHON_MEMBER(removeNeutral); ACTS_PYTHON_MEMBER(removeSecondaries); @@ -136,6 +139,8 @@ void addTruthTracking(Context& ctx) { pythonRangeProperty(c, "absEta", &Config::absEtaMin, &Config::absEtaMax); pythonRangeProperty(c, "m", &Config::mMin, &Config::mMax); pythonRangeProperty(c, "pt", &Config::ptMin, &Config::ptMax); + pythonRangeProperty(c, "measurements", &Config::measurementsMin, + &Config::measurementsMax); } { diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 93ae06d6aee..6b32b90944a 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -22,22 +22,22 @@ test_itk_seeding__particles_final.root: 0d38f1a4ede5fb86844c156c002ada340011e140 test_itk_seeding__particles_initial.root: 88315e93ed4cb5d40a8721502048a9d1fc100e0a7d504e25fd4502c8302f1578 test_propagation__propagation_steps.root: 174301b25784dbb881196b658f2d7f99c0a2ea688a0129e6110fc19aa5cf8e54 test_material_recording__geant4_material_tracks.root: e411152d370775463c22b19a351dfc7bfe40b51985e10a7c1a010aebde80715d -test_truth_tracking_kalman[generic-0.0]__trackstates_fitter.root: c8454669d2ce93296ef7be86cd9e900c0e57154c30199dcb92790983472d6847 -test_truth_tracking_kalman[generic-0.0]__tracksummary_fitter.root: e9c9c700ef3ebd6aafd3bbb95eaf6e24bc38e90c00f207f1e94a33b2271b0607 +test_truth_tracking_kalman[generic-0.0]__trackstates_fitter.root: 559321c478b63859b5d1b687e73d79655d1b3b12ebc34c2bd5ffa1049827b020 +test_truth_tracking_kalman[generic-0.0]__tracksummary_fitter.root: 74ed2ec7cb70eec8767d8023af8bdb1b488c8b80d6299f2c234b14f71f84d2e5 test_truth_tracking_kalman[generic-0.0]__performance_track_finder.root: ada9cda2ec3c648b144bdaa081d6eff923c80f3d484c4196006e847428cf67a8 -test_truth_tracking_kalman[generic-1000.0]__trackstates_fitter.root: 23bcd151eda147e9270924f76a422e807d8d38d75b6193e4ffbc442ab94e65c1 -test_truth_tracking_kalman[generic-1000.0]__tracksummary_fitter.root: 80672652eb750f117013ef22c9ff75fa09ac08a0f38dd788249fa40395e38e47 +test_truth_tracking_kalman[generic-1000.0]__trackstates_fitter.root: 0befe1ec82e461f64bac9d44f6b0d96f01263b0d6c094afdb4b1cc3d2c965095 +test_truth_tracking_kalman[generic-1000.0]__tracksummary_fitter.root: e8b1a8a7097cfc3e3973ea7eb8210126a223dd2a771c0bc2f8a2af4e152908ba test_truth_tracking_kalman[generic-1000.0]__performance_track_finder.root: ada9cda2ec3c648b144bdaa081d6eff923c80f3d484c4196006e847428cf67a8 -test_truth_tracking_kalman[odd-0.0]__trackstates_fitter.root: bcd562b4d55374e892f2bdbb742db4c7761637d523ec8879897c091901e16a95 -test_truth_tracking_kalman[odd-0.0]__tracksummary_fitter.root: 9b9bffa011bb22f54c3c988cf9b66bc84c848ca88513e0a8025de8861924866e +test_truth_tracking_kalman[odd-0.0]__trackstates_fitter.root: d197c1e2ea724cabbbe1a21673034810451cadde504b1f9e813ed95b74603cd9 +test_truth_tracking_kalman[odd-0.0]__tracksummary_fitter.root: 71e1b3d4c466df61edd758757192a6a74a4638be0114cc45bb306b994a4a8ece test_truth_tracking_kalman[odd-0.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593 -test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: ef499a4f0d0c81c9e39ec513708ff16a11003f668192fa2c7ba1fc24aba0d7c8 -test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: 6a4dda749c0b2f89aaf74e41b25938da5f8a979bac3b895b7ab15004a72e11e5 +test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: 707c0784deac8d7c63de38542f08f66dbe375f502b26f7d0aeaccf976d4bcf67 +test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: 7b15fa2d69a26ba0a9c0ca0626c36dd904e00fb31748d9697507c8e372b49b7b test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593 -test_truth_tracking_gsf[generic]__trackstates_gsf.root: 849f09ea6f45bcaad56172bd5d17cf201076f4f2760638debfbbc9b31b90b636 -test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 10c7f080b5d9de9eebd9dcd35344322a30bdb39a21ae5e5e28e3b5b99814bb12 -test_truth_tracking_gsf[odd]__trackstates_gsf.root: bf4e15ba68937c3cf5138245939ee2cd9abaa04dafd70cf09bc11b0f6edcd79f -test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 912c3e8f6fc71d7312e18e4e571ed91b2624c929b63a7d7ab407857fc3157118 +test_truth_tracking_gsf[generic]__trackstates_gsf.root: f7b0a2365519725d6836a433f2ed2d4184b83f0d455a80aae92977704abfdf8f +test_truth_tracking_gsf[generic]__tracksummary_gsf.root: e631394d26d1d42fd123229af7fce03f899c82549f878b36369745e657840231 +test_truth_tracking_gsf[odd]__trackstates_gsf.root: 16c480cf80ccedb5327cdfc534387bf8ad0a331d423c5f707a945efa8e24aa8f +test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 11c378ecb70ebb375857236d2ea3a0aaa5c71bf5ad5b4360c0d0f107c5347ed6 test_particle_gun__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec test_material_mapping__material-map_tracks.root: 2deac7a48ff1185ba3889cfd4cc0bd0a6de57293048dfd6766f3c3907232e45e test_material_mapping__propagation-material.root: 84b04ebd5721550867f0f193b5eb4e9f3f205df542cb46090d320be6bff565c5 @@ -49,22 +49,22 @@ test_digitization_example_input[smeared]__particles.root: 8549ba6e20338004ab8ba2 test_digitization_example_input[smeared]__measurements.root: 2e4fd9d3e6244e53486496e024f147e9ab2dcadb07453aa926d7ebad473089a5 test_digitization_example_input[geometric]__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec test_digitization_example_input[geometric]__measurements.root: 2dfd1889028fe3a0964fe9f91bc617987fbecd122256a9bd07433da1931323bc -test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: c95a1e289b588477cc3946ab68f0aa396ca515edc18471cf6e964e806ccd2980 -test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: 79d0c8f6a636e1ba19abcc3fb786735b1ea7b89290b01957f50927ba3f2f3481 +test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: 61c958a047d431133f332883b975ce522d3e7780e8f4cc57eaddbe8c3550cd6a +test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: 004a69290c04ea7858d8040f043ac3ca7c2be55b99a22c4ea7436c13421cd0dc test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0e0676ffafdb27112fbda50d1cf627859fa745760f98073261dcf6db3f2f991e -test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: ef8ff1d65a63d61ca1c6004960d7ff537cfd7d8fbd3920cf9411064f8d48ef9b -test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 149cb73cef1526dba115192aa2b66c4de6dcd9c7f81c73f1d54190d543b1e41f +test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: 346462706c4cd884708bc18814ec66c4f62f3fa358cf591491220863d1e7ab4f +test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 744852871167a083f10cc95b12baf51816db139d032ce839839433a4a416433b test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e -test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 10eb393cb419ff81b34c627db84bcfe89134bd58838325d894e184e00a75a9c4 -test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: e8e56e965799191c643eb4339f36044edf272807656fa73014cab4e78b1f1619 -test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 9e15151308c0983fc681cbf60b478542814235512143755c092b921c7a9c19e3 -test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: c3feeabe87c49e7e8140f420af936767f6171276cd431dfb17a04f721c1895b0 +test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 8a2e5e7751b2b1497ab77902b99ce9298b1b143bd04d156ba7639899b4f86042 +test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: c3868964a1b9197b295e47def41cfad3b120f9617f323d9a6164fcf08641a807 +test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: e8166fdd3dace260e4fb7dd71963ddd7b44960e6108c1f68d407d65c06b8f152 +test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 309fddc19969549807ef5919eb8c262264bcbd37d0791321ddf2b7b26610207e test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f -test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 7c11478d9d6ab190b3e1ff447f54f0ba680d5f9b70ee1db158d69db272627b41 -test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: d747ca518ce5a9ad01827a87052c6590a6688eb61ada213c012250cfa7e3e762 +test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 2cf76733df5d7eb29fcda9eb7525d5508b9277eff932b2cf71910824031d9402 +test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 3eb7e180bac29a82add4bd8e795e404cbe57d084f4538dc34299dce9cbc338e0 test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575 test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 37c95a3f1ad5a4b606c0d902e58ac83b814d269a8d614fe0a2a986ec58838a88 -test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: ff43def36609f25318ca901ff1dd58e78a29b8ded8fbe4a6aa0679763a23e14d +test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 0b6755e54161371b0374211c8b4f054c6158d0e44f47bcfd6602ce12bf14f6bc test_vertex_fitting_reading[Truth-False-100]__performance_vertexing.root: 76ef6084d758dfdfc0151ddec2170e12d73394424e3dac4ffe46f0f339ec8293 test_vertex_fitting_reading[Iterative-False-100]__performance_vertexing.root: 60372210c830a04f95ceb78c6c68a9b0de217746ff59e8e73053750c837b57eb test_vertex_fitting_reading[Iterative-True-100]__performance_vertexing.root: e34f217d524a5051dbb04a811d3407df3ebe2cc4bb7f54f6bda0847dbd7b52c3 diff --git a/Examples/Python/tests/test_algorithms.py b/Examples/Python/tests/test_algorithms.py index 3db13b2c9ff..b6232035819 100644 --- a/Examples/Python/tests/test_algorithms.py +++ b/Examples/Python/tests/test_algorithms.py @@ -1,7 +1,5 @@ import pytest -import acts - from acts.examples import ( TutorialVertexFinderAlgorithm, AdaptiveMultiVertexFinderAlgorithm, diff --git a/Examples/Python/tests/test_detectors.py b/Examples/Python/tests/test_detectors.py index 1a01809e5dd..ee2053a4358 100644 --- a/Examples/Python/tests/test_detectors.py +++ b/Examples/Python/tests/test_detectors.py @@ -1,5 +1,3 @@ -from pathlib import Path - import pytest from helpers import dd4hepEnabled diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py index 0f3bbb24c5c..7b7a918eaf1 100644 --- a/Examples/Python/tests/test_examples.py +++ b/Examples/Python/tests/test_examples.py @@ -20,8 +20,6 @@ exatrkxEnabled, onnxEnabled, AssertCollectionExistsAlg, - isCI, - doHashChecks, failure_threshold, ) @@ -33,7 +31,6 @@ Sequencer, GenericDetector, AlignedDetector, - RootParticleWriter, ) from acts.examples.odd import getOpenDataDetector @@ -569,7 +566,6 @@ def test_truth_tracking_kalman( root_files = [ ("trackstates_fitter.root", "trackstates", 19), ("tracksummary_fitter.root", "tracksummary", 10), - ("performance_track_finder.root", "track_finder_tracks", 19), ("performance_track_fitter.root", None, -1), ] @@ -608,11 +604,6 @@ def test_truth_tracking_gsf(tmp_path, assert_root_hash, detector_config): events=10, numThreads=1, fpeMasks=[ - ( - "Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp:177", - acts.FpeType.FLTUND, - 1, - ), ( "Core/include/Acts/TrackFitting/detail/GsfUtils.hpp:197", acts.FpeType.FLTUND, diff --git a/Examples/Python/tests/test_magnetic_field.py b/Examples/Python/tests/test_magnetic_field.py index 9f2979cccfb..3c1aae5c537 100644 --- a/Examples/Python/tests/test_magnetic_field.py +++ b/Examples/Python/tests/test_magnetic_field.py @@ -1,4 +1,3 @@ -from typing import Type import pytest import acts diff --git a/Examples/Python/tests/test_writer.py b/Examples/Python/tests/test_writer.py index e050195d34f..095b53638a7 100644 --- a/Examples/Python/tests/test_writer.py +++ b/Examples/Python/tests/test_writer.py @@ -1,4 +1,3 @@ -from typing import Type import os import inspect from pathlib import Path @@ -748,33 +747,10 @@ def test_edm4hep_tracks_writer(tmp_path): return from podio.root_io import Reader - from podio.frame import Frame import cppyy reader = Reader(str(out)) - expected = [ - (31.986961364746094, 30, 16), - (28.64777374267578, 30, 16), - (11.607606887817383, 22, 12), - (5.585886001586914, 22, 12), - (20.560943603515625, 20, 11), - (28.742727279663086, 28, 15), - (27.446802139282227, 22, 12), - (30.82959747314453, 24, 13), - (24.671127319335938, 26, 14), - (16.479907989501953, 20, 11), - (10.594233512878418, 22, 12), - (25.174715042114258, 28, 15), - (27.9674072265625, 26, 14), - (4.3012871742248535, 22, 12), - (20.492422103881836, 22, 12), - (27.92759132385254, 24, 13), - (14.514887809753418, 22, 12), - (12.876864433288574, 22, 12), - (12.951473236083984, 26, 14), - ] - actual = [] for frame in reader.get("events"): @@ -805,5 +781,5 @@ def test_edm4hep_tracks_writer(tmp_path): assert rp.x == 0.0 assert rp.y == 0.0 assert rp.z == 0.0 - assert abs(perigee.D0) < 1e-1 - assert abs(perigee.Z0) < 1 + assert abs(perigee.D0) < 1e0 + assert abs(perigee.Z0) < 1e1 diff --git a/Examples/Scripts/Python/truth_tracking_gsf.py b/Examples/Scripts/Python/truth_tracking_gsf.py index 01565bfdd3a..64328e61abf 100755 --- a/Examples/Scripts/Python/truth_tracking_gsf.py +++ b/Examples/Scripts/Python/truth_tracking_gsf.py @@ -1,36 +1,36 @@ #!/usr/bin/env python3 from pathlib import Path -from typing import Optional, Union - -from acts.examples import Sequencer, GenericDetector, RootParticleReader +from typing import Optional import acts +import acts.examples u = acts.UnitConstants def runTruthTrackingGsf( - trackingGeometry, + trackingGeometry: acts.TrackingGeometry, + field: acts.MagneticFieldProvider, digiConfigFile: Path, - field, outputDir: Path, - outputCsv=True, inputParticlePath: Optional[Path] = None, decorators=[], - s=None, + s: acts.examples.Sequencer = None, ): from acts.examples.simulation import ( addParticleGun, + ParticleConfig, EtaConfig, PhiConfig, - ParticleConfig, + MomentumConfig, addFatras, addDigitization, ) from acts.examples.reconstruction import ( addSeeding, SeedingAlgorithm, + TruthSeedRanges, addTruthTrackingGsf, ) @@ -47,10 +47,15 @@ def runTruthTrackingGsf( if inputParticlePath is None: addParticleGun( s, - EtaConfig(-2.0, 2.0), - ParticleConfig(4, acts.PdgParticle.eElectron, True), + ParticleConfig(num=1, pdg=acts.PdgParticle.eElectron, randomizeCharge=True), + EtaConfig(-3.0, 3.0, uniform=True), + MomentumConfig(1.0 * u.GeV, 100.0 * u.GeV, transverse=True), PhiConfig(0.0, 360.0 * u.degree), - multiplicity=2, + vtxGen=acts.examples.GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4(0, 0, 0, 0), + ), + multiplicity=1, rnd=rnd, ) else: @@ -59,7 +64,7 @@ def runTruthTrackingGsf( ) assert inputParticlePath.exists() s.addReader( - RootParticleReader( + acts.examples.RootParticleReader( level=acts.logging.INFO, filePath=str(inputParticlePath.resolve()), particleCollection="particles_input", @@ -87,22 +92,41 @@ def runTruthTrackingGsf( s, trackingGeometry, field, + rnd=rnd, + inputParticles="particles_input", seedingAlgorithm=SeedingAlgorithm.TruthSmeared, particleHypothesis=acts.ParticleHypothesis.electron, + truthSeedRanges=TruthSeedRanges( + pt=(1 * u.GeV, None), + nHits=(7, None), + ), ) - truthTrkFndAlg = acts.examples.TruthTrackFinder( - level=acts.logging.INFO, - inputParticles="truth_seeds_selected", - inputMeasurementParticlesMap="measurement_particles_map", - outputProtoTracks="prototracks", + addTruthTrackingGsf( + s, + trackingGeometry, + field, ) - s.addAlgorithm(truthTrkFndAlg) - - addTruthTrackingGsf(s, trackingGeometry, field) + s.addAlgorithm( + acts.examples.TrackSelectorAlgorithm( + level=acts.logging.INFO, + inputTracks="tracks", + outputTracks="selected-tracks", + selectorConfig=acts.TrackSelector.Config( + minMeasurements=7, + ), + ) + ) + s.addAlgorithm( + acts.examples.TracksToTrajectories( + level=acts.logging.INFO, + inputTracks="selected-tracks", + outputTrajectories="trajectories-from-tracks", + ) + ) + s.addWhiteboardAlias("trajectories", "trajectories-from-tracks") - # Output s.addWriter( acts.examples.RootTrajectoryStatesWriter( level=acts.logging.INFO, @@ -141,7 +165,7 @@ def runTruthTrackingGsf( if "__main__" == __name__: srcdir = Path(__file__).resolve().parent.parent.parent.parent - detector, trackingGeometry, decorators = GenericDetector.create() + detector, trackingGeometry, decorators = acts.examples.GenericDetector.create() field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T)) @@ -151,11 +175,9 @@ def runTruthTrackingGsf( runTruthTrackingGsf( trackingGeometry=trackingGeometry, - decorators=decorators, field=field, digiConfigFile=srcdir / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json", - outputCsv=True, inputParticlePath=inputParticlePath, outputDir=Path.cwd(), ).run() diff --git a/Examples/Scripts/Python/truth_tracking_kalman.py b/Examples/Scripts/Python/truth_tracking_kalman.py index 93ed96831e2..b06acfceb34 100755 --- a/Examples/Scripts/Python/truth_tracking_kalman.py +++ b/Examples/Scripts/Python/truth_tracking_kalman.py @@ -12,17 +12,20 @@ def runTruthTrackingKalman( trackingGeometry: acts.TrackingGeometry, field: acts.MagneticFieldProvider, - outputDir: Path, digiConfigFile: Path, + outputDir: Path, + inputParticlePath: Optional[Path] = None, + decorators=[], directNavigation=False, reverseFilteringMomThreshold=0 * u.GeV, s: acts.examples.Sequencer = None, - inputParticlePath: Optional[Path] = None, ): from acts.examples.simulation import ( addParticleGun, - EtaConfig, ParticleConfig, + EtaConfig, + PhiConfig, + MomentumConfig, addFatras, addDigitization, ) @@ -37,17 +40,25 @@ def runTruthTrackingKalman( events=100, numThreads=-1, logLevel=acts.logging.INFO ) - rnd = acts.examples.RandomNumbers() + for d in decorators: + s.addContextDecorator(d) + + rnd = acts.examples.RandomNumbers(seed=42) outputDir = Path(outputDir) if inputParticlePath is None: addParticleGun( s, - EtaConfig(-2.0, 2.0), - ParticleConfig(2, acts.PdgParticle.eMuon, False), + ParticleConfig(num=1, pdg=acts.PdgParticle.eMuon, randomizeCharge=True), + EtaConfig(-3.0, 3.0, uniform=True), + MomentumConfig(1.0 * u.GeV, 100.0 * u.GeV, transverse=True), + PhiConfig(0.0, 360.0 * u.degree), + vtxGen=acts.examples.GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4(0, 0, 0, 0), + ), multiplicity=1, rnd=rnd, - outputDirRoot=outputDir, ) else: acts.logging.getLogger("Truth tracking example").info( @@ -55,7 +66,7 @@ def runTruthTrackingKalman( ) assert inputParticlePath.exists() s.addReader( - RootParticleReader( + acts.examples.RootParticleReader( level=acts.logging.INFO, filePath=str(inputParticlePath.resolve()), particleCollection="particles_input", @@ -83,11 +94,13 @@ def runTruthTrackingKalman( s, trackingGeometry, field, - seedingAlgorithm=SeedingAlgorithm.TruthSmeared, rnd=rnd, + inputParticles="particles_input", + seedingAlgorithm=SeedingAlgorithm.TruthSmeared, + particleHypothesis=acts.ParticleHypothesis.muon, truthSeedRanges=TruthSeedRanges( - pt=(500 * u.MeV, None), - nHits=(9, None), + pt=(1 * u.GeV, None), + nHits=(7, None), ), ) @@ -99,7 +112,25 @@ def runTruthTrackingKalman( reverseFilteringMomThreshold, ) - # Output + s.addAlgorithm( + acts.examples.TrackSelectorAlgorithm( + level=acts.logging.INFO, + inputTracks="tracks", + outputTracks="selected-tracks", + selectorConfig=acts.TrackSelector.Config( + minMeasurements=7, + ), + ) + ) + s.addAlgorithm( + acts.examples.TracksToTrajectories( + level=acts.logging.INFO, + inputTracks="selected-tracks", + outputTrajectories="trajectories-from-tracks", + ) + ) + s.addWhiteboardAlias("trajectories", "trajectories-from-tracks") + s.addWriter( acts.examples.RootTrajectoryStatesWriter( level=acts.logging.INFO, @@ -122,18 +153,6 @@ def runTruthTrackingKalman( ) ) - s.addWriter( - acts.examples.TrackFinderPerformanceWriter( - level=acts.logging.INFO, - inputProtoTracks="sorted_truth_particle_tracks" - if directNavigation - else "truth_particle_tracks", - inputParticles="truth_seeds_selected", - inputMeasurementParticlesMap="measurement_particles_map", - filePath=str(outputDir / "performance_track_finder.root"), - ) - ) - s.addWriter( acts.examples.TrackFitterPerformanceWriter( level=acts.logging.INFO, @@ -148,7 +167,6 @@ def runTruthTrackingKalman( if "__main__" == __name__: - srcdir = Path(__file__).resolve().parent.parent.parent.parent # detector, trackingGeometry, _ = getOpenDataDetector() diff --git a/Fatras/include/ActsFatras/EventData/Particle.hpp b/Fatras/include/ActsFatras/EventData/Particle.hpp index 7eb43ab1373..f10a30ede3d 100644 --- a/Fatras/include/ActsFatras/EventData/Particle.hpp +++ b/Fatras/include/ActsFatras/EventData/Particle.hpp @@ -11,7 +11,12 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/PdgParticle.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/ParticleHypothesis.hpp" +#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Utilities/VectorHelpers.hpp" #include "ActsFatras/EventData/Barcode.hpp" #include "ActsFatras/EventData/ProcessType.hpp" @@ -19,6 +24,8 @@ #include #include #include +#include +#include namespace ActsFatras { @@ -191,6 +198,10 @@ class Particle { } /// Unit three-direction, i.e. the normalized momentum three-vector. const Vector3 &direction() const { return m_direction; } + /// Polar angle. + Scalar theta() const { return Acts::VectorHelpers::theta(direction()); } + /// Azimuthal angle. + Scalar phi() const { return Acts::VectorHelpers::phi(direction()); } /// Absolute momentum in the x-y plane. Scalar transverseMomentum() const { return m_absMomentum * m_direction.segment<2>(Acts::eMom0).norm(); @@ -236,6 +247,43 @@ class Particle { /// Accumulated path within material measured in interaction lengths. constexpr Scalar pathInL0() const { return m_pathInL0; } + /// Set the reference surface. + /// + /// @param surface reference surface + Particle &setReferenceSurface(const Acts::Surface *surface) { + m_referenceSurface = surface; + return *this; + } + + /// Reference surface. + const Acts::Surface *referenceSurface() const { return m_referenceSurface; } + + /// Check if the particle has a reference surface. + bool hasReferenceSurface() const { return m_referenceSurface != nullptr; } + + /// Bound track parameters. + Acts::Result boundParameters( + const Acts::GeometryContext &gctx) const { + if (!hasReferenceSurface()) { + return Acts::Result::failure( + std::error_code()); + } + Acts::Result localResult = + m_referenceSurface->globalToLocal(gctx, position(), direction()); + if (!localResult.ok()) { + return localResult.error(); + } + Acts::BoundVector params; + params << localResult.value(), phi(), theta(), qOverP(), time(); + return Acts::BoundTrackParameters(referenceSurface()->getSharedPtr(), + params, std::nullopt, hypothesis()); + } + + Acts::CurvilinearTrackParameters curvilinearParameters() const { + return Acts::CurvilinearTrackParameters( + fourPosition(), direction(), qOverP(), std::nullopt, hypothesis()); + } + private: // identity, i.e. things that do not change over the particle lifetime. /// Particle identifier within the event. @@ -256,6 +304,8 @@ class Particle { // accumulated material Scalar m_pathInX0 = Scalar(0); Scalar m_pathInL0 = Scalar(0); + // reference surface + const Acts::Surface *m_referenceSurface{nullptr}; }; std::ostream &operator<<(std::ostream &os, const Particle &particle); diff --git a/Fatras/include/ActsFatras/Kernel/Simulation.hpp b/Fatras/include/ActsFatras/Kernel/Simulation.hpp index 4ce5a226dc9..6c0cc402a08 100644 --- a/Fatras/include/ActsFatras/Kernel/Simulation.hpp +++ b/Fatras/include/ActsFatras/Kernel/Simulation.hpp @@ -10,6 +10,7 @@ #include "Acts/EventData/Charge.hpp" #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp" +#include "Acts/EventData/TrackParameters.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "Acts/Propagator/AbortList.hpp" @@ -89,16 +90,23 @@ struct SingleParticleSimulation { actor.interactions = interactions; actor.selectHitSurface = selectHitSurface; actor.initialParticle = particle; - // use AnyCharge to be able to handle neutral and charged parameters - Acts::CurvilinearTrackParameters start( - particle.fourPosition(), particle.direction(), particle.qOverP(), - std::nullopt, particle.hypothesis()); - auto result = propagator.propagate(start, options); + + if (particle.hasReferenceSurface()) { + auto result = propagator.propagate( + particle.boundParameters(geoCtx).value(), options); + if (not result.ok()) { + return result.error(); + } + auto &value = result.value().template get(); + return std::move(value); + } + + auto result = + propagator.propagate(particle.curvilinearParameters(), options); if (not result.ok()) { return result.error(); } auto &value = result.value().template get(); - return std::move(value); } }; diff --git a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp index 40eafe60d94..267d073a099 100644 --- a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp +++ b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp @@ -92,6 +92,13 @@ struct SimulationActor { return; } + // check if we are still on the start surface and skip if so + if ((navigator.startSurface(state.navigation) != nullptr) && + (navigator.startSurface(state.navigation) == + navigator.currentSurface(state.navigation))) { + return; + } + // update the particle state first. this also computes the proper time which // needs the particle state from the previous step for reference. that means // this must happen for every step (not just on surface) and before @@ -99,11 +106,13 @@ struct SimulationActor { if (std::isnan(result.properTimeLimit)) { // first step is special: there is no previous state and we need to arm // the decay simulation for all future steps. - result.particle = makeParticle(initialParticle, stepper, state.stepping); + result.particle = + makeParticle(initialParticle, state, stepper, navigator); result.properTimeLimit = decay.generateProperTimeLimit(*generator, initialParticle); } else { - result.particle = makeParticle(result.particle, stepper, state.stepping); + result.particle = + makeParticle(result.particle, state, stepper, navigator); } // decay check. needs to happen at every step, not just on surfaces. @@ -174,7 +183,7 @@ struct SimulationActor { // apply abs in case `normal` and `before` produce an angle > 90° slab.scaleThickness(std::abs(cosIncidenceInv)); // run the interaction simulation - interact(slab, result); + interact(slab, result); // MARK: fpeMask(FLTUND, 1, #2346) } } } @@ -199,27 +208,34 @@ struct SimulationActor { after.qOverP(), after.time()); } - /// Construct the current particle state from the stepper state. - template - Particle makeParticle(const Particle &previous, const stepper_t &stepper, - const typename stepper_t::State &state) const { + /// Construct the current particle state from the propagation state. + template + Particle makeParticle(const Particle &previous, propagator_state_t &state, + stepper_t &stepper, navigator_t &navigator) const { // a particle can lose energy and thus its gamma factor is not a constant // of motion. since the stepper provides only the lab time, we need to // compute the change in proper time for each step separately. this assumes // that the gamma factor is constant over one stepper step. - const auto deltaLabTime = stepper.time(state) - previous.time(); + const auto deltaLabTime = stepper.time(state.stepping) - previous.time(); // proper-time = time / gamma = (1/gamma) * time // beta² = p²/E² // gamma = 1 / sqrt(1 - beta²) = sqrt(m² + p²) / m // 1/gamma = m / sqrt(m² + p²) = m / E const auto gammaInv = previous.mass() / previous.energy(); const auto properTime = previous.properTime() + gammaInv * deltaLabTime; + const Acts::Surface *currentSurface = nullptr; + if (navigator.currentSurface(state.navigation) != nullptr) { + currentSurface = navigator.currentSurface(state.navigation); + } // copy all properties and update kinematic state from stepper return Particle(previous) - .setPosition4(stepper.position(state), stepper.time(state)) - .setDirection(stepper.direction(state)) - .setAbsoluteMomentum(stepper.absoluteMomentum(state)) - .setProperTime(properTime); + .setPosition4(stepper.position(state.stepping), + stepper.time(state.stepping)) + .setDirection(stepper.direction(state.stepping)) + .setAbsoluteMomentum(stepper.absoluteMomentum(state.stepping)) + .setProperTime(properTime) + .setReferenceSurface(currentSurface); } /// Prepare limits and process selection for the next point-like interaction. diff --git a/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp b/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp index 346f4f07616..4e976f69882 100644 --- a/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp +++ b/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp @@ -281,13 +281,15 @@ std::array PhotonConversion::generateChildren( .setPosition4(photon.fourPosition()) .setDirection(childDirection) .setAbsoluteMomentum(momentum1) - .setProcess(ProcessType::ePhotonConversion), + .setProcess(ProcessType::ePhotonConversion) + .setReferenceSurface(photon.referenceSurface()), Particle(photon.particleId().makeDescendant(1), Acts::ePositron, 1_e, kElectronMass) .setPosition4(photon.fourPosition()) .setDirection(childDirection) .setAbsoluteMomentum(momentum2) - .setProcess(ProcessType::ePhotonConversion), + .setProcess(ProcessType::ePhotonConversion) + .setReferenceSurface(photon.referenceSurface()), }; return children; } diff --git a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp index e6e981bbdb0..a408b1ff97f 100644 --- a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp +++ b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp @@ -528,7 +528,8 @@ std::vector NuclearInteraction::convertParametersToParticles( p.setProcess(ProcessType::eNuclearInteraction) .setPosition4(initialParticle.fourPosition()) .setAbsoluteMomentum(momentum) - .setDirection(direction); + .setDirection(direction) + .setReferenceSurface(initialParticle.referenceSurface()); // Store the particle if (i == 0 && soft) { diff --git a/Fatras/src/Physics/BetheHeitler.cpp b/Fatras/src/Physics/BetheHeitler.cpp index cf7142e9cde..9ef32c94e3f 100644 --- a/Fatras/src/Physics/BetheHeitler.cpp +++ b/Fatras/src/Physics/BetheHeitler.cpp @@ -62,6 +62,7 @@ ActsFatras::Particle ActsFatras::BetheHeitler::bremPhoton( photon.setProcess(ActsFatras::ProcessType::eBremsstrahlung) .setPosition4(particle.fourPosition()) .setDirection(photonDirection) - .setAbsoluteMomentum(gammaE); + .setAbsoluteMomentum(gammaE) + .setReferenceSurface(particle.referenceSurface()); return photon; } diff --git a/Tests/UnitTests/Core/Detector/CMakeLists.txt b/Tests/UnitTests/Core/Detector/CMakeLists.txt index 2f3c8768b60..b5733893b67 100644 --- a/Tests/UnitTests/Core/Detector/CMakeLists.txt +++ b/Tests/UnitTests/Core/Detector/CMakeLists.txt @@ -1,6 +1,7 @@ add_unittest(Blueprint BlueprintTests.cpp) add_unittest(BlueprintHelper BlueprintHelperTests.cpp) add_unittest(CylindricalContainerBuilder CylindricalContainerBuilderTests.cpp) +add_unittest(CylindricalDetectorFromBlueprint CylindricalDetectorFromBlueprintTests.cpp) add_unittest(CylindricalDetectorHelper CylindricalDetectorHelperTests.cpp) add_unittest(GridAxisGenerators GridAxisGeneratorsTests.cpp) add_unittest(Detector DetectorTests.cpp) diff --git a/Tests/UnitTests/Core/Detector/CylindricalDetectorFromBlueprintTests.cpp b/Tests/UnitTests/Core/Detector/CylindricalDetectorFromBlueprintTests.cpp new file mode 100644 index 00000000000..6f62ee892f9 --- /dev/null +++ b/Tests/UnitTests/Core/Detector/CylindricalDetectorFromBlueprintTests.cpp @@ -0,0 +1,282 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#include "Acts/Detector/Blueprint.hpp" +#include "Acts/Detector/CylindricalContainerBuilder.hpp" +#include "Acts/Detector/DetectorBuilder.hpp" +#include "Acts/Detector/DetectorComponents.hpp" +#include "Acts/Detector/DetectorVolume.hpp" +#include "Acts/Detector/GeometryIdGenerator.hpp" +#include "Acts/Detector/IndexedRootVolumeFinderBuilder.hpp" +#include "Acts/Detector/detail/BlueprintHelper.hpp" +#include "Acts/Detector/interface/IInternalStructureBuilder.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Navigation/DetectorVolumeFinders.hpp" +#include "Acts/Navigation/SurfaceCandidatesUpdators.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Utilities/BinningData.hpp" + +#include + +template +class SurfaceBuilder : public Acts::Experimental::IInternalStructureBuilder { + public: + SurfaceBuilder(const Acts::Transform3& trf, Acts::ActsScalar p0, + Acts::ActsScalar p1) + : m_surface(Acts::Surface::makeShared(trf, p0, p1)) {} + /// Conrstruct and return the internal structure creation + /// + /// @param gctx the geometry context at the creation of the internal structure + /// + /// @return a consistent set of detector volume internals + Acts::Experimental::InternalStructure construct( + [[maybe_unused]] const Acts::GeometryContext& gctx) const final { + // Trivialities first: internal volumes + std::vector> + internalVolumes = {}; + Acts::Experimental::DetectorVolumeUpdator internalVolumeUpdator = + Acts::Experimental::tryNoVolumes(); + + // Retrieve the layer surfaces + Acts::Experimental::SurfaceCandidatesUpdator internalCandidatesUpdator = + Acts::Experimental::tryAllPortalsAndSurfaces(); + + // Return the internal structure + return Acts::Experimental::InternalStructure{ + {m_surface}, + internalVolumes, + std::move(internalCandidatesUpdator), + std::move(internalVolumeUpdator)}; + } + + private: + std::shared_ptr m_surface; +}; + +BOOST_AUTO_TEST_SUITE(Detector) + +BOOST_AUTO_TEST_CASE(CylindricalDetectorFromBlueprintTest) { + Acts::GeometryContext tContext; + + // This tests shows how to careate cylindrical detector from a detector + // blueprint. + // + // In general, the blueprint (lines below) is generated through reading in + // or by parsing the geometry model (DD4heo, TGeo, Geant4, etc.). For + // testing purpose, let us create the blueprint manually. + // + + // Blueprint starts here ---------------- + + // Detector dimensions + Acts::ActsScalar detectorIr = 0.; + Acts::ActsScalar detectorOr = 120.; + Acts::ActsScalar detectorHz = 400.; + + // Beam pipe + Acts::ActsScalar beamPipeOr = 20.; + + // Pixel system + Acts::ActsScalar pixelIr = 25; + Acts::ActsScalar pixelOr = 115; + Acts::ActsScalar pixelEcHz = 50; + Acts::ActsScalar pixelEcLayerHz = 10; + + // Create root node + std::vector detectorBinning = {Acts::binR}; + std::vector detectorBoundaries = {detectorIr, detectorOr, + detectorHz}; + + // The root node - detector + auto detectorBpr = std::make_unique( + "detector", Acts::Transform3::Identity(), Acts::VolumeBounds::eCylinder, + detectorBoundaries, detectorBinning); + + // The beam pipe + std::vector beamPipeBoundaries = {detectorIr, beamPipeOr, + detectorHz}; + + auto beamPipeStructure = + std::make_shared>( + Acts::Transform3::Identity(), 18, 0.99 * detectorHz); + auto beamPipe = std::make_unique( + "beam_pipe", Acts::Transform3::Identity(), Acts::VolumeBounds::eCylinder, + beamPipeBoundaries, beamPipeStructure); + detectorBpr->add(std::move(beamPipe)); + + // A pixel system + std::vector pixelBoundaries = {pixelIr, pixelOr, + detectorHz}; + std::vector pixelBinning = {Acts::binZ}; + auto pixel = std::make_unique( + "pixel", Acts::Transform3::Identity(), Acts::VolumeBounds::eCylinder, + pixelBoundaries, pixelBinning); + + // Nec: Small differences to check if the adjustments are made + std::vector pixelEcBoundaries = {pixelIr, pixelOr - 5., + pixelEcHz}; + std::vector pixelEcBinning = {Acts::binZ}; + + Acts::Transform3 pixelNecTransform = + Acts::Transform3::Identity() * + Acts::Translation3(0., 0., -detectorHz + pixelEcHz); + + auto pixelNec = std::make_unique( + "pixel_nec", pixelNecTransform, Acts::VolumeBounds::eCylinder, + pixelEcBoundaries, pixelEcBinning); + + // Add a single encap layer + std::vector pixelNecBoundaries = {pixelIr + 2, pixelOr - 7., + pixelEcLayerHz}; + + auto pixelNecLayerStructure = + std::make_shared>( + pixelNecTransform, pixelIr + 10., pixelOr - 10.); + + auto pixelNecLayer = std::make_unique( + "pixel_nec_layer", pixelNecTransform, Acts::VolumeBounds::eCylinder, + pixelNecBoundaries, pixelNecLayerStructure); + + pixelNec->add(std::move(pixelNecLayer)); + + // Barrel + std::vector pixelBarrelBoundaries = { + pixelIr + 1, pixelOr - 1., detectorHz - 2 * pixelEcHz}; + std::vector pixelBarrelBinning = {Acts::binR}; + + auto pixelBarrel = std::make_unique( + "pixel_barrel", Acts::Transform3::Identity(), + Acts::VolumeBounds::eCylinder, pixelBarrelBoundaries, pixelBarrelBinning); + + auto pixelBarrelL0Structure = + std::make_shared>( + Acts::Transform3::Identity(), 62.5, detectorHz - 2 * pixelEcHz - 10.); + std::vector pixelBarrelL0Boundaries = { + 60, 65., detectorHz - 2 * pixelEcHz}; + auto pixelBarrelL0 = std::make_unique( + "pixel_barrel_l0", Acts::Transform3::Identity(), + Acts::VolumeBounds::eCylinder, pixelBarrelL0Boundaries, + pixelBarrelL0Structure); + + auto pixelBarrelL1Structure = + std::make_shared>( + Acts::Transform3::Identity(), 102.5, + detectorHz - 2 * pixelEcHz - 10.); + + std::vector pixelBarrelL1Boundaries = { + 100, 105., detectorHz - 2 * pixelEcHz}; + auto pixelBarrelL1 = std::make_unique( + "pixel_barrel_l1", Acts::Transform3::Identity(), + Acts::VolumeBounds::eCylinder, pixelBarrelL1Boundaries, + pixelBarrelL1Structure); + pixelBarrel->add(std::move(pixelBarrelL0)); + pixelBarrel->add(std::move(pixelBarrelL1)); + + Acts::Transform3 pixelPecTransform = + Acts::Transform3::Identity() * + Acts::Translation3(0., 0., detectorHz - pixelEcHz); + + auto pixelPec = std::make_unique( + "pixel_pec", pixelPecTransform, Acts::VolumeBounds::eCylinder, + pixelEcBoundaries, pixelEcBinning); + + std::vector pixelPecBoundaries = {pixelIr + 2, pixelOr - 7., + 10.}; + + auto pixelPecLayerStructure = + std::make_shared>( + pixelPecTransform, pixelIr + 10., pixelOr - 10.); + + auto pixelPecLayer = std::make_unique( + "pixel_pec_layer", pixelPecTransform, Acts::VolumeBounds::eCylinder, + pixelPecBoundaries, pixelPecLayerStructure); + + pixelPec->add(std::move(pixelPecLayer)); + + // Adding pixel + pixel->add(std::move(pixelNec)); + pixel->add(std::move(pixelPec)); + pixel->add(std::move(pixelBarrel)); + + detectorBpr->add(std::move(pixel)); + + // An Indexed volume finder will be attached + std::vector rootVolumeBinning = {Acts::binZ, Acts::binR}; + detectorBpr->rootVolumeFinderBuilder = + std::make_shared( + rootVolumeBinning); + + // A geo ID generator + detectorBpr->geoIdGenerator = + std::make_shared( + Acts::Experimental::GeometryIdGenerator::Config{}, + Acts::getDefaultLogger("RecursiveIdGenerator", + Acts::Logging::VERBOSE)); + + // Complete and fill gaps + Acts::Experimental::detail::BlueprintHelper::fillGaps(*detectorBpr); + + std::fstream fs("cylindrical_detector_blueprint.dot", std::ios::out); + detectorBpr->dotStream(fs, "blueprint2cylinder"); + fs.close(); + + // ----------------------------- end of blueprint + + // Create a Cylindrical detector builder from this blueprint + auto detectorBuilder = + std::make_shared( + *detectorBpr, Acts::Logging::VERBOSE); + + // Detector builder + Acts::Experimental::DetectorBuilder::Config dCfg; + dCfg.auxiliary = + "*** Test : auto generated cylindrical detector builder ***"; + dCfg.name = "Cylindrical detector from blueprint"; + dCfg.builder = detectorBuilder; + dCfg.geoIdGenerator = detectorBpr->geoIdGenerator; + + auto detector = Acts::Experimental::DetectorBuilder(dCfg).construct(tContext); + + BOOST_CHECK(detector != nullptr); + + // There should be 14 volumes, and they should be built in order + // beam_pipe + // detector_gap_0 + // pixel_nec_gap_0 + // pixel_nec_layer + // pixel_nec_gap_1 + // pixel_barrel_gap_0 + // pixel_barrel_l0 + // pixel_barrel_gap_1 + // pixel_barrel_l1 + // pixel_barrel_gap_2 + // pixel_pec_gap_0 + // pixel_pec_layer + // pixel_pec_gap_1 + // detector_gap_1 + BOOST_CHECK(detector->volumes().size() == 14u); + BOOST_CHECK(detector->volumes()[0]->name() == "beam_pipe"); + BOOST_CHECK(detector->volumes()[1]->name() == "detector_gap_0"); + BOOST_CHECK(detector->volumes()[2]->name() == "pixel_nec_gap_0"); + BOOST_CHECK(detector->volumes()[3]->name() == "pixel_nec_layer"); + BOOST_CHECK(detector->volumes()[4]->name() == "pixel_nec_gap_1"); + BOOST_CHECK(detector->volumes()[5]->name() == "pixel_barrel_gap_0"); + BOOST_CHECK(detector->volumes()[6]->name() == "pixel_barrel_l0"); + BOOST_CHECK(detector->volumes()[7]->name() == "pixel_barrel_gap_1"); + BOOST_CHECK(detector->volumes()[8]->name() == "pixel_barrel_l1"); + BOOST_CHECK(detector->volumes()[9]->name() == "pixel_barrel_gap_2"); + BOOST_CHECK(detector->volumes()[10]->name() == "pixel_pec_gap_0"); + BOOST_CHECK(detector->volumes()[11]->name() == "pixel_pec_layer"); + BOOST_CHECK(detector->volumes()[12]->name() == "pixel_pec_gap_1"); + BOOST_CHECK(detector->volumes()[13]->name() == "detector_gap_1"); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Core/EventData/TrackTests.cpp b/Tests/UnitTests/Core/EventData/TrackTests.cpp index 85d05a73060..3a086b0a1a9 100644 --- a/Tests/UnitTests/Core/EventData/TrackTests.cpp +++ b/Tests/UnitTests/Core/EventData/TrackTests.cpp @@ -537,38 +537,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(DynamicColumns, factory_t, holder_types) { BOOST_CHECK_EQUAL((t.template component()), 5.6f); } -BOOST_AUTO_TEST_CASE(ReverseTrackStates) { - VectorTrackContainer vtc{}; - VectorMultiTrajectory mtj{}; - TrackContainer tc{vtc, mtj}; - - auto t = tc.getTrack(tc.addTrack()); - - for (size_t i = 0; i < 10; i++) { - t.appendTrackState(); - } - - std::vector exp; - exp.resize(t.nTrackStates()); - std::iota(exp.rbegin(), exp.rend(), 0); - std::vector act; - std::transform(t.trackStatesReversed().begin(), t.trackStatesReversed().end(), - std::back_inserter(act), - [](const auto& ts) { return ts.index(); }); - - BOOST_CHECK_EQUAL_COLLECTIONS(exp.begin(), exp.end(), act.begin(), act.end()); - - // reverse! - t.reverseTrackStates(); - - std::iota(exp.begin(), exp.end(), 0); - act.clear(); - std::transform(t.trackStatesReversed().begin(), t.trackStatesReversed().end(), - std::back_inserter(act), - [](const auto& ts) { return ts.index(); }); - BOOST_CHECK_EQUAL_COLLECTIONS(exp.begin(), exp.end(), act.begin(), act.end()); -} - BOOST_AUTO_TEST_CASE(EnsureDynamicColumns) { TrackContainer tc{VectorTrackContainer{}, VectorMultiTrajectory{}}; tc.addColumn("counter"); diff --git a/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp b/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp index 0a6ab9b17eb..71b57a708d8 100644 --- a/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp +++ b/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp @@ -11,7 +11,11 @@ #include "Acts/EventData/VectorTrackContainer.hpp" #include "Acts/Utilities/Zip.hpp" +#include + using namespace Acts; +using namespace Acts::HashedStringLiteral; +using MultiTrajectoryTraits::IndexType; BOOST_AUTO_TEST_SUITE(EventDataTrack) @@ -99,4 +103,83 @@ BOOST_AUTO_TEST_CASE(CopyTracksIncludingDynamicColumns) { t5.template component("odd")); } } + +BOOST_AUTO_TEST_CASE(ReverseTrackStates) { + VectorTrackContainer vtc{}; + VectorMultiTrajectory mtj{}; + TrackContainer tc{vtc, mtj}; + + auto t = tc.getTrack(tc.addTrack()); + + for (size_t i = 0; i < 4; i++) { + auto ts = t.appendTrackState(); + ts.jacobian() = Acts::BoundMatrix::Identity() * i; + } + + std::vector exp; + exp.resize(t.nTrackStates()); + std::iota(exp.rbegin(), exp.rend(), 0); + std::vector act; + std::transform(t.trackStatesReversed().begin(), t.trackStatesReversed().end(), + std::back_inserter(act), + [](const auto& ts) { return ts.index(); }); + + // jacobians count up + for (const auto [e, ts] : zip(exp, t.trackStatesReversed())) { + BOOST_CHECK_EQUAL(ts.jacobian(), Acts::BoundMatrix::Identity() * e); + } + + BOOST_CHECK_EQUAL_COLLECTIONS(exp.begin(), exp.end(), act.begin(), act.end()); + + // reverse! + t.reverseTrackStates(); + + std::iota(exp.begin(), exp.end(), 0); + act.clear(); + std::transform(t.trackStatesReversed().begin(), t.trackStatesReversed().end(), + std::back_inserter(act), + [](const auto& ts) { return ts.index(); }); + BOOST_CHECK_EQUAL_COLLECTIONS(exp.begin(), exp.end(), act.begin(), act.end()); + + // jacobians stay with their track states + for (const auto [e, ts] : zip(exp, t.trackStatesReversed())) { + BOOST_CHECK_EQUAL(ts.jacobian(), Acts::BoundMatrix::Identity() * e); + } + + // back to original! + t.reverseTrackStates(); + + // jacobians stay with their track states + for (const auto [e, ts] : zip(exp, t.trackStates())) { + BOOST_CHECK_EQUAL(ts.jacobian(), Acts::BoundMatrix::Identity() * e); + } + + // reverse with jacobians + t.reverseTrackStates(true); + + std::reverse(exp.begin(), exp.end()); + std::rotate(exp.rbegin(), std::next(exp.rbegin()), exp.rend()); + + for (const auto [e, ts] : zip(exp, t.trackStates())) { + Acts::BoundMatrix expJac; + if (e == 0) { + expJac = Acts::BoundMatrix::Zero(); + } else { + expJac = (Acts::BoundMatrix::Identity() * e).inverse(); + } + + BOOST_CHECK_EQUAL(ts.jacobian(), expJac); + } + + // now back to original order, revert jacobians again + t.reverseTrackStates(true); + + // reset exp to range(0, N) + std::iota(exp.begin(), exp.end(), 0); + + for (const auto [e, ts] : zip(exp, t.trackStates())) { + BOOST_CHECK_EQUAL(ts.jacobian(), Acts::BoundMatrix::Identity() * e); + } +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp index 0b52f63d602..fb54f97c26d 100644 --- a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp @@ -29,7 +29,6 @@ #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp" #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp" #include "Acts/Tests/CommonHelpers/TestSourceLink.hpp" -#include "Acts/TrackFitting/GainMatrixUpdater.hpp" #include "Acts/TrackFitting/GlobalChiSquareFitter.hpp" #include "Acts/Utilities/Logger.hpp" @@ -66,18 +65,6 @@ Acts::CurvilinearTrackParameters makeParameters( Acts::ParticleHypothesis::pion()); } -// Construct a straight-line propagator. -auto makeStraightPropagator(std::shared_ptr geo) { - Acts::Navigator::Config cfg{std::move(geo)}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Acts::Navigator navigator(cfg); - Acts::StraightLineStepper stepper; - return Acts::Propagator( - stepper, std::move(navigator)); -} - static std::vector prepareSourceLinks( const std::vector& sourceLinks) { std::vector result; @@ -196,8 +183,8 @@ BOOST_AUTO_TEST_CASE(NoFit) { detector.geometry = makeToyDetector(geoCtx, nSurfaces); auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(0.1_m, 0.1_m, 0.1_m, 42_ns, - 10_degree, 80_degree, 1_GeV, 1_e); + auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, + 80_degree, 1_GeV, 1_e); MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; MeasurementResolutionMap resolutions = { @@ -213,7 +200,7 @@ BOOST_AUTO_TEST_CASE(NoFit) { using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter Fitter(simPropagator, gx2fLogger->clone()); + Gx2Fitter fitter(simPropagator, gx2fLogger->clone()); const Surface* rSurface = ¶metersMeasurements.referenceSurface(); @@ -232,7 +219,7 @@ BOOST_AUTO_TEST_CASE(NoFit) { Acts::VectorMultiTrajectory{}}; // Fit the track - auto res = Fitter.fit(sourceLinks.begin(), sourceLinks.end(), + auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); @@ -262,9 +249,9 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { ACTS_DEBUG("Go to propagator"); auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(10_mm, 10_mm, 10_mm, 42_ns, - 10_degree, 80_degree, 1_GeV, 1_e); - // auto startParametersFit = parametersMeasurements; + auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, + 80_degree, 1_GeV, 1_e); + // Context objects Acts::GeometryContext geoCtx; Acts::MagneticFieldContext magCtx; @@ -272,8 +259,6 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { std::default_random_engine rng(42); MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; - // MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {100_um}}; - // MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {150_um}}; MeasurementResolutionMap resolutions = { {Acts::GeometryIdentifier().setVolume(0), resPixel}}; @@ -294,21 +279,14 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { const Surface* rSurface = ¶metersMeasurements.referenceSurface(); - Navigator::Config cfg{detector.geometry}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Navigator rNavigator(cfg); - // Configure propagation with deactivated B-field - auto bField = std::make_shared(Vector3(0., 0., 0.)); using RecoStepper = EigenStepper<>; - RecoStepper rStepper(bField); - using RecoPropagator = Propagator; - RecoPropagator rPropagator(rStepper, rNavigator); + const auto recoPropagator = + makeConstantFieldPropagator(detector.geometry, 0_T); + using RecoPropagator = decltype(recoPropagator); using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter Fitter(rPropagator, gx2fLogger->clone()); + Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); Experimental::Gx2FitterExtensions extensions; extensions.calibrator @@ -328,7 +306,7 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { Acts::VectorMultiTrajectory{}}; // Fit the track - auto res = Fitter.fit(sourceLinks.begin(), sourceLinks.end(), + auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), startParametersFit, gx2fOptions, tracks); BOOST_REQUIRE(res.ok()); @@ -340,8 +318,8 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { BOOST_CHECK_EQUAL(track.nHoles(), 0u); // We need quite coarse checks here, since on different builds // the created measurements differ in the randomness - BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -10., 6e0); - BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -10., 6e0); + BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0); + BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0); BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3); BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3); BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1); @@ -365,9 +343,9 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { ACTS_DEBUG("Go to propagator"); auto parametersMeasurements = makeParameters(); - auto startParametersFit = makeParameters(10_mm, 10_mm, 10_mm, 42_ns, - 10_degree, 80_degree, 1_GeV, 1_e); - // auto startParametersFit = parametersMeasurements; + auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, + 80_degree, 1_GeV, 1_e); + // Context objects Acts::GeometryContext geoCtx; Acts::MagneticFieldContext magCtx; @@ -404,21 +382,14 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { const Surface* rSurface = ¶metersMeasurements.referenceSurface(); - Navigator::Config cfg{detector.geometry}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Navigator rNavigator(cfg); - // Configure propagation with deactivated B-field - auto bField = std::make_shared(Vector3(0., 0., 0.)); using RecoStepper = EigenStepper<>; - RecoStepper rStepper(bField); - using RecoPropagator = Propagator; - RecoPropagator rPropagator(rStepper, rNavigator); + const auto recoPropagator = + makeConstantFieldPropagator(detector.geometry, 0_T); + using RecoPropagator = decltype(recoPropagator); using Gx2Fitter = Experimental::Gx2Fitter; - Gx2Fitter fitter(rPropagator, gx2fLogger->clone()); + Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); Experimental::Gx2FitterExtensions extensions; extensions.calibrator @@ -450,8 +421,8 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { BOOST_CHECK_EQUAL(track.nHoles(), 0u); // We need quite coarse checks here, since on different builds // the created measurements differ in the randomness - BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -10., 6e0); - BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -10., 6e0); + BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0); + BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0); BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3); BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3); BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1); diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp index 4eb4ef838f8..19ecf2fb23f 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp @@ -489,7 +489,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { VertexInfo vtxInfo1; vtxInfo1.linPoint.setZero(); vtxInfo1.linPoint.head<3>() = vtxPos1; - vtxInfo1.constraintVertex = vtx1Constr; + vtxInfo1.constraint = vtx1Constr; vtxInfo1.oldPosition = vtxInfo1.linPoint; vtxInfo1.seedPosition = vtxInfo1.linPoint; @@ -516,7 +516,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { VertexInfo vtxInfo2; vtxInfo2.linPoint.setZero(); vtxInfo2.linPoint.head<3>() = vtxPos2; - vtxInfo2.constraintVertex = vtx2Constr; + vtxInfo2.constraint = vtx2Constr; vtxInfo2.oldPosition = vtxInfo2.linPoint; vtxInfo2.seedPosition = vtxInfo2.linPoint; diff --git a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp index 3606f8dba20..1c2b8586be0 100644 --- a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp +++ b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp @@ -125,6 +125,7 @@ struct MockStepper { struct MockNavigatorState { bool targetReached = false; + Acts::Surface *startSurface = nullptr; Acts::Surface *currentSurface = nullptr; }; @@ -133,6 +134,10 @@ struct MockNavigator { return state.targetReached; } + const Acts::Surface *startSurface(const MockNavigatorState &state) const { + return state.startSurface; + } + const Acts::Surface *currentSurface(const MockNavigatorState &state) const { return state.currentSurface; }