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