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