TrackingVolume.cpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackingVolume.cpp

// This file is part of the Acts project.
//
// Copyright (C) 2016-2019 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "Acts/Geometry/TrackingVolume.hpp"

#include "Acts/Geometry/GlueVolumesDescriptor.hpp"
#include "Acts/Geometry/VolumeBounds.hpp"
#include "Acts/Material/ProtoVolumeMaterial.hpp"
#include "Acts/Propagator/Navigator.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/BinUtility.hpp"

#include <functional>
#include <utility>

Acts::TrackingVolume::TrackingVolume(
    const Transform3D& transform, VolumeBoundsPtr volbounds,
    const std::shared_ptr<const TrackingVolumeArray>& containedVolumeArray,
    const std::string& volumeName)
    : Volume(transform, std::move(volbounds)),
      m_volumeMaterial(nullptr),
      m_boundarySurfaces(),
      m_confinedLayers(nullptr),
      m_confinedVolumes(containedVolumeArray),
      m_name(volumeName) {
  createBoundarySurfaces();
  interlinkLayers();
}

// constructor for arguments
Acts::TrackingVolume::TrackingVolume(
    const Transform3D& transform, VolumeBoundsPtr volumeBounds,
    std::shared_ptr<const IVolumeMaterial> volumeMaterial,
    std::unique_ptr<const LayerArray> staticLayerArray,
    std::shared_ptr<const TrackingVolumeArray> containedVolumeArray,
    MutableTrackingVolumeVector denseVolumeVector,
    const std::string& volumeName)
    : Volume(transform, std::move(volumeBounds)),
      m_volumeMaterial(std::move(volumeMaterial)),
      m_confinedLayers(std::move(staticLayerArray)),
      m_confinedVolumes(std::move(containedVolumeArray)),
      m_confinedDenseVolumes({}),
      m_name(volumeName) {
  createBoundarySurfaces();
  interlinkLayers();
  connectDenseBoundarySurfaces(denseVolumeVector);
}

// constructor for arguments
Acts::TrackingVolume::TrackingVolume(
    const Transform3D& transform, VolumeBoundsPtr volbounds,
    std::vector<std::unique_ptr<Volume::BoundingBox>> boxStore,
    std::vector<std::unique_ptr<const Volume>> descendants,
    const Volume::BoundingBox* top,
    std::shared_ptr<const IVolumeMaterial> volumeMaterial,
    const std::string& volumeName)
    : Volume(transform, std::move(volbounds)),
      m_volumeMaterial(std::move(volumeMaterial)),
      m_name(volumeName),
      m_descendantVolumes(std::move(descendants)),
      m_bvhTop(top) {
  createBoundarySurfaces();
  // we take a copy of the unique box pointers, but we want to
  // store them as consts.
  for (auto& uptr : boxStore) {
    m_boundingBoxes.push_back(
        std::unique_ptr<Volume::BoundingBox>(uptr.release()));
  }
}

Acts::TrackingVolume::~TrackingVolume() {
  delete m_glueVolumeDescriptor;
}

const Acts::TrackingVolume* Acts::TrackingVolume::lowestTrackingVolume(
    const GeometryContext& /*gctx*/, const Vector3D& position,
    const double tol) const {
  // confined static volumes - highest hierarchy
  if (m_confinedVolumes) {
    return (m_confinedVolumes->object(position).get());
  }

  // search for dense volumes
  if (!m_confinedDenseVolumes.empty())
    for (auto& denseVolume : m_confinedDenseVolumes)
      if (denseVolume->inside(position, tol))
        return denseVolume.get();

  // there is no lower sub structure
  return this;
}

const Acts::TrackingVolumeBoundaries& Acts::TrackingVolume::boundarySurfaces()
    const {
  return (m_boundarySurfaces);
}

