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   MaterialBlender(const std::vector<std::shared_ptr<Surface>>& surfaces = {})
0088       : m_surfaces(surfaces) {}
0089
0090   /// The state of the material accumulater, this is used
0091   /// to cache information across tracks/events
0092   class State final : public ISurfaceMaterialAccumulater::State {
0093    public:
0094     std::map<const Surface*, AccumulatedMaterialSlab> accumulatedMaterial;
0095   };
0096
0097   /// Factory for creating the state
0098   std::unique_ptr<ISurfaceMaterialAccumulater::State> createState()
0099       const override {
0100     auto state = std::make_unique<State>();
0101     for (auto& surface : m_surfaces) {
0102       state->accumulatedMaterial[surface.get()] = AccumulatedMaterialSlab();
0103     }
0104     return state;
0105   };
0106
0107   /// @brief Accumulate the material interaction on the surface
0108   ///
0109   /// @param state is the state of the accumulater
0110   /// @param interactions is the material interactions, with assigned surfaces
0111   /// @param surfacesWithoutAssignment are the surfaces without assignment
0112   ///
0113   /// @note this the track average over the binned material
0114   void accumulate(ISurfaceMaterialAccumulater::State& state,
0115                   const std::vector<MaterialInteraction>& interactions,
0116                   const std::vector<IAssignmentFinder::SurfaceAssignment>&
0117                   /*surfacesWithoutAssignment*/) const override {
0118     auto cState = static_cast<State*>(&state);
0119     for (const auto& mi : interactions) {
0120       // Get the surface
0121       const Surface* surface = mi.surface;
0122       // Get the accumulated material
0123       auto accMaterial = cState->accumulatedMaterial.find(surface);
0124       if (accMaterial == cState->accumulatedMaterial.end()) {
0125         throw std::invalid_argument(
0126             "Surface material is not found, inconsistent configuration.");
0127       }
0128       // Accumulate the material
0129       accMaterial->second.accumulate(mi.materialSlab);
0130     }
0131     // Average over the track
0132     for (auto& [surface, accumulatedMaterial] : cState->accumulatedMaterial) {
0133       accumulatedMaterial.trackAverage();
0134     }
0135   };
0136
0137   /// Finalize the surface material maps
0138   ///
0139   /// @param state the state of the accumulator
0140   ///
0141   /// @note this does the run average over the (binned) material
0142   std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
0143   finalizeMaterial(ISurfaceMaterialAccumulater::State& state) const override {
0144     auto cState = static_cast<State*>(&state);
0145
0146     std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
0147         materialMaps;
0148     for (auto& [surface, accumulatedMaterial] : cState->accumulatedMaterial) {
0149       materialMaps[surface->geometryId()] =
0150           std::make_shared<HomogeneousSurfaceMaterial>(
0151               accumulatedMaterial.totalAverage().first);
0152     }
0153     return materialMaps;
0154   }
0155
0156  private:
0157   std::vector<std::shared_ptr<Surface>> m_surfaces;
0158 };
0159
0160 BOOST_AUTO_TEST_SUITE(MaterialMapperTestSuite)
0161
0162 /// @brief This test checks the data flow of the material mapper, it is not
0163 /// a test of the single components, which are tested individually
0164 ///
0165 /// @note Currently only surface mapping is implemented, volume mapping/assigning
0166 /// is being added in future PRs
0167 BOOST_AUTO_TEST_CASE(MaterialMapperFlowTest) {
0168   // Create a vector of surfaces
0169   std::vector<std::shared_ptr<Surface>> surfaces = {
0170       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0),
0171       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30.0, 100.0),
0172       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 50.0,
0173                                            100.0)};
0174
0175   for (auto [is, surface] : enumerate(surfaces)) {
0176     surface->assignGeometryId(GeometryIdentifier().setSensitive(is + 1));
0177   }
0178
0179   // The assigner
0180   auto assigner = std::make_shared<IntersectSurfacesFinder>();
0181   assigner->surfaces = {surfaces[0].get(), surfaces[1].get(),
0182                         surfaces[2].get()};
0183
0184   // The accumulater - which blends all the material
0185   auto accumulator = std::make_shared<MaterialBlender>(surfaces);
0186
0187   // Create the mapper
0188   MaterialMapper::Config mmConfig;
0189   mmConfig.assignmentFinder = assigner;
0190   mmConfig.surfaceMaterialAccumulater = accumulator;
0191
0192   MaterialMapper mapper(mmConfig);
0193
0194   auto state = mapper.createState();
0195   BOOST_CHECK(state.get() != nullptr);
0196   BOOST_CHECK(state->surfaceMaterialAccumulaterState.get() != nullptr);
0197
0198   std::vector<RecordedMaterialTrack> mappedTracks;
0199   std::vector<RecordedMaterialTrack> unmappedTracks;
0200
0201   // Track loop
0202   Vector3 position(0., 0., 0.);
0203   for (unsigned int it = 0; it < 11; ++it) {
0204     Vector3 direction =
0205         Vector3(0.9 + it * 0.02, 1.1 - it * 0.02, 0.).normalized();
0206     RecordedMaterialTrack mTrack{{position, direction}, {}};
0207     for (unsigned int im = 0; im < 60; ++im) {
0208       MaterialInteraction mi;
0209       mi.materialSlab = MaterialSlab(
0210           Material::fromMassDensity(it + 1, it + 1, it + 1, it + 1, it + 1),
0211           0.1);
0212       mi.position = position + (im + 1) * direction;
0213       mi.direction = direction;
0214       mTrack.second.materialInteractions.push_back(mi);
0215     }
0216     auto [mapped, unmapped] = mapper.mapMaterial(*state, tContext, {}, mTrack);
0217     mappedTracks.push_back(mapped);
0218     unmappedTracks.push_back(unmapped);
0219   }
0220
0221   // Get the maps
0222   auto [surfaceMaps, volumeMaps] = mapper.finalizeMaps(*state);
0223
0224   BOOST_CHECK(surfaceMaps.size() == 3);
0225   BOOST_CHECK(volumeMaps.empty());
0226 }
0227
0228 BOOST_AUTO_TEST_CASE(MaterialMapperInvalidTest) {
0229   // The assigner
0230   auto assigner = std::make_shared<IntersectSurfacesFinder>();
0231
0232   // The accumulater - which blends all the material
0233   auto accumulator = std::make_shared<MaterialBlender>();
0234
0235   // Create the mapper
0236   MaterialMapper::Config mmConfigInvalid;
0237   mmConfigInvalid.assignmentFinder = nullptr;
0238   mmConfigInvalid.surfaceMaterialAccumulater = accumulator;
0239
0240   BOOST_CHECK_THROW(auto mapperIA = MaterialMapper(mmConfigInvalid),
0241                     std::invalid_argument);
0242
0243   // BOOST_CHECK_THROW(MaterialMapper(mmConfigInvalid), std::invalid_argument);
0244
0245   mmConfigInvalid.assignmentFinder = assigner;
0246   mmConfigInvalid.surfaceMaterialAccumulater = nullptr;
0247
0248   BOOST_CHECK_THROW(auto mapperIS = MaterialMapper(mmConfigInvalid),
0249                     std::invalid_argument);
0250 }
0251
0252 BOOST_AUTO_TEST_SUITE_END()
0253
0254 }  // namespace Acts::Test