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