MappingMaterialDecorator.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2017-2022 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/.

#pragma once

#include "Acts/Geometry/TrackingVolume.hpp"
#include "Acts/Material/IMaterialDecorator.hpp"
#include "Acts/Material/ISurfaceMaterial.hpp"
#include "Acts/Material/IVolumeMaterial.hpp"
#include "Acts/Material/ProtoSurfaceMaterial.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include <Acts/Surfaces/AnnulusBounds.hpp>
#include <Acts/Surfaces/CylinderBounds.hpp>
#include <Acts/Surfaces/RadialBounds.hpp>
#include <Acts/Surfaces/SurfaceBounds.hpp>
#include <Acts/Surfaces/TrapezoidBounds.hpp>

#include <fstream>
#include <map>
#include <mutex>

// Convenience shorthand

namespace Acts {

class MappingMaterialDecorator : public IMaterialDecorator {
 public:
  using BinningMap = std::map<uint64_t, std::pair<int, int>>;

  using VolumeMaterialMap =
      std::map<GeometryIdentifier, std::shared_ptr<const IVolumeMaterial>>;

  MappingMaterialDecorator(const Acts::TrackingGeometry& tGeometry,
                           Acts::Logging::Level level,
                           bool clearSurfaceMaterial = true,
                           bool clearVolumeMaterial = true)
      : m_clearSurfaceMaterial(clearSurfaceMaterial),
        m_clearVolumeMaterial(clearVolumeMaterial),
        m_logger{getDefaultLogger("MappingMaterialDecorator", level)} {
    volumeLoop(tGeometry.highestTrackingVolume());
  }

  void decorate(Surface& surface) const final {
    ACTS_VERBOSE("Processing surface: " << surface.geometryId());
    // Clear the material if registered to do so
    if (m_clearSurfaceMaterial) {
      ACTS_VERBOSE("-> Clearing surface material");
      surface.assignSurfaceMaterial(nullptr);
    }
    // Try to find the surface in the map
    auto bins = m_binningMap.find(surface.geometryId().value());
    if (bins != m_binningMap.end()) {
      ACTS_VERBOSE("-> Found material for surface, assigning");
      surface.assignSurfaceMaterial(
          binnedSurfaceMaterial(surface.getSharedPtr()));
    }
  }

  void decorate(TrackingVolume& volume) const final {
    ACTS_VERBOSE("Processing volume: " << volume.geometryId());
    // Clear the material if registered to do so
    if (m_clearVolumeMaterial) {
      ACTS_VERBOSE("-> Clearing volume material");
      volume.assignVolumeMaterial(nullptr);
    }
    // Try to find the volume in the map
    auto vMaterial = m_volumeMaterialMap.find(volume.geometryId());
    if (vMaterial != m_volumeMaterialMap.end()) {
      ACTS_VERBOSE("-> Found material for volume, assigning");
      volume.assignVolumeMaterial(vMaterial->second);
    }
  }

  void volumeLoop(const Acts::TrackingVolume* tVolume) {
    auto sameId =
        [tVolume](
            const std::pair<Acts::GeometryIdentifier,
                            std::shared_ptr<const IVolumeMaterial>>& pair) {
          return (tVolume->geometryId() == pair.first);
        };
    if (std::find_if(m_volumeMaterialMap.begin(), m_volumeMaterialMap.end(),
                     sameId) != m_volumeMaterialMap.end()) {
      // this volume was already visited
      return;
    }
    if (tVolume->volumeMaterial() != nullptr) {
      m_volumeMaterialMap.insert(
          {tVolume->geometryId(), tVolume->volumeMaterialSharedPtr()});
    }
    // there are confined volumes
    if (tVolume->confinedVolumes() != nullptr) {
      // get through the volumes
      auto& volumes = tVolume->confinedVolumes()->arrayObjects();
      // loop over the volumes
      for (auto& vol : volumes) {
        // recursive call
        volumeLoop(vol.get());
      }
    }
    // there are dense volumes
    if (!tVolume->denseVolumes().empty()) {
      // loop over the volumes
      for (auto& vol : tVolume->denseVolumes()) {
        // recursive call
        volumeLoop(vol.get());
      }
    }
    if (tVolume->confinedLayers() != nullptr) {
      // get the layers
      auto& layers = tVolume->confinedLayers()->arrayObjects();
      // loop over the layers
      for (auto& lay : layers) {
        auto& layRep = lay->surfaceRepresentation();
        if (layRep.surfaceMaterial() != nullptr &&
            layRep.geometryId() != GeometryIdentifier()) {
          m_binningMap.insert(
              {layRep.geometryId().value(), std::make_pair(1, 1)});
        }
        if (lay->approachDescriptor() != nullptr) {
          for (auto& asf : lay->approachDescriptor()->containedSurfaces()) {
            if (asf->surfaceMaterial() != nullptr) {
              m_binningMap.insert(
                  {asf->geometryId().value(), std::make_pair(1, 1)});
            }
          }
        }
        if (lay->surfaceArray() != nullptr) {
          for (auto& ssf : lay->surfaceArray()->surfaces()) {
            if (ssf->surfaceMaterial() != nullptr) {
              m_binningMap.insert(
                  {ssf->geometryId().value(), std::make_pair(1, 1)});
            }
          }
        }
      }
    }
    // Let's finally check the boundaries
    for (auto& bsurf : tVolume->boundarySurfaces()) {
      // the surface representation
      auto& bssfRep = bsurf->surfaceRepresentation();
      if (bssfRep.geometryId().volume() == tVolume->geometryId().volume()) {
        if (bssfRep.surfaceMaterial() != nullptr) {
          m_binningMap.insert(
              {bssfRep.geometryId().value(), std::make_pair(1, 1)});
        }
      }
    }
  }