void Acts::TrackingVolume::connectDenseBoundarySurfaces(
    MutableTrackingVolumeVector& confinedDenseVolumes) {
  if (!confinedDenseVolumes.empty()) {
    NavigationDirection navDir;
    // Walk over each dense volume
    for (auto& confDenseVol : confinedDenseVolumes) {
      // Walk over each boundary surface of the volume
      auto& boundSur = confDenseVol->boundarySurfaces();
      for (unsigned int i = 0; i < boundSur.size(); i++) {
        // Skip empty entries since we do not know the shape of the dense volume
        // and therewith the used indices
        if (boundSur.at(i) == nullptr) {
          continue;
        }

        // Use mother volume as the opposite direction of the already used
        // direction
        auto mutableBs =
            std::const_pointer_cast<BoundarySurfaceT<TrackingVolume>>(
                boundSur.at(i));
        if (mutableBs->m_oppositeVolume != nullptr &&
            mutableBs->m_alongVolume == nullptr) {
          navDir = forward;
          mutableBs->attachVolume(this, navDir);
        } else {
          if (mutableBs->m_oppositeVolume == nullptr &&
              mutableBs->m_alongVolume != nullptr) {
            navDir = backward;
            mutableBs->attachVolume(this, navDir);
          }
        }

        // Update the boundary
        confDenseVol->updateBoundarySurface((BoundarySurfaceFace)i, mutableBs);
      }
      // Store the volume
      m_confinedDenseVolumes.push_back(std::move(confDenseVol));
    }
  }
}

void Acts::TrackingVolume::createBoundarySurfaces() {
  using Boundary = BoundarySurfaceT<TrackingVolume>;

  // Transform Surfaces To BoundarySurfaces
  auto orientedSurfaces = Volume::volumeBounds().orientedSurfaces(m_transform);

  m_boundarySurfaces.reserve(orientedSurfaces.size());
  for (auto& osf : orientedSurfaces) {
    TrackingVolume* opposite = nullptr;
    TrackingVolume* along = nullptr;
    if (osf.second == backward) {
      opposite = this;
    } else {
      along = this;
    }
    m_boundarySurfaces.push_back(std::make_shared<const Boundary>(
        std::move(osf.first), opposite, along));
  }
}

void Acts::TrackingVolume::glueTrackingVolume(const GeometryContext& gctx,
                                              BoundarySurfaceFace bsfMine,
                                              TrackingVolume* neighbor,
                                              BoundarySurfaceFace bsfNeighbor) {
  // Find the connection of the two tracking volumes: binR returns the center
  // except for cylindrical volumes
  Vector3D bPosition(binningPosition(gctx, binR));
  Vector3D distance =
      Vector3D(neighbor->binningPosition(gctx, binR) - bPosition);
  // glue to the face
  std::shared_ptr<const BoundarySurfaceT<TrackingVolume>> bSurfaceMine =
      boundarySurfaces().at(bsfMine);
  // @todo - complex glueing could be possible with actual intersection for the
  // normal vector
  Vector3D nvector =
      bSurfaceMine->surfaceRepresentation().normal(gctx, bPosition);
  // estimate the orientation
  NavigationDirection navDir =
      (nvector.dot(distance) > 0.) ? forward : backward;
  // The easy case :
  // - no glue volume descriptors on either side
  if ((m_glueVolumeDescriptor == nullptr) ||
      m_glueVolumeDescriptor->glueVolumes(bsfMine) == nullptr) {
    // the boundary orientation
    auto mutableBSurfaceMine =
        std::const_pointer_cast<BoundarySurfaceT<TrackingVolume>>(bSurfaceMine);
    mutableBSurfaceMine->attachVolume(neighbor, navDir);
    // Make sure you keep the boundary material if there
    const Surface& neighborSurface =
        neighbor->m_boundarySurfaces.at(bsfNeighbor)->surfaceRepresentation();
    auto neighborMaterial = neighborSurface.surfaceMaterialSharedPtr();
    const Surface& mySurface = bSurfaceMine->surfaceRepresentation();
    auto myMaterial = mySurface.surfaceMaterialSharedPtr();
    // Keep the neighbor material
    if (myMaterial == nullptr and neighborMaterial != nullptr) {
      Surface* myMutbableSurface = const_cast<Surface*>(&mySurface);
      myMutbableSurface->assignSurfaceMaterial(neighborMaterial);
    }
    // Now set it to the neighbor volume
    (neighbor->m_boundarySurfaces).at(bsfNeighbor) = bSurfaceMine;
  }
}

