Detector/detail/BlueprintHelper.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/Detector/detail/BlueprintHelper.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Geometry/VolumeBounds.hpp"
0014
0015 #include <algorithm>
0016 #include <array>
0017 #include <ranges>
0018
0019 namespace {
0020
0021 std::array<Acts::Vector3, 2u> endPointsXYZ(
0022     const Acts::Experimental::Gen2Blueprint::Node& node,
0023     Acts::AxisDirection bVal) {
0024   unsigned int bIdx = 0;
0025   switch (bVal) {
0026     case Acts::AxisDirection::AxisX:
0027       bIdx = 0;
0028       break;
0029     case Acts::AxisDirection::AxisY:
0030       bIdx = 1;
0031       break;
0032     case Acts::AxisDirection::AxisZ:
0033       bIdx = 2;
0034       break;
0035     default:
0036       break;
0037   }
0038   Acts::Vector3 axis = node.transform.rotation().col(bIdx);
0039   auto halfL = node.boundaryValues[bIdx];
0040   Acts::Vector3 center = node.transform.translation();
0041   Acts::Vector3 p0 = center - halfL * axis;
0042   Acts::Vector3 p1 = center + halfL * axis;
0043   return {p0, p1};
0044 }
0045
0046 }  // namespace
0047
0048 void Acts::Experimental::detail::BlueprintHelper::sort(
0049     Gen2Blueprint::Node& node, bool recursive) {
0050   if (node.children.size() < 2u) {
0051     return;
0052   }
0053   // Sort along x, y, z
0054   if (node.binning.size() == 1) {
0055     auto bVal = node.binning.front();
0056     // x,y,z binning along the axis
0057     if (bVal == AxisDirection::AxisX || bVal == AxisDirection::AxisY ||
0058         bVal == AxisDirection::AxisZ) {
0059       Vector3 nodeCenter = node.transform.translation();
0060       Vector3 nodeSortAxis = node.transform.rotation().col(toUnderlying(bVal));
0061       std::ranges::sort(node.children, {}, [&](const auto& c) {
0062         return (c->transform.translation() - nodeCenter).dot(nodeSortAxis);
0063       });
0064     } else if (bVal == AxisDirection::AxisR &&
0065                node.boundsType == VolumeBounds::eCylinder) {
0066       std::ranges::sort(node.children, {}, [](const auto& c) {
0067         return c->boundaryValues[0] + c->boundaryValues[1];
0068       });
0069     }
0070   }
0071
0072   // Sort the children
0073   if (recursive) {
0074     for (auto& child : node.children) {
0075       sort(*child, true);
0076     }
0077   }
0078 }
0079
0080 void Acts::Experimental::detail::BlueprintHelper::fillGaps(
0081     Gen2Blueprint::Node& node, bool adjustToParent) {
0082   // Return if this is a leaf node
0083   if (node.isLeaf()) {
0084     return;
0085   }
0086   if (node.boundsType == VolumeBounds::eCylinder && node.binning.size() == 1) {
0087     fillGapsCylindrical(node, adjustToParent);
0088   } else if (node.boundsType == VolumeBounds::eCuboid &&
0089              node.binning.size() == 1) {
0090     // Doesn't look like NOT adjusting to parent
0091     // makes sense. The gaps are not going
0092     // to be filled in non-binned directions
0093     fillGapsCuboidal(node, adjustToParent);
0094   } else {
0095     throw std::runtime_error(
0096         "BlueprintHelper: gap filling is not implemented for "
0097         "this boundary type");
0098   }
0099 }
0100
0101 void Acts::Experimental::detail::BlueprintHelper::fillGapsCylindrical(
0102     Gen2Blueprint::Node& node, bool adjustToParent) {
0103   // Nodes must be sorted
0104   sort(node, false);
0105
0106   // Container values
0107   auto cInnerR = node.boundaryValues[0];
0108   auto cOuterR = node.boundaryValues[1];
0109   auto cHalfZ = node.boundaryValues[2];
0110
0111   std::vector<std::unique_ptr<Gen2Blueprint::Node>> gaps;
0112   // Only 1D binning implemented for the moment
0113   if (AxisDirection bVal = node.binning.front(); bVal == AxisDirection::AxisZ) {
0114     // adjust inner/outer radius
0115     if (adjustToParent) {
0116       std::ranges::for_each(node.children, [&](auto& child) {
0117         child->boundaryValues[0] = cInnerR;
0118         child->boundaryValues[1] = cOuterR;
0119       });
0120     }
0121     auto [negC, posC] = endPointsXYZ(node, bVal);
0122     // Assume sorted along the local z axis
0123     unsigned int igap = 0;
0124     for (auto& child : node.children) {
0125       auto [neg, pos] = endPointsXYZ(*child, bVal);
0126       double gapSpan = (neg - negC).norm();
0127       if (gapSpan > s_onSurfaceTolerance) {
0128         // Fill a gap node
0129         auto gapName = node.name + "_gap_" + std::to_string(igap);
0130         auto gapTransform = Transform3::Identity();
0131         gapTransform.rotate(node.transform.rotation());
0132         gapTransform.pretranslate(0.5 * (neg + negC));
0133         auto gap = std::make_unique<Gen2Blueprint::Node>(
0134             gapName, gapTransform, VolumeBounds::eCylinder,
0135             std::vector<double>{cInnerR, cOuterR, 0.5 * gapSpan});
0136         gaps.push_back(std::move(gap));
0137         ++igap;
0138       }
0139       // Set to new current negative value
0140       negC = pos;
0141     }
0142     // Check if a last one needs to be filled
0143     double gapSpan = (negC - posC).norm();
0144     if (gapSpan > s_onSurfaceTolerance) {
0145       // Fill a gap node
0146       auto gapName = node.name + "_gap_" + std::to_string(igap);
0147       auto gapTransform = Transform3::Identity();
0148       gapTransform.rotate(node.transform.rotation());
0149       gapTransform.pretranslate(0.5 * (negC + posC));
0150       auto gap = std::make_unique<Gen2Blueprint::Node>(
0151           gapName, gapTransform, VolumeBounds::eCylinder,
0152           std::vector<double>{cInnerR, cOuterR, 0.5 * gapSpan});
0153       gaps.push_back(std::move(gap));
0154     }
0155
0156   } else if (bVal == AxisDirection::AxisR) {
0157     // We have binning in R present
0158     if (adjustToParent) {
0159       std::ranges::for_each(node.children, [&](auto& child) {
0160         child->transform = node.transform;
0161         child->boundaryValues[2] = cHalfZ;
0162       });
0163     }
0164     // Fill the gaps in R
0165     unsigned int igap = 0;
0166     double lastR = cInnerR;
0167     for (auto& child : node.children) {
0168       double iR = child->boundaryValues[0];
0169       if (std::abs(iR - lastR) > s_onSurfaceTolerance) {
0170         auto gap = std::make_unique<Gen2Blueprint::Node>(
0171             node.name + "_gap_" + std::to_string(igap), node.transform,
0172             VolumeBounds::eCylinder, std::vector<double>{lastR, iR, cHalfZ});
0173         gaps.push_back(std::move(gap));
0174         ++igap;
0175       }
0176       // Set to new outer radius
0177       lastR = child->boundaryValues[1];
0178     }
0179     // Check if a last one needs to be filled
0180     if (std::abs(lastR - cOuterR) > s_onSurfaceTolerance) {
0181       auto gap = std::make_unique<Gen2Blueprint::Node>(
0182           node.name + "_gap_" + std::to_string(igap), node.transform,
0183           VolumeBounds::eCylinder, std::vector<double>{lastR, cOuterR, cHalfZ});
0184       gaps.push_back(std::move(gap));
0185     }
0186   } else {
0187     throw std::runtime_error(
0188         "BlueprintHelper: gap filling not implemented for "
0189         "cylinder and this binning type.");
0190   }
0191
0192   // Insert
0193   for (auto& gap : gaps) {
0194     node.add(std::move(gap));
0195   }
0196
0197   // Sort again after inserting
0198   sort(node, false);
0199   // Fill the gaps recursively
0200   for (auto& child : node.children) {
0201     fillGaps(*child, adjustToParent);
0202   }
0203 }
0204
0205 void Acts::Experimental::detail::BlueprintHelper::fillGapsCuboidal(
0206     Gen2Blueprint::Node& node, bool adjustToParent) {
0207   // Nodes must be sorted
0208   sort(node, false);
0209
0210   // Cuboidal detector binnings
0211   std::array<Acts::AxisDirection, 3u> allowedBinVals = {
0212       AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ};
0213
0214   std::vector<std::unique_ptr<Gen2Blueprint::Node>> gaps;
0215   auto binVal = node.binning.front();
0216
0217   // adjust non-binned directions
0218   if (adjustToParent) {
0219     std::ranges::for_each(node.children, [&](auto& child) {
0220       for (auto bv : allowedBinVals) {
0221         if (bv != binVal) {
0222           // Both boundary values and translation
0223           // have to be adjusted
0224           child->boundaryValues[toUnderlying(bv)] =
0225               node.boundaryValues[toUnderlying(bv)];
0226           child->transform.translation()[toUnderlying(bv)] =
0227               node.transform.translation()[toUnderlying(bv)];
0228         }
0229       }
0230     });
0231   }
0232   auto [negC, posC] = endPointsXYZ(node, binVal);
0233
0234   // Assume sorted along the local binned axis
0235   unsigned int igap = 0;
0236   for (auto& child : node.children) {
0237     auto [neg, pos] = endPointsXYZ(*child, binVal);
0238     double gapSpan = (neg - negC).norm();
0239     if (gapSpan > s_onSurfaceTolerance) {
0240       // Fill a gap node
0241       auto gapName = node.name + "_gap_" + std::to_string(igap);
0242       auto gapTransform = Transform3::Identity();
0243       gapTransform.rotate(node.transform.rotation());
0244       gapTransform.pretranslate(0.5 * (neg + negC));
0245       std::vector<double> gapBounds{0, 0, 0};
0246       gapBounds[toUnderlying(binVal)] = 0.5 * gapSpan;
0247       for (auto bv : allowedBinVals) {
0248         if (bv != binVal) {
0249           gapBounds[toUnderlying(bv)] = node.boundaryValues[toUnderlying(bv)];
0250         }
0251       }
0252       auto gap = std::make_unique<Gen2Blueprint::Node>(
0253           gapName, gapTransform, VolumeBounds::eCuboid, gapBounds);
0254       gaps.push_back(std::move(gap));
0255       ++igap;
0256     }
0257     // Set to new current negative value
0258     negC = pos;
0259   }
0260   // Check if a last one needs to be filled
0261   double gapSpan = (negC - posC).norm();
0262   if (gapSpan > s_onSurfaceTolerance) {
0263     // Fill a gap node
0264     auto gapName = node.name + "_gap_" + std::to_string(igap);
0265     auto gapTransform = Transform3::Identity();
0266     gapTransform.rotate(node.transform.rotation());
0267     gapTransform.pretranslate(0.5 * (negC + posC));
0268     std::vector<double> gapBounds{0, 0, 0};
0269     gapBounds[toUnderlying(binVal)] = 0.5 * gapSpan;
0270     for (auto bv : allowedBinVals) {
0271       if (bv != binVal) {
0272         gapBounds[toUnderlying(bv)] = node.boundaryValues[toUnderlying(bv)];
0273       }
0274     }
0275     auto gap = std::make_unique<Gen2Blueprint::Node>(
0276         gapName, gapTransform, VolumeBounds::eCuboid, gapBounds);
0277     gaps.push_back(std::move(gap));
0278   }
0279
0280   // Insert
0281   for (auto& gap : gaps) {
0282     node.add(std::move(gap));
0283   }
0284
0285   // Sort again after inserting
0286   sort(node, false);
0287   // Fill the gaps recursively
0288   for (auto& child : node.children) {
0289     fillGaps(*child, adjustToParent);
0290   }
0291 }