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