void Acts::TrackingVolume::glueTrackingVolumes(
    const GeometryContext& gctx, BoundarySurfaceFace bsfMine,
    const std::shared_ptr<TrackingVolumeArray>& neighbors,
    BoundarySurfaceFace bsfNeighbor) {
  // find the connection of the two tracking volumes : binR returns the center
  // except for cylindrical volumes
  std::shared_ptr<const TrackingVolume> nRefVolume =
      neighbors->arrayObjects().at(0);
  // get the distance
  Vector3D bPosition(binningPosition(gctx, binR));
  Vector3D distance =
      Vector3D(nRefVolume->binningPosition(gctx, binR) - bPosition);
  // take the normal at the binning positio
  std::shared_ptr<const BoundarySurfaceT<TrackingVolume>> bSurfaceMine =
      boundarySurfaces().at(bsfMine);
  // @todo - complex glueing could be possible with actual intersection for the
  // normal vector
  Vector3D nvector =
      bSurfaceMine->surfaceRepresentation().normal(gctx, bPosition);
  // estimate the orientation
  NavigationDirection navDir =
      (nvector.dot(distance) > 0.) ? forward : backward;
  // the easy case :
  // - no glue volume descriptors on either side
  if ((m_glueVolumeDescriptor == nullptr) ||
      !m_glueVolumeDescriptor->glueVolumes(bsfMine)) {
    // the boundary orientation
    auto mutableBSurfaceMine =
        std::const_pointer_cast<BoundarySurfaceT<TrackingVolume>>(bSurfaceMine);
    mutableBSurfaceMine->attachVolumeArray(neighbors, navDir);
    // now set it to the neighbor volumes - the optised way
    for (auto& nVolume : neighbors->arrayObjects()) {
      auto mutableNVolume = std::const_pointer_cast<TrackingVolume>(nVolume);
      (mutableNVolume->m_boundarySurfaces).at(bsfNeighbor) = bSurfaceMine;
    }
  }
}

void Acts::TrackingVolume::assignBoundaryMaterial(
    std::shared_ptr<const ISurfaceMaterial> surfaceMaterial,
    BoundarySurfaceFace bsFace) {
  auto bSurface = m_boundarySurfaces.at(bsFace);
  Surface* surface = const_cast<Surface*>(&bSurface->surfaceRepresentation());
  surface->assignSurfaceMaterial(std::move(surfaceMaterial));
}

void Acts::TrackingVolume::updateBoundarySurface(
    BoundarySurfaceFace bsf,
    std::shared_ptr<const BoundarySurfaceT<TrackingVolume>> bs,
    bool checkmaterial) {
  if (checkmaterial) {
    auto cMaterialPtr = m_boundarySurfaces.at(bsf)
                            ->surfaceRepresentation()
                            .surfaceMaterialSharedPtr();
    auto bsMaterial = bs->surfaceRepresentation().surfaceMaterial();
    if (cMaterialPtr != nullptr && bsMaterial == nullptr) {
      Surface* surface = const_cast<Surface*>(&bs->surfaceRepresentation());
      surface->assignSurfaceMaterial(cMaterialPtr);
    }
  }
  m_boundarySurfaces.at(bsf) = std::move(bs);
}

void Acts::TrackingVolume::registerGlueVolumeDescriptor(
    GlueVolumesDescriptor* gvd) {
  delete m_glueVolumeDescriptor;
  m_glueVolumeDescriptor = gvd;
}

Acts::GlueVolumesDescriptor& Acts::TrackingVolume::glueVolumesDescriptor() {
  if (m_glueVolumeDescriptor == nullptr) {
    m_glueVolumeDescriptor = new GlueVolumesDescriptor;
  }
  return (*m_glueVolumeDescriptor);
}

void Acts::TrackingVolume::synchronizeLayers(double envelope) const {
  // case a : Layers exist
  // msgstream << MSG::VERBOSE << "  -> synchronizing Layer dimensions of
  // TrackingVolume '" << volumeName() << "'." << endreq;

  if (m_confinedLayers) {
    // msgstream << MSG::VERBOSE << "  ---> working on " <<
    // m_confinedLayers->arrayObjects().size() << " (material+navigation)
    // layers." << endreq;
    for (auto& clayIter : m_confinedLayers->arrayObjects()) {
      if (clayIter) {
        // @todo implement syncrhonize layer
        //  if (clayIter->surfaceRepresentation().type() == Surface::Cylinder &&
        //  !(center().isApprox(clayIter->surfaceRepresentation().center())) )
        //      clayIter->resizeAndRepositionLayer(volumeBounds(),center(),envelope);
        //  else
        //      clayIter->resizeLayer(volumeBounds(),envelope);
      }  // else
      // msgstream << MSG::WARNING << "  ---> found 0 pointer to layer,
      // indicates problem." << endreq;
    }
  }

  // case b : container volume -> step down
  if (m_confinedVolumes) {
    // msgstream << MSG::VERBOSE << "  ---> no confined layers, working on " <<
    // m_confinedVolumes->arrayObjects().size() << " confined volumes." <<
    // endreq;
    for (auto& cVolumesIter : m_confinedVolumes->arrayObjects()) {
      cVolumesIter->synchronizeLayers(envelope);
    }
  }
}

