test/common/build_wire_chamber.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/builders/cylinder_portal_generator.hpp"
0013 #include "detray/builders/detector_builder.hpp"
0014 #include "detray/builders/grid_builder.hpp"
0015 #include "detray/builders/homogeneous_material_builder.hpp"
0016 #include "detray/builders/homogeneous_material_generator.hpp"
0017 #include "detray/core/detector.hpp"
0018 #include "detray/definitions/algebra.hpp"
0019 #include "detray/definitions/indexing.hpp"
0020 #include "detray/definitions/units.hpp"
0021 #include "detray/detectors/wire_chamber_metadata.hpp"
0022 #include "detray/geometry/shapes/line.hpp"
0023 #include "detray/material/predefined_materials.hpp"
0024 #include "detray/utils/consistency_checker.hpp"
0025 #include "detray/utils/print_detector.hpp"
0026
0027 // Detray test include(s)
0028 #include "detray/test/common/factories/wire_layer_generator.hpp"
0029
0030 // Vecmem include(s)
0031 #include <vecmem/memory/memory_resource.hpp>
0032
0033 namespace detray {
0034
0035 /// Configuration for building a wire chamber detector
0036 template <concepts::scalar scalar_t, typename wire_shape_t = line_square>
0037 struct wire_chamber_config {
0038   /// Number of layers
0039   unsigned int m_n_layers{10u};
0040   /// The inner radius of the first layer
0041   scalar_t m_first_layer_inner_rad{500.f * unit<scalar_t>::mm};
0042   /// Half z of cylinder chamber
0043   scalar_t m_half_z{1000.f * unit<scalar_t>::mm};
0044   /// Radius of the material rods
0045   scalar_t m_mat_radius{15.f * unit<scalar_t>::um};
0046   /// Type of material for the material rods
0047   material<scalar_t> m_wire_mat{tungsten<scalar_t>()};
0048   /// Config for the wire generation (barrel)
0049   wire_layer_generator_config<scalar_t> m_wire_factory_cfg{};
0050   /// Configuration for the homogeneous material generator
0051   hom_material_config<scalar_t> m_material_config{};
0052   /// Do a full detector consistency check after building
0053   bool m_do_check{true};
0054
0055   /// Setters
0056   /// @{
0057   constexpr wire_chamber_config &n_layers(const unsigned int n) {
0058     m_n_layers = n;
0059     return *this;
0060   }
0061   constexpr wire_chamber_config &first_layer_inner_radius(const scalar_t r) {
0062     m_first_layer_inner_rad = r;
0063     return *this;
0064   }
0065   constexpr wire_chamber_config &half_z(const scalar_t hz) {
0066     m_half_z = hz;
0067     return *this;
0068   }
0069   constexpr wire_chamber_config &cell_size(const scalar_t c) {
0070     m_wire_factory_cfg.cell_size(c);
0071     return *this;
0072   }
0073   constexpr wire_chamber_config &stereo_angle(const scalar_t s) {
0074     m_wire_factory_cfg.stereo_angle(s);
0075     return *this;
0076   }
0077   constexpr wire_chamber_config &mat_radius(const scalar_t r) {
0078     m_mat_radius = r;
0079     return *this;
0080   }
0081   constexpr wire_chamber_config &wire_material(const material<scalar_t> &m) {
0082     m_wire_mat = m;
0083     return *this;
0084   }
0085   constexpr wire_chamber_config &do_check(const bool check) {
0086     m_do_check = check;
0087     return *this;
0088   }
0089   /// @}
0090
0091   /// Getters
0092   /// @{
0093   constexpr unsigned int n_layers() const { return m_n_layers; }
0094   constexpr scalar_t first_layer_inner_radius() const {
0095     return m_first_layer_inner_rad;
0096   }
0097   constexpr scalar_t half_z() const { return m_half_z; }
0098   constexpr scalar_t cell_size() const {
0099     return m_wire_factory_cfg.cell_size();
0100   }
0101   constexpr scalar_t stereo_angle() const {
0102     return m_wire_factory_cfg.stereo_angle();
0103   }
0104   constexpr scalar_t mat_radius() const { return m_mat_radius; }
0105   constexpr const material<scalar_t> &wire_material() const {
0106     return m_wire_mat;
0107   }
0108   constexpr wire_layer_generator_config<scalar_t> &layer_config() {
0109     return m_wire_factory_cfg;
0110   }
0111   constexpr const wire_layer_generator_config<scalar_t> &layer_config() const {
0112     return m_wire_factory_cfg;
0113   }
0114   constexpr auto &material_config() { return m_material_config; }
0115   constexpr const auto &material_config() const { return m_material_config; }
0116   constexpr bool do_check() const { return m_do_check; }
0117   /// @}
0118
0119  private:
0120   /// Print the wire chamber configuration
0121   friend inline std::ostream &operator<<(std::ostream &out,
0122                                          const wire_chamber_config &cfg) {
0123     out << "\nWire Chamber\n"
0124         << "----------------------------\n"
0125         << "  No. layers            : " << cfg.n_layers() << "\n"
0126         << "  First layer inner rad.: " << cfg.first_layer_inner_radius()
0127         << " [mm]\n"
0128         << "  Half length z         : " << cfg.half_z() << " [mm]\n";
0129
0130     if constexpr (std::same_as<wire_shape_t, line_square>) {
0131       out << "  Shape                 : wire cell\n";
0132     } else {
0133       out << "  Shape                 : straw tube\n";
0134     }
0135
0136     out << "  Cell size             : " << cfg.cell_size() << " [mm]\n"
0137         << "  Stereo angle          : " << cfg.stereo_angle() << " [rad]\n"
0138         << "  Wire material         : " << cfg.wire_material() << "\n"
0139         << "  Material rad.         : " << cfg.mat_radius() << " [mm]\n";
0140
0141     return out;
0142   }
0143
0144 };  // wire chamber config
0145
0146 template <concepts::algebra algebra_t, typename wire_shape_t>
0147 inline auto build_wire_chamber(
0148     vecmem::memory_resource &resource,
0149     wire_chamber_config<dscalar<algebra_t>, wire_shape_t> &cfg) {
0150   using builder_t =
0151       detector_builder<wire_chamber_metadata<algebra_t>, volume_builder>;
0152   using detector_t = typename builder_t::detector_type;
0153   using scalar_t = dscalar<typename detector_t::algebra_type>;
0154
0155   // Wire chamber detector builder
0156   builder_t det_builder;
0157   det_builder.set_name("wire_chamber");
0158
0159   // Geometry context object
0160   typename detector_t::geometry_context gctx{};
0161
0162   // Navigation link when leaving the detector
0163   using nav_link_t = typename detector_t::surface_type::navigation_link;
0164   constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
0165   const scalar_t inner_rad{cfg.first_layer_inner_radius()};
0166   const scalar_t cell_size{cfg.cell_size()};
0167
0168   // Prepare grid building
0169   constexpr auto grid_id{
0170       detector_t::accel::id::e_surface_concentric_cylinder2D_grid};
0171   using cyl_grid_t = types::get<typename detector_t::accel, grid_id>;
0172   using loc_bin_idx_t = typename cyl_grid_t::loc_bin_index;
0173   static_assert(cyl_grid_t::dim == 2);
0174   using grid_builder_t =
0175       grid_builder<detector_t, cyl_grid_t, detray::fill_by_pos>;
0176
0177   // Binning of the grid
0178   axis::multi_bin_range<cyl_grid_t::dim> bin_range{};
0179   // Min, max bin indices per axis
0180   bin_range[static_cast<std::size_t>(axis::label::e_rphi)] = {0, 100};
0181   bin_range[static_cast<std::size_t>(axis::label::e_cyl_z)] = {0, 1};
0182
0183   // Spans of the grid axes
0184   std::vector<scalar_t> sf_grid_spans{-constant<scalar_t>::pi,
0185                                       constant<scalar_t>::pi, -cfg.half_z(),
0186                                       cfg.half_z()};
0187
0188   constexpr unsigned int bin_capacity{3u};
0189   std::vector<std::pair<loc_bin_idx_t, dindex>> capacities{};
0190   capacities.reserve(
0191       static_cast<std::size_t>(bin_range[0][1] * bin_range[1][1]));
0192
0193   //
0194   // Build empty inner volume, where silicon subdetectors would sit
0195   //
0196   auto inner_v_builder = det_builder.new_volume(volume_id::e_cylinder);
0197   inner_v_builder->add_volume_placement(/*identity*/);
0198   const dindex inner_vol_idx{inner_v_builder->vol_index()};
0199   inner_v_builder->set_name("inner_vol_" + std::to_string(inner_vol_idx));
0200
0201   // Configure the portal factory
0202   // TODO: Add material maps that model the silicon detector budget
0203   cylinder_portal_config<scalar_t> inner_pt_cfg{};
0204   inner_pt_cfg.do_autofit(false)
0205       .fixed_half_length(cfg.half_z())
0206       .fixed_inner_radius(0.f)
0207       .fixed_outer_radius(inner_rad)
0208       // No inner subdetectors present -> don't build inner portal
0209       .build_inner(false)
0210       .link_north(inner_vol_idx + 1u)
0211       .link_south(leaving_world)
0212       .link_east(leaving_world)
0213       .link_west(leaving_world);
0214
0215   auto inner_pt_factory =
0216       std::make_shared<cylinder_portal_generator<detector_t>>(inner_pt_cfg);
0217   inner_v_builder->add_surfaces(inner_pt_factory, gctx);
0218
0219   //
0220   // Build layer volumes
0221   //
0222
0223   // Configure the homogeneous material generation
0224   auto &mat_cfg = cfg.material_config();
0225   // Only build material on sensitive surfaces
0226   constexpr auto vac{vacuum<scalar_t>{}};
0227   mat_cfg.portal_material(vac).passive_material(vac);
0228   mat_cfg.sensitive_material(tungsten<scalar_t>{}).thickness(cfg.mat_radius());
0229
0230   for (unsigned int i_lay = 0; i_lay < cfg.n_layers(); i_lay++) {
0231     // New volume for layer
0232     auto v_builder = det_builder.new_volume(volume_id::e_cylinder);
0233     auto vm_builder =
0234         det_builder.template decorate<homogeneous_material_builder<detector_t>>(
0235             v_builder);
0236
0237     const dindex vol_idx{vm_builder->vol_index()};
0238     vm_builder->set_name("layer_vol_" + std::to_string(vol_idx));
0239
0240     // The barrel volumes are centered at the origin
0241     vm_builder->add_volume_placement(/*identity*/);
0242
0243     // The maximal inner and outer radius of the volume
0244     const scalar_t inner_layer_rad =
0245         inner_rad + static_cast<scalar_t>(i_lay) * 2.f * cell_size;
0246     const scalar_t outer_layer_rad =
0247         inner_rad + static_cast<scalar_t>(i_lay + 1u) * 2.f * cell_size;
0248
0249     // Configure the wire layer factory for this layer
0250     auto &layer_cfg = cfg.layer_config();
0251     layer_cfg.inner_layer_radius(inner_layer_rad).half_length(cfg.half_z());
0252
0253     const scalar_t sign = (i_lay % 2 == 0) ? 1 : -1;
0254     layer_cfg.stereo_angle(sign * math::fabs(layer_cfg.stereo_angle()));
0255
0256     // Configure the portal factory
0257     cylinder_portal_config<scalar_t> layer_portal_cfg{};
0258     // Limit to maximum valid link
0259     auto link_north{i_lay == cfg.n_layers() - 1u ? leaving_world
0260                                                  : vol_idx + 1u};
0261
0262     layer_portal_cfg.do_autofit(false)
0263         .fixed_half_length(cfg.half_z())
0264         .fixed_inner_radius(inner_layer_rad)
0265         .fixed_outer_radius(outer_layer_rad)
0266         // Link the volume portals to its neighbors
0267         .link_north(link_north)
0268         .link_south(vol_idx - 1u)
0269         .link_east(leaving_world)
0270         .link_west(leaving_world);
0271
0272     // Register sensitive wires together with material
0273     auto wire_factory =
0274         std::make_unique<wire_layer_generator<detector_t, wire_shape_t>>(
0275             layer_cfg);
0276     auto wire_mat_factory =
0277         std::make_shared<homogeneous_material_generator<detector_t>>(
0278             std::move(wire_factory), mat_cfg);
0279
0280     // Register portals
0281     auto portal_mat_factory =
0282         std::make_shared<cylinder_portal_generator<detector_t>>(
0283             layer_portal_cfg);
0284
0285     // Add surfaces and material to volume builder
0286     vm_builder->add_surfaces(portal_mat_factory);
0287     vm_builder->add_surfaces(wire_mat_factory, gctx);
0288
0289     // Add a cylinder grid to every barrel layer
0290     auto vgr_builder =
0291         det_builder.template decorate<grid_builder_t>(vm_builder);
0292
0293     // Determine bin capacities
0294     capacities.clear();
0295     auto bin_indexer2D = axis::detail::get_bin_indexer(
0296         bin_range, std::make_integer_sequence<std::size_t, cyl_grid_t::dim>{});
0297     for (const auto [bin_idx0, bin_idx1] : bin_indexer2D) {
0298       // @Todo: fine-tune capacity
0299       loc_bin_idx_t loc_bin{static_cast<unsigned int>(bin_idx0),
0300                             static_cast<unsigned int>(bin_idx1)};
0301       capacities.emplace_back(loc_bin, bin_capacity);
0302     }
0303
0304     vgr_builder->set_type(detector_t::geo_obj_ids::e_sensitive);
0305     vgr_builder->init_grid(sf_grid_spans,
0306                            {static_cast<std::size_t>(bin_range[0][1]),
0307                             static_cast<std::size_t>(bin_range[1][1])},
0308                            capacities);
0309   }
0310
0311   // Build and return the detector and fill name map
0312   typename detector_t::name_map name_map{};
0313   auto det = det_builder.build(resource, name_map);
0314
0315   if (cfg.do_check()) {
0316     const bool verbose_check{false};
0317     detray::detail::check_consistency(det, verbose_check, name_map);
0318   }
0319
0320   DETRAY_DEBUG_HOST("\n" << detray::utils::print_detector(det, name_map));
0321
0322   return std::make_pair(std::move(det), std::move(name_map));
0323 }
0324
0325 }  // namespace detray