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 "Acts/EventData/NeutralTrackParameters.hpp"
12 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
13 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
14 #include "Acts/Geometry/CylinderVolumeBuilder.hpp"
15 #include "Acts/Geometry/CylinderVolumeHelper.hpp"
16 #include "Acts/Geometry/GeometryContext.hpp"
17 #include "Acts/Geometry/TrackingGeometry.hpp"
18 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
19 #include "Acts/Geometry/TrackingVolume.hpp"
20 #include "Acts/Geometry/TrackingVolumeArrayCreator.hpp"
21 #include "Acts/MagneticField/MagneticFieldContext.hpp"
22 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
23 #include "Acts/Material/Material.hpp"
24 #include "Acts/Material/MaterialGridHelper.hpp"
25 #include "Acts/Material/MaterialMapUtils.hpp"
26 #include "Acts/Material/MaterialSlab.hpp"
27 #include "Acts/Material/ProtoVolumeMaterial.hpp"
28 #include "Acts/Material/VolumeMaterialMapper.hpp"
29 #include "Acts/Propagator/Navigator.hpp"
30 #include "Acts/Propagator/Propagator.hpp"
31 #include "Acts/Propagator/StandardAborters.hpp"
32 #include "Acts/Propagator/StraightLineStepper.hpp"
33 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
34 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
35 #include "Acts/Utilities/Definitions.hpp"
36 #include "Acts/Utilities/Helpers.hpp"
37 #include "Acts/Utilities/detail/Axis.hpp"
38 #include "Acts/Utilities/detail/Grid.hpp"
39 
40 #include <limits>
41 #include <random>
42 #include <vector>
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 =
169  std::make_shared<HomogeneousVolumeMaterial>(makeSilicon());
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  RecordedMaterialVolumePoint matRecord;
210  for (unsigned int i = 0; i < 1e4; i++) {
211  Vector3D pos(disX(gen), disYZ(gen), disYZ(gen));
212  std::vector<Vector3D> volPos;
213  volPos.push_back(pos);
214  Material tv =
215  (detector->lowestTrackingVolume(gc, pos)->volumeMaterial() != nullptr)
216  ? (detector->lowestTrackingVolume(gc, pos)->volumeMaterial())
217  ->material(pos)
218  : Material();
219  MaterialSlab matProp(tv, 1);
220  matRecord.push_back(std::make_pair(matProp, volPos));
221  }
222 
223  // Build the material grid
224  Grid3D Grid = createGrid(xAxis, yAxis, zAxis);
225  std::function<Vector3D(Vector3D)> transfoGlobalToLocal =
226  [](Vector3D pos) -> Vector3D {
227  return {pos.x(), pos.y(), pos.z()};
228  };
229  MaterialGrid3D matGrid =
230  mapMaterialPoints(Grid, matRecord, transfoGlobalToLocal);
231 
232  // Construct a simple propagation through the detector
233  StraightLineStepper sls;
234  Navigator nav(std::move(detector));
235  Propagator<StraightLineStepper, Navigator> prop(sls, nav);
236 
237  // Set some start parameters
238  Vector4D pos4(0., 0., 0., 42_ns);
239  Vector3D dir(1., 0., 0.);
240  NeutralCurvilinearTrackParameters sctp(pos4, dir, 1 / 1_GeV);
241 
242  MagneticFieldContext mc;
243  // Launch propagation and gather result
244  PropagatorOptions<ActionList<MaterialCollector>, AbortList<EndOfWorldReached>>
245  po(gc, mc, getDummyLogger());
246  po.maxStepSize = 1._mm;
247  po.maxSteps = 1e6;
248 
249  const auto& result = prop.propagate(sctp, po).value();
250  const MaterialCollector::this_result& stepResult =
251  result.get<typename MaterialCollector::result_type>();
252 
253  // Collect the material as given by the grid and test it
254  std::vector<Material> matvector;
255  double gridX0 = 0., gridL0 = 0., trueX0 = 0., trueL0 = 0.;
256  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
257  matvector.push_back(matGrid.atPosition(stepResult.position[i]));
258  gridX0 += 1 / matvector[i].X0();
259  gridL0 += 1 / matvector[i].L0();
260  trueX0 += 1 / stepResult.matTrue[i].X0();
261  trueL0 += 1 / stepResult.matTrue[i].L0();
262  }
263  CHECK_CLOSE_REL(gridX0, trueX0, 1e-1);
264  CHECK_CLOSE_REL(gridL0, trueL0, 1e-1);
265 }
266 
267 } // namespace Test
268 } // namespace Acts