Digitization/src/DigitizationAlgorithm.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 "ActsExamples/Digitization/DigitizationAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "ActsExamples/Digitization/ModuleClusters.hpp"
0015 #include "ActsExamples/EventData/GeometryContainers.hpp"
0016 #include "ActsExamples/EventData/Index.hpp"
0017 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0018
0019 #include <algorithm>
0020 #include <array>
0021 #include <limits>
0022 #include <numeric>
0023 #include <ostream>
0024 #include <stdexcept>
0025 #include <string>
0026 #include <utility>
0027
0028 namespace ActsExamples {
0029
0030 DigitizationAlgorithm::DigitizationAlgorithm(
0031     Config config, std::unique_ptr<const Acts::Logger> logger)
0032     : IAlgorithm("DigitizationAlgorithm", std::move(logger)),
0033       m_cfg(std::move(config)) {
0034   if (m_cfg.inputSimHits.empty()) {
0035     throw std::invalid_argument("Missing simulated hits input collection");
0036   }
0037   if (m_cfg.surfaceByIdentifier.empty()) {
0038     throw std::invalid_argument("Missing Surface-GeometryID association map");
0039   }
0040   if (!m_cfg.randomNumbers) {
0041     throw std::invalid_argument("Missing random numbers tool");
0042   }
0043   if (m_cfg.digitizationConfigs.empty()) {
0044     throw std::invalid_argument("Missing digitization configuration");
0045   }
0046
0047   if (m_cfg.doClusterization) {
0048     if (m_cfg.outputMeasurements.empty()) {
0049       throw std::invalid_argument("Missing measurements output collection");
0050     }
0051     if (m_cfg.outputClusters.empty()) {
0052       throw std::invalid_argument("Missing cluster output collection");
0053     }
0054     if (m_cfg.outputMeasurementParticlesMap.empty()) {
0055       throw std::invalid_argument(
0056           "Missing hit-to-particles map output collection");
0057     }
0058     if (m_cfg.outputMeasurementSimHitsMap.empty()) {
0059       throw std::invalid_argument(
0060           "Missing hit-to-simulated-hits map output collection");
0061     }
0062     if (m_cfg.outputParticleMeasurementsMap.empty()) {
0063       throw std::invalid_argument(
0064           "Missing particle-to-measurements map output collection");
0065     }
0066     if (m_cfg.outputSimHitMeasurementsMap.empty()) {
0067       throw std::invalid_argument(
0068           "Missing particle-to-simulated-hits map output collection");
0069     }
0070
0071     m_outputMeasurements.initialize(m_cfg.outputMeasurements);
0072     m_outputMeasurementSubset.initialize(m_cfg.outputMeasurementSubset);
0073     m_outputClusters.initialize(m_cfg.outputClusters);
0074     m_outputMeasurementParticlesMap.initialize(
0075         m_cfg.outputMeasurementParticlesMap);
0076     m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
0077     m_outputParticleMeasurementsMap.initialize(
0078         m_cfg.outputParticleMeasurementsMap);
0079     m_outputSimHitMeasurementsMap.initialize(m_cfg.outputSimHitMeasurementsMap);
0080   }
0081
0082   if (m_cfg.doOutputCells) {
0083     if (m_cfg.outputCells.empty()) {
0084       throw std::invalid_argument("Missing cell output collection");
0085     }
0086
0087     m_outputCells.initialize(m_cfg.outputCells);
0088   }
0089
0090   m_inputHits.initialize(m_cfg.inputSimHits);
0091
0092   // Create the digitizers from the configuration
0093   std::vector<std::pair<Acts::GeometryIdentifier, Digitizer>> digitizerInput;
0094
0095   for (std::size_t i = 0; i < m_cfg.digitizationConfigs.size(); ++i) {
0096     GeometricConfig geoCfg;
0097     Acts::GeometryIdentifier geoId = m_cfg.digitizationConfigs.idAt(i);
0098
0099     const auto& digiCfg = m_cfg.digitizationConfigs.valueAt(i);
0100     geoCfg = digiCfg.geometricDigiConfig;
0101
0102     if (!geoCfg.indices.empty() && geoCfg.segmentation == nullptr) {
0103       throw std::invalid_argument(
0104           "Geometric digitization requires a segmentation");
0105     }
0106     // Copy so we can sort in-place
0107     SmearingConfig smCfg = digiCfg.smearingDigiConfig;
0108
0109     std::vector<Acts::BoundIndices> indices;
0110     for (auto& gcf : smCfg.params) {
0111       indices.push_back(gcf.index);
0112     }
0113     indices.insert(indices.begin(), geoCfg.indices.begin(),
0114                    geoCfg.indices.end());
0115
0116     // Make sure the configured input parameter indices are sorted and unique
0117     std::ranges::sort(indices);
0118
0119     auto dup = std::adjacent_find(indices.begin(), indices.end());
0120     if (dup != indices.end()) {
0121       throw std::invalid_argument(
0122           "Digitization configuration contains duplicate parameter indices");
0123     }
0124
0125     switch (smCfg.params.size()) {
0126       case 0u:
0127         digitizerInput.emplace_back(geoId, makeDigitizer<0u>(digiCfg));
0128         break;
0129       case 1u:
0130         digitizerInput.emplace_back(geoId, makeDigitizer<1u>(digiCfg));
0131         break;
0132       case 2u:
0133         digitizerInput.emplace_back(geoId, makeDigitizer<2u>(digiCfg));
0134         break;
0135       case 3u:
0136         digitizerInput.emplace_back(geoId, makeDigitizer<3u>(digiCfg));
0137         break;
0138       case 4u:
0139         digitizerInput.emplace_back(geoId, makeDigitizer<4u>(digiCfg));
0140         break;
0141       default:
0142         throw std::invalid_argument("Unsupported smearer size");
0143     }
0144   }
0145
0146   m_digitizers = Acts::GeometryHierarchyMap<Digitizer>(digitizerInput);
0147 }
0148
0149 ProcessCode DigitizationAlgorithm::execute(const AlgorithmContext& ctx) const {
0150   // Retrieve input
0151   const auto& simHits = m_inputHits(ctx);
0152   ACTS_DEBUG("Loaded " << simHits.size() << " sim hits");
0153
0154   // Prepare output containers
0155   // need list here for stable addresses
0156   MeasurementContainer measurements;
0157   ClusterContainer clusters;
0158
0159   MeasurementParticlesMap measurementParticlesMap;
0160   MeasurementSimHitsMap measurementSimHitsMap;
0161   measurements.reserve(simHits.size());
0162   measurementParticlesMap.reserve(simHits.size());
0163   measurementSimHitsMap.reserve(simHits.size());
0164
0165   // Setup random number generator
0166   auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
0167
0168   // Some statistics
0169   std::size_t skippedHits = 0;
0170
0171   // Some algorithms do the clusterization themselves such as the traccc chain.
0172   // Thus we need to store the cell data from the simulation.
0173   CellsMap cellsMap;
0174
0175   ACTS_DEBUG("Starting loop over modules ...");
0176   for (const auto& simHitsGroup : groupByModule(simHits)) {
0177     // Manual pair unpacking instead of using
0178     //   auto [moduleGeoId, moduleSimHits] : ...
0179     // otherwise clang on macos complains that it is unable to capture the local
0180     // binding in the lambda used for visiting the smearer below.
0181     Acts::GeometryIdentifier moduleGeoId = simHitsGroup.first;
0182     const auto& moduleSimHits = simHitsGroup.second;
0183
0184     auto surfaceItr = m_cfg.surfaceByIdentifier.find(moduleGeoId);
0185
0186     if (surfaceItr == m_cfg.surfaceByIdentifier.end()) {
0187       // this is either an invalid geometry id or a misconfigured smearer
0188       // setup; both cases can not be handled and should be fatal.
0189       ACTS_ERROR("Could not find surface " << moduleGeoId
0190                                            << " for configured smearer");
0191       return ProcessCode::ABORT;
0192     }
0193
0194     const Acts::Surface* surfacePtr = surfaceItr->second;
0195
0196     auto digitizerItr = m_digitizers.find(moduleGeoId);
0197     if (digitizerItr == m_digitizers.end()) {
0198       ACTS_VERBOSE("No digitizer present for module " << moduleGeoId);
0199       continue;
0200     } else {
0201       ACTS_VERBOSE("Digitizer found for module " << moduleGeoId);
0202     }
0203
0204     // Run the digitizer. Iterate over the hits for this surface inside the
0205     // visitor so we do not need to lookup the variant object per-hit.
0206     std::visit(
0207         [&](const auto& digitizer) {
0208           ModuleClusters moduleClusters(
0209               digitizer.geometric.segmentation, digitizer.geometric.indices,
0210               m_cfg.doMerge, m_cfg.mergeNsigma, m_cfg.mergeCommonCorner);
0211
0212           for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
0213             const auto& simHit = *h;
0214             const auto simHitIdx = simHits.index_of(h);
0215
0216             DigitizedParameters dParameters;
0217
0218             if (simHit.depositedEnergy() < m_cfg.minEnergyDeposit) {
0219               ACTS_VERBOSE("Skip hit because energy deposit to small");
0220               continue;
0221             }
0222
0223             // Geometric part - 0, 1, 2 local parameters are possible
0224             if (!digitizer.geometric.indices.empty()) {
0225               ACTS_VERBOSE("Configured to geometric digitize "
0226                            << digitizer.geometric.indices.size()
0227                            << " parameters.");
0228               const auto& cfg = digitizer.geometric;
0229               Acts::Vector3 driftDir = cfg.drift(simHit.position(), rng);
0230               auto channelsRes = m_channelizer.channelize(
0231                   simHit, *surfacePtr, ctx.geoContext, driftDir,
0232                   *cfg.segmentation, cfg.thickness);
0233               if (!channelsRes.ok() || channelsRes->empty()) {
0234                 ACTS_DEBUG(
0235                     "Geometric channelization did not work, skipping this "
0236                     "hit.");
0237                 continue;
0238               }
0239               ACTS_VERBOSE("Activated " << channelsRes->size()
0240                                         << " channels for this hit.");
0241               dParameters =
0242                   localParameters(digitizer.geometric, *channelsRes, rng);
0243               if (dParameters.cluster.channels.empty()) {
0244                 ACTS_DEBUG("All channels below threshold, skipping this hit.");
0245                 continue;
0246               }
0247             }
0248
0249             // Smearing part - (optionally) rest
0250             if (!digitizer.smearing.indices.empty()) {
0251               ACTS_VERBOSE("Configured to smear "
0252                            << digitizer.smearing.indices.size()
0253                            << " parameters.");
0254               auto res =
0255                   digitizer.smearing(rng, simHit, *surfacePtr, ctx.geoContext);
0256               if (!res.ok()) {
0257                 ++skippedHits;
0258                 ACTS_DEBUG("Problem in hit smearing, skip hit ("
0259                            << res.error().message() << ")");
0260                 continue;
0261               }
0262               const auto& [par, cov] = res.value();
0263               for (Eigen::Index ip = 0; ip < par.rows(); ++ip) {
0264                 dParameters.indices.push_back(digitizer.smearing.indices[ip]);
0265                 dParameters.values.push_back(par[ip]);
0266                 dParameters.variances.push_back(cov(ip, ip));
0267               }
0268             }
0269
0270             // Check on success - threshold could have eliminated all channels
0271             if (dParameters.values.empty()) {
0272               ACTS_VERBOSE(
0273                   "Parameter digitization did not yield a measurement.");
0274               continue;
0275             }
0276
0277             moduleClusters.add(std::move(dParameters), simHitIdx);
0278           }
0279
0280           auto digitizeParametersResult = moduleClusters.digitizedParameters();
0281
0282           // Store the cell data into a map.
0283           if (m_cfg.doOutputCells) {
0284             std::vector<Cluster::Cell> cells;
0285             for (const auto& [dParameters, simHitsIdxs] :
0286                  digitizeParametersResult) {
0287               for (const auto& cell : dParameters.cluster.channels) {
0288                 cells.push_back(cell);
0289               }
0290             }
0291             cellsMap.insert({moduleGeoId, std::move(cells)});
0292           }
0293
0294           if (m_cfg.doClusterization) {
0295             for (auto& [dParameters, simHitsIdxs] : digitizeParametersResult) {
0296               auto measurement =
0297                   createMeasurement(measurements, moduleGeoId, dParameters);
0298
0299               dParameters.cluster.globalPosition = measurementGlobalPosition(
0300                   dParameters, *surfacePtr, ctx.geoContext);
0301               clusters.emplace_back(std::move(dParameters.cluster));
0302
0303               for (auto simHitIdx : simHitsIdxs) {
0304                 measurementParticlesMap.emplace_hint(
0305                     measurementParticlesMap.end(), measurement.index(),
0306                     simHits.nth(simHitIdx)->particleId());
0307                 measurementSimHitsMap.emplace_hint(measurementSimHitsMap.end(),
0308                                                    measurement.index(),
0309                                                    simHitIdx);
0310               }
0311             }
0312           }
0313         },
0314         *digitizerItr);
0315   }
0316
0317   if (skippedHits > 0) {
0318     ACTS_WARNING(
0319         skippedHits
0320         << " skipped in Digitization. Enable DEBUG mode to see more details.");
0321   }
0322
0323   if (m_cfg.doClusterization) {
0324     ACTS_DEBUG("Created " << measurements.size() << " measurements, "
0325                           << clusters.size() << " clusters" << " from "
0326                           << simHits.size() << " sim hits.");
0327
0328     const auto& storedMeasurements =
0329         m_outputMeasurements(ctx, std::move(measurements));
0330
0331     // Build initial full subset: all measurements, indices in original space.
0332     std::vector<MeasurementContainer::Index> allIndices(
0333         storedMeasurements.size());
0334     std::iota(allIndices.begin(), allIndices.end(), Index{0});
0335     m_outputMeasurementSubset(
0336         ctx, MeasurementSubset(storedMeasurements, std::move(allIndices)));
0337
0338     m_outputClusters(ctx, std::move(clusters));
0339
0340     // invert them before they are moved
0341     m_outputParticleMeasurementsMap(
0342         ctx, invertIndexMultimap(measurementParticlesMap));
0343     m_outputSimHitMeasurementsMap(ctx,
0344                                   invertIndexMultimap(measurementSimHitsMap));
0345
0346     m_outputMeasurementParticlesMap(ctx, std::move(measurementParticlesMap));
0347     m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
0348   }
0349
0350   if (m_cfg.doOutputCells) {
0351     m_outputCells(ctx, std::move(cellsMap));
0352   }
0353
0354   return ProcessCode::SUCCESS;
0355 }
0356
0357 DigitizedParameters DigitizationAlgorithm::localParameters(
0358     const GeometricConfig& geoCfg,
0359     const std::vector<ActsFatras::Segmentizer::ChannelSegment>& channels,
0360     RandomEngine& rng) const {
0361   DigitizedParameters dParameters;
0362
0363   // For digital readout, the weight needs to be split in x and y
0364   std::array<double, 2u> pos = {0., 0.};
0365   std::array<double, 2u> totalWeight = {0., 0.};
0366   std::array<std::size_t, 2u> bmin = {std::numeric_limits<std::size_t>::max(),
0367                                       std::numeric_limits<std::size_t>::max()};
0368   std::array<std::size_t, 2u> bmax = {0, 0};
0369
0370   // The component digital store
0371   std::array<std::set<std::size_t>, 2u> componentChannels;
0372
0373   // Combine the channels
0374   for (const auto& ch : channels) {
0375     auto bin = ch.bin;
0376     double charge = geoCfg.charge(ch.activation, rng);
0377     // Loop and check
0378     if (charge > geoCfg.threshold) {
0379       double weight = geoCfg.digital ? 1. : charge;
0380       for (std::size_t ib = 0; ib < 2; ++ib) {
0381         // Cell bins are zero-based while the axis bin indices start at one
0382         const double binCenter =
0383             geoCfg.segmentation->getAxis(ib).getBinCenter(bin[ib] + 1);
0384         if (geoCfg.digital && geoCfg.componentDigital) {
0385           // only fill component of this row/column if not yet filled
0386           if (!componentChannels[ib].contains(bin[ib])) {
0387             totalWeight[ib] += weight;
0388             pos[ib] += weight * binCenter;
0389             componentChannels[ib].insert(bin[ib]);
0390           }
0391         } else {
0392           totalWeight[ib] += weight;
0393           pos[ib] += weight * binCenter;
0394         }
0395         // min max channels
0396         bmin[ib] = std::min(bmin[ib], static_cast<std::size_t>(bin[ib]));
0397         bmax[ib] = std::max(bmax[ib], static_cast<std::size_t>(bin[ib]));
0398       }
0399       // Create a copy of the channel, as activation may change
0400       auto chdig = ch;
0401       chdig.bin = ch.bin;
0402       chdig.activation = charge;
0403       dParameters.cluster.channels.push_back(chdig);
0404     }
0405   }
0406   if (totalWeight[0] > 0. && totalWeight[1] > 0.) {
0407     pos[0] /= totalWeight[0];
0408     pos[1] /= totalWeight[1];
0409     dParameters.indices = geoCfg.indices;
0410     for (auto idx : dParameters.indices) {
0411       dParameters.values.push_back(pos[idx]);
0412     }
0413     std::size_t size0 = (bmax[0] - bmin[0] + 1);
0414     std::size_t size1 = (bmax[1] - bmin[1] + 1);
0415
0416     dParameters.variances = geoCfg.variances({size0, size1}, bmin);
0417     dParameters.cluster.sizeLoc0 = size0;
0418     dParameters.cluster.sizeLoc1 = size1;
0419   }
0420
0421   return dParameters;
0422 }
0423
0424 }  // namespace ActsExamples