clusterization/impl/measurement_creation.ipp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2021-2026 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 #pragma once
0009
0010 // Library include(s).
0011 #include "traccc/utils/detray_conversion.hpp"
0012
0013 // System include(s).
0014 #include <cassert>
0015
0016 namespace traccc::details {
0017
0018 template <typename TCell, typename TDesign>
0019 TRACCC_HOST_DEVICE inline vector2 position_from_cell(
0020     const edm::silicon_cell<TCell>& cell,
0021     const traccc::detector_design_description_interface<TDesign>& module_dd) {
0022   // Calculate / construct the local cell position.
0023   vector2 cell_lower_position = {(module_dd.bin_edges_x()).at(cell.channel0()),
0024                                  (module_dd.bin_edges_y()).at(cell.channel1())};
0025
0026   vector2 cell_upper_position = {
0027       (module_dd.bin_edges_x()).at(cell.channel0() + 1),
0028       (module_dd.bin_edges_y()).at(cell.channel1() + 1)};
0029
0030   vector2 cell_middle_position = {
0031       scalar{0.5f} * (cell_upper_position[0] + cell_lower_position[0]),
0032       scalar{0.5f} * (cell_upper_position[1] + cell_lower_position[1])};
0033
0034   return cell_middle_position;
0035 }
0036
0037 template <typename T, typename TDesign>
0038 TRACCC_HOST_DEVICE inline void calc_cluster_properties(
0039     const edm::silicon_cluster<T>& cluster,
0040     const edm::silicon_cell_collection::const_device& cells,
0041     const traccc::detector_design_description_interface<TDesign>& module_dd,
0042     point2& mean, point2& var, scalar& totalWeight) {
0043   point2 offset{0.f, 0.f}, width{0.f, 0.f};
0044   bool first_processed = false;
0045
0046   unsigned int min_channel0 = std::numeric_limits<unsigned int>::max();
0047   unsigned int max_channel0 = std::numeric_limits<unsigned int>::lowest();
0048   unsigned int min_channel1 = std::numeric_limits<unsigned int>::max();
0049   unsigned int max_channel1 = std::numeric_limits<unsigned int>::lowest();
0050
0051   // Loop over the cell indices of the cluster.
0052   for (const unsigned int cell_idx : cluster.cell_indices()) {
0053     // The cell object.
0054     const edm::silicon_cell cell = cells.at(cell_idx);
0055
0056     // Translate the cell readout value into a weight.
0057     const scalar weight = cell.activation();
0058
0059     // Update all output properties with this cell.
0060     totalWeight += weight;
0061     scalar weight_factor = weight / totalWeight;
0062
0063     point2 cell_position = position_from_cell(cell, module_dd);
0064
0065     min_channel0 = std::min(min_channel0, cell.channel0());
0066     min_channel1 = std::min(min_channel1, cell.channel1());
0067     max_channel0 = std::max(max_channel0, cell.channel0());
0068     max_channel1 = std::max(max_channel1, cell.channel1());
0069
0070     if (!first_processed) {
0071       offset = cell_position;
0072       first_processed = true;
0073     }
0074
0075     cell_position = cell_position - offset;
0076
0077     const point2 diff_old = cell_position - mean;
0078     mean = mean + diff_old * weight_factor;
0079     const point2 diff_new = cell_position - mean;
0080
0081     var[0] = (1.f - weight_factor) * var[0] +
0082              weight_factor * (diff_old[0] * diff_new[0]);
0083     var[1] = (1.f - weight_factor) * var[1] +
0084              weight_factor * (diff_old[1] * diff_new[1]);
0085   }
0086
0087   // cluster width in the number of cells
0088   unsigned int delta0 = (max_channel0 - min_channel0) + 1;
0089   unsigned int delta1 = (max_channel1 - min_channel1) + 1;
0090
0091   vector2 cluster_lower_position = {(module_dd.bin_edges_x()).at(min_channel0),
0092                                     (module_dd.bin_edges_y()).at(min_channel1)};
0093
0094   vector2 cluster_upper_position = {
0095       (module_dd.bin_edges_x()).at(max_channel0 + 1),
0096       (module_dd.bin_edges_y()).at(max_channel1 + 1)};
0097
0098   width[0] = cluster_upper_position[0] - cluster_lower_position[0];
0099   width[1] = cluster_upper_position[1] - cluster_lower_position[1];
0100
0101   point2 pitch = {width[0] / static_cast<float>(delta0),
0102                   width[1] / static_cast<float>(delta1)};
0103
0104   var = var + point2{pitch[0] * pitch[0] / static_cast<scalar>(12.),
0105                      pitch[1] * pitch[1] / static_cast<scalar>(12.)};
0106
0107   mean = mean + offset;
0108 }
0109
0110 template <typename T1, typename T2>
0111 TRACCC_HOST_DEVICE inline void fill_measurement(
0112     edm::measurement<T1>& measurement, const edm::silicon_cluster<T2>& cluster,
0113     const unsigned int index,
0114     const edm::silicon_cell_collection::const_device& cells,
0115     const detector_design_description::const_device& det_descr,
0116     const detector_conditions_description::const_device& det_cond) {
0117   // To calculate the mean and variance with high numerical stability
0118   // we use a weighted variant of Welford's algorithm. This is a
0119   // single-pass online algorithm that works well for large numbers
0120   // of samples, as well as samples with very high values.
0121   //
0122   // To learn more about this algorithm please refer to:
0123   // [1] https://doi.org/10.1080/00401706.1962.10490022
0124   // [2] The Art of Computer Programming, Donald E. Knuth, second
0125   //     edition, chapter 4.2.2.
0126
0127   // Security checks.
0128   assert(cluster.cell_indices().empty() == false);
0129   assert([&]() {
0130     const unsigned int module_idx =
0131         cells.module_index().at(cluster.cell_indices().front());
0132     for (const unsigned int cell_idx : cluster.cell_indices()) {
0133       if (cells.module_index().at(cell_idx) != module_idx) {
0134         return false;
0135       }
0136     }
0137     return true;
0138   }() == true);
0139
0140   // The index of the module the cluster is on.
0141   const unsigned int module_idx =
0142       cells.module_index().at(cluster.cell_indices().front());
0143   const auto module_cd = det_cond.at(module_idx);
0144   const unsigned int design_idx = module_cd.module_to_design_id();
0145   const auto module_dd = det_descr.at(design_idx);
0146
0147   // Calculate the cluster properties
0148   scalar totalWeight = 0.f;
0149   point2 mean{0.f, 0.f}, var{0.f, 0.f};
0150   calc_cluster_properties(cluster, cells, module_dd, mean, var, totalWeight);
0151   assert(totalWeight > 0.f);
0152
0153   // Fill the measurement object.
0154   measurement.surface_link() = module_cd.geometry_id();
0155
0156   // apply lorentz shift to the cell position
0157   measurement.local_position() = utils::to_float_array<default_algebra>(
0158       mean + module_cd.measurement_translation());
0159
0160   // plus pitch^2 / 12
0161   measurement.local_variance() = utils::to_float_array<default_algebra>(var);
0162
0163   // For the ambiguity resolution algorithm, give a unique measurement ID
0164   measurement.identifier() = index;
0165   measurement.cluster_index() = index;
0166
0167   // Set the measurement dimensionality.
0168   measurement.dimensions() = module_dd.dimensions();
0169
0170   // Set the measurement's subspace.
0171   measurement.set_subspace(module_dd.subspace());
0172
0173   // Save the index of the cluster that produced this measurement
0174   measurement.cluster_index() = static_cast<unsigned int>(index);
0175 }
0176
0177 }  // namespace traccc::details