diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index e575a4f1bdc..16da26e3d37 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -20,6 +20,7 @@ #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/EventData/VectorTrackContainer.hpp" #include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "Acts/Material/MaterialSlab.hpp" #include "Acts/Propagator/AbortList.hpp" @@ -221,6 +222,10 @@ struct Gx2FitterResult { // Count how many surfaces have been hit std::size_t surfaceCount = 0; + + // Monitor which volume we start in. We do not allow to switch the start of a + // following iteration in a different volume. + const TrackingVolume* startVolume = nullptr; }; /// addToGx2fSums Function @@ -397,6 +402,10 @@ class Gx2Fitter { /// Calibration context for the fit const CalibrationContext* calibrationContext{nullptr}; + /// Monitor which volume we start in. We do not allow to switch the start of + /// a following iteration in a different volume. + const TrackingVolume* startVolume = nullptr; + /// @brief Gx2f actor operation /// /// @tparam propagator_state_t is the type of Propagator state @@ -424,7 +433,7 @@ class Gx2Fitter { } // End the propagation and return to the fitter - if (result.finished) { + if (result.finished || !result.result.ok()) { // Remove the missing surfaces that occur after the last measurement if (result.measurementStates > 0) { result.missedActiveSurfaces.resize(result.measurementHoles); @@ -433,6 +442,19 @@ class Gx2Fitter { return; } + if (startVolume != nullptr && + startVolume != state.navigation.startVolume) { + ACTS_INFO("The update pushed us to a new volume from '" + << startVolume->volumeName() << "' to '" + << state.navigation.startVolume->volumeName() + << "'. Starting to abort."); + result.result = Result( + Experimental::GlobalChiSquareFitterError::UpdatePushedToNewVolume); + + return; + } + result.startVolume = state.navigation.startVolume; + // Update: // - Waiting for a current surface auto surface = navigator.currentSurface(state.navigation); @@ -716,6 +738,10 @@ class Gx2Fitter { // want to fit e.g. q/p and adjusts itself later. std::size_t ndfSystem = std::numeric_limits::max(); + // Monitor which volume we start in. We do not allow to switch the start of + // a following iteration in a different volume. + const TrackingVolume* startVolume = nullptr; + ACTS_VERBOSE("params:\n" << params); /// Actual Fitting ///////////////////////////////////////////////////////// @@ -750,6 +776,7 @@ class Gx2Fitter { gx2fActor.extensions = gx2fOptions.extensions; gx2fActor.calibrationContext = &gx2fOptions.calibrationContext.get(); gx2fActor.actorLogger = m_actorLogger.get(); + gx2fActor.startVolume = startVolume; auto propagatorState = m_propagator.makeState(params, propagatorOptions); @@ -762,6 +789,7 @@ class Gx2Fitter { auto propagationResult = m_propagator.template propagate(propagatorState); + // Run the fitter auto result = m_propagator.template makeResult(std::move(propagatorState), propagationResult, propagatorOptions, false); @@ -776,6 +804,13 @@ class Gx2Fitter { auto& propRes = *result; GX2FResult gx2fResult = std::move(propRes.template get()); + if (!gx2fResult.result.ok()) { + ACTS_INFO("GlobalChiSquareFitter failed in actor: " + << gx2fResult.result.error() << ", " + << gx2fResult.result.error().message()); + return gx2fResult.result.error(); + } + auto track = trackContainerTemp.makeTrack(); tipIndex = gx2fResult.lastMeasurementIndex; @@ -903,6 +938,7 @@ class Gx2Fitter { } oldChi2sum = chi2sum; + startVolume = gx2fResult.startVolume; } ACTS_DEBUG("Finished to iterate"); ACTS_VERBOSE("final params:\n" << params); diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitterError.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitterError.hpp index f84505769a6..ac47142b07d 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitterError.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitterError.hpp @@ -18,6 +18,7 @@ enum class GlobalChiSquareFitterError { AIsNotInvertible = 1, DidNotConverge = 2, NotEnoughMeasurements = 3, + UpdatePushedToNewVolume = 4, }; std::error_code make_error_code( diff --git a/Core/src/TrackFitting/GlobalChiSquareFitterError.cpp b/Core/src/TrackFitting/GlobalChiSquareFitterError.cpp index dff5795c932..f1fd9accd3a 100644 --- a/Core/src/TrackFitting/GlobalChiSquareFitterError.cpp +++ b/Core/src/TrackFitting/GlobalChiSquareFitterError.cpp @@ -30,6 +30,8 @@ class GlobalChiSquareFitterErrorCategory : public std::error_category { return "Gx2f: Did not converge in 'nUpdateMax' updates."; case GlobalChiSquareFitterError::NotEnoughMeasurements: return "Gx2f: Not enough measurements."; + case GlobalChiSquareFitterError::UpdatePushedToNewVolume: + return "Gx2f: Update pushed the parameters to a new volume."; default: return "unknown"; } diff --git a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp index caa15e9c6d8..44ebc963829 100644 --- a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp @@ -31,6 +31,9 @@ #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp" #include "Acts/TrackFitting/GlobalChiSquareFitter.hpp" #include "Acts/Utilities/Logger.hpp" +#include "Acts/Visualization/EventDataView3D.hpp" +#include "Acts/Visualization/GeometryView3D.hpp" +#include "Acts/Visualization/ObjVisualization3D.hpp" #include @@ -44,6 +47,37 @@ const auto gx2fLogger = Acts::getDefaultLogger("Gx2f", logLevel); namespace Acts::Test { +/// @brief Helper function to visualise measurements in a 3D environment. +/// +/// This function iterates through the provided measurements and visualises each +/// one using the specified 3D visualisation helper. The visualisation takes +/// into account the surface transformations and localisation errors. +/// +/// @param helper The 3D visualisation helper used to draw the measurements. +/// @param measurements A collection of measurements to be visualised, containing source links with parameters and covariance information. +/// @param geometry A shared pointer to the constant tracking geometry used to find surfaces associated with measurements. +/// @param geoCtx The geometry context used for transformations and accessing geometry-related information. +/// @param locErrorScale Scaling factor for localisation errors. Default value is 1.0. +/// @param viewConfig Configuration settings for the visualisation view. Default value is s_viewMeasurement. +static void drawMeasurements( + IVisualization3D& helper, const Measurements& measurements, + const std::shared_ptr& geometry, + const Acts::GeometryContext& geoCtx, double locErrorScale = 1., + const ViewConfig& viewConfig = s_viewMeasurement) { + std::cout << "\n*** Draw measurements ***\n" << std::endl; + + for (auto& singleMeasurement : measurements.sourceLinks) { + auto cov = singleMeasurement.covariance; + auto lposition = singleMeasurement.parameters; + + auto surf = geometry->findSurface(singleMeasurement.m_geometryId); + auto transf = surf->transform(geoCtx); + + EventDataView3D::drawMeasurement(helper, lposition, cov, transf, + locErrorScale, viewConfig); + } +} + //// Construct initial track parameters. Acts::CurvilinearTrackParameters makeParameters( const ActsScalar x = 0.0_m, const ActsScalar y = 0.0_m, @@ -78,6 +112,7 @@ static std::vector prepareSourceLinks( /// /// @param geoCtx /// @param nSurfaces Number of surfaces +/// @param surfaceIndexWithMaterial A set of index of the material surfaces std::shared_ptr makeToyDetector( const Acts::GeometryContext& geoCtx, const std::size_t nSurfaces = 5, const std::set& surfaceIndexWithMaterial = {}) { @@ -88,7 +123,7 @@ std::shared_ptr makeToyDetector( // Define the dimensions of the square surfaces const double halfSizeSurface = 1_m; - // Rotation of the surfaces + // Rotation of the surfaces around the y-axis const double rotationAngle = M_PI * 0.5; const Vector3 xPos(cos(rotationAngle), 0., sin(rotationAngle)); const Vector3 yPos(0., 1., 0.); @@ -150,7 +185,7 @@ std::shared_ptr makeToyDetector( 2 * halfSizeSurface}; volumeConfig.position = {volumeConfig.length.x() / 2, 0., 0.}; volumeConfig.layerCfg = layerConfig; - volumeConfig.name = "Test volume"; + volumeConfig.name = "TestVolume"; // Outer volume - Build TrackingGeometry configuration CuboidVolumeBuilder::Config config; @@ -175,6 +210,117 @@ std::shared_ptr makeToyDetector( return detector; } +/// @brief Create a simple telescope detector in the Y direction. +/// +/// We cannot reuse the previous detector, since the cuboid volume builder only +/// allows merging of YZ-faces. +/// +/// @param geoCtx +/// @param nSurfaces Number of surfaces +std::shared_ptr makeToyDetectorYdirection( + const Acts::GeometryContext& geoCtx, const std::size_t nSurfaces = 5) { + if (nSurfaces < 1) { + throw std::invalid_argument("At least 1 surfaces needs to be created."); + } + + // Define the dimensions of the square surfaces + const double halfSizeSurface = 1_m; + + // Rotation of the surfaces around the x-axis + const double rotationAngle = M_PI * 0.5; + const Vector3 xPos(1., 0., 0.); + const Vector3 yPos(0., cos(rotationAngle), sin(rotationAngle)); + const Vector3 zPos(0., -sin(rotationAngle), cos(rotationAngle)); + + // Construct builder + CuboidVolumeBuilder cvb; + + // Create configurations for surfaces + std::vector surfaceConfig; + for (std::size_t surfPos = 1; surfPos <= nSurfaces; surfPos++) { + // Position of the surfaces + CuboidVolumeBuilder::SurfaceConfig cfg; + cfg.position = {0., surfPos * UnitConstants::m, 0.}; + + // Rotation of the surfaces + cfg.rotation.col(0) = xPos; + cfg.rotation.col(1) = yPos; + cfg.rotation.col(2) = zPos; + + // Boundaries of the surfaces (shape) + cfg.rBounds = std::make_shared( + RectangleBounds(halfSizeSurface, halfSizeSurface)); + + // Thickness of the detector element + cfg.thickness = 1_um; + + cfg.detElementConstructor = + [](const Transform3& trans, + const std::shared_ptr& bounds, + double thickness) { + return new DetectorElementStub(trans, bounds, thickness); + }; + surfaceConfig.push_back(cfg); + } + + // Build layer configurations + std::vector layerConfig; + for (auto& sCfg : surfaceConfig) { + CuboidVolumeBuilder::LayerConfig cfg; + cfg.surfaceCfg = {sCfg}; + cfg.active = true; + cfg.envelopeX = {-0.1_mm, 0.1_mm}; + cfg.envelopeY = {-0.1_mm, 0.1_mm}; + cfg.envelopeZ = {-0.1_mm, 0.1_mm}; + cfg.binningDimension = Acts::BinningValue::binY; + layerConfig.push_back(cfg); + } + + // Inner Volume - Build volume configuration + CuboidVolumeBuilder::VolumeConfig volumeConfig; + volumeConfig.length = {2 * halfSizeSurface, (nSurfaces + 1) * 1_m, + 2 * halfSizeSurface}; + volumeConfig.position = {0., volumeConfig.length.y() / 2, 0.}; + volumeConfig.layerCfg = layerConfig; + volumeConfig.name = "TestVolume"; + volumeConfig.binningDimension = Acts::BinningValue::binY; + + // This basically adds an empty volume in y-direction + // Second inner Volume - Build volume configuration + CuboidVolumeBuilder::VolumeConfig volumeConfig2; + // volumeConfig2.length = volumeConfig.length; + volumeConfig2.length = {2 * halfSizeSurface, (nSurfaces + 1) * 1_m, + 2 * halfSizeSurface}; + ; + volumeConfig2.position = {volumeConfig2.length.x(), + volumeConfig2.length.y() / 2, 0.}; + volumeConfig2.name = "AdditionalVolume"; + volumeConfig2.binningDimension = Acts::BinningValue::binY; + + // Outer volume - Build TrackingGeometry configuration and fill + CuboidVolumeBuilder::Config config; + config.length = {4 * halfSizeSurface, (nSurfaces + 1) * 1_m, + 2 * halfSizeSurface}; + config.position = {volumeConfig.length.x() / 2, volumeConfig.length.y() / 2, + 0.}; + config.volumeCfg = {volumeConfig, volumeConfig2}; + + cvb.setConfig(config); + + TrackingGeometryBuilder::Config tgbCfg; + + tgbCfg.trackingVolumeBuilders.push_back( + [=](const auto& context, const auto& inner, const auto&) { + return cvb.trackingVolume(context, inner, nullptr); + }); + + TrackingGeometryBuilder tgb(tgbCfg); + + std::unique_ptr detector = + tgb.trackingGeometry(geoCtx); + return detector; +} + struct Detector { // geometry std::shared_ptr geometry; @@ -377,6 +523,120 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { ACTS_INFO("*** Test: Fit5Iterations -- Finish"); } +BOOST_AUTO_TEST_CASE(UpdatePushedToNewVolume) { + ACTS_INFO("*** Test: UpdatePushedToNewVolume -- Start"); + + std::default_random_engine rng(42); + + ACTS_DEBUG("Create the detector"); + const std::size_t nSurfaces = 5; + Detector detector; + detector.geometry = makeToyDetectorYdirection(geoCtx, nSurfaces); + + ACTS_DEBUG("Set the start parameters for measurement creation and fit"); + const auto parametersMeasurements = + makeParameters(0_mm, 0_mm, 0_mm, 42_ns, 90_degree, 90_degree, 1_GeV, 1_e); + const auto startParametersFit = makeParameters( + 1500_mm, 0_mm, 0_mm, 42_ns, -45_degree, -90_degree, 1_GeV, 1_e); + + ACTS_DEBUG("Create the measurements"); + using SimPropagator = + Acts::Propagator; + const SimPropagator simPropagator = makeStraightPropagator(detector.geometry); + const auto measurements = + createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements, + resMapAllPixel, rng); + const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size()); + + BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces); + + ACTS_DEBUG("Set up the fitter"); + const Surface* rSurface = ¶metersMeasurements.referenceSurface(); + + using RecoStepper = EigenStepper<>; + const auto recoPropagator = + makeConstantFieldPropagator(detector.geometry, 0_T); + + using RecoPropagator = decltype(recoPropagator); + using Gx2Fitter = + Experimental::Gx2Fitter; + const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); + + Experimental::Gx2FitterExtensions extensions; + extensions.calibrator + .connect<&testSourceLinkCalibrator>(); + TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry}; + extensions.surfaceAccessor + .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); + + const Experimental::Gx2FitterOptions gx2fOptions( + geoCtx, magCtx, calCtx, extensions, + PropagatorPlainOptions(geoCtx, magCtx), rSurface, false, false, + FreeToBoundCorrection(false), 5, 0); + + Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, + Acts::VectorMultiTrajectory{}}; + + ACTS_DEBUG("Fit the track"); + ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements); + ACTS_VERBOSE("startParameter fit:\n" << startParametersFit); + const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), + startParametersFit, gx2fOptions, tracks); + + // Helper to visualise the detector + { + std::cout << "\n*** Create .obj of Detector ***\n" << std::endl; + // Only need for obj + ObjVisualization3D obj; + + bool triangulate = true; + ViewConfig viewSensitive = ViewConfig({0, 180, 240}); + viewSensitive.triangulate = triangulate; + ViewConfig viewPassive = ViewConfig({240, 280, 0}); + viewPassive.triangulate = triangulate; + ViewConfig viewVolume = ViewConfig({220, 220, 0}); + viewVolume.triangulate = triangulate; + ViewConfig viewContainer = ViewConfig({220, 220, 0}); + viewContainer.triangulate = triangulate; + ViewConfig viewGrid = ViewConfig({220, 0, 0}); + viewGrid.nSegments = 8; + viewGrid.offset = 3.; + viewGrid.triangulate = triangulate; + + std::string tag = "gx2f_toydet"; + + const Acts::TrackingVolume& tgVolume = + *(detector.geometry->highestTrackingVolume()); + + GeometryView3D::drawTrackingVolume(obj, tgVolume, geoCtx, viewContainer, + viewVolume, viewPassive, viewSensitive, + viewGrid, true, tag); + } + // Helper to visualise the measurements + { + std::cout << "\n*** Create .obj of measurements ***\n" << std::endl; + ObjVisualization3D obj; + + double localErrorScale = 10000000.; + ViewConfig mcolor({255, 145, 48}); + mcolor.offset = 2; + // mcolor.visible = true; + + drawMeasurements(obj, measurements, detector.geometry, geoCtx, + localErrorScale, mcolor); + + obj.write("meas"); + } + + BOOST_REQUIRE(!res.ok()); + BOOST_CHECK_EQUAL( + res.error(), + Acts::Experimental::GlobalChiSquareFitterError::UpdatePushedToNewVolume); + + ACTS_INFO("*** Test: UpdatePushedToNewVolume -- Finish"); +} + BOOST_AUTO_TEST_CASE(MixedDetector) { ACTS_INFO("*** Test: MixedDetector -- Start"); @@ -853,7 +1113,7 @@ BOOST_AUTO_TEST_CASE(FindHoles) { sourceLinks.pop_back(); ACTS_VERBOSE("sourceLinks.size() [after pop] = " << sourceLinks.size()); - // We remove the second to last measurement in the list. This effectivly + // We remove the second to last measurement in the list. This effectively // creates a hole on that surface. const std::size_t indexHole = sourceLinks.size() - 2; ACTS_VERBOSE("Remove measurement " << indexHole);