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