void Acts::TrackingVolume::interlinkLayers() {
  if (m_confinedLayers) {
    auto& layers = m_confinedLayers->arrayObjects();

    // forward register the last one as the previous one
    //  first <- | -> second, first <- | -> second, first <- | -> second
    const Layer* lastLayer = nullptr;
    for (auto& layerPtr : layers) {
      // we'll need to mutate our confined layers to perform this operation
      Layer& mutableLayer = *(std::const_pointer_cast<Layer>(layerPtr));
      // register the layers
      mutableLayer.m_nextLayerUtility = m_confinedLayers->binUtility();
      mutableLayer.m_nextLayers.first = lastLayer;
      // register the volume
      mutableLayer.encloseTrackingVolume(*this);
      // remember the last layer
      lastLayer = &mutableLayer;
    }
    // backward loop
    lastLayer = nullptr;
    for (auto layerIter = layers.rbegin(); layerIter != layers.rend();
         ++layerIter) {
      // set the other next volume
      Layer& mutableLayer = *(std::const_pointer_cast<Layer>(*layerIter));
      mutableLayer.m_nextLayers.second = lastLayer;
      lastLayer = &mutableLayer;
    }
  }
}

void Acts::TrackingVolume::closeGeometry(
    const IMaterialDecorator* materialDecorator,
    std::map<std::string, const TrackingVolume*>& volumeMap, size_t& vol) {
  // insert the volume into the map
  volumeMap[volumeName()] = this;

  // we can construct the volume ID from this
  auto volumeID = GeometryIdentifier().setVolume(++vol);
  // assign the Volume ID to the volume itself
  auto thisVolume = const_cast<TrackingVolume*>(this);
  thisVolume->assignGeometryId(volumeID);

  // assign the material if you have a decorator
  if (materialDecorator != nullptr) {
    materialDecorator->decorate(*thisVolume);
  }
  if (thisVolume->volumeMaterial() == nullptr && thisVolume->motherVolume() &&
      thisVolume->motherVolume()->volumeMaterial() != nullptr) {
    auto protoMaterial = dynamic_cast<const Acts::ProtoVolumeMaterial*>(
        thisVolume->motherVolume()->volumeMaterial());
    if (protoMaterial == nullptr) {
      thisVolume->assignVolumeMaterial(
          thisVolume->motherVolume()->volumeMaterialSharedPtr());
    }
  }

  this->assignGeometryId(volumeID);
  // loop over the boundary surfaces
  GeometryIdentifier::Value iboundary = 0;
  // loop over the boundary surfaces
  for (auto& bSurfIter : boundarySurfaces()) {
    // get the intersection soltuion
    auto& bSurface = bSurfIter->surfaceRepresentation();
    // create the boundary surface id
    auto boundaryID = GeometryIdentifier(volumeID).setBoundary(++iboundary);
    // now assign to the boundary surface
    auto& mutableBSurface = *(const_cast<Surface*>(&bSurface));
    mutableBSurface.assignGeometryId(boundaryID);
    // assign the material if you have a decorator
    if (materialDecorator != nullptr) {
      materialDecorator->decorate(mutableBSurface);
    }
  }

  // A) this is NOT a container volume, volumeID is already incremented
  if (!m_confinedVolumes) {
    // loop over the confined layers
    if (m_confinedLayers) {
      GeometryIdentifier::Value ilayer = 0;
      // loop over the layers
      for (auto& layerPtr : m_confinedLayers->arrayObjects()) {
        // create the layer identification
        auto layerID = GeometryIdentifier(volumeID).setLayer(++ilayer);
        // now close the geometry
        auto mutableLayerPtr = std::const_pointer_cast<Layer>(layerPtr);
        mutableLayerPtr->closeGeometry(materialDecorator, layerID);
      }
    } else if (m_bvhTop != nullptr) {
      GeometryIdentifier::Value isurface = 0;
      for (const auto& descVol : m_descendantVolumes) {
        // Attempt to cast to AbstractVolume: only one we'll handle
        const AbstractVolume* avol =
            dynamic_cast<const AbstractVolume*>(descVol.get());
        if (avol != nullptr) {
          const auto& bndSrf = avol->boundarySurfaces();
          for (const auto& bnd : bndSrf) {
            const auto& srf = bnd->surfaceRepresentation();
            Surface* mutableSurfcePtr = const_cast<Surface*>(&srf);
            auto geoID = GeometryIdentifier(volumeID).setSensitive(++isurface);
            mutableSurfcePtr->assignGeometryId(geoID);
          }
        }
      }
    }
  } else {
    // B) this is a container volume, go through sub volume
    // do the loop
    for (auto& volumesIter : m_confinedVolumes->arrayObjects()) {
      auto mutableVolumesIter =
          std::const_pointer_cast<TrackingVolume>(volumesIter);
      mutableVolumesIter->setMotherVolume(this);
      mutableVolumesIter->closeGeometry(materialDecorator, volumeMap, vol);
    }
  }

  if (!m_confinedDenseVolumes.empty()) {
    for (auto& volumesIter : m_confinedDenseVolumes) {
      auto mutableVolumesIter =
          std::const_pointer_cast<TrackingVolume>(volumesIter);
      mutableVolumesIter->setMotherVolume(this);
      mutableVolumesIter->closeGeometry(materialDecorator, volumeMap, vol);
    }
  }
}

