Core/Material/VolumeMaterialMapperTests.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/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/Geometry/TrackingGeometry.hpp"
0018 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
0019 #include "Acts/Geometry/TrackingVolume.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Material/AccumulatedVolumeMaterial.hpp"
0022 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
0023 #include "Acts/Material/IVolumeMaterial.hpp"
0024 #include "Acts/Material/Material.hpp"
0025 #include "Acts/Material/MaterialGridHelper.hpp"
0026 #include "Acts/Material/MaterialSlab.hpp"
0027 #include "Acts/Material/ProtoVolumeMaterial.hpp"
0028 #include "Acts/Material/VolumeMaterialMapper.hpp"
0029 #include "Acts/Propagator/ActorList.hpp"
0030 #include "Acts/Propagator/Navigator.hpp"
0031 #include "Acts/Propagator/Propagator.hpp"
0032 #include "Acts/Propagator/StandardAborters.hpp"
0033 #include "Acts/Propagator/StraightLineStepper.hpp"
0034 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0035 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
0036 #include "Acts/Utilities/BinUtility.hpp"
0037 #include "Acts/Utilities/BinningType.hpp"
0038 #include "Acts/Utilities/Logger.hpp"
0039 #include "Acts/Utilities/Result.hpp"
0040
0041 #include <functional>
0042 #include <map>
0043 #include <memory>
0044 #include <random>
0045 #include <string>
0046 #include <tuple>
0047 #include <utility>
0048 #include <vector>
0049
0050 namespace Acts {
0051
0052 /// @brief Collector of material and position along propagation
0053 struct MaterialCollector {
0054   struct this_result {
0055     std::vector<Material> matTrue;
0056     std::vector<Vector3> position;
0057   };
0058   using result_type = this_result;
0059
0060   template <typename propagator_state_t, typename stepper_t,
0061             typename navigator_t>
0062   void act(propagator_state_t& state, const stepper_t& stepper,
0063            const navigator_t& navigator, result_type& result,
0064            const Logger& /*logger*/) const {
0065     if (navigator.currentVolume(state.navigation) != nullptr) {
0066       auto position = stepper.position(state.stepping);
0067       result.matTrue.push_back(
0068           (navigator.currentVolume(state.navigation)->volumeMaterial() !=
0069            nullptr)
0070               ? navigator.currentVolume(state.navigation)
0071                     ->volumeMaterial()
0072                     ->material(position)
0073               : Material());
0074
0075       result.position.push_back(position);
0076     }
0077   }
0078 };
0079
0080 }  // namespace Acts
0081
0082 namespace Acts::Test {
0083
0084 /// Test the filling and conversion
0085 BOOST_AUTO_TEST_CASE(SurfaceMaterialMapper_tests) {
0086   using namespace Acts::UnitLiterals;
0087
0088   BinUtility bu1(4, 0_m, 1_m, open, AxisDirection::AxisX);
0089   bu1 += BinUtility(2, -0.5_m, 0.5_m, open, AxisDirection::AxisY);
0090   bu1 += BinUtility(2, -0.5_m, 0.5_m, open, AxisDirection::AxisZ);
0091
0092   BinUtility bu2(4, 1_m, 2_m, open, AxisDirection::AxisX);
0093   bu2 += BinUtility(2, -0.5_m, 0.5_m, open, AxisDirection::AxisY);
0094   bu2 += BinUtility(2, -0.5_m, 0.5_m, open, AxisDirection::AxisZ);
0095
0096   BinUtility bu3(4, 2_m, 3_m, open, AxisDirection::AxisX);
0097   bu3 += BinUtility(2, -0.5_m, 0.5_m, open, AxisDirection::AxisY);
0098   bu3 += BinUtility(2, -0.5_m, 0.5_m, open, AxisDirection::AxisZ);
0099
0100   // Build a vacuum volume
0101   CuboidVolumeBuilder::VolumeConfig vCfg1;
0102   vCfg1.position = Vector3(0.5_m, 0., 0.);
0103   vCfg1.length = Vector3(1_m, 1_m, 1_m);
0104   vCfg1.name = "Vacuum volume";
0105   vCfg1.volumeMaterial = std::make_shared<const ProtoVolumeMaterial>(bu1);
0106
0107   // Build a material volume
0108   CuboidVolumeBuilder::VolumeConfig vCfg2;
0109   vCfg2.position = Vector3(1.5_m, 0., 0.);
0110   vCfg2.length = Vector3(1_m, 1_m, 1_m);
0111   vCfg2.name = "First material volume";
0112   vCfg2.volumeMaterial = std::make_shared<const ProtoVolumeMaterial>(bu2);
0113
0114   // Build another material volume with different material
0115   CuboidVolumeBuilder::VolumeConfig vCfg3;
0116   vCfg3.position = Vector3(2.5_m, 0., 0.);
0117   vCfg3.length = Vector3(1_m, 1_m, 1_m);
0118   vCfg3.name = "Second material volume";
0119   vCfg3.volumeMaterial = std::make_shared<const ProtoVolumeMaterial>(bu3);
0120
0121   // Configure world
0122   CuboidVolumeBuilder::Config cfg;
0123   cfg.position = Vector3(1.5_m, 0., 0.);
0124   cfg.length = Vector3(3_m, 1_m, 1_m);
0125   cfg.volumeCfg = {vCfg1, vCfg2, vCfg3};
0126
0127   GeometryContext gc;
0128
0129   // Build a detector
0130   CuboidVolumeBuilder cvb(cfg);
0131   TrackingGeometryBuilder::Config tgbCfg;
0132   tgbCfg.trackingVolumeBuilders.push_back(
0133       [=](const auto& context, const auto& inner, const auto&) {
0134         return cvb.trackingVolume(context, inner, nullptr);
0135       });
0136   TrackingGeometryBuilder tgb(tgbCfg);
0137   std::shared_ptr<const TrackingGeometry> tGeometry = tgb.trackingGeometry(gc);
0138
0139   /// We need a Navigator, Stepper to build a Propagator
0140   Navigator navigator({tGeometry});
0141   StraightLineStepper stepper;
0142   VolumeMaterialMapper::StraightLinePropagator propagator(stepper,
0143                                                           std::move(navigator));
0144
0145   /// The config object
0146   Acts::VolumeMaterialMapper::Config vmmConfig;
0147   Acts::VolumeMaterialMapper vmMapper(
0148       vmmConfig, std::move(propagator),
0149       getDefaultLogger("VolumeMaterialMapper", Logging::VERBOSE));
0150
0151   /// Create some contexts
0152   GeometryContext gCtx;
0153   MagneticFieldContext mfCtx;
0154
0155   /// Now create the mapper state
0156   auto mState = vmMapper.createState(gCtx, mfCtx, *tGeometry);
0157
0158   /// Test if this is not null
0159   BOOST_CHECK_EQUAL(mState.materialBin.size(), 3u);
0160 }
0161
0162 /// @brief Test case for comparison between the mapped material and the
0163 /// associated material by propagation
0164 BOOST_AUTO_TEST_CASE(VolumeMaterialMapper_comparison_tests) {
0165   using namespace Acts::UnitLiterals;
0166
0167   // Build a vacuum volume
0168   CuboidVolumeBuilder::VolumeConfig vCfg1;
0169   vCfg1.position = Vector3(0.5_m, 0., 0.);
0170   vCfg1.length = Vector3(1_m, 1_m, 1_m);
0171   vCfg1.name = "Vacuum volume";
0172   vCfg1.volumeMaterial =
0173       std::make_shared<const HomogeneousVolumeMaterial>(Material());
0174
0175   // Build a material volume
0176   CuboidVolumeBuilder::VolumeConfig vCfg2;
0177   vCfg2.position = Vector3(1.5_m, 0., 0.);
0178   vCfg2.length = Vector3(1_m, 1_m, 1_m);
0179   vCfg2.name = "First material volume";
0180   vCfg2.volumeMaterial =
0181       std::make_shared<HomogeneousVolumeMaterial>(makeSilicon());
0182
0183   // Build another material volume with different material
0184   CuboidVolumeBuilder::VolumeConfig vCfg3;
0185   vCfg3.position = Vector3(2.5_m, 0., 0.);
0186   vCfg3.length = Vector3(1_m, 1_m, 1_m);
0187   vCfg3.name = "Second material volume";
0188   vCfg3.volumeMaterial =
0189       std::make_shared<const HomogeneousVolumeMaterial>(Material());
0190
0191   // Configure world
0192   CuboidVolumeBuilder::Config cfg;
0193   cfg.position = Vector3(1.5_m, 0., 0.);
0194   cfg.length = Vector3(3_m, 1_m, 1_m);
0195   cfg.volumeCfg = {vCfg1, vCfg2, vCfg3};
0196
0197   GeometryContext gc;
0198
0199   // Build a detector
0200   CuboidVolumeBuilder cvb(cfg);
0201   TrackingGeometryBuilder::Config tgbCfg;
0202   tgbCfg.trackingVolumeBuilders.push_back(
0203       [=](const auto& context, const auto& inner, const auto&) {
0204         return cvb.trackingVolume(context, inner, nullptr);
0205       });
0206   TrackingGeometryBuilder tgb(tgbCfg);
0207   std::unique_ptr<const TrackingGeometry> detector = tgb.trackingGeometry(gc);
0208
0209   // Set up the grid axes
0210   Acts::MaterialGridAxisData xAxis{0_m, 3_m, 7};
0211   Acts::MaterialGridAxisData yAxis{-0.5_m, 0.5_m, 7};
0212   Acts::MaterialGridAxisData zAxis{-0.5_m, 0.5_m, 7};
0213
0214   // Set up a random engine for sampling material
0215   std::random_device rd;
0216   std::mt19937 gen(42);
0217   std::uniform_real_distribution<double> disX(0., 3_m);
0218   std::uniform_real_distribution<double> disYZ(-0.5_m, 0.5_m);
0219
0220   // Sample the Material in the detector
0221   RecordedMaterialVolumePoint matRecord;
0222   for (unsigned int i = 0; i < 1e4; i++) {
0223     Vector3 pos(disX(gen), disYZ(gen), disYZ(gen));
0224     std::vector<Vector3> volPos;
0225     volPos.push_back(pos);
0226     Material tv =
0227         (detector->lowestTrackingVolume(gc, pos)->volumeMaterial() != nullptr)
0228             ? (detector->lowestTrackingVolume(gc, pos)->volumeMaterial())
0229                   ->material(pos)
0230             : Material();
0231     MaterialSlab matProp(tv, 1);
0232     matRecord.push_back(std::make_pair(matProp, volPos));
0233   }
0234
0235   // Build the material grid
0236   Grid3D Grid = createGrid(xAxis, yAxis, zAxis);
0237   std::function<Vector3(Vector3)> transfoGlobalToLocal =
0238       [](Vector3 pos) -> Vector3 { return {pos.x(), pos.y(), pos.z()}; };
0239
0240   // Walk over each property
0241   for (const auto& rm : matRecord) {
0242     // Walk over each point associated with the properties
0243     for (const auto& point : rm.second) {
0244       // Search for fitting grid point and accumulate
0245       Acts::Grid3D::index_t index =
0246           Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
0247       Grid.atLocalBins(index).accumulate(rm.first);
0248     }
0249   }
0250
0251   MaterialGrid3D matGrid = mapMaterialPoints(Grid);
0252
0253   // Construct a simple propagation through the detector
0254   StraightLineStepper sls;
0255   Navigator::Config navCfg;
0256   navCfg.trackingGeometry = std::move(detector);
0257   Navigator nav(navCfg);
0258   Propagator<StraightLineStepper, Navigator> prop(sls, nav);
0259
0260   // Set some start parameters
0261   Vector4 pos4(0., 0., 0., 42_ns);
0262   Vector3 dir(1., 0., 0.);
0263   CurvilinearTrackParameters sctp(pos4, dir, 1 / 1_GeV, std::nullopt,
0264                                   ParticleHypothesis::pion0());
0265
0266   MagneticFieldContext mc;
0267   // Launch propagation and gather result
0268   using PropagatorOptions = Propagator<StraightLineStepper, Navigator>::Options<
0269       ActorList<MaterialCollector, EndOfWorldReached>>;
0270   PropagatorOptions po(gc, mc);
0271   po.stepping.maxStepSize = 1._mm;
0272   po.maxSteps = 1e6;
0273
0274   const auto& result = prop.propagate(sctp, po).value();
0275   const MaterialCollector::this_result& stepResult =
0276       result.get<typename MaterialCollector::result_type>();
0277
0278   // Collect the material as given by the grid and test it
0279   std::vector<Material> matvector;
0280   double gridX0 = 0., gridL0 = 0., trueX0 = 0., trueL0 = 0.;
0281   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
0282     matvector.push_back(
0283         Acts::Material{matGrid.atPosition(stepResult.position[i])});
0284     gridX0 += 1 / matvector[i].X0();
0285     gridL0 += 1 / matvector[i].L0();
0286     trueX0 += 1 / stepResult.matTrue[i].X0();
0287     trueL0 += 1 / stepResult.matTrue[i].L0();
0288   }
0289   CHECK_CLOSE_REL(gridX0, trueX0, 1e-1);
0290   CHECK_CLOSE_REL(gridL0, trueL0, 1e-1);
0291 }
0292
0293 }  // namespace Acts::Test