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