Merge branch 'main' into ex-rename-vtx-perf-writer

CI/physmon/phys_perf_mon.sh

@@ -225,11 +225,19 @@ function trackfinding() {
     run_histcmp \
         $outdir/data/$path/performance_ckf.root \
         $refdir/$path/performance_ckf.root \
-        "CKF | ${name}" \
+        "CKF finding performance | ${name}" \
         $path/performance_ckf.html \
         $path/performance_ckf_plots \
         --config $default_config
 
+    run_histcmp \
+        $outdir/data/$path/performance_fitting_ckf.root \
+        $refdir/$path/performance_fitting_ckf.root \
+        "CKF fitting performance | ${name}" \
+        $path/performance_fitting_ckf.html \
+        $path/performance_fitting_ckf_plots \
+        --config $default_config
+
     run Examples/Scripts/generic_plotter.py \
         $outdir/data/$path/tracksummary_ckf.root \
         tracksummary \
@@ -252,7 +260,7 @@ function trackfinding() {
         run_histcmp \
             $outdir/data/$path/performance_ckf_ambi.root \
             $refdir/$path/performance_ckf_ambi.root \
-            "Ambisolver | ${name}" \
+            "Ambisolver finding performance | ${name}" \
             $path/performance_ckf_ambi.html \
             $path/performance_ckf_ambi
     fi

CI/physmon/workflows/physmon_trackfinding_1muon.py

@@ -112,7 +112,7 @@ def run_ckf_tracking(label, seeding):
                 maxSeedsPerSpM=1,
                 sigmaScattering=5,
                 radLengthPerSeed=0.1,
-                minPt=500 * u.MeV,
+                minPt=0.5 * u.GeV,
                 impactMax=3 * u.mm,
             ),
             SeedFinderOptionsArg(bFieldInZ=2 * u.T),
@@ -163,8 +163,9 @@ def run_ckf_tracking(label, seeding):
             if seeding != SeedingAlgorithm.TruthSmeared
             else []
         ) + [
-            "performance_ckf.root",
             "tracksummary_ckf.root",
+            "performance_fitting_ckf.root",
+            "performance_ckf.root",
         ]:
             perf_file = tp / file
             assert perf_file.exists(), f"Performance file not found {perf_file}"

CI/physmon/workflows/physmon_trackfinding_4muon_50vertices.py

@@ -96,7 +96,7 @@
             maxSeedsPerSpM=1,
             sigmaScattering=5,
             radLengthPerSeed=0.1,
-            minPt=500 * u.MeV,
+            minPt=0.5 * u.GeV,
             impactMax=3 * u.mm,
         ),
         SeedFinderOptionsArg(bFieldInZ=2 * u.T, beamPos=(0.0, 0.0)),
@@ -195,6 +195,10 @@
 
     s.run()
 
