src/detail/GeoUnionDoubleTrdConverter.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 "ActsPlugins/GeoModel/detail/GeoUnionDoubleTrdConverter.hpp"
0010
0011 #include "Acts/Surfaces/PlaneSurface.hpp"
0012 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0013 #include "Acts/Utilities/TransformHelpers.hpp"
0014 #include "ActsPlugins/GeoModel/GeoModelConversionError.hpp"
0015 #include "ActsPlugins/GeoModel/detail/GeoShiftConverter.hpp"
0016
0017 #include <GeoModelHelpers/GeoShapeUtils.h>
0018
0019 namespace {
0020
0021 double distanceLinePoint(const Acts::Vector3 &lineA, const Acts::Vector3 &lineB,
0022                          const Acts::Vector3 &p) {
0023   auto dir = lineB - lineA;
0024   auto ap = p - lineA;
0025   return ap.cross(dir).norm() / dir.norm();
0026 }
0027
0028 /// Checks with the following properties if the trapezoids are mergeable
0029 bool trapezoidsAreMergeable(const std::vector<Acts::Vector3> &vtxsa,
0030                             const std::vector<Acts::Vector3> &vtxsb) {
0031   // Compute the distance of the lines connecting A3 and B0 and the midpoint of
0032   // the gap (resp. for other trapezoid side) These should be close to zero,
0033   // otherwise we cannot merge the trapezoids
0034   auto P1 = vtxsa[0] + 0.5 * (vtxsb[3] - vtxsa[0]);
0035   auto dist1 = distanceLinePoint(vtxsa[3], vtxsb[0], P1);
0036
0037   auto P2 = vtxsa[1] + 0.5 * (vtxsb[2] - vtxsa[1]);
0038   auto dist2 = distanceLinePoint(vtxsa[2], vtxsb[1], P2);
0039
0040   if (dist1 > 1.e-3 || dist2 > 1.e-3) {
0041     return false;
0042   }
0043
0044   // We could verify other properties, but this seems sufficient for know
0045   return true;
0046 }
0047 /// Extract the shifted trapezoidal vertices
0048 /// @param bounds: Surface bounds providing the 2D vertices
0049 /// @param shiftTrf: Shift to be applied on each vertex
0050 std::vector<Acts::Vector3> extactVertices(const Acts::TrapezoidBounds &bounds,
0051                                           const Acts::Transform3 &shiftTrf) {
0052   std::vector<Acts::Vector3> vertices{};
0053   std::ranges::transform(bounds.vertices(0), std::back_inserter(vertices),
0054                          [&shiftTrf](const Acts::Vector2 &vertex) {
0055                            return shiftTrf *
0056                                   Acts::Vector3{vertex.x(), vertex.y(), 0.};
0057                          });
0058   return vertices;
0059 }
0060 }  // namespace
0061
0062 using namespace Acts;
0063
0064 namespace ActsPlugins::detail {
0065
0066 Result<GeoModelSensitiveSurface> GeoUnionDoubleTrdConverter::operator()(
0067     const PVConstLink &geoPV, const GeoShapeUnion &geoUnion,
0068     const Transform3 &absTransform, SurfaceBoundFactory &boundFactory,
0069     bool sensitive) const {
0070   const auto shiftA = dynamic_pointer_cast<const GeoShapeShift>(
0071       compressShift(geoUnion.getOpA()));
0072   const auto shiftB = dynamic_pointer_cast<const GeoShapeShift>(
0073       compressShift(geoUnion.getOpB()));
0074
0075   if (shiftA == nullptr || shiftB == nullptr) {
0076     return GeoModelConversionError::WrongShapeForConverter;
0077   }
0078
0079   auto shiftARes = detail::GeoShiftConverter{}(geoPV, *shiftA, absTransform,
0080                                                boundFactory, sensitive);
0081   if (!shiftARes.ok()) {
0082     return shiftARes.error();
0083   }
0084   auto shiftBRes = detail::GeoShiftConverter{}(geoPV, *shiftB, absTransform,
0085                                                boundFactory, sensitive);
0086   if (!shiftBRes.ok()) {
0087     return shiftBRes.error();
0088   }
0089
0090   const auto &[elA, surfaceA] = shiftARes.value();
0091   const auto &[elB, surfaceB] = shiftBRes.value();
0092
0093   if (!(surfaceA && surfaceB)) {
0094     return GeoModelConversionError::WrongShapeForConverter;
0095   }
0096
0097   if (surfaceA->bounds().type() != Acts::SurfaceBounds::eTrapezoid ||
0098       surfaceB->bounds().type() != Acts::SurfaceBounds::eTrapezoid) {
0099     return GeoModelConversionError::WrongShapeForConverter;
0100   }
0101
0102   // At this point we know, that we have two trapezoids
0103   // We assume the following Layout for this:
0104   //
0105   //  0 ______________________ 1
0106   //    \                    /
0107   //     \        B         /
0108   //      \________________/
0109   //     3                 2
0110   //      0 _____________ 1
0111   //        \           /
0112   //         \    A    /
0113   //          \_______/
0114   //          3       2
0115
0116   // Compute the gap between trapezoids
0117   const auto &boundsA =
0118       static_cast<const TrapezoidBounds &>(surfaceA->bounds());
0119   const auto &boundsB =
0120       static_cast<const TrapezoidBounds &>(surfaceB->bounds());
0121
0122   // First check now, if this actually is correct
0123   const auto vtxsa = extactVertices(boundsA, makeTransform3(shiftA->getX()));
0124   const auto vtxsb = extactVertices(boundsB, makeTransform3(shiftB->getX()));
0125
0126   if (!trapezoidsAreMergeable(vtxsa, vtxsb)) {
0127     return GeoModelConversionError::WrongShapeForConverter;
0128   }
0129
0130   // Compute half length y as the distance between the middle points of the
0131   // outer lines of the trapezoids
0132   Acts::Vector3 mpA = vtxsa[3] + 0.5 * (vtxsa[2] - vtxsa[3]);
0133   Acts::Vector3 mpB = vtxsb[0] + 0.5 * (vtxsb[1] - vtxsb[0]);
0134
0135   auto halfLengthY = 0.5 * (mpB - mpA).norm();
0136
0137   const auto gap =
0138       halfLengthY - (boundsA.values()[TrapezoidBounds::eHalfLengthY] +
0139                      boundsB.values()[TrapezoidBounds::eHalfLengthY]);
0140
0141   if (gap > gapTolerance) {
0142     return GeoModelConversionError::WrongShapeForConverter;
0143   }
0144
0145   // For the x half lengths, we can take the values from the 2 trapezoids
0146   auto hlxny = boundsA.values()[TrapezoidBounds::eHalfLengthXposY];
0147   auto hlxpy = boundsB.values()[TrapezoidBounds::eHalfLengthXnegY];
0148
0149   auto trapezoidBounds =
0150       boundFactory.makeBounds<TrapezoidBounds>(hlxpy, hlxny, halfLengthY);
0151   // Create transform from the transform of surfaceA and translate it in y
0152   // direction using the half length
0153   auto transform = absTransform * makeTransform3(shiftA->getX());
0154   transform.translate(Vector3{
0155       0.f, boundsA.values()[TrapezoidBounds::eHalfLengthY] - halfLengthY, 0.f});
0156
0157   // TODO extract this code snipped since it is copied quite a lot
0158   if (!sensitive) {
0159     auto surface =
0160         Surface::makeShared<PlaneSurface>(transform, trapezoidBounds);
0161     return std::make_tuple(nullptr, surface);
0162   }
0163
0164   // Create the element and the surface (we assume both have equal thickness)
0165   auto detectorElement =
0166       GeoModelDetectorElement::createDetectorElement<PlaneSurface>(
0167           geoPV, trapezoidBounds, transform, elA->thickness());
0168   auto surface = detectorElement->surface().getSharedPtr();
0169
0170   return std::make_tuple(detectorElement, surface);
0171 }
0172
0173 }  // namespace ActsPlugins::detail