void Acts::TrackingVolume::visitSurfaces(
    const std::function<void(const Acts::Surface*)>& visitor) const {
  if (!m_confinedVolumes) {
    // no sub volumes => loop over the confined layers
    if (m_confinedLayers) {
      for (const auto& layer : m_confinedLayers->arrayObjects()) {
        if (layer->surfaceArray() == nullptr) {
          // no surface array (?)
          continue;
        }
        for (const auto& srf : layer->surfaceArray()->surfaces()) {
          visitor(srf);
        }
      }
    }
  } else {
    // contains sub volumes
    for (const auto& volume : m_confinedVolumes->arrayObjects()) {
      volume->visitSurfaces(visitor);
    }
  }
}

// Returns the boundary surfaces ordered in probability to hit them based on
std::vector<Acts::BoundaryIntersection>
Acts::TrackingVolume::compatibleBoundaries(
    const GeometryContext& gctx, const Vector3D& position,
    const Vector3D& direction, const NavigationOptions<Surface>& options,
    LoggerWrapper logger) const {
  ACTS_VERBOSE("Finding compatibleBoundaries");
  // Loop over boundarySurfaces and calculate the intersection
  auto excludeObject = options.startObject;
  std::vector<BoundaryIntersection> bIntersections;

  // The signed direction: solution (except overstepping) is positive
  auto sDirection = options.navDir * direction;

  // The Limits: current, path & overstepping
  double pLimit = options.pathLimit;
  double oLimit = options.overstepLimit;

  // Helper function to test intersection
  auto checkIntersection =
      [&](SurfaceIntersection& sIntersection,
          const BoundarySurface* bSurface) -> BoundaryIntersection {
    // Avoid doing anything if that's a rotten apple already
    if (!sIntersection) {
      return BoundaryIntersection();
    }

    ACTS_VERBOSE("Check intersection with surface "
                 << &bSurface->surfaceRepresentation());
    double cLimit = sIntersection.intersection.pathLength;
    ACTS_VERBOSE(" -> pLimit, oLimit, cLimit: " << pLimit << ", " << oLimit
                                                << ", " << cLimit);

    // Check if the surface is within limit
    bool withinLimit =
        (cLimit > oLimit and
         cLimit * cLimit <= pLimit * pLimit + s_onSurfaceTolerance);
    if (withinLimit) {
      ACTS_VERBOSE("Intersection is WITHIN limit");
      sIntersection.intersection.pathLength *=
          std::copysign(1., options.navDir);
      return BoundaryIntersection(sIntersection.intersection, bSurface,
                                  sIntersection.object);
    } else {
      ACTS_VERBOSE("Intersection is OUTSIDE limit");
    }

    // Check the alternative
    if (sIntersection.alternative) {
      ACTS_VERBOSE("Consider alternative");
      // Test the alternative
      cLimit = sIntersection.alternative.pathLength;
      ACTS_VERBOSE(" -> pLimit, oLimit, cLimit: " << pLimit << ", " << oLimit
                                                  << ", " << cLimit);
      withinLimit = (cLimit > oLimit and
                     cLimit * cLimit <= pLimit * pLimit + s_onSurfaceTolerance);
      if (sIntersection.alternative and withinLimit) {
        ACTS_VERBOSE("Intersection is WITHIN limit");
        sIntersection.alternative.pathLength *=
            std::copysign(1., options.navDir);
        return BoundaryIntersection(sIntersection.alternative, bSurface,
                                    sIntersection.object);
      } else {
        ACTS_VERBOSE("Intersection is OUTSIDE limit");
      }
    } else {
      ACTS_VERBOSE("No alternative for intersection");
    }
    // Return an invalid one
    ACTS_VERBOSE("No intersection accepted");
    return BoundaryIntersection();
  };

  auto processBoundaries =
      [&](const TrackingVolumeBoundaries& bSurfaces) -> void {
    ACTS_VERBOSE("Processing boundaries");
    // Loop over the boundary surfaces
    for (auto& bsIter : bSurfaces) {
      // Get the boundary surface pointer
      const auto& bSurfaceRep = bsIter->surfaceRepresentation();
      if (logger().doPrint(Logging::VERBOSE)) {
        auto os = logger().log(Logging::VERBOSE);
        os << "Consider boundary surface " << &bSurfaceRep << " :\n";
        std::stringstream strm;
        bSurfaceRep.toStream(gctx, strm);
        os << strm.str();
      }

      // Exclude the boundary where you are on
      if (excludeObject != &bSurfaceRep) {
        auto bCandidate = bSurfaceRep.intersect(gctx, position, sDirection,
                                                options.boundaryCheck);
        // Intersect and continue
        auto bIntersection = checkIntersection(bCandidate, bsIter.get());
        if (bIntersection) {
          ACTS_VERBOSE(" - Proceed with surface");
          bIntersections.push_back(bIntersection);
        } else {
          ACTS_VERBOSE(" - Surface intersecion invalid");
        }
      } else {
        ACTS_VERBOSE(" - Surface is excluded surface");
      }
    }
  };

  // Process the boundaries of the current volume
  auto& bSurfaces = boundarySurfaces();
  ACTS_VERBOSE("Volume reports " << bSurfaces.size() << " boundary surfaces");
  processBoundaries(bSurfaces);

  // Process potential boundaries of contained volumes
  auto confinedDenseVolumes = denseVolumes();
  ACTS_VERBOSE("Volume reports " << confinedDenseVolumes.size()
                                 << " confined dense volumes");
  for (const auto& dv : confinedDenseVolumes) {
    auto& bSurfacesConfined = dv->boundarySurfaces();
    ACTS_VERBOSE(" -> " << bSurfacesConfined.size() << " boundary surfaces");
    processBoundaries(bSurfacesConfined);
  }

  // Sort them accordingly to the navigation direction
  if (options.navDir == forward) {
    std::sort(bIntersections.begin(), bIntersections.end());
  } else {
    std::sort(bIntersections.begin(), bIntersections.end(), std::greater<>());
  }
  return bIntersections;
}

