detray/geometry/surface.hpp

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 #pragma once
0010
0011 // Project include(s)
0012 #include "detray/definitions/detail/qualifiers.hpp"
0013 #include "detray/definitions/geometry.hpp"
0014 #include "detray/definitions/indexing.hpp"
0015 #include "detray/definitions/math.hpp"
0016 #include "detray/geometry/detail/surface_kernels.hpp"
0017 #include "detray/geometry/identifier.hpp"
0018 #include "detray/material/material.hpp"
0019 #include "detray/utils/ranges/ranges.hpp"
0020
0021 // System include(s)
0022 #include <ostream>
0023 #include <type_traits>
0024
0025 namespace detray::geometry {
0026
0027 /// @brief Facade for a detray detector surface.
0028 ///
0029 /// Provides an interface to geometry specific functionality like
0030 /// local-to-global coordinate transforms or mask and material visitors. It
0031 /// wraps a detector instance that contains the data and a surface descriptor
0032 /// that contains the indices into the detector data containers for the
0033 /// specific surface instance.
0034 template <typename det_t>  // @TODO: This needs a concept
0035 class surface {
0036   /// Make sure the detector is always evaluated as constant type
0037   using detector_t = std::add_const_t<det_t>;
0038   /// Surface descriptor type
0039   using descr_t = typename detector_t::surface_type;
0040   /// Implementation
0041   using kernels = detail::surface_kernels<typename detector_t::algebra_type>;
0042
0043  public:
0044   using algebra_type = typename detector_t::algebra_type;
0045   using scalar_type = dscalar<algebra_type>;
0046   using point2_type = dpoint2D<algebra_type>;
0047   using point3_type = dpoint3D<algebra_type>;
0048   using vector3_type = dvector3D<algebra_type>;
0049   using transform3_type = dtransform3D<algebra_type>;
0050   using context = typename detector_t::geometry_context;
0051
0052   /// Not allowed: always needs a detector and a descriptor.
0053   surface() = delete;
0054
0055   /// Constructor from detector @param det and surface descriptor
0056   /// @param desc from that detector.
0057   DETRAY_HOST_DEVICE
0058   constexpr surface(const detector_t &det, const descr_t &desc)
0059       : m_detector{det}, m_desc{desc} {
0060     assert(!m_desc.identifier().is_invalid());
0061     assert(m_desc.index() < det.surfaces().size());
0062     assert(m_desc.transform() < det.transform_store().size());
0063   }
0064
0065   /// Constructor from detector @param det and identifier @param geo_id in
0066   /// that detector.
0067   DETRAY_HOST_DEVICE
0068   constexpr surface(const detector_t &det, const geometry::identifier geo_id)
0069       : surface(det, det.surface(geo_id)) {}
0070
0071   /// Constructor from detector @param det and surface index @param sf_idx
0072   DETRAY_HOST_DEVICE
0073   constexpr surface(const detector_t &det, const dindex sf_idx)
0074       : surface(det, det.surface(sf_idx)) {}
0075
0076   /// Allow conversion from surface<detector_t> to surface<const detector_t>,
0077   /// which has no semantic difference (the detector is always const internally)
0078   template <typename detector_type = detector_t>
0079     requires(!std::is_const_v<detector_type>)
0080   // NOLINTNEXTLINE
0081   DETRAY_HOST_DEVICE constexpr operator surface<const detector_type>() const {
0082     return surface<const detector_type>{this->m_detector, this->m_desc};
0083   }
0084
0085   /// Equality operator
0086   ///
0087   /// @param rhs is the right hand side to be compared to
0088   DETRAY_HOST_DEVICE
0089   constexpr auto operator==(const surface &rhs) const -> bool {
0090     return (&m_detector == &(rhs.m_detector) && m_desc == rhs.m_desc);
0091   }
0092
0093   /// @returns the surface identifier
0094   DETRAY_HOST_DEVICE
0095   constexpr auto identifier() const -> geometry::identifier {
0096     assert(!m_desc.identifier().is_invalid());
0097     return m_desc.identifier();
0098   }
0099
0100   /// @returns the index of the mother volume
0101   DETRAY_HOST_DEVICE
0102   constexpr auto volume() const -> dindex {
0103     assert(identifier().volume() < m_detector.volumes().size());
0104     return identifier().volume();
0105   }
0106
0107   /// @returns the index of the surface in the detector surface lookup
0108   DETRAY_HOST_DEVICE
0109   constexpr auto index() const -> dindex {
0110     assert(identifier().index() < m_detector.surfaces().size());
0111     return identifier().index();
0112   }
0113
0114   /// @returns the surface id (sensitive, passive or portal)
0115   DETRAY_HOST_DEVICE
0116   constexpr auto id() const -> surface_id {
0117     assert(identifier().id() != surface_id::e_unknown);
0118     return identifier().id();
0119   }
0120
0121   /// @returns the extra bits in the identifier
0122   DETRAY_HOST_DEVICE
0123   constexpr auto extra() const -> dindex { return identifier().extra(); }
0124
0125   /// @returns an id for the surface type (e.g. 'rectangle')
0126   DETRAY_HOST_DEVICE
0127   constexpr auto shape_id() const { return m_desc.mask().id(); }
0128
0129   /// @returns the surface source link
0130   DETRAY_HOST_DEVICE
0131   constexpr auto source() const {
0132     return m_detector.surface(identifier()).source;
0133   }
0134
0135   /// @returns true if the surface is a sensitive detector module.
0136   DETRAY_HOST_DEVICE
0137   constexpr auto is_sensitive() const -> bool {
0138     return identifier().id() == surface_id::e_sensitive;
0139   }
0140
0141   /// @returns true if the surface is a portal.
0142   DETRAY_HOST_DEVICE
0143   constexpr auto is_portal() const -> bool {
0144     return identifier().id() == surface_id::e_portal;
0145   }
0146
0147   /// @returns true if the surface is a passive detector element.
0148   DETRAY_HOST_DEVICE
0149   constexpr auto is_passive() const -> bool {
0150     return identifier().id() == surface_id::e_passive;
0151   }
0152
0153   /// @returns true if the surface carries material.
0154   DETRAY_HOST_DEVICE
0155   constexpr bool has_material() const { return m_desc.has_material(); }
0156
0157   /// @returns the mask volume link
0158   DETRAY_HOST_DEVICE
0159   constexpr auto volume_links() const {
0160     return visit_mask<typename kernels::get_volume_links>();
0161   }
0162
0163   /// @returns the mask shape name
0164   DETRAY_HOST
0165   std::string shape_name() const {
0166     return visit_mask<typename kernels::get_shape_name>();
0167   }
0168
0169   /// @returns the number of masks associated with this surface
0170   DETRAY_HOST_DEVICE
0171   std::size_t n_masks() const {
0172     if constexpr (concepts::interval<decltype(m_desc.mask().index())>) {
0173       return m_desc.mask().index().size();
0174     } else {
0175       return 1u;
0176     }
0177   }
0178
0179   /// @returns the coordinate transform matrix of the surface
0180   DETRAY_HOST_DEVICE
0181   constexpr auto transform(const context &ctx) const
0182       -> const transform3_type & {
0183     assert(m_desc.transform() < m_detector.transform_store().size());
0184     return m_detector.transform_store().at(m_desc.transform(), ctx);
0185   }
0186
0187   /// @returns the mask volume link
0188   template <concepts::point point_t = point2_type>
0189   DETRAY_HOST_DEVICE constexpr bool is_inside(const point_t &loc_p,
0190                                               const scalar_type tol) const {
0191     return visit_mask<typename kernels::is_inside>(loc_p, tol);
0192   }
0193
0194   /// @returns a boundary value of the surface, according to @param index
0195   DETRAY_HOST_DEVICE
0196   constexpr auto boundary(std::size_t index) const {
0197     return visit_mask<typename kernels::get_mask_value>(index);
0198   }
0199
0200   /// @returns the centroid of the surface mask in local cartesian coordinates
0201   DETRAY_HOST_DEVICE
0202   constexpr auto centroid() const -> point3_type {
0203     return visit_mask<typename kernels::centroid>();
0204   }
0205
0206   /// @returns the center position of the surface in global coordinates
0207   /// @note for shapes like the annulus this is not synonymous to the centroid
0208   /// but instead the focal point of the strip system
0209   DETRAY_HOST_DEVICE
0210   constexpr auto center(const context &ctx) const -> point3_type {
0211     return transform(ctx).translation();
0212   }
0213
0214   /// @returns the surface normal in global coordinates at a given bound/local
0215   /// position @param p
0216   template <concepts::point point_t = point2_type>
0217   DETRAY_HOST_DEVICE constexpr auto normal(const context &ctx,
0218                                            const point_t &p) const
0219       -> vector3_type {
0220     return visit_mask<typename kernels::normal>(transform(ctx), p);
0221   }
0222
0223   /// @returns a pointer to the material parameters at the local position
0224   /// @param loc_p
0225   DETRAY_HOST_DEVICE constexpr const material<scalar_type> *material_parameters(
0226       const point2_type &loc_p) const {
0227     return visit_material<typename kernels::get_material_params>(loc_p);
0228   }
0229
0230   /// @returns the local position to the global point @param global for
0231   /// a given geometry context @param ctx and track direction @param dir
0232   DETRAY_HOST_DEVICE
0233   constexpr point3_type global_to_local(const context &ctx,
0234                                         const point3_type &global,
0235                                         const vector3_type &dir) const {
0236     return visit_mask<typename kernels::global_to_local>(transform(ctx), global,
0237                                                          dir);
0238   }
0239
0240   /// @returns the bound (2D) position to the global point @param global for
0241   /// a given geometry context @param ctx and track direction @param dir
0242   DETRAY_HOST_DEVICE
0243   constexpr point2_type global_to_bound(const context &ctx,
0244                                         const point3_type &global,
0245                                         const vector3_type &dir) const {
0246     const point3_type local{global_to_local(ctx, global, dir)};
0247
0248     return {local[0], local[1]};
0249   }
0250
0251   /// @returns the global position to the given local position @param local
0252   /// for a given geometry context @param ctx
0253   template <concepts::point point_t = point2_type>
0254   DETRAY_HOST_DEVICE constexpr point3_type local_to_global(
0255       const context &ctx, const point_t &local, const vector3_type &dir) const {
0256     return visit_mask<typename kernels::local_to_global>(transform(ctx), local,
0257                                                          dir);
0258   }
0259
0260   /// Call a functor on the surfaces mask with additional arguments.
0261   ///
0262   /// @tparam functor_t the prescription to be applied to the mask
0263   /// @tparam Args      types of additional arguments to the functor
0264   template <typename functor_t, typename... Args>
0265   DETRAY_HOST_DEVICE constexpr auto visit_mask(Args &&...args) const {
0266     assert(!m_desc.mask().is_invalid());
0267     const auto &masks = m_detector.mask_store();
0268     return masks.template visit<functor_t>(m_desc.mask(),
0269                                            std::forward<Args>(args)...);
0270   }
0271
0272   /// Call a functor on the surfaces material with additional arguments.
0273   ///
0274   /// @tparam functor_t the prescription to be applied to the material
0275   /// @tparam Args      types of additional arguments to the functor
0276   template <typename functor_t, typename... Args>
0277   DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&...args) const {
0278     assert(has_material());
0279     const auto &materials = m_detector.material_store();
0280     return materials.template visit<functor_t>(m_desc.material(),
0281                                                std::forward<Args>(args)...);
0282   }
0283
0284   /// Do a consistency check on the surface after building the detector.
0285   ///
0286   /// @param os output stream for error messages.
0287   ///
0288   /// @returns true if the surface is consistent
0289   DETRAY_HOST bool self_check(std::ostream &os) const {
0290     if (identifier().is_invalid()) {
0291       os << "DETRAY ERROR (HOST): Invalid identifier for surface:\n"
0292          << *this << std::endl;
0293       return false;
0294     }
0295     if (index() >= m_detector.surfaces().size()) {
0296       os << "DETRAY ERROR (HOST): Surface index out of bounds for "
0297             "surface:\n"
0298          << *this << std::endl;
0299       return false;
0300     }
0301     if (volume() >= m_detector.volumes().size()) {
0302       os << "DETRAY ERROR (HOST): Surface volume index out of bounds for "
0303             "surface:\n"
0304          << *this << std::endl;
0305       return false;
0306     }
0307     if (detail::is_invalid_value(m_desc.transform())) {
0308       os << "DETRAY ERROR (HOST): Surface transform undefined for "
0309             "surface:\n"
0310          << *this << std::endl;
0311       return false;
0312     }
0313     if (m_desc.transform() >= m_detector.transform_store().size()) {
0314       os << "DETRAY ERROR (HOST): Surface transform index out of bounds "
0315             "for surface:\n"
0316          << *this << std::endl;
0317       return false;
0318     }
0319     if (detail::is_invalid_value(m_desc.mask())) {
0320       os << "DETRAY ERROR (HOST): Surface does not have a valid mask "
0321             "link:\n"
0322          << *this << std::endl;
0323       return false;
0324     }
0325     // Only check, if there is material in the detector
0326     if (!m_detector.material_store().all_empty() && has_material() &&
0327         m_desc.material().is_invalid_index()) {
0328       os << "DETRAY ERROR (HOST): Surface does not have valid material "
0329             "link:\n"
0330          << *this << std::endl;
0331       return false;
0332     }
0333     // Check the mask boundaries
0334     if (!visit_mask<typename kernels::mask_self_check>(os)) {
0335       os << "\nSurface: " << *this << std::endl;
0336       return false;
0337     }
0338
0339     // Check the mask volume link
0340     const auto vol_links = visit_mask<typename kernels::get_volume_links>();
0341     for (const auto vol_link : vol_links) {
0342       if (is_portal()) {
0343         if (vol_link == volume()) {
0344           os << "DETRAY ERROR (HOST): Portal surface links to mother "
0345                 "volume:\n"
0346              << *this << std::endl;
0347           return false;
0348         }
0349       } else if (vol_link != volume()) {
0350         os << "DETRAY ERROR (HOST): Passive/sensitive surface does not "
0351               "link to mother volume: Mask volume link : "
0352            << vol_link << "\n"
0353            << *this << std::endl;
0354         return false;
0355       }
0356     }
0357
0358     return true;
0359   }
0360
0361  protected:
0362   /// @returns a string stream that prints the surface details
0363   DETRAY_HOST
0364   friend std::ostream &operator<<(std::ostream &os, const surface &sf) {
0365     os << sf.m_desc;
0366     return os;
0367   }
0368
0369   /// @returns access to the underlying detector object
0370   DETRAY_HOST_DEVICE
0371   const detector_t &detector() const { return m_detector; }
0372
0373   /// @returns access to the underlying surface descriptor object
0374   DETRAY_HOST_DEVICE
0375   descr_t descriptor() const { return m_desc; }
0376
0377  private:
0378   /// Access to the detector stores
0379   const detector_t &m_detector;
0380   /// Access to the descriptor
0381   const descr_t m_desc;
0382 };
0383
0384 template <typename detector_t, typename descr_t>
0385 DETRAY_HOST_DEVICE surface(const detector_t &, const descr_t &)
0386     -> surface<detector_t>;
0387
0388 template <typename detector_t>
0389 DETRAY_HOST_DEVICE surface(const detector_t &, const geometry::identifier)
0390     -> surface<detector_t>;
0391
0392 }  // namespace detray::geometry