Podio/src/PodioUtil.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/Podio//PodioUtil.hpp"
0010
0011 #include "Acts/Geometry/GeometryContext.hpp"
0012 #include "Acts/Geometry/GeometryIdentifier.hpp"
0013 #include "Acts/Surfaces/AnnulusBounds.hpp"
0014 #include "Acts/Surfaces/ConeSurface.hpp"
0015 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
0016 #include "Acts/Surfaces/CylinderBounds.hpp"
0017 #include "Acts/Surfaces/CylinderSurface.hpp"
0018 #include "Acts/Surfaces/DiamondBounds.hpp"
0019 #include "Acts/Surfaces/DiscBounds.hpp"
0020 #include "Acts/Surfaces/DiscSurface.hpp"
0021 #include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
0022 #include "Acts/Surfaces/EllipseBounds.hpp"
0023 #include "Acts/Surfaces/PerigeeSurface.hpp"
0024 #include "Acts/Surfaces/PlanarBounds.hpp"
0025 #include "Acts/Surfaces/PlaneSurface.hpp"
0026 #include "Acts/Surfaces/RadialBounds.hpp"
0027 #include "Acts/Surfaces/StrawSurface.hpp"
0028 #include "Acts/Surfaces/Surface.hpp"
0029 #include "Acts/Utilities/ThrowAssert.hpp"
0030 #include "Acts/Utilities/TypeList.hpp"
0031 #include "ActsPodioEdm/Surface.h"
0032
0033 #include <limits>
0034 #include <memory>
0035
0036 using namespace Acts;
0037
0038 namespace ActsPlugins {
0039 namespace PodioUtil {
0040
0041 namespace {
0042 template <typename bounds_t>
0043 std::shared_ptr<const bounds_t> createBounds(
0044     const ActsPodioEdm::Surface& surface) {
0045   constexpr std::size_t S = bounds_t::eSize;
0046   throw_assert(surface.boundValuesSize == S,
0047                "Unexpected number of bound values");
0048
0049   std::array<double, S> values{};
0050   for (std::size_t i = 0; i < S; i++) {
0051     values.at(i) = surface.boundValues.at(i);
0052   }
0053   return std::make_shared<bounds_t>(values);
0054 }
0055 }  // namespace
0056
0057 ActsPodioEdm::Surface convertSurfaceToPodio(const ConversionHelper& helper,
0058                                             const Surface& surface) {
0059   ActsPodioEdm::Surface result;
0060
0061   std::optional<Identifier> identifier = helper.surfaceToIdentifier(surface);
0062   if (identifier.has_value()) {
0063     result.identifier = identifier.value();
0064   } else {
0065     result.identifier = kNoIdentifier;
0066     assert(surface.associatedDetectorElement() == nullptr &&
0067            "Unidentified surface does not have detector element");
0068     // @TODO: Surface type is not well-defined for curvilinear surface: looks like any plane surface
0069     result.surfaceType = surface.type();
0070     // @TODO: Test line bounds, does not have bounds, so nullptr
0071     result.boundsType = surface.bounds().type();
0072     result.geometryId = surface.geometryId().value();
0073     auto values = surface.bounds().values();
0074
0075     if (values.size() > result.boundValues.size()) {
0076       throw std::runtime_error{"Too many bound values to store"};
0077     }
0078
0079     for (std::size_t i = 0; i < values.size(); i++) {
0080       result.boundValues.at(i) = values.at(i);
0081     }
0082     result.boundValuesSize = values.size();
0083
0084     Eigen::Map<ActsSquareMatrix<4>> trf{result.transform.data()};
0085
0086     // This is safe ONLY(!) if there is no associated detector element, since
0087     // the surface will not inspect the geometry context at all by itself.
0088     GeometryContext gctx;
0089     trf = surface.transform(gctx).matrix();
0090   }
0091
0092   return result;
0093 }
0094
0095 std::shared_ptr<const Surface> convertSurfaceFromPodio(
0096     const ConversionHelper& helper, const ActsPodioEdm::Surface& surface) {
0097   if (surface.surfaceType == kNoSurface) {
0098     return nullptr;
0099   }
0100
0101   Eigen::Map<const ActsSquareMatrix<4>> mat{surface.transform.data()};
0102   Transform3 transform{mat};
0103
0104   using T = Surface::SurfaceType;
0105   using B = SurfaceBounds;
0106
0107   std::shared_ptr<const Surface> result;
0108
0109   if (const Surface* srf = helper.identifierToSurface(surface.identifier);
0110       srf != nullptr) {
0111     result = srf->getSharedPtr();
0112   }
0113
0114   if (result) {
0115     return result;
0116   }
0117
0118   switch (surface.surfaceType) {
0119     default:
0120       throw std::runtime_error{"Invalid surface type encountered"};
0121
0122     case T::Cone:
0123       throw_assert(surface.boundsType == B::eCone, "Unexpected bounds type");
0124       result = Surface::makeShared<ConeSurface>(
0125           transform, createBounds<ConeBounds>(surface));
0126       break;
0127
0128     case T::Cylinder:
0129       throw_assert(surface.boundsType == B::eCylinder,
0130                    "Unexpected bounds type");
0131       result = Surface::makeShared<CylinderSurface>(
0132           transform, createBounds<CylinderBounds>(surface));
0133       break;
0134
0135     case T::Disc: {
0136       std::shared_ptr<const DiscBounds> dBounds;
0137       switch (surface.boundsType) {
0138         default:
0139           throw std::runtime_error{"Invalid bounds type encountered"};
0140
0141         case B::eDisc:
0142           dBounds = createBounds<RadialBounds>(surface);
0143           break;
0144
0145         case B::eAnnulus:
0146           dBounds = createBounds<AnnulusBounds>(surface);
0147           break;
0148
0149         case B::eDiscTrapezoid:
0150           dBounds = createBounds<DiscTrapezoidBounds>(surface);
0151           break;
0152       }
0153       result = Surface::makeShared<DiscSurface>(transform, dBounds);
0154       break;
0155     }
0156
0157     case T::Perigee:
0158       throw_assert(surface.boundsType == B::eBoundless,
0159                    "Unexpected bounds type");
0160       result = Surface::makeShared<PerigeeSurface>(transform);
0161       break;
0162
0163     case T::Plane: {
0164       std::shared_ptr<const PlanarBounds> pBounds;
0165       switch (surface.boundsType) {
0166         default:
0167           throw std::runtime_error{"Invalid bounds type encountered"};
0168
0169         case B::eDiamond:
0170           pBounds = createBounds<DiamondBounds>(surface);
0171           break;
0172         case B::eEllipse:
0173           pBounds = createBounds<EllipseBounds>(surface);
0174           break;
0175         case B::eRectangle:
0176           pBounds = createBounds<RectangleBounds>(surface);
0177           break;
0178         case B::eConvexPolygon:
0179           template_switch_lambda<6, 32>(surface.boundValuesSize, [&](auto N) {
0180             constexpr std::size_t nValues = decltype(N)::value;
0181             constexpr std::size_t nVertices = nValues / 2;
0182             pBounds = createBounds<ConvexPolygonBounds<nVertices>>(surface);
0183           });
0184           // @TODO: Maybe handle dynamic convex polygons?
0185           break;
0186       }
0187       assert(pBounds && "No PlanarBounds");
0188       result = Surface::makeShared<PlaneSurface>(transform, pBounds);
0189
0190       break;
0191     }
0192
0193     case T::Straw:
0194       throw_assert(surface.boundsType == B::eLine, "Unexpected bounds type");
0195       result = Surface::makeShared<StrawSurface>(
0196           transform, createBounds<LineBounds>(surface));
0197       break;
0198
0199     case T::Curvilinear:
0200       throw_assert(surface.boundsType == B::eBoundless,
0201                    "Unexpected bounds type");
0202       result = Surface::makeShared<PlaneSurface>(transform);
0203       break;
0204   }
0205
0206   return result;
0207 }
0208
0209 }  // namespace PodioUtil
0210 namespace podio_detail {
0211
0212 template <typename F, typename... Args>
0213 void apply(F&& f, TypeList<Args...> /*unused*/) {
0214   f(Args{}...);
0215 }
0216
0217 void recoverDynamicColumns(
0218     const podio::Frame& frame, const std::string& stem,
0219     std::unordered_map<HashedString,
0220                        std::unique_ptr<podio_detail::ConstDynamicColumnBase>>&
0221         dynamic) {
0222   // See
0223   // https://github.com/AIDASoft/podio/blob/858c0ff0b841705d1b18aafd57569fcbd1beda91/include/podio/UserDataCollection.h#L30-L31
0224   using types = TypeList<float, double, std::int8_t, std::int16_t, std::int32_t,
0225                          std::int64_t, std::uint8_t, std::uint16_t,
0226                          std::uint32_t, std::uint64_t>;
0227
0228   std::vector<std::string> available = frame.getAvailableCollections();
0229
0230   for (const auto& col : available) {
0231     std::string prefix = stem + "_extra__";
0232     std::size_t p = col.find(prefix);
0233     if (p == std::string::npos) {
0234       continue;
0235     }
0236     std::string dynName = col.substr(prefix.size());
0237     const podio::CollectionBase* coll = frame.get(col);
0238
0239     std::unique_ptr<podio_detail::ConstDynamicColumnBase> up;
0240
0241     apply(
0242         [&](auto... args) {
0243           auto inner = [&](auto arg) {
0244             if (up) {
0245               return;
0246             }
0247             using T = decltype(arg);
0248             const auto* dyn =
0249                 dynamic_cast<const podio::UserDataCollection<T>*>(coll);
0250             if (dyn == nullptr) {
0251               return;
0252             }
0253             up = std::make_unique<podio_detail::ConstDynamicColumn<T>>(dynName,
0254                                                                        *dyn);
0255           };
0256
0257           ((inner(args)), ...);
0258         },
0259         types{});
0260
0261     if (!up) {
0262       throw std::runtime_error{"Dynamic column '" + dynName +
0263                                "' is not of allowed type"};
0264     }
0265
0266     HashedString hashedKey = hashStringDynamic(dynName);
0267     dynamic.insert({hashedKey, std::move(up)});
0268   }
0269 }
0270
0271 }  // namespace podio_detail
0272 }  // namespace ActsPlugins