  std::shared_ptr<const Acts::ISurfaceMaterial> binnedSurfaceMaterial(
      const std::shared_ptr<const Acts::Surface>& surface) const {
    auto bin = m_binningMap.find(surface->geometryId().value());
    Acts::BinUtility bUtility;
    if (bin == m_binningMap.end()) {
      ACTS_ERROR("No corresponding binning in the map to surface "
                 << surface->geometryId());
    } else {
      auto binning = bin->second;
      // Check which type of bounds is associated to the surface
      const Acts::SurfaceBounds& surfaceBounds = surface->bounds();
      const Acts::RadialBounds* radialBounds =
          dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds);
      const Acts::CylinderBounds* cylinderBounds =
          dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds);
      const Acts::AnnulusBounds* annulusBounds =
          dynamic_cast<const Acts::AnnulusBounds*>(&surfaceBounds);
      const Acts::RectangleBounds* rectangleBounds =
          dynamic_cast<const Acts::RectangleBounds*>(&surfaceBounds);
      const Acts::TrapezoidBounds* trapezoidBounds =
          dynamic_cast<const Acts::TrapezoidBounds*>(&surfaceBounds);

      if (radialBounds != nullptr) {
        bUtility += Acts::BinUtility(
            binning.first,
            radialBounds->get(Acts::RadialBounds::eAveragePhi) -
                radialBounds->get(Acts::RadialBounds::eHalfPhiSector),
            radialBounds->get(Acts::RadialBounds::eAveragePhi) +
                radialBounds->get(Acts::RadialBounds::eHalfPhiSector),
            (radialBounds->get(Acts::RadialBounds::eHalfPhiSector) - M_PI) <
                    Acts::s_epsilon
                ? Acts::closed
                : Acts::open,
            Acts::binPhi);
        bUtility +=
            Acts::BinUtility(binning.second, radialBounds->rMin(),
                             radialBounds->rMax(), Acts::open, Acts::binR);
      }
      if (cylinderBounds != nullptr) {
        bUtility += Acts::BinUtility(
            binning.first,
            cylinderBounds->get(Acts::CylinderBounds::eAveragePhi) -
                cylinderBounds->get(Acts::CylinderBounds::eHalfPhiSector),
            cylinderBounds->get(Acts::CylinderBounds::eAveragePhi) +
                cylinderBounds->get(Acts::CylinderBounds::eHalfPhiSector),
            (cylinderBounds->get(Acts::CylinderBounds::eHalfPhiSector) - M_PI) <
                    Acts::s_epsilon
                ? Acts::closed
                : Acts::open,
            Acts::binPhi);
        bUtility += Acts::BinUtility(
            binning.second,
            -1 * cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ),
            cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ), Acts::open,
            Acts::binZ);
      }
      if (annulusBounds != nullptr) {
        bUtility += Acts::BinUtility(
            binning.first, annulusBounds->get(Acts::AnnulusBounds::eMinPhiRel),
            annulusBounds->get(Acts::AnnulusBounds::eMaxPhiRel), Acts::open,
            Acts::binPhi);
        bUtility +=
            Acts::BinUtility(binning.second, annulusBounds->rMin(),
                             annulusBounds->rMax(), Acts::open, Acts::binR);
      }
      if (rectangleBounds != nullptr) {
        bUtility += Acts::BinUtility(
            binning.first, rectangleBounds->get(Acts::RectangleBounds::eMinX),
            rectangleBounds->get(Acts::RectangleBounds::eMaxX), Acts::open,
            Acts::binX);
        bUtility += Acts::BinUtility(
            binning.second, rectangleBounds->get(Acts::RectangleBounds::eMinY),
            rectangleBounds->get(Acts::RectangleBounds::eMaxY), Acts::open,
            Acts::binY);
      }
      if (trapezoidBounds != nullptr) {
        double halfLengthX = std::max(
            trapezoidBounds->get(Acts::TrapezoidBounds::eHalfLengthXnegY),
            trapezoidBounds->get(Acts::TrapezoidBounds::eHalfLengthXposY));
        bUtility += Acts::BinUtility(binning.first, -1 * halfLengthX,
                                     halfLengthX, Acts::open, Acts::binX);
        bUtility += Acts::BinUtility(
            binning.second,
            -1 * trapezoidBounds->get(Acts::TrapezoidBounds::eHalfLengthY),
            trapezoidBounds->get(Acts::TrapezoidBounds::eHalfLengthY),
            Acts::open, Acts::binY);
      }
    }
    return std::make_shared<Acts::ProtoSurfaceMaterial>(bUtility);
  }

  BinningMap& binningMap() { return m_binningMap; }

  void setBinningMap(BinningMap binningMap) {
    m_binningMap = std::move(binningMap);
  }

 private:
  BinningMap m_binningMap;

  VolumeMaterialMap m_volumeMaterialMap;

  bool m_clearSurfaceMaterial{true};
  bool m_clearVolumeMaterial{true};

  std::unique_ptr<const Logger> m_logger;

  const Logger& logger() const { return *m_logger; }
};
}  // namespace Acts