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