ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VolumeMaterialMapperTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file VolumeMaterialMapperTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <boost/test/unit_test.hpp>
10 
11 #include <limits>
12 #include <random>
13 #include <vector>
14 
15 #include "Acts/EventData/SingleCurvilinearTrackParameters.hpp"
16 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
17 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
18 #include "Acts/Geometry/CylinderVolumeBuilder.hpp"
19 #include "Acts/Geometry/CylinderVolumeHelper.hpp"
20 #include "Acts/Geometry/GeometryContext.hpp"
21 #include "Acts/Geometry/TrackingGeometry.hpp"
22 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
23 #include "Acts/Geometry/TrackingVolume.hpp"
24 #include "Acts/Geometry/TrackingVolumeArrayCreator.hpp"
25 #include "Acts/MagneticField/MagneticFieldContext.hpp"
26 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
27 #include "Acts/Material/Material.hpp"
28 #include "Acts/Material/MaterialGridHelper.hpp"
29 #include "Acts/Material/MaterialMapUtils.hpp"
30 #include "Acts/Material/MaterialProperties.hpp"
31 #include "Acts/Material/ProtoVolumeMaterial.hpp"
32 #include "Acts/Material/VolumeMaterialMapper.hpp"
33 #include "Acts/Propagator/DebugOutputActor.hpp"
34 #include "Acts/Propagator/Navigator.hpp"
35 #include "Acts/Propagator/Propagator.hpp"
36 #include "Acts/Propagator/StandardAborters.hpp"
37 #include "Acts/Propagator/StraightLineStepper.hpp"
38 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
39 #include "Acts/Utilities/Definitions.hpp"
40 #include "Acts/Utilities/Helpers.hpp"
41 #include "Acts/Utilities/detail/Axis.hpp"
42 #include "Acts/Utilities/detail/Grid.hpp"
43 
44 namespace Acts {
45 
47 struct MaterialCollector {
48  struct this_result {
49  std::vector<Material> matTrue;
50  std::vector<Vector3D> position;
51  };
52  using result_type = this_result;
53 
54  template <typename propagator_state_t, typename stepper_t>
55  void operator()(propagator_state_t& state, const stepper_t& stepper,
56  result_type& result) const {
57  if (state.navigation.currentVolume != nullptr) {
58  auto position = stepper.position(state.stepping);
59  result.matTrue.push_back(
60  (state.navigation.currentVolume->volumeMaterial() != nullptr)
61  ? state.navigation.currentVolume->volumeMaterial()->material(
62  position)
63  : Material());
64 
65  result.position.push_back(position);
66  }
67  }
68 };
69 
70 namespace Test {
71 
73 BOOST_AUTO_TEST_CASE(SurfaceMaterialMapper_tests) {
74  using namespace Acts::UnitLiterals;
75 
76  BinUtility bu1(4, 0_m, 1_m, open, binX);
77  bu1 += BinUtility(2, -0.5_m, 0.5_m, open, binY);
78  bu1 += BinUtility(2, -0.5_m, 0.5_m, open, binZ);
79 
80  BinUtility bu2(4, 1_m, 2_m, open, binX);
81  bu2 += BinUtility(2, -0.5_m, 0.5_m, open, binY);
82  bu2 += BinUtility(2, -0.5_m, 0.5_m, open, binZ);
83 
84  BinUtility bu3(4, 2_m, 3_m, open, binX);
85  bu3 += BinUtility(2, -0.5_m, 0.5_m, open, binY);
86  bu3 += BinUtility(2, -0.5_m, 0.5_m, open, binZ);
87 
88  // Build a vacuum volume
89  CuboidVolumeBuilder::VolumeConfig vCfg1;
90  vCfg1.position = Vector3D(0.5_m, 0., 0.);
91  vCfg1.length = Vector3D(1_m, 1_m, 1_m);
92  vCfg1.name = "Vacuum volume";
93  vCfg1.volumeMaterial = std::make_shared<const ProtoVolumeMaterial>(bu1);
94 
95  // Build a material volume
96  CuboidVolumeBuilder::VolumeConfig vCfg2;
97  vCfg2.position = Vector3D(1.5_m, 0., 0.);
98  vCfg2.length = Vector3D(1_m, 1_m, 1_m);
99  vCfg2.name = "First material volume";
100  vCfg2.volumeMaterial = std::make_shared<const ProtoVolumeMaterial>(bu2);
101 
102  // Build another material volume with different material
103  CuboidVolumeBuilder::VolumeConfig vCfg3;
104  vCfg3.position = Vector3D(2.5_m, 0., 0.);
105  vCfg3.length = Vector3D(1_m, 1_m, 1_m);
106  vCfg3.name = "Second material volume";
107  vCfg3.volumeMaterial = std::make_shared<const ProtoVolumeMaterial>(bu3);
108 
109  // Configure world
110  CuboidVolumeBuilder::Config cfg;
111  cfg.position = Vector3D(1.5_m, 0., 0.);
112  cfg.length = Vector3D(3_m, 1_m, 1_m);
113  cfg.volumeCfg = {vCfg1, vCfg2, vCfg3};
114 
115  GeometryContext gc;
116 
117  // Build a detector
118  CuboidVolumeBuilder cvb(cfg);
119  TrackingGeometryBuilder::Config tgbCfg;
120  tgbCfg.trackingVolumeBuilders.push_back(
121  [=](const auto& context, const auto& inner, const auto&) {
122  return cvb.trackingVolume(context, inner, nullptr);
123  });
124  TrackingGeometryBuilder tgb(tgbCfg);
125  std::shared_ptr<const TrackingGeometry> tGeometry = tgb.trackingGeometry(gc);
126 
128  Navigator navigator(tGeometry);
129  StraightLineStepper stepper;
130  VolumeMaterialMapper::StraightLinePropagator propagator(std::move(stepper),
131  std::move(navigator));
132 
134  Acts::VolumeMaterialMapper::Config vmmConfig;
135  Acts::VolumeMaterialMapper vmMapper(
136  vmmConfig, std::move(propagator),
137  getDefaultLogger("VolumeMaterialMapper", Logging::VERBOSE));
138 
140  GeometryContext gCtx;
141  MagneticFieldContext mfCtx;
142 
144  auto mState = vmMapper.createState(gCtx, mfCtx, *tGeometry);
145 
147  BOOST_CHECK_EQUAL(mState.recordedMaterial.size(), 3u);
148 }
149 
152 BOOST_AUTO_TEST_CASE(VolumeMaterialMapper_comparison_tests) {
153  using namespace Acts::UnitLiterals;
154 
155  // Build a vacuum volume
156  CuboidVolumeBuilder::VolumeConfig vCfg1;
157  vCfg1.position = Vector3D(0.5_m, 0., 0.);
158  vCfg1.length = Vector3D(1_m, 1_m, 1_m);
159  vCfg1.name = "Vacuum volume";
160  vCfg1.volumeMaterial =
161  std::make_shared<const HomogeneousVolumeMaterial>(Material());
162 
163  // Build a material volume
164  CuboidVolumeBuilder::VolumeConfig vCfg2;
165  vCfg2.position = Vector3D(1.5_m, 0., 0.);
166  vCfg2.length = Vector3D(1_m, 1_m, 1_m);
167  vCfg2.name = "First material volume";
168  vCfg2.volumeMaterial = std::make_shared<const HomogeneousVolumeMaterial>(
169  Material(95.7, 465.2, 28.03, 14., 2.32e-3));
170 
171  // Build another material volume with different material
172  CuboidVolumeBuilder::VolumeConfig vCfg3;
173  vCfg3.position = Vector3D(2.5_m, 0., 0.);
174  vCfg3.length = Vector3D(1_m, 1_m, 1_m);
175  vCfg3.name = "Second material volume";
176  vCfg3.volumeMaterial =
177  std::make_shared<const HomogeneousVolumeMaterial>(Material());
178 
179  // Configure world
180  CuboidVolumeBuilder::Config cfg;
181  cfg.position = Vector3D(1.5_m, 0., 0.);
182  cfg.length = Vector3D(3_m, 1_m, 1_m);
183  cfg.volumeCfg = {vCfg1, vCfg2, vCfg3};
184 
185  GeometryContext gc;
186 
187  // Build a detector
188  CuboidVolumeBuilder cvb(cfg);
189  TrackingGeometryBuilder::Config tgbCfg;
190  tgbCfg.trackingVolumeBuilders.push_back(
191  [=](const auto& context, const auto& inner, const auto&) {
192  return cvb.trackingVolume(context, inner, nullptr);
193  });
194  TrackingGeometryBuilder tgb(tgbCfg);
195  std::unique_ptr<const TrackingGeometry> detector = tgb.trackingGeometry(gc);
196 
197  // Set up the grid axes
198  std::array<double, 3> xAxis{0_m, 3_m, 7};
199  std::array<double, 3> yAxis{-0.5_m, 0.5_m, 7};
200  std::array<double, 3> zAxis{-0.5_m, 0.5_m, 7};
201 
202  // Set up a random engine for sampling material
203  std::random_device rd;
204  std::mt19937 gen(42);
205  std::uniform_real_distribution<> disX(0., 3_m);
206  std::uniform_real_distribution<> disYZ(-0.5_m, 0.5_m);
207 
208  // Sample the Material in the detector
209  RecordedMaterialPoint matRecord;
210  for (unsigned int i = 0; i < 1e4; i++) {
211  Vector3D pos(disX(gen), disYZ(gen), disYZ(gen));
212  Material tv =
213  (detector->lowestTrackingVolume(gc, pos)->volumeMaterial() != nullptr)
214  ? (detector->lowestTrackingVolume(gc, pos)->volumeMaterial())
215  ->material(pos)
216  : Material();
217  MaterialProperties matProp(tv, 1);
218  matRecord.push_back(std::make_pair(matProp, pos));
219  }
220 
221  // Build the material grid
222  Grid3D Grid = createGrid(xAxis, yAxis, zAxis);
223  std::function<Vector3D(Vector3D)> transfoGlobalToLocal =
224  [](Vector3D pos) -> Vector3D {
225  return {pos.x(), pos.y(), pos.z()};
226  };
227  MaterialGrid3D matGrid =
228  mapMaterialPoints(Grid, matRecord, transfoGlobalToLocal);
229 
230  // Construct a simple propagation through the detector
231  StraightLineStepper sls;
232  Navigator nav(std::move(detector));
233  Propagator<StraightLineStepper, Navigator> prop(sls, nav);
234 
235  // Set some start parameters
236  Vector3D pos(0., 0., 0.);
237  Vector3D mom(1_GeV, 0., 0.);
238  SingleCurvilinearTrackParameters<NeutralPolicy> sctp(std::nullopt, pos, mom,
239  42_ns);
240 
241  MagneticFieldContext mc;
242  // Launch propagation and gather result
243  PropagatorOptions<ActionList<MaterialCollector, DebugOutputActor>,
244  AbortList<EndOfWorldReached>>
245  po(gc, mc);
246  po.maxStepSize = 1._mm;
247  po.maxSteps = 1e6;
248  po.debug = true;
249 
250  const auto& result = prop.propagate(sctp, po).value();
251  const MaterialCollector::this_result& stepResult =
252  result.get<typename MaterialCollector::result_type>();
253 
254  if (po.debug) {
255  auto screenOutput = result.get<DebugOutputActor::result_type>();
256  std::cout << screenOutput.debugString << std::endl;
257  }
258 
259  // Collect the material as given by the grid and test it
260  std::vector<Material> matvector;
261  double gridX0 = 0., gridL0 = 0., trueX0 = 0., trueL0 = 0.;
262  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
263  matvector.push_back(matGrid.atPosition(stepResult.position[i]));
264  gridX0 += 1 / matvector[i].X0();
265  gridL0 += 1 / matvector[i].L0();
266  trueX0 += 1 / stepResult.matTrue[i].X0();
267  trueL0 += 1 / stepResult.matTrue[i].L0();
268  }
269  CHECK_CLOSE_REL(gridX0, trueX0, 1e-1);
270  CHECK_CLOSE_REL(gridL0, trueL0, 1e-1);
271 }
272 
273 } // namespace Test
274 } // namespace Acts