Merge branch 'main' of github.com:acts-project/acts into refactor-mat…

…erial-interactor

Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp

@@ -98,6 +98,9 @@ struct BetheHeitlerApproxSingleCmp {
 /// mixture. To enable an approximation for continuous input variables, the
 /// weights, means and variances are internally parametrized as a Nth order
 /// polynomial.
+/// @todo This class is rather inflexible: It forces two data representations,
+/// making it a bit awkward to add a single parameterization. It would be good
+/// to generalize this at some point.
 template <int NComponents, int PolyDegree>
 class AtlasBetheHeitlerApprox {
   static_assert(NComponents > 0);
@@ -112,40 +115,47 @@ class AtlasBetheHeitlerApprox {
 
   using Data = std::array<PolyData, NComponents>;
 
-  constexpr static double noChangeLimit = 0.0001;
-  constexpr static double singleGaussianLimit = 0.002;
-  constexpr static double lowerLimit = 0.10;
-  constexpr static double higherLimit = 0.20;
-
  private:
-  Data m_low_data;
-  Data m_high_data;
-  bool m_low_transform;
-  bool m_high_transform;
+  Data m_lowData;
+  Data m_highData;
+  bool m_lowTransform;
+  bool m_highTransform;
+
+  constexpr static double m_noChangeLimit = 0.0001;
+  constexpr static double m_singleGaussianLimit = 0.002;
+  double m_lowLimit = 0.10;
+  double m_highLimit = 0.20;
 
  public:
-  /// Construct the Bethe-Heitler approximation description. Additional to the
-  /// coefficients of the polynomials, the information whether these values need
-  /// to be transformed beforehand must be given (see ATLAS code).
+  /// Construct the Bethe-Heitler approximation description with two
+  /// parameterizations, one for lower ranges, one for higher ranges.
+  /// Is it assumed that the lower limit of the high-x/x0 data is equal
+  /// to the upper limit of the low-x/x0 data.
   ///
-  /// @param low_data data for the lower x/x0 range
-  /// @param high_data data for the higher x/x0 range
-  /// @param low_transform whether the low data need to be transformed
-  /// @param high_transform whether the high data need to be transformed
-  constexpr AtlasBetheHeitlerApprox(const Data &low_data, const Data &high_data,
-                                    bool low_transform, bool high_transform)
-      : m_low_data(low_data),
-        m_high_data(high_data),
-        m_low_transform(low_transform),
-        m_high_transform(high_transform) {}
+  /// @param lowData data for the lower x/x0 range
+  /// @param highData data for the higher x/x0 range
+  /// @param lowTransform whether the low data need to be transformed
+  /// @param highTransform whether the high data need to be transformed
+  /// @param lowLimit the upper limit for the low data
+  /// @param highLimit the upper limit for the high data
+  constexpr AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData,
+                                    bool lowTransform, bool highTransform,
+                                    double lowLimit = 0.1,
+                                    double highLimit = 0.2)
+      : m_lowData(lowData),
+        m_highData(highData),
+        m_lowTransform(lowTransform),
+        m_highTransform(highTransform),
+        m_lowLimit(lowLimit),
+        m_highLimit(highLimit) {}
 
   /// Returns the number of components the returned mixture will have
   constexpr auto numComponents() const { return NComponents; }
 
   /// Checks if an input is valid for the parameterization
   ///
   /// @param x pathlength in terms of the radiation length
-  constexpr bool validXOverX0(ActsScalar x) const { return x < higherLimit; }
+  constexpr bool validXOverX0(ActsScalar x) const { return x < m_highLimit; }
 
   /// Generates the mixture from the polynomials and reweights them, so
   /// that the sum of all weights is 1
@@ -194,7 +204,7 @@ class AtlasBetheHeitlerApprox {
     };
 
     // Return no change
-    if (x < noChangeLimit) {
+    if (x < m_noChangeLimit) {
       Array ret(1);
 
       ret[0].weight = 1.0;
@@ -204,19 +214,19 @@ class AtlasBetheHeitlerApprox {
       return ret;
     }
     // Return single gaussian approximation
-    if (x < singleGaussianLimit) {
+    if (x < m_singleGaussianLimit) {
       Array ret(1);
       ret[0] = BetheHeitlerApproxSingleCmp::mixture(x)[0];
       return ret;
     }
     // Return a component representation for lower x0
-    if (x < lowerLimit) {
-      return make_mixture(m_low_data, x, m_low_transform);
+    if (x < m_lowLimit) {
+      return make_mixture(m_lowData, x, m_lowTransform);
     }
     // Return a component representation for higher x0
     // Cap the x because beyond the parameterization goes wild
-    const auto high_x = std::min(higherLimit, x);
-    return make_mixture(m_high_data, high_x, m_high_transform);
+    const auto high_x = std::min(m_highLimit, x);
+    return make_mixture(m_highData, high_x, m_highTransform);
   }
 
   /// Loads a parameterization from a file according to the Atlas file
@@ -226,8 +236,11 @@ class AtlasBetheHeitlerApprox {
   /// the parameterization for low x/x0
   /// @param high_parameters_path Path to the foo.par file that stores
   /// the parameterization for high x/x0
+  /// @param lowLimit the upper limit for the low x/x0-data
+  /// @param highLimit the upper limit for the high x/x0-data
   static auto loadFromFiles(const std::string &low_parameters_path,
-                            const std::string &high_parameters_path) {
+                            const std::string &high_parameters_path,
+                            double lowLimit = 0.1, double highLimit = 0.2) {
     auto read_file = [](const std::string &filepath) {
       std::ifstream file(filepath);
 
@@ -267,11 +280,11 @@ class AtlasBetheHeitlerApprox {
       return std::make_tuple(data, transform_code);
     };
 
-    const auto [low_data, low_transform] = read_file(low_parameters_path);
-    const auto [high_data, high_transform] = read_file(high_parameters_path);
+    const auto [lowData, lowTransform] = read_file(low_parameters_path);
+    const auto [highData, highTransform] = read_file(high_parameters_path);
 
-    return AtlasBetheHeitlerApprox(low_data, high_data, low_transform,
-                                   high_transform);
+    return AtlasBetheHeitlerApprox(lowData, highData, lowTransform,
+                                   highTransform, lowLimit, highLimit);
   }
 };
 

Core/include/Acts/TrackFitting/GaussianSumFitter.hpp

@@ -322,6 +322,7 @@ struct GaussianSumFitter {
     ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates);
 
     std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler;
+    double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0;
 
     //////////////////
     // Backward pass
@@ -408,13 +409,16 @@ struct GaussianSumFitter {
     }
 
     nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler;
+    maxPathXOverX0 = std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0);
 
     if (nInvalidBetheHeitler > 0) {
-      ACTS_WARNING("Encountered "
-                   << nInvalidBetheHeitler
-                   << " cases where the material thickness exceeds the range "
-                      "of the Bethe-Heitler-Approximation. Enable DEBUG output "
-                      "for more information.");
+      ACTS_WARNING("Encountered " << nInvalidBetheHeitler
+                                  << " cases where x/X0 exceeds the range "
+                                     "of the Bethe-Heitler-Approximation. The "
+                                     "maximum x/X0 encountered was "
+                                  << maxPathXOverX0
+                                  << ". Enable DEBUG output "
+                                     "for more information.");
     }
 
     ////////////////////////////////////

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

@@ -56,7 +56,9 @@ struct GsfResult {
   std::vector<const Acts::Surface*> visitedSurfaces;
   std::vector<const Acts::Surface*> surfacesVisitedBwdAgain;
 
+  /// Statistics about encounterings of invalid bethe heitler approximations
   std::size_t nInvalidBetheHeitler = 0;
+  double maxPathXOverX0 = 0.0;
 
   // Propagate potential errors to the outside
   Result<void> result{Result<void>::success()};
@@ -287,7 +289,7 @@ struct GsfActor {
       componentCache.clear();
 
       convoluteComponents(state, stepper, navigator, tmpStates, componentCache,
-                          result.nInvalidBetheHeitler);
+                          result);
 
       if (componentCache.empty()) {
         ACTS_WARNING(
@@ -328,7 +330,7 @@ struct GsfActor {
                            const navigator_t& navigator,
                            const TemporaryStates& tmpStates,
                            std::vector<ComponentCache>& componentCache,
-                           std::size_t& nInvalidBetheHeitler) const {
+                           result_type& result) const {
     auto cmps = stepper.componentIterable(state.stepping);
     for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
       auto proxy = tmpStates.traj.getTrackState(idx);
@@ -338,7 +340,7 @@ struct GsfActor {
                                  stepper.particleHypothesis(state.stepping));
 
       applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx),
-                        componentCache, nInvalidBetheHeitler);
+                        componentCache, result);
     }
   }
 
@@ -348,7 +350,7 @@ struct GsfActor {
                          const BoundTrackParameters& old_bound,
                          const double old_weight,
                          std::vector<ComponentCache>& componentCaches,
-                         std::size_t& nInvalidBetheHeitler) const {
+                         result_type& result) const {
     const auto& surface = *navigator.currentSurface(state.navigation);
     const auto p_prev = old_bound.absoluteMomentum();
 
@@ -357,22 +359,24 @@ struct GsfActor {
         old_bound.position(state.stepping.geoContext), state.options.direction,
         MaterialUpdateStage::FullUpdate);
 
-    auto pathCorrection = surface.pathCorrection(
+    const auto pathCorrection = surface.pathCorrection(
         state.stepping.geoContext,
         old_bound.position(state.stepping.geoContext), old_bound.direction());
     slab.scaleThickness(pathCorrection);
 
+    const double pathXOverX0 = slab.thicknessInX0();
+
     // Emit a warning if the approximation is not valid for this x/x0
-    if (!m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0())) {
-      ++nInvalidBetheHeitler;
+    if (!m_cfg.bethe_heitler_approx->validXOverX0(pathXOverX0)) {
+      ++result.nInvalidBetheHeitler;
+      result.maxPathXOverX0 = std::max(result.maxPathXOverX0, pathXOverX0);
       ACTS_DEBUG(
           "Bethe-Heitler approximation encountered invalid value for x/x0="
-          << slab.thicknessInX0() << " at surface " << surface.geometryId());
+          << pathXOverX0 << " at surface " << surface.geometryId());
     }
 
     // Get the mixture
-    const auto mixture =
-        m_cfg.bethe_heitler_approx->mixture(slab.thicknessInX0());
+    const auto mixture = m_cfg.bethe_heitler_approx->mixture(pathXOverX0);
 
     // Create all possible new components
     for (const auto& gaussian : mixture) {

Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp

@@ -58,12 +58,12 @@ class AdaptiveMultiVertexFinder {
     /// @param ipEst ImpactPointEstimator
     /// @param lin Track linearizer
     /// @param bIn Input magnetic field
-    Config(vfitter_t fitter, const sfinder_t& sfinder,
-           const ImpactPointEstimator<InputTrack_t, Propagator_t>& ipEst,
+    Config(vfitter_t fitter, sfinder_t sfinder,
+           ImpactPointEstimator<InputTrack_t, Propagator_t> ipEst,
            Linearizer_t lin, std::shared_ptr<const MagneticFieldProvider> bIn)
         : vertexFitter(std::move(fitter)),
-          seedFinder(sfinder),
-          ipEstimator(ipEst),
+          seedFinder(std::move(sfinder)),
+          ipEstimator(std::move(ipEst)),
           linearizer(std::move(lin)),
           bField{std::move(bIn)} {}
 

Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp

@@ -56,7 +56,7 @@ auto Acts::AdaptiveMultiVertexFinder<vfitter_t, sfinder_t>::find(
     Vertex<InputTrack_t>& vtxCandidate = *allVertices.back();
     allVerticesPtr.push_back(&vtxCandidate);
 
-    ACTS_DEBUG("Position of current vertex candidate after seeding: "
+    ACTS_DEBUG("Position of vertex candidate after seeding: "
                << vtxCandidate.fullPosition().transpose());
     if (vtxCandidate.position().z() ==
         vertexingOptions.constraint.position().z()) {
@@ -79,7 +79,8 @@ auto Acts::AdaptiveMultiVertexFinder<vfitter_t, sfinder_t>::find(
       return prepResult.error();
     }
     if (!(*prepResult)) {
-      ACTS_DEBUG("Could not prepare for fit anymore. Break.");
+      ACTS_DEBUG(
+          "Could not prepare for fit. Discarding the vertex candindate.");
       allVertices.pop_back();
       allVerticesPtr.pop_back();
       break;
@@ -93,7 +94,7 @@ auto Acts::AdaptiveMultiVertexFinder<vfitter_t, sfinder_t>::find(
     if (!fitResult.ok()) {
       return fitResult.error();
     }
-    ACTS_DEBUG("New position of current vertex candidate after fit: "
+    ACTS_DEBUG("Position of vertex candidate after the fit: "
                << vtxCandidate.fullPosition().transpose());
     // Check if vertex is good vertex
     auto [nCompatibleTracks, isGoodVertex] =

Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp

@@ -119,7 +119,7 @@ class AdaptiveMultiVertexFitter {
     /// @brief Config constructor
     ///
     /// @param est ImpactPointEstimator
-    Config(const IPEstimator& est) : ipEst(est) {}
+    Config(IPEstimator est) : ipEst(std::move(est)) {}
 
     // ImpactPointEstimator
     IPEstimator ipEst;
@@ -297,6 +297,13 @@ class AdaptiveMultiVertexFitter {
   ///
   /// @param state Fitter state
   void doVertexSmoothing(State& state) const;
+
+  /// @brief Logs vertices in state.vertexCollection and associated tracks
+  ///
+  /// @param state Fitter state
+  /// @param geoContext Geometry context
+  void logDebugData(const State& state,
+                    const GeometryContext& geoContext) const;
 };
 
 }  // namespace Acts