Podio/src/PodioUtil.cpp

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