std::vector<Acts::LayerIntersection> Acts::TrackingVolume::compatibleLayers(
    const GeometryContext& gctx, const Vector3D& position,
    const Vector3D& direction, const NavigationOptions<Layer>& options) const {
  // the layer intersections which are valid
  std::vector<LayerIntersection> lIntersections;

  // the confinedLayers
  if (m_confinedLayers != nullptr) {
    // start layer given or not - test layer
    const Layer* tLayer = options.startObject != nullptr
                              ? options.startObject
                              : associatedLayer(gctx, position);
    while (tLayer != nullptr) {
      // check if the layer needs resolving
      // - resolveSensitive -> always take layer if it has a surface array
      // - resolveMaterial -> always take layer if it has material
      // - resolvePassive -> always take, unless it's a navigation layer
      // skip the start object
      if (tLayer != options.startObject && tLayer->resolve(options)) {
        // if it's a resolveable start layer, you are by definition on it
        // layer on approach intersection
        auto atIntersection =
            tLayer->surfaceOnApproach(gctx, position, direction, options);
        auto path = atIntersection.intersection.pathLength;
        bool withinLimit =
            (path * path <= options.pathLimit * options.pathLimit);
        // Intersection is ok - take it (move to surface on appraoch)
        if (atIntersection &&
            (atIntersection.object != options.targetSurface) && withinLimit) {
          // create a layer intersection
          lIntersections.push_back(LayerIntersection(
              atIntersection.intersection, tLayer, atIntersection.object));
        }
      }
      // move to next one or break because you reached the end layer
      tLayer =
          (tLayer == options.endObject)
              ? nullptr
              : tLayer->nextLayer(gctx, position, options.navDir * direction);
    }
    // sort them accordingly to the navigation direction
    if (options.navDir == forward) {
      std::sort(lIntersections.begin(), lIntersections.end());
    } else {
      std::sort(lIntersections.begin(), lIntersections.end(), std::greater<>());
    }
  }
  // and return
  return lIntersections;
}