+    shutil.move(
+        tp / "performance_fitting_ambi.root",
+        tp / "performance_fitting_ckf_ambi.root",
+    )
     shutil.move(
         tp / "performance_ambi.root",
         tp / "performance_ckf_ambi.root",
@@ -208,7 +212,9 @@
     for file in [
         "performance_seeding.root",
         "tracksummary_ckf.root",
+        "performance_fitting_ckf.root",
         "performance_ckf.root",
+        "performance_fitting_ckf_ambi.root",
         "performance_ckf_ambi.root",
         "performance_vertexing_ivf_notime.root",
         "performance_vertexing_amvf_gauss_notime.root",

CI/physmon/workflows/physmon_trackfinding_ttbar_pu200.py

@@ -94,7 +94,7 @@
             maxSeedsPerSpM=1,
             sigmaScattering=5,
             radLengthPerSeed=0.1,
-            minPt=500 * u.MeV,
+            minPt=0.5 * u.GeV,
             impactMax=3 * u.mm,
         ),
         SeedFinderOptionsArg(bFieldInZ=2 * u.T, beamPos=(0.0, 0.0)),
@@ -118,7 +118,7 @@
         setup.trackingGeometry,
         setup.field,
         TrackSelectorConfig(
-            pt=(500 * u.MeV, None),
+            pt=(0.5 * u.GeV, None),
             loc0=(-4.0 * u.mm, 4.0 * u.mm),
             nMeasurementsMin=6,
             maxHoles=2,
@@ -179,6 +179,10 @@
 
     s.run()
 
+    shutil.move(
+        tp / "performance_fitting_ambi.root",
+        tp / "performance_fitting_ckf_ambi.root",
+    )
     shutil.move(
         tp / "performance_ambi.root",
         tp / "performance_ckf_ambi.root",
@@ -192,7 +196,9 @@
     for file in [
         "performance_seeding.root",
         "tracksummary_ckf.root",
+        "performance_fitting_ckf.root",
         "performance_ckf.root",
+        "performance_fitting_ckf_ambi.root",
         "performance_ckf_ambi.root",
         "performance_vertexing_amvf_gauss_notime.root",
         "performance_vertexing_amvf_grid_time.root",

Core/include/Acts/EventData/TrackParameterHelpers.hpp

@@ -0,0 +1,44 @@
+// 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/Definitions/Algebra.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Utilities/detail/periodic.hpp"
+
+namespace Acts {
+
+/// Normalize the bound parameter angles
+///
+/// @param boundParams The bound parameters to normalize
+///
+/// @return The normalized bound parameters
+inline BoundVector normalizeBoundParameters(const BoundVector& boundParams) {
+  BoundVector result = boundParams;
+  std::tie(result[eBoundPhi], result[eBoundTheta]) =
+      detail::normalizePhiTheta(result[eBoundPhi], result[eBoundTheta]);
+  return result;
+}
+
+/// Subtract bound parameters and take care of angle periodicity for phi and
+/// theta.
+///
+/// @param lhs The left hand side bound parameters
+/// @param rhs The right hand side bound parameters
+///
+/// @return The difference of the bound parameters
+inline BoundVector subtractBoundParameters(const BoundVector& lhs,
+                                           const BoundVector& rhs) {
+  BoundVector result = lhs - rhs;
+  result[eBoundPhi] =
+      detail::difference_periodic(lhs[eBoundPhi], rhs[eBoundPhi], 2 * M_PI);
+  return result;
+}
+
+}  // namespace Acts

Core/include/Acts/TrackFitting/detail/GainMatrixUpdaterImpl.hpp

@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include "Acts/EventData/TrackParameterHelpers.hpp"
 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
@@ -60,6 +61,8 @@ std::tuple<double, std::error_code> GainMatrixUpdater::visitMeasurementImpl(
 
   trackState.filtered =
       trackState.predicted + K * (calibrated - H * trackState.predicted);
+  // Normalize phi and theta
+  trackState.filtered = normalizeBoundParameters(trackState.filtered);
   trackState.filteredCovariance =
       (BoundSquareMatrix::Identity() - K * H) * trackState.predictedCovariance;
   ACTS_VERBOSE("Filtered parameters: " << trackState.filtered.transpose());

Core/src/TrackFitting/GainMatrixSmoother.cpp

@@ -9,6 +9,7 @@
 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
 
 #include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/EventData/TrackParameterHelpers.hpp"
 #include "Acts/EventData/detail/CovarianceHelper.hpp"
 #include "Acts/TrackFitting/KalmanFitterError.hpp"
 
@@ -57,7 +58,10 @@ Result<void> GainMatrixSmoother::calculate(
       "Prev. predicted parameters: " << predicted(prev_ts).transpose());
 
   // Calculate the smoothed parameters
-  smoothed(ts) = filtered(ts) + G * (smoothed(prev_ts) - predicted(prev_ts));
+  smoothed(ts) = filtered(ts) + G * subtractBoundParameters(smoothed(prev_ts),
+                                                            predicted(prev_ts));
+  // Normalize phi and theta
+  smoothed(ts) = normalizeBoundParameters(smoothed(ts));
 
   ACTS_VERBOSE("Smoothed parameters are: " << smoothed(ts).transpose());
   ACTS_VERBOSE("Calculate smoothed covariance:");

Core/src/TrackFitting/MbfSmoother.cpp

@@ -8,6 +8,8 @@
 
 #include "Acts/TrackFitting/MbfSmoother.hpp"
 
+#include "Acts/EventData/TrackParameterHelpers.hpp"
+
 namespace Acts {
 
 void MbfSmoother::calculateSmoothed(InternalTrackState& ts,
@@ -17,6 +19,8 @@ void MbfSmoother::calculateSmoothed(InternalTrackState& ts,
                                                       bigLambdaHat *
                                                       ts.filteredCovariance;
   ts.smoothed = ts.filtered - ts.filteredCovariance * smallLambdaHat;
+  // Normalize phi and theta
+  ts.smoothed = normalizeBoundParameters(ts.smoothed);
 }
 
 void MbfSmoother::visitNonMeasurement(const InternalTrackState& ts,

Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp

@@ -179,6 +179,7 @@ class TrackFindingAlgorithm final : public IAlgorithm {
   mutable std::atomic<std::size_t> m_nFoundTracks{0};
   mutable std::atomic<std::size_t> m_nSelectedTracks{0};
   mutable std::atomic<std::size_t> m_nStoppedBranches{0};
+  mutable std::atomic<std::size_t> m_nSkippedSecondPass{0};
 
   mutable tbb::combinable<Acts::VectorMultiTrajectory::Statistics>
       m_memoryStatistics{[]() {

Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp

@@ -532,6 +532,15 @@ ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
           Acts::BoundTrackParameters secondInitialParameters =
               trackCandidate.createParametersFromState(firstMeasurement);
 
+          if (!secondInitialParameters.referenceSurface().insideBounds(
+                  secondInitialParameters.localPosition())) {
+            m_nSkippedSecondPass++;
+            ACTS_DEBUG(
+                "Smoothing of first pass fit produced out-of-bounds parameters "
+                "relative to the surface. Skipping second pass.");
+            continue;
+          }
+
           auto secondRootBranch = tracksTemp.makeTrack();
           secondRootBranch.copyFrom(trackCandidate, false);
           auto secondResult =
@@ -681,6 +690,7 @@ ProcessCode TrackFindingAlgorithm::finalize() {
   ACTS_INFO("- found tracks: " << m_nFoundTracks);
   ACTS_INFO("- selected tracks: " << m_nSelectedTracks);
   ACTS_INFO("- stopped branches: " << m_nStoppedBranches);
+  ACTS_INFO("- skipped second pass: " << m_nSkippedSecondPass);
 
   auto memoryStatistics =
       m_memoryStatistics.combine([](const auto& a, const auto& b) {

Examples/Framework/src/Validation/ResPlotTool.cpp

@@ -158,21 +158,27 @@ void ActsExamples::ResPlotTool::fill(
   auto trackParameter = fittedParamters.parameters();
 
   // get the perigee surface
-  auto pSurface = &fittedParamters.referenceSurface();
+  const auto& pSurface = fittedParamters.referenceSurface();
 
   // get the truth position and momentum
   ParametersVector truthParameter = ParametersVector::Zero();
 
   // get the truth perigee parameter
-  auto lpResult = pSurface->globalToLocal(gctx, truthParticle.position(),
-                                          truthParticle.direction());
-  if (lpResult.ok()) {
+  auto intersection =
+      pSurface
+          .intersect(gctx, truthParticle.position(), truthParticle.direction())
+          .closest();
+  if (intersection.isValid()) {
+    auto lpResult = pSurface.globalToLocal(gctx, intersection.position(),
+                                           truthParticle.direction());
+    assert(lpResult.ok());
+
     truthParameter[Acts::BoundIndices::eBoundLoc0] =
         lpResult.value()[Acts::BoundIndices::eBoundLoc0];
     truthParameter[Acts::BoundIndices::eBoundLoc1] =
         lpResult.value()[Acts::BoundIndices::eBoundLoc1];
   } else {
-    ACTS_ERROR("Global to local transformation did not succeed.");
+    ACTS_ERROR("Cannot get the truth perigee parameter");
   }
   truthParameter[Acts::BoundIndices::eBoundPhi] =
       phi(truthParticle.direction());

Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp

@@ -173,6 +173,7 @@ class RootAthenaDumpReader : public IReader {
   std::vector<std::tuple<std::uint32_t, std::size_t, std::size_t>> m_eventMap;
   std::shared_ptr<TChain> m_inputchain;
   long unsigned int m_events;
+  bool m_haveStripFeatures = true;
 
   static constexpr unsigned int maxCL = 1500000;
   static constexpr unsigned int maxSP = 1500000;
@@ -189,7 +190,7 @@ class RootAthenaDumpReader : public IReader {
   int CLindex[maxCL] = {};  //[nCL]
 
   // Clusters
-  std::vector<std::string> *CLhardware;
+  std::vector<std::string> *CLhardware{};
   Double_t CLx[maxCL] = {};           //[nCL]
   Double_t CLy[maxCL] = {};           //[nCL]
   Double_t CLz[maxCL] = {};           //[nCL]
@@ -199,13 +200,13 @@ class RootAthenaDumpReader : public IReader {
   Int_t CLphi_module[maxCL] = {};     //[nCL]
   Int_t CLside[maxCL] = {};           //[nCL]
   ULong64_t CLmoduleID[maxCL] = {};   //[nCL]
-  std::vector<std::vector<int>> *CLparticleLink_eventIndex;
-  std::vector<std::vector<int>> *CLparticleLink_barcode;
-  std::vector<std::vector<bool>> *CLbarcodesLinked;
-  std::vector<std::vector<float>> *CLparticle_charge;
-  std::vector<std::vector<int>> *CLphis;
-  std::vector<std::vector<int>> *CLetas;
-  std::vector<std::vector<int>> *CLtots;
+  std::vector<std::vector<int>> *CLparticleLink_eventIndex{};
+  std::vector<std::vector<int>> *CLparticleLink_barcode{};
+  std::vector<std::vector<bool>> *CLbarcodesLinked{};
+  std::vector<std::vector<float>> *CLparticle_charge{};
+  std::vector<std::vector<int>> *CLphis{};
+  std::vector<std::vector<int>> *CLetas{};
+  std::vector<std::vector<int>> *CLtots{};
   Double_t CLloc_direction1[maxCL] = {};      //[nCL]
   Double_t CLloc_direction2[maxCL] = {};      //[nCL]
   Double_t CLloc_direction3[maxCL] = {};      //[nCL]
@@ -223,7 +224,7 @@ class RootAthenaDumpReader : public IReader {
   Float_t CLnorm_x[maxCL] = {};               //[nCL]
   Float_t CLnorm_y[maxCL] = {};               //[nCL]
   Float_t CLnorm_z[maxCL] = {};               //[nCL]
-  std::vector<std::vector<double>> *CLlocal_cov;
+  std::vector<std::vector<double>> *CLlocal_cov{};
 
   // Particles
   Int_t nPartEVT = 0;
@@ -246,8 +247,8 @@ class RootAthenaDumpReader : public IReader {
   Int_t Part_vProdNout[maxPart] = {};     //[nPartEVT]
   Int_t Part_vProdStatus[maxPart] = {};   //[nPartEVT]
   Int_t Part_vProdBarcode[maxPart] = {};  //[nPartEVT]
-  std::vector<std::vector<int>> *Part_vParentID;
-  std::vector<std::vector<int>> *Part_vParentBarcode;
+  std::vector<std::vector<int>> *Part_vParentID{};
+  std::vector<std::vector<int>> *Part_vParentBarcode{};
 
   // Spacepoints
   Int_t nSP = 0;
@@ -263,27 +264,27 @@ class RootAthenaDumpReader : public IReader {
   double SPcovz[maxSP] = {};        //[nSP]
   float SPhl_topstrip[maxSP] = {};  //[nSP]
   float SPhl_botstrip[maxSP] = {};  //[nSP]
-  std::vector<std::vector<float>> *SPtopStripDirection;
-  std::vector<std::vector<float>> *SPbottomStripDirection;
-  std::vector<std::vector<float>> *SPstripCenterDistance;
-  std::vector<std::vector<float>> *SPtopStripCenterPosition;
+  std::vector<std::vector<float>> *SPtopStripDirection{};
+  std::vector<std::vector<float>> *SPbottomStripDirection{};
+  std::vector<std::vector<float>> *SPstripCenterDistance{};
+  std::vector<std::vector<float>> *SPtopStripCenterPosition{};
 
   // Tracks
   Int_t nTRK = 0;
   Int_t TRKindex[maxTRK] = {};                //[nTRK]
   Int_t TRKtrack_fitter[maxTRK] = {};         //[nTRK]
   Int_t TRKparticle_hypothesis[maxTRK] = {};  //[nTRK]
-  std::vector<std::vector<int>> *TRKproperties;
-  std::vector<std::vector<int>> *TRKpattern;
+  std::vector<std::vector<int>> *TRKproperties{};
+  std::vector<std::vector<int>> *TRKpattern{};
   Int_t TRKndof[maxTRK] = {};     //[nTRK]
   Int_t TRKmot[maxTRK] = {};      //[nTRK]
   Int_t TRKoot[maxTRK] = {};      //[nTRK]
   Float_t TRKchiSq[maxTRK] = {};  //[nTRK]
-  std::vector<std::vector<int>> *TRKmeasurementsOnTrack_pixcl_sctcl_index;
-  std::vector<std::vector<int>> *TRKoutliersOnTrack_pixcl_sctcl_index;
+  std::vector<std::vector<int>> *TRKmeasurementsOnTrack_pixcl_sctcl_index{};
+  std::vector<std::vector<int>> *TRKoutliersOnTrack_pixcl_sctcl_index{};
   Int_t TRKcharge[maxTRK] = {};  //[nTRK]
-  std::vector<std::vector<double>> *TRKperigee_position;
-  std::vector<std::vector<double>> *TRKperigee_momentum;
+  std::vector<std::vector<double>> *TRKperigee_position{};
+  std::vector<std::vector<double>> *TRKperigee_momentum{};
   Int_t TTCindex[maxTRK] = {};          //[nTRK]
   Int_t TTCevent_index[maxTRK] = {};    //[nTRK]
   Int_t TTCparticle_link[maxTRK] = {};  //[nTRK]
@@ -293,10 +294,10 @@ class RootAthenaDumpReader : public IReader {
   Int_t nDTT = 0;
   Int_t DTTindex[maxDTT] = {};  //[nDTT]
   Int_t DTTsize[maxDTT] = {};   //[nDTT]
-  std::vector<std::vector<int>> *DTTtrajectory_eventindex;
-  std::vector<std::vector<int>> *DTTtrajectory_barcode;
-  std::vector<std::vector<int>> *DTTstTruth_subDetType;
-  std::vector<std::vector<int>> *DTTstTrack_subDetType;
-  std::vector<std::vector<int>> *DTTstCommon_subDetType;
+  std::vector<std::vector<int>> *DTTtrajectory_eventindex{};
+  std::vector<std::vector<int>> *DTTtrajectory_barcode{};
+  std::vector<std::vector<int>> *DTTstTruth_subDetType{};
+  std::vector<std::vector<int>> *DTTstTrack_subDetType{};
+  std::vector<std::vector<int>> *DTTstCommon_subDetType{};
 };
 }  // namespace ActsExamples