src/Material/VolumeMaterialMapper.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 "Acts/Material/VolumeMaterialMapper.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/EventData/ParticleHypothesis.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Geometry/ApproachDescriptor.hpp"
0016 #include "Acts/Geometry/BoundarySurfaceT.hpp"
0017 #include "Acts/Geometry/Layer.hpp"
0018 #include "Acts/Geometry/TrackingGeometry.hpp"
0019 #include "Acts/Material/AccumulatedVolumeMaterial.hpp"
0020 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
0021 #include "Acts/Material/IVolumeMaterial.hpp"
0022 #include "Acts/Material/InterpolatedMaterialMap.hpp"
0023 #include "Acts/Material/Material.hpp"
0024 #include "Acts/Material/MaterialGridHelper.hpp"
0025 #include "Acts/Material/MaterialInteraction.hpp"
0026 #include "Acts/Material/ProtoVolumeMaterial.hpp"
0027 #include "Acts/Propagator/ActorList.hpp"
0028 #include "Acts/Propagator/SurfaceCollector.hpp"
0029 #include "Acts/Propagator/VolumeCollector.hpp"
0030 #include "Acts/Surfaces/SurfaceArray.hpp"
0031 #include "Acts/Utilities/BinAdjustmentVolume.hpp"
0032 #include "Acts/Utilities/BinnedArray.hpp"
0033 #include "Acts/Utilities/Grid.hpp"
0034 #include "Acts/Utilities/Result.hpp"
0035
0036 #include <cmath>
0037 #include <cstddef>
0038 #include <iosfwd>
0039 #include <ostream>
0040 #include <stdexcept>
0041 #include <string>
0042 #include <tuple>
0043 #include <type_traits>
0044 #include <vector>
0045
0046 namespace Acts {
0047 struct EndOfWorldReached;
0048 }  // namespace Acts
0049
0050 Acts::VolumeMaterialMapper::VolumeMaterialMapper(
0051     const Config& cfg, StraightLinePropagator propagator,
0052     std::unique_ptr<const Logger> slogger)
0053     : m_cfg(cfg),
0054       m_propagator(std::move(propagator)),
0055       m_logger(std::move(slogger)) {}
0056
0057 Acts::VolumeMaterialMapper::State Acts::VolumeMaterialMapper::createState(
0058     const GeometryContext& gctx, const MagneticFieldContext& mctx,
0059     const TrackingGeometry& tGeometry) const {
0060   // Parse the geometry and find all surfaces with material proxies
0061   auto world = tGeometry.highestTrackingVolume();
0062
0063   // The Surface material mapping state
0064   State mState(gctx, mctx);
0065   resolveMaterialVolume(mState, *world);
0066   collectMaterialSurfaces(mState, *world);
0067   return mState;
0068 }
0069
0070 void Acts::VolumeMaterialMapper::resolveMaterialVolume(
0071     State& mState, const TrackingVolume& tVolume) const {
0072   ACTS_VERBOSE("Checking volume '" << tVolume.volumeName()
0073                                    << "' for material surfaces.");
0074
0075   ACTS_VERBOSE("- Insert Volume ...");
0076   checkAndInsert(mState, tVolume);
0077
0078   // Step down into the sub volume
0079   if (tVolume.confinedVolumes()) {
0080     ACTS_VERBOSE("- Check children volume ...");
0081     for (auto& sVolume : tVolume.confinedVolumes()->arrayObjects()) {
0082       // Recursive call
0083       resolveMaterialVolume(mState, *sVolume);
0084     }
0085   }
0086   if (!tVolume.denseVolumes().empty()) {
0087     for (auto& sVolume : tVolume.denseVolumes()) {
0088       // Recursive call
0089       resolveMaterialVolume(mState, *sVolume);
0090     }
0091   }
0092 }
0093
0094 void Acts::VolumeMaterialMapper::checkAndInsert(
0095     State& mState, const TrackingVolume& volume) const {
0096   auto volumeMaterial = volume.volumeMaterial();
0097   // Check if the volume has a proxy
0098   if (volumeMaterial != nullptr) {
0099     auto geoID = volume.geometryId();
0100     std::size_t volumeID = geoID.volume();
0101     ACTS_DEBUG("Material volume found with volumeID " << volumeID);
0102     ACTS_DEBUG("       - ID is " << geoID);
0103
0104     // We need a dynamic_cast to either a volume material proxy or
0105     // proper surface material
0106     auto psm = dynamic_cast<const ProtoVolumeMaterial*>(volumeMaterial);
0107     // Get the bin utility: try proxy material first
0108     const BinUtility* bu = (psm != nullptr) ? (&psm->binUtility()) : nullptr;
0109     if (bu != nullptr) {
0110       // Screen output for Binned Surface material
0111       ACTS_DEBUG("       - (proto) binning is " << *bu);
0112       // Now update
0113       BinUtility buAdjusted = adjustBinUtility(*bu, volume);
0114       // Screen output for Binned Surface material
0115       ACTS_DEBUG("       - adjusted binning is " << buAdjusted);
0116       mState.materialBin[geoID] = buAdjusted;
0117       if (bu->dimensions() == 0) {
0118         ACTS_DEBUG("Binning of dimension 0 create AccumulatedVolumeMaterial");
0119         Acts::AccumulatedVolumeMaterial homogeneousAccumulation;
0120         mState.homogeneousGrid[geoID] = homogeneousAccumulation;
0121       } else if (bu->dimensions() == 2) {
0122         ACTS_DEBUG("Binning of dimension 2 create 2D Grid");
0123         std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal;
0124         Acts::Grid2D Grid = createGrid2D(buAdjusted, transfoGlobalToLocal);
0125         mState.grid2D.insert(std::make_pair(geoID, Grid));
0126         mState.transform2D.insert(std::make_pair(geoID, transfoGlobalToLocal));
0127       } else if (bu->dimensions() == 3) {
0128         ACTS_DEBUG("Binning of dimension 3 create 3D Grid");
0129         std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal;
0130         Acts::Grid3D Grid = createGrid3D(buAdjusted, transfoGlobalToLocal);
0131         mState.grid3D.insert(std::make_pair(geoID, Grid));
0132         mState.transform3D.insert(std::make_pair(geoID, transfoGlobalToLocal));
0133       } else {
0134         throw std::invalid_argument(
0135             "Incorrect bin dimension, only 0, 2 and 3 are accepted");
0136       }
0137       return;
0138     }
0139     // Second attempt: 2D binned material
0140     auto bmp2 = dynamic_cast<
0141         const InterpolatedMaterialMap<MaterialMapper<Acts::MaterialGrid2D>>*>(
0142         volumeMaterial);
0143     bu = (bmp2 != nullptr) ? (&bmp2->binUtility()) : nullptr;
0144     if (bu != nullptr) {
0145       // Screen output for Binned Surface material
0146       ACTS_DEBUG("       - (2D grid) binning is " << *bu);
0147       mState.materialBin[geoID] = *bu;
0148       std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal;
0149       Acts::Grid2D Grid = createGrid2D(*bu, transfoGlobalToLocal);
0150       mState.grid2D.insert(std::make_pair(geoID, Grid));
0151       mState.transform2D.insert(std::make_pair(geoID, transfoGlobalToLocal));
0152       return;
0153     }
0154     // Third attempt: 3D binned material
0155     auto bmp3 = dynamic_cast<
0156         const InterpolatedMaterialMap<MaterialMapper<Acts::MaterialGrid3D>>*>(
0157         volumeMaterial);
0158     bu = (bmp3 != nullptr) ? (&bmp3->binUtility()) : nullptr;
0159     if (bu != nullptr) {
0160       // Screen output for Binned Surface material
0161       ACTS_DEBUG("       - (3D grid) binning is " << *bu);
0162       mState.materialBin[geoID] = *bu;
0163       std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal;
0164       Acts::Grid3D Grid = createGrid3D(*bu, transfoGlobalToLocal);
0165       mState.grid3D.insert(std::make_pair(geoID, Grid));
0166       mState.transform3D.insert(std::make_pair(geoID, transfoGlobalToLocal));
0167       return;
0168     } else {
0169       // Create a homogeneous type of material
0170       ACTS_DEBUG("       - this is homogeneous material.");
0171       BinUtility buHomogeneous;
0172       mState.materialBin[geoID] = buHomogeneous;
0173       Acts::AccumulatedVolumeMaterial homogeneousAccumulation;
0174       mState.homogeneousGrid[geoID] = homogeneousAccumulation;
0175       return;
0176     }
0177   }
0178 }
0179
0180 void Acts::VolumeMaterialMapper::collectMaterialSurfaces(
0181     State& mState, const TrackingVolume& tVolume) const {
0182   ACTS_VERBOSE("Checking volume '" << tVolume.volumeName()
0183                                    << "' for material surfaces.");
0184
0185   ACTS_VERBOSE("- boundary surfaces ...");
0186   // Check the boundary surfaces
0187   for (auto& bSurface : tVolume.boundarySurfaces()) {
0188     if (bSurface->surfaceRepresentation().surfaceMaterial() != nullptr) {
0189       mState.surfaceMaterial[bSurface->surfaceRepresentation().geometryId()] =
0190           bSurface->surfaceRepresentation().surfaceMaterialSharedPtr();
0191     }
0192   }
0193
0194   ACTS_VERBOSE("- confined layers ...");
0195   // Check the confined layers
0196   if (tVolume.confinedLayers() != nullptr) {
0197     for (auto& cLayer : tVolume.confinedLayers()->arrayObjects()) {
0198       // Take only layers that are not navigation layers
0199       if (cLayer->layerType() != navigation) {
0200         // Check the representing surface
0201         if (cLayer->surfaceRepresentation().surfaceMaterial() != nullptr) {
0202           mState.surfaceMaterial[cLayer->surfaceRepresentation().geometryId()] =
0203               cLayer->surfaceRepresentation().surfaceMaterialSharedPtr();
0204         }
0205         // Get the approach surfaces if present
0206         if (cLayer->approachDescriptor() != nullptr) {
0207           for (auto& aSurface :
0208                cLayer->approachDescriptor()->containedSurfaces()) {
0209             if (aSurface != nullptr) {
0210               if (aSurface->surfaceMaterial() != nullptr) {
0211                 mState.surfaceMaterial[aSurface->geometryId()] =
0212                     aSurface->surfaceMaterialSharedPtr();
0213               }
0214             }
0215           }
0216         }
0217         // Get the sensitive surface is present
0218         if (cLayer->surfaceArray() != nullptr) {
0219           // Sensitive surface loop
0220           for (auto& sSurface : cLayer->surfaceArray()->surfaces()) {
0221             if (sSurface != nullptr) {
0222               if (sSurface->surfaceMaterial() != nullptr) {
0223                 mState.surfaceMaterial[sSurface->geometryId()] =
0224                     sSurface->surfaceMaterialSharedPtr();
0225               }
0226             }
0227           }
0228         }
0229       }
0230     }
0231   }
0232   // Step down into the sub volume
0233   if (tVolume.confinedVolumes()) {
0234     for (auto& sVolume : tVolume.confinedVolumes()->arrayObjects()) {
0235       // Recursive call
0236       collectMaterialSurfaces(mState, *sVolume);
0237     }
0238   }
0239 }
0240
0241 void Acts::VolumeMaterialMapper::createExtraHits(
0242     State& mState,
0243     std::pair<const GeometryIdentifier, BinUtility>& currentBinning,
0244     Acts::MaterialSlab properties, const Vector3& position,
0245     Vector3 direction) const {
0246   if (currentBinning.second.dimensions() == 0) {
0247     // Writing homogeneous material for the current volumes no need to create
0248     // extra hits. We directly accumulate the material
0249     mState.homogeneousGrid[currentBinning.first].accumulate(properties);
0250     return;
0251   }
0252
0253   // Computing the extra hits properties based on the mappingStep length
0254   int volumeStep =
0255       static_cast<int>(std::floor(properties.thickness() / m_cfg.mappingStep));
0256   float remainder = properties.thickness() - m_cfg.mappingStep * volumeStep;
0257   properties.scaleThickness(m_cfg.mappingStep / properties.thickness());
0258   direction = direction * (m_cfg.mappingStep / direction.norm());
0259
0260   for (int extraStep = 0; extraStep < volumeStep; extraStep++) {
0261     Vector3 extraPosition = position + extraStep * direction;
0262     // Create additional extrapolated points for the grid mapping
0263
0264     if (currentBinning.second.dimensions() == 2) {
0265       auto grid = mState.grid2D.find(currentBinning.first);
0266       if (grid != mState.grid2D.end()) {
0267         // Find which grid bin the material fall into then accumulate
0268         Acts::Grid2D::index_t index = grid->second.localBinsFromLowerLeftEdge(
0269             mState.transform2D[currentBinning.first](extraPosition));
0270         grid->second.atLocalBins(index).accumulate(properties);
0271       } else {
0272         throw std::domain_error("No grid 2D was found");
0273       }
0274     } else if (currentBinning.second.dimensions() == 3) {
0275       auto grid = mState.grid3D.find(currentBinning.first);
0276       if (grid != mState.grid3D.end()) {
0277         // Find which grid bin the material fall into then accumulate
0278         Acts::Grid3D::index_t index = grid->second.localBinsFromLowerLeftEdge(
0279             mState.transform3D[currentBinning.first](extraPosition));
0280         grid->second.atLocalBins(index).accumulate(properties);
0281       } else {
0282         throw std::domain_error("No grid 3D was found");
0283       }
0284     }
0285   }
0286
0287   if (remainder > 0) {
0288     // We need to have an additional extra hit with the remainder length. Adjust
0289     // the thickness of the last extrapolated step
0290     properties.scaleThickness(remainder / properties.thickness());
0291     Vector3 extraPosition = position + volumeStep * direction;
0292     if (currentBinning.second.dimensions() == 2) {
0293       auto grid = mState.grid2D.find(currentBinning.first);
0294       if (grid != mState.grid2D.end()) {
0295         // Find which grid bin the material fall into then accumulate
0296         Acts::Grid2D::index_t index = grid->second.localBinsFromLowerLeftEdge(
0297             mState.transform2D[currentBinning.first](extraPosition));
0298         grid->second.atLocalBins(index).accumulate(properties);
0299       } else {
0300         throw std::domain_error("No grid 2D was found");
0301       }
0302     } else if (currentBinning.second.dimensions() == 3) {
0303       auto grid = mState.grid3D.find(currentBinning.first);
0304       if (grid != mState.grid3D.end()) {
0305         // Find which grid bin the material fall into then accumulate
0306         Acts::Grid3D::index_t index = grid->second.localBinsFromLowerLeftEdge(
0307             mState.transform3D[currentBinning.first](extraPosition));
0308         grid->second.atLocalBins(index).accumulate(properties);
0309       } else {
0310         throw std::domain_error("No grid 3D was found");
0311       }
0312     }
0313   }
0314 }
0315
0316 void Acts::VolumeMaterialMapper::finalizeMaps(State& mState) const {
0317   // iterate over the volumes
0318   for (auto& matBin : mState.materialBin) {
0319     ACTS_DEBUG("Create the material for volume  " << matBin.first);
0320     if (matBin.second.dimensions() == 0) {
0321       // Average the homogeneous volume material then store it
0322       ACTS_DEBUG("Homogeneous material volume");
0323       Acts::Material mat = mState.homogeneousGrid[matBin.first].average();
0324       mState.volumeMaterial[matBin.first] =
0325           std::make_unique<HomogeneousVolumeMaterial>(mat);
0326     } else if (matBin.second.dimensions() == 2) {
0327       // Average the material in the 2D grid then create an
0328       // InterpolatedMaterialMap
0329       ACTS_DEBUG("Grid material volume");
0330       auto grid = mState.grid2D.find(matBin.first);
0331       if (grid != mState.grid2D.end()) {
0332         Acts::MaterialGrid2D matGrid = mapMaterialPoints(grid->second);
0333         MaterialMapper<Acts::MaterialGrid2D> matMap(
0334             mState.transform2D[matBin.first], matGrid);
0335         mState.volumeMaterial[matBin.first] = std::make_unique<
0336             InterpolatedMaterialMap<MaterialMapper<Acts::MaterialGrid2D>>>(
0337             std::move(matMap), matBin.second);
0338       } else {
0339         throw std::domain_error("No grid 2D was found");
0340       }
0341     } else if (matBin.second.dimensions() == 3) {
0342       // Average the material in the 3D grid then create an
0343       // InterpolatedMaterialMap
0344       ACTS_DEBUG("Grid material volume");
0345       auto grid = mState.grid3D.find(matBin.first);
0346       if (grid != mState.grid3D.end()) {
0347         Acts::MaterialGrid3D matGrid = mapMaterialPoints(grid->second);
0348         MaterialMapper<Acts::MaterialGrid3D> matMap(
0349             mState.transform3D[matBin.first], matGrid);
0350         mState.volumeMaterial[matBin.first] = std::make_unique<
0351             InterpolatedMaterialMap<MaterialMapper<Acts::MaterialGrid3D>>>(
0352             std::move(matMap), matBin.second);
0353       } else {
0354         throw std::domain_error("No grid 3D was found");
0355       }
0356     } else {
0357       throw std::invalid_argument(
0358           "Incorrect bin dimension, only 0, 2 and 3 are accepted");
0359     }
0360   }
0361 }
0362
0363 void Acts::VolumeMaterialMapper::mapMaterialTrack(
0364     State& mState, RecordedMaterialTrack& mTrack) const {
0365   using VectorHelpers::makeVector4;
0366
0367   // Neutral curvilinear parameters
0368   NeutralCurvilinearTrackParameters start(
0369       makeVector4(mTrack.first.first, 0), mTrack.first.second,
0370       1 / mTrack.first.second.norm(), std::nullopt,
0371       NeutralParticleHypothesis::geantino());
0372
0373   // Prepare Action list and abort list
0374   using BoundSurfaceCollector = SurfaceCollector<BoundSurfaceSelector>;
0375   using MaterialVolumeCollector = VolumeCollector<MaterialVolumeSelector>;
0376   using ActionList = ActorList<BoundSurfaceCollector, MaterialVolumeCollector,
0377                                EndOfWorldReached>;
0378
0379   StraightLinePropagator::Options<ActionList> options(mState.geoContext,
0380                                                       mState.magFieldContext);
0381
0382   // Now collect the material volume by using the straight line propagator
0383   const auto& result = m_propagator.propagate(start, options).value();
0384   auto mcResult = result.get<BoundSurfaceCollector::result_type>();
0385   auto mvcResult = result.get<MaterialVolumeCollector::result_type>();
0386
0387   auto mappingSurfaces = mcResult.collected;
0388   auto mappingVolumes = mvcResult.collected;
0389
0390   // Retrieve the recorded material from the recorded material track
0391   auto& rMaterial = mTrack.second.materialInteractions;
0392   ACTS_VERBOSE("Retrieved " << rMaterial.size()
0393                             << " recorded material steps to map.");
0394
0395   // These should be mapped onto the mapping surfaces found
0396   ACTS_VERBOSE("Found     " << mappingVolumes.size()
0397                             << " mapping volumes for this track.");
0398   ACTS_VERBOSE("Mapping volumes are :");
0399   for (auto& mVolumes : mappingVolumes) {
0400     ACTS_VERBOSE(" - Volume : " << mVolumes.volume->geometryId()
0401                                 << " at position = (" << mVolumes.position.x()
0402                                 << ", " << mVolumes.position.y() << ", "
0403                                 << mVolumes.position.z() << ")");
0404   }
0405   // Run the mapping process, i.e. take the recorded material and map it
0406   // onto the mapping volume:
0407   auto rmIter = rMaterial.begin();
0408   auto sfIter = mappingSurfaces.begin();
0409   auto volIter = mappingVolumes.begin();
0410
0411   // Use those to minimize the lookup
0412   GeometryIdentifier lastID = GeometryIdentifier();
0413   GeometryIdentifier currentID = GeometryIdentifier();
0414   auto currentBinning = mState.materialBin.end();
0415
0416   // store end position of the last material slab
0417   Acts::Vector3 lastPositionEnd = {0, 0, 0};
0418   Acts::Vector3 direction = {0, 0, 0};
0419
0420   if (volIter != mappingVolumes.end()) {
0421     lastPositionEnd = volIter->position;
0422   }
0423
0424   // loop over all the material hit in the track or until there no more volume
0425   // to map onto
0426   while (rmIter != rMaterial.end() && volIter != mappingVolumes.end()) {
0427     if (volIter != mappingVolumes.end() &&
0428         !volIter->volume->inside(rmIter->position)) {
0429       // Check if the material point is past the entry point to the current
0430       // volume (this prevents switching volume before the first volume has been
0431       // reached)
0432       double distVol = (volIter->position - mTrack.first.first).norm();
0433       double distMat = (rmIter->position - mTrack.first.first).norm();
0434       if (distMat - distVol > s_epsilon) {
0435         // Switch to next material volume
0436         ++volIter;
0437         continue;
0438       }
0439     }
0440     if (volIter != mappingVolumes.end() &&
0441         volIter->volume->inside(rmIter->position, s_epsilon)) {
0442       currentID = volIter->volume->geometryId();
0443       direction = rmIter->direction;
0444       if (!(currentID == lastID)) {
0445         // Let's (re-)assess the information
0446         lastID = currentID;
0447         lastPositionEnd = volIter->position;
0448         currentBinning = mState.materialBin.find(currentID);
0449       }
0450       // If the current volume has a ProtoVolumeMaterial
0451       // and the material hit has a non 0 thickness
0452       if (currentBinning != mState.materialBin.end() &&
0453           rmIter->materialSlab.thickness() > 0) {
0454         // check if there is vacuum between this material point and the last one
0455         float vacuumThickness = (rmIter->position - lastPositionEnd).norm();
0456         if (vacuumThickness > s_epsilon) {
0457           auto properties = Acts::MaterialSlab(vacuumThickness);
0458           // creat vacuum hits
0459           createExtraHits(mState, *currentBinning, properties, lastPositionEnd,
0460                           direction);
0461         }
0462         // determine the position of the last material slab using the track
0463         // direction
0464         direction =
0465             direction * (rmIter->materialSlab.thickness() / direction.norm());
0466         lastPositionEnd = rmIter->position + direction;
0467         // create additional material point
0468         createExtraHits(mState, *currentBinning, rmIter->materialSlab,
0469                         rmIter->position, direction);
0470       }
0471
0472       // check if we have reached the end of the volume or the last hit of the
0473       // track.
0474       if ((rmIter + 1) == rMaterial.end() ||
0475           !volIter->volume->inside((rmIter + 1)->position, s_epsilon)) {
0476         // find the boundary surface corresponding to the end of the volume
0477         while (sfIter != mappingSurfaces.end()) {
0478           if (sfIter->surface->geometryId().volume() == lastID.volume() ||
0479               ((volIter + 1) != mappingVolumes.end() &&
0480                sfIter->surface->geometryId().volume() ==
0481                    (volIter + 1)->volume->geometryId().volume())) {
0482             double distVol = (volIter->position - mTrack.first.first).norm();
0483             double distSur = (sfIter->position - mTrack.first.first).norm();
0484             if (distSur - distVol > s_epsilon) {
0485               float vacuumThickness =
0486                   (sfIter->position - lastPositionEnd).norm();
0487               // if the last material slab stop before the boundary surface
0488               // create vacuum hits
0489               if (vacuumThickness > s_epsilon) {
0490                 auto properties = Acts::MaterialSlab(vacuumThickness);
0491                 createExtraHits(mState, *currentBinning, properties,
0492                                 lastPositionEnd, direction);
0493                 lastPositionEnd = sfIter->position;
0494               }
0495               break;
0496             }
0497           }
0498           sfIter++;
0499         }
0500       }
0501       rmIter->volume = volIter->volume;
0502       rmIter->intersectionID = currentID;
0503       rmIter->intersection = rmIter->position;
0504     }
0505     ++rmIter;
0506   }
0507 }