src/Geometry/LayerArrayCreator.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/Geometry/LayerArrayCreator.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryObjectSorter.hpp"
0013 #include "Acts/Geometry/Layer.hpp"
0014 #include "Acts/Geometry/NavigationLayer.hpp"
0015 #include "Acts/Surfaces/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/CylinderSurface.hpp"
0017 #include "Acts/Surfaces/DiscSurface.hpp"
0018 #include "Acts/Surfaces/PlaneSurface.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Surfaces/SurfaceBounds.hpp"
0021 #include "Acts/Utilities/BinUtility.hpp"
0022 #include "Acts/Utilities/BinnedArrayXD.hpp"
0023 #include "Acts/Utilities/BinningType.hpp"
0024
0025 #include <algorithm>
0026 #include <memory>
0027 #include <ostream>
0028 #include <utility>
0029 #include <vector>
0030
0031 std::unique_ptr<const Acts::LayerArray> Acts::LayerArrayCreator::layerArray(
0032     const GeometryContext& gctx, const LayerVector& layersInput, double min,
0033     double max, BinningType bType, AxisDirection aDir) const {
0034   ACTS_VERBOSE("Build LayerArray with " << layersInput.size()
0035                                         << " layers at input.");
0036   ACTS_VERBOSE("       min/max provided : " << min << " / " << max);
0037   ACTS_VERBOSE("       binning type     : " << bType);
0038   ACTS_VERBOSE("       binning value    : " << aDir);
0039
0040   // create a local copy of the layer vector
0041   LayerVector layers(layersInput);
0042
0043   // sort it accordingly to the binning value
0044   GeometryObjectSorterT<std::shared_ptr<const Layer>> layerSorter(gctx, aDir);
0045   std::ranges::sort(layers, layerSorter);
0046   // useful typedef
0047   using LayerOrderPosition = std::pair<std::shared_ptr<const Layer>, Vector3>;
0048   // needed for all cases
0049   std::shared_ptr<const Layer> layer = nullptr;
0050   std::unique_ptr<const BinUtility> binUtility = nullptr;
0051   std::vector<LayerOrderPosition> layerOrderVector;
0052
0053   // switch the binning type
0054   switch (bType) {
0055     // equidistant binning - no navigation layers built - only equdistant layers
0056     case equidistant: {
0057       // loop over layers and put them in
0058       for (auto& layIter : layers) {
0059         ACTS_VERBOSE("equidistant : registering a Layer at binning position : "
0060                      << (layIter->referencePosition(gctx, aDir)));
0061         layerOrderVector.push_back(LayerOrderPosition(
0062             layIter, layIter->referencePosition(gctx, aDir)));
0063       }
0064       // create the binUitlity
0065       binUtility = std::make_unique<const BinUtility>(layers.size(), min, max,
0066                                                       open, aDir);
0067       ACTS_VERBOSE("equidistant : created a BinUtility as " << *binUtility);
0068     } break;
0069
0070     // arbitrary binning
0071     case arbitrary: {
0072       std::vector<float> boundaries;
0073       // initial step
0074       boundaries.push_back(min);
0075       double layerValue = 0.;
0076       double layerThickness = 0.;
0077       std::shared_ptr<const Layer> navLayer = nullptr;
0078       std::shared_ptr<const Layer> lastLayer = nullptr;
0079       // loop over layers
0080       for (auto& layIter : layers) {
0081         // estimate the offset
0082         layerThickness = layIter->thickness();
0083         layerValue = layIter->referencePositionValue(gctx, aDir);
0084         // register the new boundaries in the step vector
0085         boundaries.push_back(layerValue - 0.5 * layerThickness);
0086         boundaries.push_back(layerValue + 0.5 * layerThickness);
0087         // calculate the layer value for the offset
0088         double navigationValue = 0.5 * ((layerValue - 0.5 * layerThickness) +
0089                                         boundaries.at(boundaries.size() - 3));
0090         // if layers are attached to each other bail out - navigation will not
0091         // work anymore
0092         if (navigationValue == (layerValue - 0.5 * layerThickness)) {
0093           ACTS_ERROR(
0094               "Layers are attached to each other at: "
0095               << layerValue - 0.5 * layerThickness
0096               << ", which corrupts "
0097                  "navigation. This should never happen. Please detach the "
0098                  "layers in your geometry description.");
0099         }
0100         // if layers are overlapping bail out
0101         if (navigationValue > (layerValue - 0.5 * layerThickness)) {
0102           ACTS_ERROR("Layers are overlapping at: "
0103                      << layerValue - 0.5 * layerThickness
0104                      << ". This should never happen. "
0105                         "Please check your geometry description.");
0106         }
0107
0108         // create the navigation layer surface from the layer
0109         std::shared_ptr<const Surface> navLayerSurface =
0110             createNavigationSurface(gctx, *layIter, aDir,
0111                                     -std::abs(layerValue - navigationValue));
0112         ACTS_VERBOSE(
0113             "arbitrary : creating a  NavigationLayer at "
0114             << (navLayerSurface->referencePosition(gctx, aDir)).x() << ", "
0115             << (navLayerSurface->referencePosition(gctx, aDir)).y() << ", "
0116             << (navLayerSurface->referencePosition(gctx, aDir)).z());
0117         navLayer = NavigationLayer::create(std::move(navLayerSurface));
0118         // push the navigation layer in
0119         layerOrderVector.push_back(LayerOrderPosition(
0120             navLayer, navLayer->referencePosition(gctx, aDir)));
0121
0122         // push the original layer in
0123         layerOrderVector.push_back(LayerOrderPosition(
0124             layIter, layIter->referencePosition(gctx, aDir)));
0125         ACTS_VERBOSE("arbitrary : registering MaterialLayer at  "
0126                      << (layIter->referencePosition(gctx, aDir)).x() << ", "
0127                      << (layIter->referencePosition(gctx, aDir)).y() << ", "
0128                      << (layIter->referencePosition(gctx, aDir)).z());
0129         // remember the last
0130         lastLayer = layIter;
0131       }
0132       // a final navigation layer
0133       // calculate the layer value for the offset
0134       double navigationValue =
0135           0.5 * (boundaries.at(boundaries.size() - 1) + max);
0136       // create navigation layer only when necessary
0137       if (navigationValue != max && lastLayer != nullptr) {
0138         // create the navigation layer surface from the layer
0139         std::shared_ptr<const Surface> navLayerSurface =
0140             createNavigationSurface(gctx, *lastLayer, aDir,
0141                                     navigationValue - layerValue);
0142         ACTS_VERBOSE(
0143             "arbitrary : creating a  NavigationLayer at "
0144             << (navLayerSurface->referencePosition(gctx, aDir)).x() << ", "
0145             << (navLayerSurface->referencePosition(gctx, aDir)).y() << ", "
0146             << (navLayerSurface->referencePosition(gctx, aDir)).z());
0147         navLayer = NavigationLayer::create(std::move(navLayerSurface));
0148         // push the navigation layer in
0149         layerOrderVector.push_back(LayerOrderPosition(
0150             navLayer, navLayer->referencePosition(gctx, aDir)));
0151       }
0152       // now close the boundaries
0153       boundaries.push_back(max);
0154       // some screen output
0155       ACTS_VERBOSE(layerOrderVector.size()
0156                    << " Layers (material + navigation) built. ");
0157       // create the BinUtility
0158       binUtility = std::make_unique<const BinUtility>(boundaries, open, aDir);
0159       ACTS_VERBOSE("arbitrary : created a BinUtility as " << *binUtility);
0160
0161     } break;
0162     // default return nullptr
0163     default: {
0164       return nullptr;
0165     }
0166   }
0167   // return the binned array
0168   return std::make_unique<const BinnedArrayXD<LayerPtr>>(layerOrderVector,
0169                                                          std::move(binUtility));
0170 }
0171
0172 std::shared_ptr<Acts::Surface> Acts::LayerArrayCreator::createNavigationSurface(
0173     const GeometryContext& gctx, const Layer& layer, AxisDirection aDir,
0174     double offset) const {
0175   // surface reference
0176   const Surface& layerSurface = layer.surfaceRepresentation();
0177   // translation to be applied
0178   Vector3 translation(0., 0., 0.);
0179   // switching he binnig values
0180   switch (aDir) {
0181     // case x
0182     case AxisDirection::AxisX: {
0183       translation = Vector3(offset, 0., 0.);
0184     } break;
0185     // case y
0186     case AxisDirection::AxisY: {
0187       translation = Vector3(0., offset, 0.);
0188     } break;
0189     // case z
0190     case AxisDirection::AxisZ: {
0191       translation = Vector3(0., 0., offset);
0192     } break;
0193     // case R
0194     case AxisDirection::AxisR: {
0195       // binning in R and cylinder surface means something different
0196       if (layerSurface.type() == Surface::Cylinder) {
0197         break;
0198       }
0199       translation = Vector3(offset, 0., 0.);
0200     } break;
0201     // do nothing for the default
0202     default: {
0203       ACTS_WARNING("Not yet implemented.");
0204     }
0205   }
0206   // navigation surface
0207   std::shared_ptr<Surface> navigationSurface;
0208   // for everything else than a cylinder it's a copy with shift
0209   if (layerSurface.type() == Surface::Plane) {
0210     // create a transform that does the shift
0211     Transform3 shift = Transform3(Translation3(translation));
0212     const PlaneSurface* plane =
0213         dynamic_cast<const PlaneSurface*>(&layerSurface);
0214     navigationSurface = Surface::makeShared<PlaneSurface>(gctx, *plane, shift);
0215   } else if (layerSurface.type() == Surface::Disc) {
0216     // create a transform that does the shift
0217     Transform3 shift = Transform3(Translation3(translation));
0218     const DiscSurface* disc = dynamic_cast<const DiscSurface*>(&layerSurface);
0219     navigationSurface = Surface::makeShared<DiscSurface>(gctx, *disc, shift);
0220   } else if (layerSurface.type() == Surface::Cylinder) {
0221     // get the bounds
0222     const CylinderBounds* cBounds =
0223         dynamic_cast<const CylinderBounds*>(&(layerSurface.bounds()));
0224     double navigationR = cBounds->get(CylinderBounds::eR) + offset;
0225     double halflengthZ = cBounds->get(CylinderBounds::eHalfLengthZ);
0226     // new navigation layer
0227     auto cylinderBounds =
0228         std::make_shared<CylinderBounds>(navigationR, halflengthZ);
0229     navigationSurface = Surface::makeShared<CylinderSurface>(
0230         layerSurface.transform(gctx), cylinderBounds);
0231   } else {
0232     ACTS_WARNING("Not implemented.");
0233   }
0234   return navigationSurface;
0235 }