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/ISurfaceMaterialAccumulator.hpp"
0021 #include "Acts/Surfaces/CylinderSurface.hpp"
0022 #include "Acts/Utilities/Enumerate.hpp"
0023 #include "Acts/Utilities/Intersection.hpp"
0024
0025 using namespace Acts;
0026
0027 namespace ActsTests {
0028
0029 auto tContext = GeometryContext::dangerouslyDefaultConstruct();
0030
0031 /// @brief Interface for the material mapping that seeks the possible
0032 /// assignment candidates for the material interactiosn
0033 class IntersectSurfacesFinder : public IAssignmentFinder {
0034  public:
0035   std::vector<const Surface*> surfaces;
0036
0037   /// @brief Interface method for generating assignment candidates for the
0038   /// material interaction assignment to surfaces or volumes
0039   ///
0040   /// @param gctx is the geometry context
0041   /// @param mctx is the magnetic field context
0042   /// @param position is the position of the initial ray
0043   /// @param direction is the direction of initial ray
0044   ///
0045   /// @return a vector of Surface Assignments and Volume Assignments
0046   std::pair<std::vector<IAssignmentFinder::SurfaceAssignment>,
0047             std::vector<IAssignmentFinder::VolumeAssignment>>
0048   assignmentCandidates(const GeometryContext& gctx,
0049                        const MagneticFieldContext& /*ignored*/,
0050                        const Vector3& position,
0051                        const Vector3& direction) const override {
0052     std::vector<IAssignmentFinder::SurfaceAssignment> surfaceAssignments;
0053     std::vector<IAssignmentFinder::VolumeAssignment> volumeAssignments;
0054     // Intersect the surfaces
0055     for (auto& surface : surfaces) {
0056       // Get the intersection
0057       MultiIntersection3D multiIntersection = surface->intersect(
0058           gctx, position, direction, BoundaryTolerance::None());
0059       // One solution, take it
0060       if (multiIntersection.size() == 1u &&
0061           multiIntersection.at(0).status() >= IntersectionStatus::reachable &&
0062           multiIntersection.at(0).pathLength() >= 0.0) {
0063         surfaceAssignments.push_back(
0064             {surface, multiIntersection.at(0).position(), direction});
0065         continue;
0066       }
0067       if (multiIntersection.size() > 1u) {
0068         // Multiple intersections, take the closest
0069         Intersection3D closestForward = multiIntersection.closestForward();
0070         if (closestForward.status() >= IntersectionStatus::reachable &&
0071             closestForward.pathLength() > 0.0) {
0072           surfaceAssignments.push_back(
0073               {surface, closestForward.position(), direction});
0074           continue;
0075         }
0076       }
0077     }
0078     return {surfaceAssignments, volumeAssignments};
0079   }
0080 };
0081
0082 /// @brief Interface for the material mapping, this is the accumulation step
0083 class MaterialBlender : public ISurfaceMaterialAccumulator {
0084  public:
0085   explicit MaterialBlender(
0086       const std::vector<std::shared_ptr<Surface>>& surfaces = {})
0087       : m_surfaces(surfaces) {}
0088
0089   /// The state of the material accumulator, this is used
0090   /// to cache information across tracks/events
0091   class State final : public ISurfaceMaterialAccumulator::State {
0092    public:
0093     std::map<const Surface*, AccumulatedMaterialSlab> accumulatedMaterial;
0094   };
0095
0096   /// Factory for creating the state
0097   std::unique_ptr<ISurfaceMaterialAccumulator::State> createState(
0098       const GeometryContext& /*gctx*/) const override {
0099     auto state = std::make_unique<State>();
0100     for (auto& surface : m_surfaces) {
0101       state->accumulatedMaterial[surface.get()] = AccumulatedMaterialSlab();
0102     }
0103     return state;
0104   };
0105
0106   /// @brief Accumulate the material interaction on the surface
0107   ///
0108   /// @param state is the state of the accumulator
0109   /// @param interactions is the material interactions, with assigned surfaces
0110   /// @param surfacesWithoutAssignment are the surfaces without assignment
0111   ///
0112   /// @note this the track average over the binned material
0113   void accumulate(ISurfaceMaterialAccumulator::State& state,
0114                   const GeometryContext& /*gctx*/,
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(ISurfaceMaterialAccumulator::State& state,
0144                    const GeometryContext& /*gctx*/) 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(MaterialSuite)
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 accumulator - 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.surfaceMaterialAccumulator = accumulator;
0192
0193   MaterialMapper mapper(mmConfig);
0194
0195   auto state = mapper.createState(tContext);
0196   BOOST_CHECK(state.get() != nullptr);
0197   BOOST_CHECK(state->surfaceMaterialAccumulatorState.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, tContext);
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 accumulator - 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.surfaceMaterialAccumulator = 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.surfaceMaterialAccumulator = 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 ActsTests