Core/Material/MaterialMapperTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Material/AccumulatedMaterialSlab.hpp"
0015 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0016 #include "Acts/Material/MaterialInteraction.hpp"
0017 #include "Acts/Material/MaterialMapper.hpp"
0018 #include "Acts/Material/MaterialSlab.hpp"
0019 #include "Acts/Material/interface/IAssignmentFinder.hpp"
0020 #include "Acts/Material/interface/ISurfaceMaterialAccumulater.hpp"
0021 #include "Acts/Propagator/SurfaceCollector.hpp"
0022 #include "Acts/Surfaces/CylinderSurface.hpp"
0023 #include "Acts/Utilities/Enumerate.hpp"
0024 #include "Acts/Utilities/VectorHelpers.hpp"
0025
0026 #include <limits>
0027
0028 namespace Acts::Test {
0029
0030 auto tContext = GeometryContext();
0031
0032 /// @brief Interface for the material mapping that seeks the possible
0033 /// assignment candidates for the material interactiosn
0034 class IntersectSurfacesFinder : public IAssignmentFinder {
0035  public:
0036   std::vector<const Surface*> surfaces;
0037
0038   /// @brief Interface method for generating assignment candidates for the
0039   /// material interaction assignment to surfaces or volumes
0040   ///
0041   /// @param gctx is the geometry context
0042   /// @param mctx is the magnetic field context
0043   /// @param position is the position of the initial ray
0044   /// @param direction is the direction of initial ray
0045   ///
0046   /// @return a vector of Surface Assignments and Volume Assignments
0047   std::pair<std::vector<IAssignmentFinder::SurfaceAssignment>,
0048             std::vector<IAssignmentFinder::VolumeAssignment>>
0049   assignmentCandidates(const GeometryContext& gctx,
0050                        const MagneticFieldContext& /*ignored*/,
0051                        const Vector3& position,
0052                        const Vector3& direction) const override {
0053     std::vector<IAssignmentFinder::SurfaceAssignment> surfaceAssignments;
0054     std::vector<IAssignmentFinder::VolumeAssignment> volumeAssignments;
0055     // Intersect the surfaces
0056     for (auto& surface : surfaces) {
0057       // Get the intersection
0058       auto sMultiIntersection = surface->intersect(gctx, position, direction,
0059                                                    BoundaryTolerance::None());
0060       // One solution, take it
0061       if (sMultiIntersection.size() == 1u &&
0062           sMultiIntersection[0u].status() >=
0063               Acts::IntersectionStatus::reachable &&
0064           sMultiIntersection[0u].pathLength() >= 0.0) {
0065         surfaceAssignments.push_back(
0066             {surface, sMultiIntersection[0u].position(), direction});
0067         continue;
0068       }
0069       if (sMultiIntersection.size() > 1u) {
0070         // Multiple intersections, take the closest
0071         auto closestForward = sMultiIntersection.closestForward();
0072         if (closestForward.status() >= Acts::IntersectionStatus::reachable &&
0073             closestForward.pathLength() > 0.0) {
0074           surfaceAssignments.push_back(
0075               {surface, closestForward.position(), direction});
0076           continue;
0077         }
0078       }
0079     }
0080     return {surfaceAssignments, volumeAssignments};
0081   }
0082 };
0083
0084 /// @brief Interface for the material mapping, this is the accumulation step
0085 class MaterialBlender : public ISurfaceMaterialAccumulater {
0086  public:
0087   explicit MaterialBlender(
0088       const std::vector<std::shared_ptr<Surface>>& surfaces = {})
0089       : m_surfaces(surfaces) {}
0090
0091   /// The state of the material accumulater, this is used
0092   /// to cache information across tracks/events
0093   class State final : public ISurfaceMaterialAccumulater::State {
0094    public:
0095     std::map<const Surface*, AccumulatedMaterialSlab> accumulatedMaterial;
0096   };
0097
0098   /// Factory for creating the state
0099   std::unique_ptr<ISurfaceMaterialAccumulater::State> createState()
0100       const override {
0101     auto state = std::make_unique<State>();
0102     for (auto& surface : m_surfaces) {
0103       state->accumulatedMaterial[surface.get()] = AccumulatedMaterialSlab();
0104     }
0105     return state;
0106   };
0107
0108   /// @brief Accumulate the material interaction on the surface
0109   ///
0110   /// @param state is the state of the accumulater
0111   /// @param interactions is the material interactions, with assigned surfaces
0112   /// @param surfacesWithoutAssignment are the surfaces without assignment
0113   ///
0114   /// @note this the track average over the binned material
0115   void accumulate(ISurfaceMaterialAccumulater::State& state,
0116                   const std::vector<MaterialInteraction>& interactions,
0117                   const std::vector<IAssignmentFinder::SurfaceAssignment>&
0118                   /*surfacesWithoutAssignment*/) const override {
0119     auto cState = static_cast<State*>(&state);
0120     for (const auto& mi : interactions) {
0121       // Get the surface
0122       const Surface* surface = mi.surface;
0123       // Get the accumulated material
0124       auto accMaterial = cState->accumulatedMaterial.find(surface);
0125       if (accMaterial == cState->accumulatedMaterial.end()) {
0126         throw std::invalid_argument(
0127             "Surface material is not found, inconsistent configuration.");
0128       }
0129       // Accumulate the material
0130       accMaterial->second.accumulate(mi.materialSlab);
0131     }
0132     // Average over the track
0133     for (auto& [surface, accumulatedMaterial] : cState->accumulatedMaterial) {
0134       accumulatedMaterial.trackAverage();
0135     }
0136   };
0137
0138   /// Finalize the surface material maps
0139   ///
0140   /// @param state the state of the accumulator
0141   ///
0142   /// @note this does the run average over the (binned) material
0143   std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
0144   finalizeMaterial(ISurfaceMaterialAccumulater::State& state) const override {
0145     auto cState = static_cast<State*>(&state);
0146
0147     std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
0148         materialMaps;
0149     for (auto& [surface, accumulatedMaterial] : cState->accumulatedMaterial) {
0150       materialMaps[surface->geometryId()] =
0151           std::make_shared<HomogeneousSurfaceMaterial>(
0152               accumulatedMaterial.totalAverage().first);
0153     }
0154     return materialMaps;
0155   }
0156
0157  private:
0158   std::vector<std::shared_ptr<Surface>> m_surfaces;
0159 };
0160
0161 BOOST_AUTO_TEST_SUITE(MaterialMapperTestSuite)
0162
0163 /// @brief This test checks the data flow of the material mapper, it is not
0164 /// a test of the single components, which are tested individually
0165 ///
0166 /// @note Currently only surface mapping is implemented, volume mapping/assigning
0167 /// is being added in future PRs
0168 BOOST_AUTO_TEST_CASE(MaterialMapperFlowTest) {
0169   // Create a vector of surfaces
0170   std::vector<std::shared_ptr<Surface>> surfaces = {
0171       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0),
0172       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30.0, 100.0),
0173       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 50.0,
0174                                            100.0)};
0175
0176   for (auto [is, surface] : enumerate(surfaces)) {
0177     surface->assignGeometryId(GeometryIdentifier().withSensitive(is + 1));
0178   }
0179
0180   // The assigner
0181   auto assigner = std::make_shared<IntersectSurfacesFinder>();
0182   assigner->surfaces = {surfaces[0].get(), surfaces[1].get(),
0183                         surfaces[2].get()};
0184
0185   // The accumulater - which blends all the material
0186   auto accumulator = std::make_shared<MaterialBlender>(surfaces);
0187
0188   // Create the mapper
0189   MaterialMapper::Config mmConfig;
0190   mmConfig.assignmentFinder = assigner;
0191   mmConfig.surfaceMaterialAccumulater = accumulator;
0192
0193   MaterialMapper mapper(mmConfig);
0194
0195   auto state = mapper.createState();
0196   BOOST_CHECK(state.get() != nullptr);
0197   BOOST_CHECK(state->surfaceMaterialAccumulaterState.get() != nullptr);
0198
0199   std::vector<RecordedMaterialTrack> mappedTracks;
0200   std::vector<RecordedMaterialTrack> unmappedTracks;
0201
0202   // Track loop
0203   Vector3 position(0., 0., 0.);
0204   for (unsigned int it = 0; it < 11; ++it) {
0205     Vector3 direction =
0206         Vector3(0.9 + it * 0.02, 1.1 - it * 0.02, 0.).normalized();
0207     RecordedMaterialTrack mTrack{{position, direction}, {}};
0208     for (unsigned int im = 0; im < 60; ++im) {
0209       MaterialInteraction mi;
0210       mi.materialSlab = MaterialSlab(
0211           Material::fromMassDensity(it + 1, it + 1, it + 1, it + 1, it + 1),
0212           0.1);
0213       mi.position = position + (im + 1) * direction;
0214       mi.direction = direction;
0215       mTrack.second.materialInteractions.push_back(mi);
0216     }
0217     auto [mapped, unmapped] = mapper.mapMaterial(*state, tContext, {}, mTrack);
0218     mappedTracks.push_back(mapped);
0219     unmappedTracks.push_back(unmapped);
0220   }
0221
0222   // Get the maps
0223   auto [surfaceMaps, volumeMaps] = mapper.finalizeMaps(*state);
0224
0225   BOOST_CHECK(surfaceMaps.size() == 3);
0226   BOOST_CHECK(volumeMaps.empty());
0227 }
0228
0229 BOOST_AUTO_TEST_CASE(MaterialMapperInvalidTest) {
0230   // The assigner
0231   auto assigner = std::make_shared<IntersectSurfacesFinder>();
0232
0233   // The accumulater - which blends all the material
0234   auto accumulator = std::make_shared<MaterialBlender>();
0235
0236   // Create the mapper
0237   MaterialMapper::Config mmConfigInvalid;
0238   mmConfigInvalid.assignmentFinder = nullptr;
0239   mmConfigInvalid.surfaceMaterialAccumulater = accumulator;
0240
0241   BOOST_CHECK_THROW(auto mapperIA = MaterialMapper(mmConfigInvalid),
0242                     std::invalid_argument);
0243
0244   // BOOST_CHECK_THROW(MaterialMapper(mmConfigInvalid), std::invalid_argument);
0245
0246   mmConfigInvalid.assignmentFinder = assigner;
0247   mmConfigInvalid.surfaceMaterialAccumulater = nullptr;
0248
0249   BOOST_CHECK_THROW(auto mapperIS = MaterialMapper(mmConfigInvalid),
0250                     std::invalid_argument);
0251 }
0252
0253 BOOST_AUTO_TEST_SUITE_END()
0254
0255 }  // namespace Acts::Test