ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SurfaceMaterialMapper.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SurfaceMaterialMapper.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-2020 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 
21 
23  const Config& cfg, StraightLinePropagator propagator,
24  std::unique_ptr<const Logger> slogger)
25  : m_cfg(cfg),
26  m_propagator(std::move(propagator)),
27  m_logger(std::move(slogger)) {}
28 
30  const GeometryContext& gctx, const MagneticFieldContext& mctx,
31  const TrackingGeometry& tGeometry) const {
32  // Parse the geometry and find all surfaces with material proxies
33  auto world = tGeometry.highestTrackingVolume();
34 
35  // The Surface material mapping state
36  State mState(gctx, mctx);
37  resolveMaterialSurfaces(mState, *world);
38 
39  ACTS_DEBUG(mState.accumulatedMaterial.size()
40  << " Surfaces with PROXIES collected ... ");
41  for (auto& smg : mState.accumulatedMaterial) {
42  ACTS_VERBOSE(" -> Surface in with id " << smg.first);
43  }
44  return mState;
45 }
46 
48  State& mState, const TrackingVolume& tVolume) const {
49  ACTS_VERBOSE("Checking volume '" << tVolume.volumeName()
50  << "' for material surfaces.")
51 
52  ACTS_VERBOSE("- boundary surfaces ...");
53  // Check the boundary surfaces
54  for (auto& bSurface : tVolume.boundarySurfaces()) {
55  checkAndInsert(mState, bSurface->surfaceRepresentation());
56  }
57 
58  ACTS_VERBOSE("- confined layers ...");
59  // Check the confined layers
60  if (tVolume.confinedLayers() != nullptr) {
61  for (auto& cLayer : tVolume.confinedLayers()->arrayObjects()) {
62  // Take only layers that are not navigation layers
63  if (cLayer->layerType() != navigation) {
64  // Check the representing surface
65  checkAndInsert(mState, cLayer->surfaceRepresentation());
66  // Get the approach surfaces if present
67  if (cLayer->approachDescriptor() != nullptr) {
68  for (auto& aSurface :
69  cLayer->approachDescriptor()->containedSurfaces()) {
70  if (aSurface != nullptr) {
71  checkAndInsert(mState, *aSurface);
72  }
73  }
74  }
75  // Get the sensitive surface is present
76  if (cLayer->surfaceArray() != nullptr) {
77  // Sensitive surface loop
78  for (auto& sSurface : cLayer->surfaceArray()->surfaces()) {
79  if (sSurface != nullptr) {
80  checkAndInsert(mState, *sSurface);
81  }
82  }
83  }
84  }
85  }
86  }
87  // Step down into the sub volume
88  if (tVolume.confinedVolumes()) {
89  for (auto& sVolume : tVolume.confinedVolumes()->arrayObjects()) {
90  // Recursive call
91  resolveMaterialSurfaces(mState, *sVolume);
92  }
93  }
94 }
95 
97  const Surface& surface) const {
98  auto surfaceMaterial = surface.surfaceMaterial();
99  // Check if the surface has a proxy
100  if (surfaceMaterial != nullptr) {
101  auto geoID = surface.geoID();
102  size_t volumeID = geoID.volume();
103  ACTS_DEBUG("Material surface found with volumeID " << volumeID);
104  ACTS_DEBUG(" - surfaceID is " << geoID);
105 
106  // We need a dynamic_cast to either a surface material proxy or
107  // proper surface material
108  auto psm = dynamic_cast<const ProtoSurfaceMaterial*>(surfaceMaterial);
109 
110  // Get the bin utility: try proxy material first
111  const BinUtility* bu = (psm != nullptr) ? (&psm->binUtility()) : nullptr;
112  if (bu != nullptr) {
113  // Screen output for Binned Surface material
114  ACTS_DEBUG(" - (proto) binning is " << *bu);
115  // Now update
116  BinUtility buAdjusted = adjustBinUtility(*bu, surface);
117  // Screen output for Binned Surface material
118  ACTS_DEBUG(" - adjusted binning is " << buAdjusted);
119  mState.accumulatedMaterial[geoID] =
120  AccumulatedSurfaceMaterial(buAdjusted);
121  return;
122  }
123 
124  // Second attempt: binned material
125  auto bmp = dynamic_cast<const BinnedSurfaceMaterial*>(surfaceMaterial);
126  bu = (bmp != nullptr) ? (&bmp->binUtility()) : nullptr;
127  // Creaete a binned type of material
128  if (bu != nullptr) {
129  // Screen output for Binned Surface material
130  ACTS_DEBUG(" - binning is " << *bu);
131  mState.accumulatedMaterial[geoID] = AccumulatedSurfaceMaterial(*bu);
132  } else {
133  // Create a homogeneous type of material
134  ACTS_DEBUG(" - this is homogeneous material.");
136  }
137  }
138 }
139 
141  // iterate over the map to call the total average
142  for (auto& accMaterial : mState.accumulatedMaterial) {
143  ACTS_DEBUG("Finalizing map for Surface " << accMaterial.first);
144  mState.surfaceMaterial[accMaterial.first] =
145  accMaterial.second.totalAverage();
146  }
147 }
148 
150  State& mState, RecordedMaterialTrack& mTrack) const {
151  // Neutral curvilinear parameters
152  NeutralCurvilinearParameters start(std::nullopt, mTrack.first.first,
153  mTrack.first.second, 0.);
154 
155  // Prepare Action list and abort list
157  using MaterialSurfaceCollector = SurfaceCollector<MaterialSurface>;
158  using MaterialVolumeCollector = VolumeCollector<MaterialVolume>;
159  using ActionList = ActionList<MaterialSurfaceCollector,
160  MaterialVolumeCollector, DebugOutput>;
162 
164  mState.magFieldContext);
165  options.debug = m_cfg.mapperDebugOutput;
166 
167  // Now collect the material layers by using the straight line propagator
168  const auto& result = m_propagator.propagate(start, options).value();
169  auto mcResult = result.get<MaterialSurfaceCollector::result_type>();
170  auto mvcResult = result.get<MaterialVolumeCollector::result_type>();
171  // Massive screen output
172  if (m_cfg.mapperDebugOutput) {
173  auto debugOutput = result.get<DebugOutput::result_type>();
174  ACTS_VERBOSE("Debug propagation output.");
175  ACTS_VERBOSE(debugOutput.debugString);
176  }
177 
178  auto mappingSurfaces = mcResult.collected;
179  auto mappingVolumes = mvcResult.collected;
180 
181  // Retrieve the recorded material from the recorded material track
182  auto& rMaterial = mTrack.second.materialInteractions;
183  std::map<GeometryID, unsigned int> assignedMaterial;
184  ACTS_VERBOSE("Retrieved " << rMaterial.size()
185  << " recorded material steps to map.")
186 
187  // These should be mapped onto the mapping surfaces found
188  ACTS_VERBOSE("Found " << mappingSurfaces.size()
189  << " mapping surfaces for this track.");
190  ACTS_VERBOSE("Mapping surfaces are :")
191  for (auto& mSurface : mappingSurfaces) {
192  ACTS_VERBOSE(" - Surface : " << mSurface.surface->geoID()
193  << " at position = (" << mSurface.position.x()
194  << ", " << mSurface.position.y() << ", "
195  << mSurface.position.z() << ")");
196  assignedMaterial[mSurface.surface->geoID()] = 0;
197  }
198 
199  // Run the mapping process, i.e. take the recorded material and map it
200  // onto the mapping surfaces:
201  // - material steps and surfaces are assumed to be ordered along the
202  // mapping ray
203  //- do not record the material inside a volume with material
204  auto rmIter = rMaterial.begin();
205  auto sfIter = mappingSurfaces.begin();
206  auto volIter = mappingVolumes.begin();
207  bool encounterVolume = false;
208 
209  // Use those to minimize the lookup
210  GeometryID lastID = GeometryID();
211  GeometryID currentID = GeometryID();
212  Vector3D currentPos(0., 0., 0);
213  double currentPathCorrection = 0.;
214  auto currentAccMaterial = mState.accumulatedMaterial.end();
215 
216  // To remember the bins of this event
217  using MapBin = std::pair<AccumulatedSurfaceMaterial*, std::array<size_t, 3>>;
218  std::multimap<AccumulatedSurfaceMaterial*, std::array<size_t, 3>>
219  touchedMapBins;
220 
221  // Assign the recorded ones, break if you hit an end
222  while (rmIter != rMaterial.end() && sfIter != mappingSurfaces.end()) {
223  if (volIter != mappingVolumes.end() && encounterVolume == true &&
224  !volIter->volume->inside(rmIter->position)) {
225  encounterVolume = false;
226  // Switch to next material volume
227  ++volIter;
228  }
230  if (volIter != mappingVolumes.end() &&
231  volIter->volume->inside(rmIter->position)) {
232  encounterVolume = true;
233  ++rmIter;
234  continue;
235  }
236  if (sfIter != mappingSurfaces.end() - 1 &&
237  (rmIter->position - sfIter->position).norm() >
238  (rmIter->position - (sfIter + 1)->position).norm()) {
239  // Switch to next assignment surface
240  ++sfIter;
241  }
242  // get the current Surface ID
243  currentID = sfIter->surface->geoID();
244  // We have work to do: the assignemnt surface has changed
245  if (not(currentID == lastID)) {
246  // Let's (re-)assess the information
247  lastID = currentID;
248  currentPos = (sfIter)->position;
249  currentPathCorrection = sfIter->surface->pathCorrection(
250  mState.geoContext, currentPos, sfIter->direction);
251  currentAccMaterial = mState.accumulatedMaterial.find(currentID);
252  }
253  // Now assign the material for the accumulation process
254  auto tBin = currentAccMaterial->second.accumulate(
255  currentPos, rmIter->materialProperties, currentPathCorrection);
256  touchedMapBins.insert(MapBin(&(currentAccMaterial->second), tBin));
257  ++assignedMaterial[currentID];
258  // Update the material interaction with the associated surface and direction
259  rmIter->direction = mTrack.first.second.normalized();
260  rmIter->surface = sfIter->surface;
261  // Switch to next material
262  ++rmIter;
263  }
264 
265  ACTS_VERBOSE("Surfaces have following number of assigned hits :")
266  for (auto& [key, value] : assignedMaterial) {
267  ACTS_VERBOSE(" + Surface : " << key << " has " << value << " hits.");
268  }
269 
270  // After mapping this track, average the touched bins
271  for (auto tmapBin : touchedMapBins) {
272  std::vector<std::array<size_t, 3>> trackBins = {tmapBin.second};
273  tmapBin.first->trackAverage(trackBins);
274  }
275 
276  // After mapping this track, average the untouched but intersected bins
277  if (m_cfg.emptyBinCorrection) {
278  // Use the assignedMaterial map to account for empty hits, i.e.
279  // the material surface has been intersected by the mapping ray
280  // but no material step was assigned to this surface
281  for (auto& mSurface : mappingSurfaces) {
282  auto mgID = mSurface.surface->geoID();
283  // Count an empty hit only if the surface does not appear in the
284  // list of assigned surfaces
285  if (assignedMaterial[mgID] == 0) {
286  auto missedMaterial = mState.accumulatedMaterial.find(mgID);
287  missedMaterial->second.trackAverage(mSurface.position, true);
288  }
289  }
290  }
291 }