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