namespace {
template <typename T>
std::vector<const Acts::Volume*> intersectSearchHierarchy(
    const T obj, const Acts::Volume::BoundingBox* lnode) {
  std::vector<const Acts::Volume*> hits;
  hits.reserve(20);  // arbitrary
  do {
    if (lnode->intersect(obj)) {
      if (lnode->hasEntity()) {
        // found primitive
        // check obb to limit false positivies
        const Acts::Volume* vol = lnode->entity();
        const auto& obb = vol->orientedBoundingBox();
        if (obb.intersect(obj.transformed(vol->itransform()))) {
          hits.push_back(vol);
        }
        // we skip in any case, whether we actually hit the OBB or not
        lnode = lnode->getSkip();
      } else {
        // go over children
        lnode = lnode->getLeftChild();
      }
    } else {
      lnode = lnode->getSkip();
    }
  } while (lnode != nullptr);

  return hits;
}
}  // namespace

std::vector<Acts::SurfaceIntersection>
Acts::TrackingVolume::compatibleSurfacesFromHierarchy(
    const GeometryContext& gctx, const Vector3D& position,
    const Vector3D& direction, double angle,
    const NavigationOptions<Surface>& options) const {
  std::vector<SurfaceIntersection> sIntersections;
  sIntersections.reserve(20);  // arbitrary

  // The limits for this navigation step
  double pLimit = options.pathLimit;
  double oLimit = options.overstepLimit;

  if (m_bvhTop == nullptr || !options.navDir) {
    return sIntersections;
  }

  // The signed direction
  Vector3D sdir = options.navDir * direction;

  std::vector<const Volume*> hits;
  if (angle == 0) {
    // use ray
    Ray3D obj(position, sdir);
    hits = intersectSearchHierarchy(std::move(obj), m_bvhTop);
  } else {
    Acts::Frustum<double, 3, 4> obj(position, sdir, angle);
    hits = intersectSearchHierarchy(std::move(obj), m_bvhTop);
  }

  // have cells, decompose to surfaces
  for (const Volume* vol : hits) {
    const AbstractVolume* avol = dynamic_cast<const AbstractVolume*>(vol);
    const std::vector<std::shared_ptr<const BoundarySurfaceT<AbstractVolume>>>&
        boundarySurfaces = avol->boundarySurfaces();
    for (const auto& bs : boundarySurfaces) {
      const Surface& srf = bs->surfaceRepresentation();
      auto sfi = srf.intersect(gctx, position, sdir, false);
      if (sfi and sfi.intersection.pathLength > oLimit and
          sfi.intersection.pathLength <= pLimit) {
        sIntersections.push_back(std::move(sfi));
      }
    }
  }

  // Sort according to the path length
  if (options.navDir == forward) {
    std::sort(sIntersections.begin(), sIntersections.end());
  } else {
    std::sort(sIntersections.begin(), sIntersections.end(), std::greater<>());
  }

  return sIntersections;
}