Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-27 07:24:12

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/builders/surface_factory_interface.hpp"
0013 #include "detray/core/detail/data_context.hpp"
0014 #include "detray/definitions/algebra.hpp"
0015 #include "detray/definitions/detail/qualifiers.hpp"
0016 #include "detray/definitions/geometry.hpp"
0017 #include "detray/definitions/indexing.hpp"
0018 #include "detray/geometry/mask.hpp"
0019 #include "detray/geometry/shapes/rectangle2D.hpp"
0020 #include "detray/material/material.hpp"
0021 #include "detray/utils/logging.hpp"
0022 
0023 // System include(s)
0024 #include <cassert>
0025 #include <iostream>
0026 #include <limits>
0027 
0028 namespace detray {
0029 
0030 /// @brief configuration for the barrel module generator
0031 ///
0032 /// The default values correspond to the toy detector inner pixel layer
0033 template <concepts::scalar scalar_t>
0034 struct endcap_generator_config {
0035   /// Construct positive or negative endcap
0036   int m_side{-1};
0037   /// Center position (|z|)
0038   scalar_t m_center_z{0.f};
0039   /// Inner layer radius
0040   scalar_t m_inner_radius{25.f * unit<scalar_t>::mm};
0041   /// Outer layer radius
0042   scalar_t m_outer_radius{180.f * unit<scalar_t>::mm};
0043   /// Boundary values for the module masks (per ring)
0044   std::vector<std::vector<scalar_t>> m_mask_values{
0045       {3.f * unit<scalar_t>::mm, 9.5f * unit<scalar_t>::mm,
0046        39.f * unit<scalar_t>::mm},
0047       {6.f * unit<scalar_t>::mm, 10.f * unit<scalar_t>::mm,
0048        39.f * unit<scalar_t>::mm}};
0049   /// Stagger between the two rings
0050   scalar_t m_ring_stagger{2.f * unit<scalar_t>::mm};
0051   /// Stagger in phi (per ring)
0052   std::vector<scalar_t> m_phi_stagger = {4.f * unit<scalar_t>::mm,
0053                                          4.f * unit<scalar_t>::mm};
0054   /// Substagger in phi (per ring)
0055   std::vector<scalar_t> m_phi_sub_stagger = {0.5f * unit<scalar_t>::mm,
0056                                              0.5f * unit<scalar_t>::mm};
0057   /// Module tilt (per ring)
0058   std::vector<scalar_t> m_tilt = {0.f, 0.f};
0059   /// Number of modules in phi (per ring)
0060   std::vector<unsigned int> m_binning = {40u, 68u};
0061 
0062   /// Setters
0063   /// @{
0064   constexpr endcap_generator_config &side(const int s) {
0065     m_side = s;
0066     return *this;
0067   }
0068   constexpr endcap_generator_config &center(const scalar_t z) {
0069     m_center_z = std::abs(z);
0070     return *this;
0071   }
0072   constexpr endcap_generator_config &inner_radius(const scalar_t r) {
0073     m_inner_radius = r;
0074     return *this;
0075   }
0076   constexpr endcap_generator_config &outer_radius(const scalar_t r) {
0077     m_outer_radius = r;
0078     return *this;
0079   }
0080   endcap_generator_config &module_bounds(
0081       const std::vector<std::vector<scalar_t>> &bnds) {
0082     m_mask_values.clear();
0083     m_mask_values.push_back(bnds.at(0));
0084     m_mask_values.push_back(bnds.at(1));
0085     return *this;
0086   }
0087   constexpr endcap_generator_config &ring_stagger(const scalar_t rs) {
0088     m_ring_stagger = rs;
0089     return *this;
0090   }
0091   constexpr endcap_generator_config &phi_stagger(
0092       const std::vector<scalar_t> &ps) {
0093     m_phi_stagger = ps;
0094     return *this;
0095   }
0096   constexpr endcap_generator_config &phi_sub_stagger(
0097       const std::vector<scalar_t> &ps) {
0098     m_phi_sub_stagger = ps;
0099     return *this;
0100   }
0101   constexpr endcap_generator_config &module_tilt(
0102       const std::vector<scalar_t> &t) {
0103     m_tilt = t;
0104     return *this;
0105   }
0106   constexpr endcap_generator_config &binning(
0107       const std::vector<unsigned int> &b) {
0108     m_binning = b;
0109     return *this;
0110   }
0111   /// @}
0112 
0113   /// Getters
0114   /// @{
0115   constexpr int side() const { return m_side; }
0116   constexpr scalar_t center() const { return m_center_z; }
0117   constexpr scalar_t inner_radius() const { return m_inner_radius; }
0118   constexpr scalar_t outer_radius() const { return m_outer_radius; }
0119   const auto &module_bounds() const { return m_mask_values; }
0120   constexpr scalar_t ring_stagger() const { return m_ring_stagger; }
0121   constexpr const auto &phi_stagger() const { return m_phi_stagger; }
0122   constexpr const auto &phi_sub_stagger() const { return m_phi_sub_stagger; }
0123   constexpr const auto &module_tilt() const { return m_tilt; }
0124   constexpr const auto &binning() const { return m_binning; }
0125   /// @}
0126 };
0127 
0128 /// @brief Generates surfaces for an endcap layer consisting of two rings
0129 ///
0130 /// @tparam detector_t the type of detector the layer should be added to
0131 template <typename detector_t, typename mask_shape_t = trapezoid2D>
0132 class endcap_generator final : public surface_factory_interface<detector_t> {
0133   using algebra_t = typename detector_t::algebra_type;
0134   using scalar_t = dscalar<algebra_t>;
0135   using transform3_t = dtransform3D<algebra_t>;
0136   using point3_t = dpoint3D<algebra_t>;
0137   using vector3_t = dvector3D<algebra_t>;
0138 
0139  public:
0140   /// Build an endcap layer according to the parameters given in @param cfg
0141   DETRAY_HOST
0142   explicit endcap_generator(const endcap_generator_config<scalar_t> &cfg)
0143       : m_cfg{cfg} {}
0144 
0145   /// @returns the number of surfaces this factory will produce
0146   DETRAY_HOST
0147   auto size() const -> dindex override {
0148     dindex n_modules{0u};
0149     for (const unsigned int n_bins : m_cfg.binning()) {
0150       n_modules += n_bins;
0151     }
0152     return n_modules;
0153   }
0154 
0155   /// This is a surface generator, no external surface data needed
0156   /// @{
0157   DETRAY_HOST
0158   void clear() override { /*Do nothing*/ };
0159   DETRAY_HOST
0160   void push_back(
0161       surface_data<detector_t> && /*unused*/) override { /*Do nothing*/ }
0162   DETRAY_HOST
0163   auto push_back(std::vector<surface_data<detector_t>> && /*unused*/)
0164       -> void override { /*Do nothing*/ }
0165   /// @}
0166 
0167   /// Create a pixel tracker endcap layer with two rings.
0168   ///
0169   /// @param volume the volume the portals need to be added to.
0170   /// @param surfaces the surface collection to wrap and to add the portals to
0171   /// @param transforms the transforms of the surfaces.
0172   /// @param masks the masks of the surfaces.
0173   /// @param ctx the geometry context (not needed for portals).
0174   DETRAY_HOST
0175   auto operator()(typename detector_t::volume_type &volume,
0176                   typename detector_t::surface_lookup_container &surfaces,
0177                   typename detector_t::transform_container &transforms,
0178                   typename detector_t::mask_container &masks,
0179                   typename detector_t::geometry_context ctx = {})
0180       -> dindex_range override {
0181     DETRAY_VERBOSE_HOST("Generate silicon tracker endcap modules...");
0182     DETRAY_VERBOSE_HOST("-> Generate " << size() << " surfaces");
0183 
0184     using surface_t = typename detector_t::surface_type;
0185     using nav_link_t = typename surface_t::navigation_link;
0186     using mask_link_t = typename surface_t::mask_link;
0187     using material_link_t = typename surface_t::material_link;
0188 
0189     auto volume_idx{volume.index()};
0190     const dindex surfaces_offset{static_cast<dindex>(surfaces.size())};
0191     constexpr auto invalid_src_link{detail::invalid_value<std::uint64_t>()};
0192 
0193     // The type id of the surface mask shape
0194     constexpr auto mask_id{
0195         types::id<typename detector_t::masks, mask<mask_shape_t, algebra_t>>};
0196     // The material will be added in a later step
0197     constexpr auto no_material{surface_t::material_id::e_none};
0198     // Modules link back to mother volume in navigation
0199     const auto mask_volume_link{static_cast<nav_link_t>(volume_idx)};
0200 
0201     // The radii of the rings
0202     std::vector<scalar_t> radii{};
0203     // The radial span of the disc
0204     scalar_t delta_r{m_cfg.outer_radius() - m_cfg.inner_radius()};
0205 
0206     //
0207     // Calculate the radial positions of the rings
0208     //
0209 
0210     // Only one ring
0211     if (m_cfg.binning().size() == 1u) {
0212       radii.push_back(0.5f * (m_cfg.inner_radius() + m_cfg.outer_radius()));
0213     } else {
0214       // Sum up the total length of the modules along r
0215       scalar_t tot_length{0.f};
0216       for (const auto &bounds : m_cfg.module_bounds()) {
0217         // Add an extra buffer, so that the trapezoid corners don't poke
0218         // out of the cylinder portals
0219         tot_length += 2.f * bounds.at(trapezoid2D::e_half_length_2) +
0220                       1.f * unit<scalar_t>::mm;
0221       }
0222 
0223       // Calculate the overlap (equal pay)
0224       scalar_t r_overlap{(tot_length - delta_r) /
0225                          static_cast<scalar_t>(m_cfg.module_bounds().size())};
0226 
0227       // Fill the radii and gaps
0228       scalar_t prev_r{m_cfg.inner_radius() + r_overlap};
0229       scalar_t prev_hl{0.f};
0230       scalar_t prev_ol{0.f};
0231 
0232       for (const auto &bounds : m_cfg.module_bounds()) {
0233         const scalar_t mod_hlength{bounds.at(trapezoid2D::e_half_length_2)};
0234         // Calculate the radius
0235         radii.push_back(prev_r + prev_hl - prev_ol + mod_hlength);
0236         prev_r = radii.back();
0237         prev_ol = 2.f * r_overlap;
0238         prev_hl = mod_hlength;
0239       }
0240     }
0241 
0242     //
0243     // Build the modules in every ring
0244     //
0245     for (unsigned int ir = 0u; ir < radii.size(); ++ir) {
0246       // Generate the position in z (observe ring stagger)
0247       // Convention: inner ring is closer to origin
0248       const scalar_t center{m_cfg.center()};
0249       const scalar_t staggered_center{
0250           (ir % 2u != 0 ? center + 0.5f * m_cfg.ring_stagger()
0251                         : center - 0.5f * m_cfg.ring_stagger())};
0252       const scalar_t rz{radii.size() == 1u ? center : staggered_center};
0253 
0254       // Generate the ring module positions (observe phi stagger)
0255       const scalar_t sub_stagger{m_cfg.phi_sub_stagger().size() > 1u
0256                                      ? m_cfg.phi_sub_stagger().at(ir)
0257                                      : 0.f};
0258 
0259       std::vector<point3_t> module_positions =
0260           module_positions_ring(rz, radii.at(ir), m_cfg.phi_stagger().at(ir),
0261                                 sub_stagger, m_cfg.binning().at(ir));
0262 
0263       // Build the modules
0264       for (const point3_t &mod_position : module_positions) {
0265         // Module mask
0266         mask_link_t mask_link{mask_id, {masks.template size<mask_id>(), 1u}};
0267         material_link_t material_link{no_material, dindex_invalid};
0268 
0269         // Surface descriptor
0270         surfaces.push_back({transforms.size(ctx), mask_link, material_link,
0271                             volume_idx, surface_id::e_sensitive},
0272                            invalid_src_link);
0273 
0274         // Build the module transform
0275 
0276         // The center position of the module
0277         point3_t mod_center{mod_position};
0278         mod_center[2] *= static_cast<scalar_t>(m_cfg.side());
0279         // The rotation matrix of the module
0280         const scalar_t mod_phi{vector::phi(mod_position)};
0281         const vector3_t mod_loc_y{math::cos(mod_phi), math::sin(mod_phi),
0282                                   static_cast<scalar_t>(0)};
0283         // Take different axis to have the same readout direction
0284         const vector3_t mod_loc_z{static_cast<scalar_t>(0),
0285                                   static_cast<scalar_t>(0),
0286                                   static_cast<scalar_t>(m_cfg.side())};
0287         const vector3_t mod_loc_x{vector::cross(mod_loc_y, mod_loc_z)};
0288 
0289         // Create the module transform object
0290         transforms.emplace_back(ctx, mod_center, mod_loc_z, mod_loc_x);
0291       }
0292 
0293       // Build the mask for this ring
0294       std::vector<scalar_t> mask_values{m_cfg.module_bounds().at(ir)};
0295 
0296       // Precompute trapezoid divisor
0297       if constexpr (std::is_same_v<mask_shape_t, trapezoid2D>) {
0298         const scalar_t div{
0299             1.f / (2.f * mask_values.at(trapezoid2D::e_half_length_2))};
0300 
0301         mask_values.insert(mask_values.begin() + trapezoid2D::e_divisor, div);
0302       }
0303 
0304       masks.template emplace_back<mask_id>(empty_context{}, mask_values,
0305                                            mask_volume_link);
0306     }
0307 
0308     return {surfaces_offset, static_cast<dindex>(surfaces.size())};
0309   }
0310 
0311  private:
0312   /// Helper method for positioning of modules in an endcap ring
0313   ///
0314   /// @param z is the z position of the ring
0315   /// @param radius is the ring radius
0316   /// @param phi_stagger is the radial staggering along phi
0317   /// @param phi_sub_stagger is the overlap of the modules
0318   /// @param n_phi_bins is the number of bins in phi
0319   ///
0320   /// @return a vector of the module positions in a ring
0321   auto module_positions_ring(const scalar_t z, const scalar_t radius,
0322                              const scalar_t phi_stagger,
0323                              const scalar_t phi_sub_stagger,
0324                              const unsigned int n_phi_bins) const {
0325     // create and fill the positions
0326     std::vector<point3_t> mod_positions;
0327     mod_positions.reserve(n_phi_bins);
0328 
0329     // prep work
0330     const scalar_t phi_step{2.f * constant<scalar_t>::pi /
0331                             static_cast<scalar_t>(n_phi_bins)};
0332     const scalar_t min_phi{-constant<scalar_t>::pi + 0.5f * phi_step};
0333 
0334     for (unsigned int iphi = 0u; iphi < n_phi_bins; ++iphi) {
0335       // If we have a phi sub stagger presents
0336       scalar_t rzs{0.f};
0337       // Phi stagger affects 0 vs 1, 2 vs 3 ... etc
0338       // -> only works if it is a %4
0339       // Phi sub stagger affects 2 vs 4, 1 vs 3 etc.
0340       if (phi_sub_stagger != 0.f && (n_phi_bins % 4u == 0u)) {
0341         // switch sides
0342         if ((iphi % 4u) == 0u) {
0343           rzs = phi_sub_stagger;
0344         } else if (((iphi + 1u) % 4u) == 0u) {
0345           rzs = -phi_sub_stagger;
0346         }
0347       }
0348       // The module phi
0349       const scalar_t phi{min_phi + static_cast<scalar_t>(iphi) * phi_step};
0350       // Main z position depending on phi bin
0351       const scalar_t rz{iphi % 2u != 0u ? z - 0.5f * phi_stagger
0352                                         : z + 0.5f * phi_stagger};
0353       mod_positions.push_back(
0354           {radius * math::cos(phi), radius * math::sin(phi), rz + rzs});
0355     }
0356     return mod_positions;
0357   }
0358 
0359   /// The generator configuration
0360   endcap_generator_config<scalar_t> m_cfg{};
0361 };
0362 
0363 }  // namespace detray