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