LayerCreator.cpp

// This file is part of the ACTS project.
//
// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.

#include "Acts/Geometry/LayerCreator.hpp"

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Geometry/CylinderLayer.hpp"
#include "Acts/Geometry/DiscLayer.hpp"
#include "Acts/Geometry/Extent.hpp"
#include "Acts/Geometry/Layer.hpp"
#include "Acts/Geometry/PlaneLayer.hpp"
#include "Acts/Geometry/ProtoLayer.hpp"
#include "Acts/Geometry/SurfaceArrayCreator.hpp"
#include "Acts/Surfaces/CylinderBounds.hpp"
#include "Acts/Surfaces/RadialBounds.hpp"
#include "Acts/Surfaces/RectangleBounds.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/Diagnostics.hpp"

#include <algorithm>
#include <iterator>
#include <ostream>
#include <set>
#include <utility>

namespace Acts {

using VectorHelpers::perp;
using VectorHelpers::phi;

LayerCreator::LayerCreator(const LayerCreator::Config& lcConfig,
                           std::unique_ptr<const Logger> logger)
    : m_cfg(lcConfig), m_logger(std::move(logger)) {}

void LayerCreator::setConfiguration(const LayerCreator::Config& lcConfig) {
  // @todo check consistency
  // copy the configuration
  m_cfg = lcConfig;
}

void LayerCreator::setLogger(std::unique_ptr<const Logger> newLogger) {
  m_logger = std::move(newLogger);
}

MutableLayerPtr LayerCreator::cylinderLayer(
    const GeometryContext& gctx,
    std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsPhi,
    std::size_t binsZ, std::optional<ProtoLayer> _protoLayer,
    const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad,
    std::uint8_t maxNeighborDistance) const {
  ProtoLayer protoLayer =
      _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
  if (!_protoLayer) {
    protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
    protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
  }

  // Remaining layer parameters - they include the envelopes
  double layerR = protoLayer.medium(AxisDirection::AxisR);
  double layerZ = protoLayer.medium(AxisDirection::AxisZ);
  double layerHalfZ = 0.5 * protoLayer.range(AxisDirection::AxisZ);
  double layerThickness = protoLayer.range(AxisDirection::AxisR);

  ACTS_VERBOSE("Creating a cylindrical Layer:");
  ACTS_VERBOSE(" - with layer R     = " << layerR);
  ACTS_VERBOSE(" - from R min/max   = "
               << protoLayer.min(AxisDirection::AxisR, false) << " / "
               << protoLayer.max(AxisDirection::AxisR, false));
  ACTS_VERBOSE(" - with R thickness = " << layerThickness);
  ACTS_VERBOSE("   - incl envelope  = "
               << protoLayer.envelope[AxisDirection::AxisR][0u] << " / "
               << protoLayer.envelope[AxisDirection::AxisR][1u]);

  ACTS_VERBOSE(" - with z min/max   = "
               << protoLayer.min(AxisDirection::AxisZ, false) << " (-"
               << protoLayer.envelope[AxisDirection::AxisZ][0u] << ") / "
               << protoLayer.max(AxisDirection::AxisZ, false) << " (+"
               << protoLayer.envelope[AxisDirection::AxisZ][1u] << ")");

  ACTS_VERBOSE(" - z center         = " << layerZ);
  ACTS_VERBOSE(" - halflength z     = " << layerHalfZ);

  // create the layer transforms if not given
  // we need to transform in case layerZ != 0, so that the layer will be
  // correctly defined using the halflength
  Transform3 fullTransform = transform;
  if (fullTransform.isApprox(Transform3::Identity())) {
    fullTransform = Translation3(0, 0, layerZ) * fullTransform;
    ACTS_VERBOSE(" - layer z shift  = " << -layerZ);
  }

  ACTS_VERBOSE(" - with phi min/max = "
               << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
               << protoLayer.max(AxisDirection::AxisPhi, false));
  ACTS_VERBOSE(" - # of modules     = " << surfaces.size() << " ordered in ( "
                                        << binsPhi << " x " << binsZ << ")");
  std::unique_ptr<SurfaceArray> sArray;
  if (!surfaces.empty()) {
    sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnCylinder(
        gctx, std::move(surfaces), binsPhi, binsZ, protoLayer, fullTransform,
        maxNeighborDistance);
  }

  // create the layer and push it back
  std::shared_ptr<const CylinderBounds> cBounds(
      new CylinderBounds(layerR, layerHalfZ));

  // create the layer
  MutableLayerPtr cLayer =
      CylinderLayer::create(fullTransform, cBounds, std::move(sArray),
                            layerThickness, std::move(ad), active);

  if (!cLayer) {
    ACTS_ERROR("Creation of cylinder layer did not succeed!");
  }
  associateSurfacesToLayer(*cLayer);

  // now return
  return cLayer;
}

MutableLayerPtr LayerCreator::cylinderLayer(
    const GeometryContext& gctx,
    std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypePhi,
    BinningType bTypeZ, std::optional<ProtoLayer> _protoLayer,
    const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad,
    std::uint8_t maxNeighborDistance) const {
  ProtoLayer protoLayer =
      _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
  if (!_protoLayer) {
    protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
    protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
  }

  // remaining layer parameters
  double layerR = protoLayer.medium(AxisDirection::AxisR);
  double layerZ = protoLayer.medium(AxisDirection::AxisZ);
  double layerHalfZ = 0.5 * protoLayer.range(AxisDirection::AxisZ);
  double layerThickness = protoLayer.range(AxisDirection::AxisR);

  // adjust the layer radius
  ACTS_VERBOSE("Creating a cylindrical Layer:");
  ACTS_VERBOSE(" - with layer R     = " << layerR);
  ACTS_VERBOSE(" - from R min/max   = "
               << protoLayer.min(AxisDirection::AxisR, false) << " / "
               << protoLayer.max(AxisDirection::AxisR, false));
  ACTS_VERBOSE(" - with R thickness = " << layerThickness);
  ACTS_VERBOSE("   - incl envelope  = "
               << protoLayer.envelope[AxisDirection::AxisR][0u] << " / "
               << protoLayer.envelope[AxisDirection::AxisR][1u]);
  ACTS_VERBOSE(" - with z min/max   = "
               << protoLayer.min(AxisDirection::AxisZ, false) << " (-"
               << protoLayer.envelope[AxisDirection::AxisZ][0u] << ") / "
               << protoLayer.max(AxisDirection::AxisZ, false) << " (+"
               << protoLayer.envelope[AxisDirection::AxisZ][1u] << ")");
  ACTS_VERBOSE(" - z center         = " << layerZ);
  ACTS_VERBOSE(" - halflength z     = " << layerHalfZ);

  // create the layer transforms if not given
  // we need to transform in case layerZ != 0, so that the layer will be
  // correctly defined using the halflength
  // create the layer transforms if not given
  Transform3 fullTransform = transform;
  if (fullTransform.isApprox(Transform3::Identity()) && bTypeZ == equidistant) {
    fullTransform = Translation3(0, 0, layerZ) * fullTransform;
    ACTS_VERBOSE(" - layer z shift    = " << -layerZ);
  }

  ACTS_VERBOSE(" - with phi min/max = "
               << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
               << protoLayer.max(AxisDirection::AxisPhi, false));
  ACTS_VERBOSE(" - # of modules     = " << surfaces.size() << "");

  // create the surface array
  std::unique_ptr<SurfaceArray> sArray;
  if (!surfaces.empty()) {
    sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnCylinder(
        gctx, std::move(surfaces), bTypePhi, bTypeZ, protoLayer, fullTransform,
        maxNeighborDistance);
  }

  // create the layer and push it back
  std::shared_ptr<const CylinderBounds> cBounds(
      new CylinderBounds(layerR, layerHalfZ));

  // create the layer
  MutableLayerPtr cLayer =
      CylinderLayer::create(fullTransform, cBounds, std::move(sArray),
                            layerThickness, std::move(ad), active);

  if (!cLayer) {
    ACTS_ERROR("Creation of cylinder layer did not succeed!");
  }
  associateSurfacesToLayer(*cLayer);

  // now return
  return cLayer;
}

MutableLayerPtr LayerCreator::discLayer(
    const GeometryContext& gctx,
    std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsR,
    std::size_t binsPhi, std::optional<ProtoLayer> _protoLayer,
    const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad,
    std::uint8_t maxNeighborDistance) const {
  ProtoLayer protoLayer =
      _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
  if (!_protoLayer) {
    protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
    protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
  }

  double layerZ = protoLayer.medium(AxisDirection::AxisZ);
  double layerThickness = protoLayer.range(AxisDirection::AxisZ);

  // adjust the layer radius
  ACTS_VERBOSE("Creating a disk Layer:");
  ACTS_VERBOSE(" - at Z position    = " << layerZ);
  ACTS_VERBOSE(" - from Z min/max   = "
               << protoLayer.min(AxisDirection::AxisZ, false) << " / "
               << protoLayer.max(AxisDirection::AxisZ, false));
  ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
  ACTS_VERBOSE("   - incl envelope  = "
               << protoLayer.envelope[AxisDirection::AxisZ][0u] << " / "
               << protoLayer.envelope[AxisDirection::AxisZ][1u]);
  ACTS_VERBOSE(" - with R min/max   = "
               << protoLayer.min(AxisDirection::AxisR, false) << " (-"
               << protoLayer.envelope[AxisDirection::AxisR][0u] << ") / "
               << protoLayer.max(AxisDirection::AxisR, false) << " (+"
               << protoLayer.envelope[AxisDirection::AxisR][1u] << ")");
  ACTS_VERBOSE(" - with phi min/max = "
               << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
               << protoLayer.max(AxisDirection::AxisPhi, false));
  ACTS_VERBOSE(" - # of modules    = " << surfaces.size() << " ordered in ( "
                                       << binsR << " x " << binsPhi << ")");

  // create the layer transforms if not given
  Transform3 fullTransform = transform;
  if (fullTransform.isApprox(Transform3::Identity())) {
    fullTransform = Translation3(0, 0, layerZ) * fullTransform;
  }
  // create the surface array
  std::unique_ptr<SurfaceArray> sArray;
  if (!surfaces.empty()) {
    sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnDisc(
        gctx, std::move(surfaces), binsR, binsPhi, protoLayer, fullTransform,
        maxNeighborDistance);
  }

  // create the share disc bounds
  auto dBounds = std::make_shared<const RadialBounds>(
      protoLayer.min(AxisDirection::AxisR),
      protoLayer.max(AxisDirection::AxisR));

  // create the layers
  // we use the same transform here as for the layer itself
  // for disk this is fine since we don't bin in Z, so does not matter
  MutableLayerPtr dLayer =
      DiscLayer::create(fullTransform, dBounds, std::move(sArray),
                        layerThickness, std::move(ad), active);

  if (!dLayer) {
    ACTS_ERROR("Creation of disc layer did not succeed!");
  }
  associateSurfacesToLayer(*dLayer);
  // return the layer
  return dLayer;
}

MutableLayerPtr LayerCreator::discLayer(
    const GeometryContext& gctx,
    std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypeR,
    BinningType bTypePhi, std::optional<ProtoLayer> _protoLayer,
    const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad,
    std::uint8_t maxNeighborDistance) const {
  ProtoLayer protoLayer =
      _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
  if (!_protoLayer) {
    protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
    protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
  }

  const double layerZ = protoLayer.medium(AxisDirection::AxisZ);
  const double layerThickness = protoLayer.range(AxisDirection::AxisZ);

  // adjust the layer radius
  ACTS_VERBOSE("Creating a disk Layer:");
  ACTS_VERBOSE(" - at Z position    = " << layerZ);
  ACTS_VERBOSE(" - from Z min/max   = "
               << protoLayer.min(AxisDirection::AxisZ, false) << " / "
               << protoLayer.max(AxisDirection::AxisZ, false));
  ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
  ACTS_VERBOSE("   - incl envelope  = "
               << protoLayer.envelope[AxisDirection::AxisZ][0u] << " / "
               << protoLayer.envelope[AxisDirection::AxisZ][1u]);
  ACTS_VERBOSE(" - with R min/max   = "
               << protoLayer.min(AxisDirection::AxisR, false) << " (-"
               << protoLayer.envelope[AxisDirection::AxisR][0u] << ") / "
               << protoLayer.max(AxisDirection::AxisR, false) << " (+"
               << protoLayer.envelope[AxisDirection::AxisR][1u] << ")");
  ACTS_VERBOSE(" - with phi min/max = "
               << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
               << protoLayer.max(AxisDirection::AxisPhi, false));
  ACTS_VERBOSE(" - # of modules     = " << surfaces.size());

  // create the layer transforms if not given
  Transform3 fullTransform = transform;
  if (fullTransform.isApprox(Transform3::Identity())) {
    fullTransform = Translation3(0, 0, layerZ) * fullTransform;
  }

  // create the surface array
  std::unique_ptr<SurfaceArray> sArray;
  if (!surfaces.empty()) {
    sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnDisc(
        gctx, std::move(surfaces), bTypeR, bTypePhi, protoLayer, fullTransform,
        maxNeighborDistance);
  }

  // create the shared disc bounds
  auto dBounds = std::make_shared<const RadialBounds>(
      protoLayer.min(AxisDirection::AxisR),
      protoLayer.max(AxisDirection::AxisR));

  // create the layers
  MutableLayerPtr dLayer =
      DiscLayer::create(fullTransform, dBounds, std::move(sArray),
                        layerThickness, std::move(ad), active);
  if (!dLayer) {
    ACTS_ERROR("Creation of disc layer did not succeed!");
  }
  associateSurfacesToLayer(*dLayer);
  // return the layer
  return dLayer;
}

MutableLayerPtr LayerCreator::planeLayer(
    const GeometryContext& gctx,
    std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t bins1,
    std::size_t bins2, AxisDirection aDir,
    std::optional<ProtoLayer> _protoLayer, const Transform3& transform,
    std::unique_ptr<ApproachDescriptor> ad,
    std::uint8_t maxNeighborDistance) const {
  ProtoLayer protoLayer =
      _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
  if (!_protoLayer) {
    protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
    protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
  }

  // remaining layer parameters
  double layerHalf1 = 0, layerHalf2 = 0, layerThickness = 0;
  switch (aDir) {
    case AxisDirection::AxisX: {
      layerHalf1 = 0.5 * protoLayer.range(AxisDirection::AxisY);
      layerHalf2 = 0.5 * protoLayer.range(AxisDirection::AxisZ);
      layerThickness = protoLayer.range(AxisDirection::AxisX);
      break;
    }
    case AxisDirection::AxisY: {
      layerHalf1 = 0.5 * protoLayer.range(AxisDirection::AxisX);
      layerHalf2 = 0.5 * protoLayer.range(AxisDirection::AxisZ);
      layerThickness = protoLayer.range(AxisDirection::AxisY);
      break;
    }
    case AxisDirection::AxisZ: {
      layerHalf1 = 0.5 * protoLayer.range(AxisDirection::AxisX);
      layerHalf2 = 0.5 * protoLayer.range(AxisDirection::AxisY);
      layerThickness = protoLayer.range(AxisDirection::AxisZ);
      break;
    }
    default:
      throw std::invalid_argument("Invalid binning value");
  }

  const double centerX = protoLayer.medium(AxisDirection::AxisX);
  const double centerY = protoLayer.medium(AxisDirection::AxisY);
  const double centerZ = protoLayer.medium(AxisDirection::AxisZ);

  ACTS_VERBOSE("Creating a plane Layer:");
  ACTS_VERBOSE(" - with layer center     = "
               << "(" << centerX << ", " << centerY << ", " << centerZ << ")");
  ACTS_VERBOSE(" - from X min/max   = "
               << protoLayer.min(AxisDirection::AxisX) << " / "
               << protoLayer.max(AxisDirection::AxisX));
  ACTS_VERBOSE(" - from Y min/max   = "
               << protoLayer.min(AxisDirection::AxisY) << " / "
               << protoLayer.max(AxisDirection::AxisY));
  ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
  ACTS_VERBOSE("   - incl envelope  = " << protoLayer.envelope[aDir][0u]
                                        << " / "
                                        << protoLayer.envelope[aDir][1u]);

  // create the layer transforms if not given
  // we need to transform in case centerX/centerY/centerZ != 0, so that the
  // layer will be correctly defined
  Transform3 fullTransform = transform;
  if (fullTransform.isApprox(Transform3::Identity())) {
    fullTransform = Translation3(centerX, centerY, centerZ) * fullTransform;
    ACTS_VERBOSE(" - layer shift  = " << "(" << centerX << ", " << centerY
                                      << ", " << centerZ << ")");
  }

  std::unique_ptr<SurfaceArray> sArray;
  if (!surfaces.empty()) {
    sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnPlane(
        gctx, std::move(surfaces), bins1, bins2, aDir, protoLayer,
        fullTransform, maxNeighborDistance);
  }

  // create the layer and push it back
  auto pBounds = std::make_shared<RectangleBounds>(layerHalf1, layerHalf2);

  // create the layer
  MutableLayerPtr pLayer =
      PlaneLayer::create(fullTransform, pBounds, std::move(sArray),
                         layerThickness, std::move(ad), active);

  if (!pLayer) {
    ACTS_ERROR("Creation of plane layer did not succeed!");
  }
  associateSurfacesToLayer(*pLayer);

  // now return
  return pLayer;
}

void LayerCreator::associateSurfacesToLayer(Layer& layer) const {
  if (layer.surfaceArray() != nullptr) {
    auto surfaces = layer.surfaceArray()->surfaces();

    for (auto& surface : surfaces) {
      auto mutableSurface = const_cast<Surface*>(surface);
      mutableSurface->associateLayer(layer);
    }
  }
}

}  // namespace Acts