Core/Material/VolumeMaterialMapperTests.cpp

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