src/Digitization/Segmentizer.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 "ActsFatras/Digitization/Segmentizer.hpp"
0010
0011 #include "Acts/Surfaces/Surface.hpp"
0012 #include "Acts/Surfaces/detail/IntersectionHelper2D.hpp"
0013 #include "Acts/Utilities/BinUtility.hpp"
0014 #include "Acts/Utilities/BinningType.hpp"
0015 #include "Acts/Utilities/Helpers.hpp"
0016 #include "Acts/Utilities/Intersection.hpp"
0017
0018 #include <algorithm>
0019 #include <cmath>
0020 #include <memory>
0021
0022 std::vector<ActsFatras::Segmentizer::ChannelSegment>
0023 ActsFatras::Segmentizer::segments(const Acts::GeometryContext& geoCtx,
0024                                   const Acts::Surface& surface,
0025                                   const Acts::BinUtility& segmentation,
0026                                   const Segment2D& segment) const {
0027   // Return if the segmentation is not two-dimensional
0028   // (strips need to have one bin along the strip)
0029   if (segmentation.dimensions() != 2) {
0030     return {};
0031   }
0032
0033   // Start and end point
0034   const auto& start = segment[0];
0035   const auto& end = segment[1];
0036
0037   // Full path length - the full channel
0038   auto segment2d = (end - start);
0039   std::vector<ChannelStep> cSteps;
0040   Bin2D bstart = {0, 0};
0041   Bin2D bend = {0, 0};
0042
0043   if (surface.type() == Acts::Surface::SurfaceType::Plane) {
0044     // Get the segmentation and convert it to lines & arcs
0045     bstart = {static_cast<unsigned int>(segmentation.bin(start, 0)),
0046               static_cast<unsigned int>(segmentation.bin(start, 1))};
0047     bend = {static_cast<unsigned int>(segmentation.bin(end, 0)),
0048             static_cast<unsigned int>(segmentation.bin(end, 1))};
0049     // Fast single channel exit
0050     if (bstart == bend) {
0051       return {ChannelSegment(bstart, {start, end}, segment2d.norm())};
0052     }
0053     // The lines channel segment lines along x
0054     if (bstart[0] != bend[0]) {
0055       double k = segment2d.y() / segment2d.x();
0056       double d = start.y() - k * start.x();
0057
0058       const auto& xboundaries = segmentation.binningData()[0].boundaries();
0059       std::vector<double> xbbounds = {
0060           xboundaries.begin() + std::min(bstart[0], bend[0]) + 1,
0061           xboundaries.begin() + std::max(bstart[0], bend[0]) + 1};
0062       for (const auto x : xbbounds) {
0063         cSteps.push_back(ChannelStep{
0064             {(bstart[0] < bend[0] ? 1 : -1), 0}, {x, k * x + d}, start});
0065       }
0066     }
0067     // The lines channel segment lines along y
0068     if (bstart[1] != bend[1]) {
0069       double k = segment2d.x() / segment2d.y();
0070       double d = start.x() - k * start.y();
0071       const auto& yboundaries = segmentation.binningData()[1].boundaries();
0072       std::vector<double> ybbounds = {
0073           yboundaries.begin() + std::min(bstart[1], bend[1]) + 1,
0074           yboundaries.begin() + std::max(bstart[1], bend[1]) + 1};
0075       for (const auto y : ybbounds) {
0076         cSteps.push_back(ChannelStep{
0077             {0, (bstart[1] < bend[1] ? 1 : -1)}, {k * y + d, y}, start});
0078       }
0079     }
0080
0081   } else if (surface.type() == Acts::Surface::SurfaceType::Disc) {
0082     Acts::Vector2 pstart(Acts::VectorHelpers::perp(start),
0083                          Acts::VectorHelpers::phi(start));
0084     Acts::Vector2 pend(Acts::VectorHelpers::perp(end),
0085                        Acts::VectorHelpers::phi(end));
0086
0087     // Get the segmentation and convert it to lines & arcs
0088     bstart = {static_cast<unsigned int>(segmentation.bin(pstart, 0)),
0089               static_cast<unsigned int>(segmentation.bin(pstart, 1))};
0090     bend = {static_cast<unsigned int>(segmentation.bin(pend, 0)),
0091             static_cast<unsigned int>(segmentation.bin(pend, 1))};
0092
0093     // Fast single channel exit
0094     if (bstart == bend) {
0095       return {ChannelSegment(bstart, {start, end}, segment2d.norm())};
0096     }
0097
0098     double phistart = pstart[1];
0099     double phiend = pend[1];
0100
0101     // The radial boundaries
0102     if (bstart[0] != bend[0]) {
0103       const auto& rboundaries = segmentation.binningData()[0].boundaries();
0104       std::vector<double> rbbounds = {
0105           rboundaries.begin() + std::min(bstart[0], bend[0]) + 1,
0106           rboundaries.begin() + std::max(bstart[0], bend[0]) + 1};
0107       for (const auto& r : rbbounds) {
0108         auto radIntersection =
0109             Acts::detail::IntersectionHelper2D::intersectCircleSegment(
0110                 r, std::min(phistart, phiend), std::max(phistart, phiend),
0111                 start, (end - start).normalized());
0112         cSteps.push_back(ChannelStep{{(bstart[0] < bend[0] ? 1 : -1), 0},
0113                                      radIntersection.position(),
0114                                      start});
0115       }
0116     }
0117     // The phi boundaries
0118     if (bstart[1] != bend[1]) {
0119       double referenceR =
0120           surface.referencePositionValue(geoCtx, Acts::AxisDirection::AxisR);
0121       Acts::Vector2 origin = {0., 0.};
0122       const auto& phiboundaries = segmentation.binningData()[1].boundaries();
0123       std::vector<double> phibbounds = {
0124           phiboundaries.begin() + std::min(bstart[1], bend[1]) + 1,
0125           phiboundaries.begin() + std::max(bstart[1], bend[1]) + 1};
0126
0127       for (const auto& phi : phibbounds) {
0128         Acts::Vector2 philine(referenceR * std::cos(phi),
0129                               referenceR * std::sin(phi));
0130         auto phiIntersection =
0131             Acts::detail::IntersectionHelper2D::intersectSegment(
0132                 origin, philine, start, (end - start).normalized());
0133         cSteps.push_back(ChannelStep{{0, (bstart[1] < bend[1] ? 1 : -1)},
0134                                      phiIntersection.position(),
0135                                      start});
0136       }
0137     }
0138   }
0139
0140   // Register the last step if successful
0141   if (!cSteps.empty()) {
0142     cSteps.push_back(ChannelStep({0, 0}, end, start));
0143     std::ranges::sort(cSteps, std::less<ChannelStep>{});
0144   }
0145
0146   std::vector<ChannelSegment> cSegments;
0147   cSegments.reserve(cSteps.size());
0148
0149   Bin2D currentBin = {bstart[0], bstart[1]};
0150   BinDelta2D lastDelta = {0, 0};
0151   Acts::Vector2 lastIntersect = start;
0152   double lastPath = 0.;
0153   for (auto& cStep : cSteps) {
0154     currentBin[0] += lastDelta[0];
0155     currentBin[1] += lastDelta[1];
0156     double path = cStep.path - lastPath;
0157     cSegments.push_back(
0158         ChannelSegment(currentBin, {lastIntersect, cStep.intersect}, path));
0159     lastPath = cStep.path;
0160     lastDelta = cStep.delta;
0161     lastIntersect = cStep.intersect;
0162   }
0163
0164   return cSegments;
0165 }