Core/Material/MaterialInteractionAssignmentTests.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/Material.hpp"
0015 #include "Acts/Material/MaterialInteraction.hpp"
0016 #include "Acts/Material/MaterialInteractionAssignment.hpp"
0017 #include "Acts/Material/MaterialSlab.hpp"
0018 #include "Acts/Material/interface/IAssignmentFinder.hpp"
0019 #include "Acts/Propagator/SurfaceCollector.hpp"
0020 #include "Acts/Surfaces/CylinderSurface.hpp"
0021 #include "Acts/Utilities/Enumerate.hpp"
0022 #include "Acts/Utilities/VectorHelpers.hpp"
0023
0024 #include <limits>
0025
0026 using namespace Acts;
0027
0028 namespace ActsTests {
0029
0030 auto tContext = GeometryContext();
0031
0032 BOOST_AUTO_TEST_SUITE(MaterialSuite)
0033
0034 BOOST_AUTO_TEST_CASE(AssignToClosest) {
0035   // Create a vector of surfaces
0036   std::vector<std::shared_ptr<Surface>> surfaces = {
0037       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0),
0038       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30.0, 100.0),
0039       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 50.0,
0040                                            100.0)};
0041
0042   for (auto [is, surface] : enumerate(surfaces)) {
0043     surface->assignGeometryId(GeometryIdentifier().withSensitive(is + 1));
0044   }
0045
0046   std::vector<IAssignmentFinder::SurfaceAssignment> intersectedSurfaces = {
0047       {surfaces[0].get(), {20., 0., 0.}, {1., 0., 0.}},
0048       {surfaces[1].get(), {30., 0., 0.}, {1., 0., 0.}},
0049       {surfaces[2].get(), {50., 0., 0.}, {1., 0., 0.}}};
0050
0051   // Create a material
0052   Material material = Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0);
0053
0054   std::vector<MaterialInteraction> materialInteractions;
0055   materialInteractions.reserve(50);
0056   // Generate some material interactions
0057   for (unsigned int i = 1; i < 50; ++i) {
0058     MaterialInteraction materialInteraction;
0059     materialInteraction.materialSlab = MaterialSlab(material, 0.1);
0060     materialInteraction.position = Vector3{i * 1.0, 0, 0};
0061     materialInteraction.direction = Vector3{1.0, 0., 0.};
0062     materialInteractions.push_back(materialInteraction);
0063   }
0064
0065   MaterialInteractionAssignment::Options options;
0066
0067   // Assign the material interactions to the surface hits
0068   auto [assigned, unassigned, surfacesLeft] =
0069       MaterialInteractionAssignment::assign(tContext, materialInteractions,
0070                                             intersectedSurfaces, options);
0071
0072   // Check that the material interaction was assigned
0073   BOOST_CHECK_EQUAL(assigned.size(), materialInteractions.size());
0074   BOOST_CHECK_EQUAL(unassigned.size(), 0u);
0075   BOOST_CHECK_EQUAL(surfacesLeft.size(), 0u);
0076
0077   // Check that it is assigned to the closest surface always
0078   for (const auto& mi : assigned) {
0079     double minDistance = std::numeric_limits<double>::max();
0080     const Surface* closestSurface = nullptr;
0081     for (const auto& [surface, position, direction] : intersectedSurfaces) {
0082       double distance = (mi.position - position).norm();
0083       if (distance < minDistance) {
0084         minDistance = distance;
0085         closestSurface = surface;
0086       }
0087     }
0088     BOOST_CHECK_EQUAL(mi.surface, closestSurface);
0089   }
0090 }
0091
0092 BOOST_AUTO_TEST_CASE(AssignToClosest_withGlobalVeto) {
0093   // Create a vector of surfaces
0094   std::vector<std::shared_ptr<Surface>> surfaces = {
0095       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0),
0096       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30.0, 100.0),
0097       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 50.0,
0098                                            100.0)};
0099
0100   for (auto [is, surface] : enumerate(surfaces)) {
0101     surface->assignGeometryId(GeometryIdentifier().withSensitive(is + 1));
0102   }
0103
0104   std::vector<IAssignmentFinder::SurfaceAssignment> intersectedSurfaces = {
0105       {surfaces[0].get(), {20., 0., 0.}, {1., 0., 0.}},
0106       {surfaces[1].get(), {30., 0., 0.}, {1., 0., 0.}},
0107       {surfaces[2].get(), {50., 0., 0.}, {1., 0., 0.}}};
0108
0109   // Create a material
0110   Material material = Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0);
0111
0112   std::vector<MaterialInteraction> materialInteractions;
0113   materialInteractions.reserve(50);
0114   // Generate some material interactions
0115   for (unsigned int i = 1; i < 50; ++i) {
0116     MaterialInteraction materialInteraction;
0117     materialInteraction.materialSlab = MaterialSlab(material, 0.1);
0118     materialInteraction.position = Vector3{i * 1.0, 0, 0};
0119     materialInteraction.direction = Vector3{1.0, 0., 0.};
0120     materialInteractions.push_back(materialInteraction);
0121   }
0122
0123   // Veto everything above 40 mm
0124   struct RadialVeto {
0125     double rMax = 40.0;
0126     bool operator()(const MaterialInteraction& mi) const {
0127       return VectorHelpers::perp(mi.position) > rMax;
0128     }
0129   };
0130   MaterialInteractionAssignment::Options options;
0131   options.globalVetos.push_back(RadialVeto{40});
0132
0133   // Assign the material interactions to the surface hits
0134   auto [assigned, unassigned, surfacesLeft] =
0135       MaterialInteractionAssignment::assign(tContext, materialInteractions,
0136                                             intersectedSurfaces, options);
0137
0138   // Check that the material interaction was assigned
0139   BOOST_CHECK_EQUAL(assigned.size(), 40u);
0140   BOOST_CHECK_EQUAL(unassigned.size(), 9u);
0141   BOOST_CHECK_EQUAL(surfacesLeft.size(), 1u);
0142 }
0143
0144 BOOST_AUTO_TEST_CASE(AssignToClosest_withLocalVeto) {
0145   // Create a vector of surfaces
0146   std::vector<std::shared_ptr<Surface>> surfaces = {
0147       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0),
0148       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30.0, 100.0),
0149       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 50.0,
0150                                            100.0)};
0151
0152   for (auto [is, surface] : enumerate(surfaces)) {
0153     surface->assignGeometryId(GeometryIdentifier().withSensitive(is + 1));
0154   }
0155
0156   std::vector<IAssignmentFinder::SurfaceAssignment> intersectedSurfaces = {
0157       {surfaces[0].get(), {20., 0., 0.}, {1., 0., 0.}},
0158       {surfaces[1].get(), {30., 0., 0.}, {1., 0., 0.}},
0159       {surfaces[2].get(), {50., 0., 0.}, {1., 0., 0.}}};
0160
0161   // Create a material
0162   Material material = Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0);
0163
0164   std::vector<MaterialInteraction> materialInteractions;
0165   materialInteractions.reserve(50);
0166   // Generate some material interactions
0167   for (unsigned int i = 1; i < 50; ++i) {
0168     MaterialInteraction materialInteraction;
0169     materialInteraction.materialSlab = MaterialSlab(material, 0.1);
0170     materialInteraction.position = Vector3{i * 1.0, 0, 0};
0171     materialInteraction.direction = Vector3{1.0, 0., 0.};
0172     materialInteractions.push_back(materialInteraction);
0173   }
0174
0175   // Veto in a specific one
0176   struct VetoThisOne {
0177     bool operator()(
0178         const MaterialInteraction& /*m*/,
0179         const IAssignmentFinder::SurfaceAssignment& /*suggestedAssignment*/)
0180         const {
0181       return true;
0182     }
0183   };
0184
0185   // We assign this to
0186   std::vector<
0187       std::pair<GeometryIdentifier, MaterialInteractionAssignment::LocalVeto>>
0188       localVetoVector = {
0189           {GeometryIdentifier().withSensitive(2), VetoThisOne{}}};
0190   GeometryHierarchyMap<MaterialInteractionAssignment::LocalVeto> localVetos(
0191       localVetoVector);
0192   MaterialInteractionAssignment::Options options;
0193   options.localVetos = localVetos;
0194
0195   // Assign the material interactions to the surface hits
0196   auto [assigned, unassigned, surfacesLeft] =
0197       MaterialInteractionAssignment::assign(tContext, materialInteractions,
0198                                             intersectedSurfaces, options);
0199
0200   // Check that the material interaction was assigned
0201   BOOST_CHECK_EQUAL(assigned.size(), 34u);
0202   BOOST_CHECK_EQUAL(unassigned.size(), 15u);
0203   BOOST_CHECK_EQUAL(surfacesLeft.size(), 1u);
0204 }
0205
0206 BOOST_AUTO_TEST_CASE(AssignToClosest_withReassignment) {
0207   // Create a vector of surfaces
0208   std::vector<std::shared_ptr<Surface>> surfaces = {
0209       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0),
0210       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30.0, 100.0),
0211       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 50.0,
0212                                            100.0)};
0213
0214   for (auto [is, surface] : enumerate(surfaces)) {
0215     surface->assignGeometryId(GeometryIdentifier().withSensitive(is + 1));
0216   }
0217
0218   std::vector<IAssignmentFinder::SurfaceAssignment> intersectedSurfaces = {
0219       {surfaces[0].get(), {20., 0., 0.}, {1., 0., 0.}},
0220       {surfaces[1].get(), {30., 0., 0.}, {1., 0., 0.}},
0221       {surfaces[2].get(), {50., 0., 0.}, {1., 0., 0.}}};
0222
0223   // Create a material
0224   Material material = Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0);
0225
0226   std::vector<MaterialInteraction> materialInteractions;
0227   materialInteractions.reserve(50);
0228   // Generate some material interactions
0229   for (unsigned int i = 1; i < 50; ++i) {
0230     MaterialInteraction materialInteraction;
0231     materialInteraction.materialSlab = MaterialSlab(material, 0.1);
0232     materialInteraction.position = Vector3{i * 1.0, 0, 0};
0233     materialInteraction.direction = Vector3{1.0, 0., 0.};
0234     materialInteractions.push_back(materialInteraction);
0235   }
0236
0237   // Veto in a specific one
0238   struct ReAssignToNeighbor {
0239     void operator()(
0240         MaterialInteraction& m,
0241         const IAssignmentFinder::SurfaceAssignment& /*suggestedAssignment*/,
0242         const IAssignmentFinder::SurfaceAssignment& n) const {
0243       auto [surface, position, direction] = n;
0244       m.surface = surface;
0245       m.position = position;
0246       m.direction = direction;
0247       m.intersectionID = surface->geometryId();
0248       return;
0249     }
0250   };
0251
0252   // We assign this to
0253   std::vector<std::pair<GeometryIdentifier,
0254                         MaterialInteractionAssignment::ReAssignment>>
0255       reassignmentVector = {
0256           {GeometryIdentifier().withSensitive(2), ReAssignToNeighbor{}}};
0257   GeometryHierarchyMap<MaterialInteractionAssignment::ReAssignment>
0258       reassignments(reassignmentVector);
0259   MaterialInteractionAssignment::Options options;
0260   options.reAssignments = reassignments;
0261
0262   // Assign the material interactions to the surface hits
0263   auto [assigned, unassigned, surfaceLeft] =
0264       MaterialInteractionAssignment::assign(tContext, materialInteractions,
0265                                             intersectedSurfaces, options);
0266
0267   // Check that the material interaction was assigned
0268   BOOST_CHECK_EQUAL(assigned.size(), 49u);
0269   BOOST_CHECK_EQUAL(unassigned.size(), 0u);
0270   BOOST_CHECK_EQUAL(surfaceLeft.size(), 1u);
0271
0272   // Check that the geoid with number 2 never shows up
0273   for (const auto& mi : assigned) {
0274     BOOST_CHECK_NE(mi.intersectionID, GeometryIdentifier().withSensitive(2));
0275   }
0276 }
0277
0278 BOOST_AUTO_TEST_CASE(AssignWithPathLength) {
0279   auto surface =
0280       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 20.0, 100.0);
0281   surface->assignGeometryId(GeometryIdentifier().withSensitive(1));
0282
0283   Vector3 position = {20., 10., 0.};
0284   Vector3 direction = position.normalized();
0285
0286   IAssignmentFinder::SurfaceAssignment surfaceHit{surface.get(), position,
0287                                                   direction};
0288
0289   Material material = Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0);
0290
0291   MaterialInteraction materialInteraction;
0292   materialInteraction.materialSlab = MaterialSlab(material, 0.1);
0293   materialInteraction.position = position + 0.5 * direction;
0294   materialInteraction.direction = direction;
0295
0296   MaterialInteractionAssignment::Options options;
0297
0298   auto [assigned, unassigned, surfaceLeft] =
0299       MaterialInteractionAssignment::assign(tContext, {materialInteraction},
0300                                             {surfaceHit}, options);
0301
0302   // Check that the material interaction was assigned
0303   BOOST_CHECK_EQUAL(assigned.size(), 1u);
0304   BOOST_CHECK_EQUAL(unassigned.size(), 0u);
0305   BOOST_CHECK_EQUAL(surfaceLeft.size(), 0u);
0306
0307   // Check that the path correction is set
0308   BOOST_CHECK_NE(assigned[0].pathCorrection, 0.);
0309 }
0310
0311 BOOST_AUTO_TEST_SUITE_END()
0312
0313 }  // namespace ActsTests