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