EIC code displayed by LXR master/​EICrecon/​src/​algorithms/​digi/​RandomNoisePixel.cc	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-08-18 08:30:23

 // SPDX-License-Identifier: LGPL-3.0-or-later
 // Copyright (C) 2025 Minjung Kim, Joshua Sobaljic, Shujie Li
 //
 // RandomNoisePixel caches the discrete addressable pixels of each sensitive
 // component in compact ranges, then samples noise from one global per-pixel rate.
 //
 // If changing this, please also make necessary updates to the
 // docs/design/svt_pixel_noise_injection.md
 
 #include "RandomNoisePixel.h"
 
 #include <DD4hep/Alignments.h>
 #include <DD4hep/IDDescriptor.h>
 #include <DD4hep/Objects.h>
 #include <DD4hep/Segmentations.h>
 #include <DD4hep/Volumes.h>
 #include <DDSegmentation/BitFieldCoder.h>
 #include <DDSegmentation/CartesianGridXY.h>
 #include <DDSegmentation/CartesianGridXZ.h>
 #include <DDSegmentation/CylindricalGridPhiZ.h>
 #include <DDSegmentation/MultiSegmentation.h>
 #include <DDSegmentation/Segmentation.h>
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
 #include <RtypesCore.h>
 #include <TGeoBBox.h>
 #include <TGeoNode.h>
 #include <TGeoShape.h>
 #include <TGeoTube.h>
 #include <algorithms/geo.h>
 #include <algorithm>
 #include <array>
 #include <cctype>
 #include <cmath>
 #include <gsl/pointers>
 #include <initializer_list>
 #include <iterator>
 #include <limits>
 #include <mutex>
 #include <numbers>
 #include <optional>
 #include <stdexcept>
 #include <tuple>
 #include <utility>
 
 #include "algorithms/digi/RandomNoisePixelConfig.h"
 
 namespace eicrecon {
 namespace {
 
   using RawSegmentation = dd4hep::DDSegmentation::Segmentation;
 
   // ROOT's geometry navigator is shared. This mutex protects only the short
   // initialization phase in which global positions are converted to cell IDs.
   std::mutex& geometryInitializationMutex() {
     static std::mutex mutex;
     return mutex;
   }
 
   // Make DD4hep volume-ID field matching insensitive to capitalization.
   std::string lower(std::string value) {
     std::transform(value.begin(), value.end(), value.begin(),
                    [](unsigned char c) { return static_cast<char>(std::tolower(c)); });
     return value;
   }
 
   // Read the layer number from a placed volume's explicit DD4hep volume IDs.
   std::optional<int> layerFromPlacement(const dd4hep::PlacedVolume& placement) {
     if (!placement.isValid()) {
       return std::nullopt;
     }
     for (const auto& [field, value] : placement.volIDs()) {
       const auto name = lower(field);
       if (name == "layer" || name.find("lay") != std::string::npos) {
         return value;
       }
     }
     return std::nullopt;
   }
 
   // Prefer the layer placement ID, then try the module placement for geometries
   // that attach the layer field one level lower in the hierarchy.
   std::optional<int> layerFromDetElement(const dd4hep::DetElement& layer,
                                          const dd4hep::DetElement& module) {
     if (auto id = layerFromPlacement(layer.placement())) {
       return id;
     }
     return layerFromPlacement(module.placement());
   }
 
   // Copy a placed volume's local-to-parent transformation into an owned matrix.
   TGeoHMatrix placementTransform(const dd4hep::PlacedVolume& placement) {
     if (auto* node = placement.ptr()) {
       if (auto* matrix = node->GetMatrix()) {
         return TGeoHMatrix{*matrix};
       }
     }
     return TGeoHMatrix{};
   }
 
   // Compose child-local -> parent-local with parent-local -> world.
   TGeoHMatrix composeTransforms(const TGeoHMatrix& parentToWorld,
                                 const TGeoHMatrix& childToParent) {
     TGeoHMatrix childToWorld{parentToWorld};
     childToWorld.Multiply(childToParent);
     return childToWorld;
   }
 
   // Accept a sensitive volume only when its readout exactly matches this
   // RandomNoisePixel instance (for example, SiBarrelHits but not MPGD hits).
   bool hasRequestedReadout(const dd4hep::PlacedVolume& placement, std::string_view readoutName) {
     if (!placement.isValid()) {
       return false;
     }
     dd4hep::SensitiveDetector sensitive{placement.volume().sensitiveDetector()};
     if (!sensitive.isValid()) {
       return false;
     }
     const auto readout = sensitive.readout();
     return readout.isValid() && readout.name() == readoutName;
   }
 
   // Record the shape and transform needed to identify one sensitive placement.
   // This initialization-only information is discarded after the cache is built.
   RandomNoisePixel::SensitiveComponentGeometry
   makeComponentGeometry(const std::string& detectorName, int layer,
                         const dd4hep::PlacedVolume& placement, const TGeoHMatrix& localToWorld) {
     auto volume = placement.volume();
     if (!volume.isValid() || !volume.ptr() || !volume->GetShape()) {
       throw std::runtime_error("RandomNoisePixel sensitive volume has no TGeo shape");
     }
     return {.detectorName          = detectorName,
             .layer                 = layer,
             .volume                = volume.ptr(),
             .localToWorldTransform = localToWorld};
   }
 
   // Recursively search a module assembly for sensitive descendants while
   // accumulating the complete sensitive-local to world transformation.
   void
   collectSensitiveDescendants(const std::string& detectorName, int layer,
                               const dd4hep::PlacedVolume& parent, const TGeoHMatrix& parentToWorld,
                               std::string_view readoutName,
                               std::vector<RandomNoisePixel::SensitiveComponentGeometry>& output) {
     auto* node = parent.ptr();
     if (!node) {
       return;
     }
 
     for (int i = 0; i < node->GetNdaughters(); ++i) {
       auto* daughterNode = node->GetDaughter(i);
       if (!daughterNode) {
         continue;
       }
       dd4hep::PlacedVolume daughter{daughterNode};
       const auto daughterToWorld = composeTransforms(parentToWorld, placementTransform(daughter));
       if (hasRequestedReadout(daughter, readoutName)) {
         output.push_back(makeComponentGeometry(detectorName, layer, daughter, daughterToWorld));
       } else {
         collectSensitiveDescendants(detectorName, layer, daughter, daughterToWorld, readoutName,
                                     output);
       }
     }
   }
 
   // Return all matching sensitive pieces belonging to one module. Some modules
   // are directly sensitive; others contain two or more sensitive daughters.
   std::vector<RandomNoisePixel::SensitiveComponentGeometry>
   findSensitiveComponents(const std::string& detectorName, int layer,
                           const dd4hep::DetElement& module, std::string_view readoutName) {
     const auto modulePlacement = module.placement();
     if (!modulePlacement.isValid()) {
       return {};
     }
 
     std::vector<RandomNoisePixel::SensitiveComponentGeometry> components;
     const auto moduleToWorld = module.nominal().worldTransformation();
     if (hasRequestedReadout(modulePlacement, readoutName)) {
       components.push_back(
           makeComponentGeometry(detectorName, layer, modulePlacement, moduleToWorld));
     } else {
       collectSensitiveDescendants(detectorName, layer, modulePlacement, moduleToWorld, readoutName,
                                   components);
     }
     return components;
   }
 
   // Access the axis-aligned local bounds supplied by ROOT for the sensor shape.
   const TGeoBBox& boundingBox(const RandomNoisePixel::SensitiveComponentGeometry& component) {
     const auto* box = dynamic_cast<const TGeoBBox*>(component.volume->GetShape());
     if (!box) {
       throw std::runtime_error("RandomNoisePixel requires a TGeo shape with bounding-box data");
     }
     return *box;
   }
 
   // Choose a point safely inside the sensor. Tube segments require their radial
   // and angular midpoint because the center of their bounding box can be empty.
   std::array<double, 3>
   referenceLocalPoint(const RandomNoisePixel::SensitiveComponentGeometry& component) {
     if (const auto* tube = dynamic_cast<const TGeoTubeSeg*>(component.volume->GetShape())) {
       const double radius = 0.5 * (tube->GetRmin() + tube->GetRmax());
       const double phi    = 0.5 * (tube->GetPhi1() + tube->GetPhi2()) * std::numbers::pi / 180.0;
       return {radius * std::cos(phi), radius * std::sin(phi), 0.0};
     }
 
     const auto& box    = boundingBox(component);
     const auto* origin = box.GetOrigin();
     return {origin[0], origin[1], origin[2]};
   }
 
   // Resolve a MultiSegmentation to the concrete grid selected by fields already
   // present in the sensor's base volume ID (the BVTX discriminator is layer).
   const RawSegmentation& selectedSegmentation(const dd4hep::Segmentation& segmentation,
                                               std::uint64_t baseVolumeID) {
     const RawSegmentation* selected = segmentation.segmentation();
     while (selected && selected->type() == "MultiSegmentation") {
       const auto* multi = dynamic_cast<const dd4hep::DDSegmentation::MultiSegmentation*>(selected);
       if (!multi) {
         throw std::runtime_error("RandomNoisePixel cannot access MultiSegmentation");
       }
       selected = &multi->subsegmentation(baseVolumeID);
     }
     if (!selected) {
       throw std::runtime_error("RandomNoisePixel readout has no segmentation");
     }
     return *selected;
   }
 
   // Convert a physical coordinate interval into the inclusive integer indices
   // of pixel centers inside it, then clamp to the cell-ID field's legal range.
   std::pair<std::int64_t, std::int64_t>
   indexRange(double minimum, double maximum, double pitch, double offset,
              const dd4hep::DDSegmentation::BitFieldElement& field) {
     if (!(pitch > 0.0) || !std::isfinite(pitch)) {
       throw std::runtime_error("RandomNoisePixel found an invalid segmentation pitch");
     }
     constexpr double indexTolerance = 1.0e-9;
     auto first = static_cast<std::int64_t>(std::ceil((minimum - offset) / pitch - indexTolerance));
     auto last  = static_cast<std::int64_t>(std::floor((maximum - offset) / pitch + indexTolerance));
     first      = std::max(first, static_cast<std::int64_t>(field.minValue()));
     last       = std::min(last, static_cast<std::int64_t>(field.maxValue()));
     return {first, last};
   }
 
   // A true box can be represented by one rectangular range without row tests.
   bool isPureBox(const TGeoShape& shape) {
     return std::string_view{shape.ClassName()} == "TGeoBBox";
   }
 
   // ROOT may classify a mathematically boundary-centered pixel as just outside
   // after a local -> global -> local floating-point round trip. Accept only
   // points within ROOT's own surface tolerance; larger misses remain errors.
   bool containsOrTouches(const TGeoShape& shape, const double point[3]) {
     return shape.Contains(point) || shape.Safety(point, false) <= 10.0 * TGeoShape::Tolerance();
   }
 
   // Build the exact discrete pixel layout for Cartesian XY or XZ grids. Boxes
   // become rectangles; trapezoids and clipped sensors become compact row spans.
   std::shared_ptr<RandomNoisePixel::PixelLayout>
   makeCartesianLayout(const RandomNoisePixel::SensitiveComponentGeometry& component,
                       RandomNoisePixel::GridKind kind, std::string firstField,
                       std::string secondField, double firstPitch, double secondPitch,
                       double firstOffset, double secondOffset,
                       const dd4hep::DDSegmentation::BitFieldCoder& decoder) {
     const auto& box           = boundingBox(component);
     const auto* origin        = box.GetOrigin();
     const bool xy             = kind == RandomNoisePixel::GridKind::CartesianXY;
     const double firstCenter  = origin[0];
     const double firstHalf    = box.GetDX();
     const double secondCenter = xy ? origin[1] : origin[2];
     const double secondHalf   = xy ? box.GetDY() : box.GetDZ();
     // Step 1: turn the sensor's physical bounds into candidate pixel-index bounds.
     const auto [firstMin, firstMax] = indexRange(firstCenter - firstHalf, firstCenter + firstHalf,
                                                  firstPitch, firstOffset, decoder[firstField]);
     const auto [secondMin, secondMax] =
         indexRange(secondCenter - secondHalf, secondCenter + secondHalf, secondPitch, secondOffset,
                    decoder[secondField]);
 
     if (firstMax < firstMin || secondMax < secondMin) {
       throw std::runtime_error("RandomNoisePixel sensitive component contains no pixel centers");
     }
 
     auto layout         = std::make_shared<RandomNoisePixel::PixelLayout>();
     layout->kind        = kind;
     layout->firstField  = std::move(firstField);
     layout->secondField = std::move(secondField);
 
     // Step 2a: a box contains every candidate pixel center, so store four limits.
     if (isPureBox(*component.volume->GetShape())) {
       layout->rectangular    = true;
       layout->firstMin       = firstMin;
       layout->firstMax       = firstMax;
       layout->secondMin      = secondMin;
       layout->secondMax      = secondMax;
       const auto firstCount  = static_cast<std::uint64_t>(firstMax - firstMin + 1);
       const auto secondCount = static_cast<std::uint64_t>(secondMax - secondMin + 1);
       if (firstCount > std::numeric_limits<std::uint64_t>::max() / secondCount) {
         throw std::overflow_error("RandomNoisePixel pixel count overflow");
       }
       layout->totalPixels = firstCount * secondCount;
       return layout;
     }
 
     // Step 2b: for a non-box shape, test pixel centers against TGeo::Contains().
     const auto* shape = component.volume->GetShape();
     auto contains     = [&](std::int64_t firstIndex, std::int64_t secondIndex) {
       const double first  = firstOffset + firstPitch * static_cast<double>(firstIndex);
       const double second = secondOffset + secondPitch * static_cast<double>(secondIndex);
       double local[3]     = {origin[0], origin[1], origin[2]};
       local[0]            = first;
       if (xy) {
         local[1] = second;
       } else {
         local[2] = second;
       }
       return containsOrTouches(*shape, local);
     };
 
     // Step 3: each supported non-box sensor is convex. Find one interior pixel
     // near the row center, then binary-search its left and right boundaries.
     const auto middleGuess = std::clamp(
         static_cast<std::int64_t>(std::llround((firstCenter - firstOffset) / firstPitch)), firstMin,
         firstMax);
     std::uint64_t cumulative = 0;
     for (std::int64_t second = secondMin; second <= secondMax; ++second) {
       if (!contains(middleGuess, second)) {
         continue;
       }
 
       std::int64_t low  = firstMin;
       std::int64_t high = middleGuess;
       while (low < high) {
         const auto middle = low + (high - low) / 2;
         if (contains(middle, second)) {
           high = middle;
         } else {
           low = middle + 1;
         }
       }
       const auto rowFirst = low;
 
       low  = middleGuess;
       high = firstMax;
       while (low < high) {
         const auto middle = low + (high - low + 1) / 2;
         if (contains(middle, second)) {
           low = middle;
         } else {
           high = middle - 1;
         }
       }
       const auto rowLast  = low;
       const auto rowCount = static_cast<std::uint64_t>(rowLast - rowFirst + 1);
       if (cumulative > std::numeric_limits<std::uint64_t>::max() - rowCount) {
         throw std::overflow_error("RandomNoisePixel pixel row count overflow");
       }
       // The cumulative count maps a flat random integer back to this row.
       cumulative += rowCount;
       layout->rows.push_back({second, rowFirst, rowLast, cumulative});
     }
 
     if (layout->rows.empty()) {
       throw std::runtime_error("RandomNoisePixel found no segmentation centers inside TGeo shape");
     }
     layout->totalPixels = cumulative;
     return layout;
   }
 
   // Build a phi-z rectangle in the sensitive volume's local coordinate system.
   // DD4hep later applies the placement transform when converting a cell ID to
   // a global position, so using world bounds here would apply the placement twice.
   std::shared_ptr<RandomNoisePixel::PixelLayout>
   makeCylindricalLayout(const RandomNoisePixel::SensitiveComponentGeometry& component,
                         const dd4hep::DDSegmentation::CylindricalGridPhiZ& grid,
                         const dd4hep::DDSegmentation::BitFieldCoder& decoder) {
     // Step 1: establish a local central phi so angles remain continuous across +/-pi.
     const auto center      = referenceLocalPoint(component);
     const double centerPhi = std::atan2(center[1], center[0]);
 
     double phiMin          = std::numeric_limits<double>::max();
     double phiMax          = std::numeric_limits<double>::lowest();
     double zMin            = std::numeric_limits<double>::max();
     double zMax            = std::numeric_limits<double>::lowest();
     auto includeLocalPoint = [&](double x, double y, double z) {
       const double phi =
           centerPhi + std::remainder(std::atan2(y, x) - centerPhi, 2.0 * std::numbers::pi);
       phiMin = std::min(phiMin, phi);
       phiMax = std::max(phiMax, phi);
       zMin   = std::min(zMin, z);
       zMax   = std::max(zMax, z);
     };
 
     // Step 2: inspect the true local TubeSeg boundaries. For other supported
     // shapes, use local bounding-box corners; initialization validation below
     // rejects a layout if the segmentation and solid frames are incompatible.
     if (const auto* tube = dynamic_cast<const TGeoTubeSeg*>(component.volume->GetShape())) {
       for (double radius : {tube->GetRmin(), tube->GetRmax()}) {
         for (double phiDegrees : {tube->GetPhi1(), tube->GetPhi2()}) {
           const double phi = phiDegrees * std::numbers::pi / 180.0;
           for (double z : {-tube->GetDz(), tube->GetDz()}) {
             includeLocalPoint(radius * std::cos(phi), radius * std::sin(phi), z);
           }
         }
       }
     } else {
       const auto& box    = boundingBox(component);
       const auto* origin = box.GetOrigin();
       for (double sx : {-1.0, 1.0}) {
         for (double sy : {-1.0, 1.0}) {
           for (double sz : {-1.0, 1.0}) {
             includeLocalPoint(origin[0] + sx * box.GetDX(), origin[1] + sy * box.GetDY(),
                               origin[2] + sz * box.GetDZ());
           }
         }
       }
     }
 
     // Step 3: convert angular and longitudinal bounds to discrete pixel indices.
     const auto [phiIndexMin, phiIndexMax] = indexRange(
         phiMin, phiMax, grid.gridSizePhi(), grid.offsetPhi(), decoder[grid.fieldNamePhi()]);
     const auto [zIndexMin, zIndexMax] =
         indexRange(zMin, zMax, grid.gridSizeZ(), grid.offsetZ(), decoder[grid.fieldNameZ()]);
     if (phiIndexMax < phiIndexMin || zIndexMax < zIndexMin) {
       throw std::runtime_error("RandomNoisePixel cylindrical component contains no pixels");
     }
 
     auto layout         = std::make_shared<RandomNoisePixel::PixelLayout>();
     layout->kind        = RandomNoisePixel::GridKind::CylindricalPhiZ;
     layout->firstField  = grid.fieldNamePhi();
     layout->secondField = grid.fieldNameZ();
     layout->rectangular = true;
     layout->firstMin    = phiIndexMin;
     layout->firstMax    = phiIndexMax;
     layout->secondMin   = zIndexMin;
     layout->secondMax   = zIndexMax;
     const auto phiCount = static_cast<std::uint64_t>(phiIndexMax - phiIndexMin + 1);
     const auto zCount   = static_cast<std::uint64_t>(zIndexMax - zIndexMin + 1);
     if (phiCount > std::numeric_limits<std::uint64_t>::max() / zCount) {
       throw std::overflow_error("RandomNoisePixel cylindrical pixel count overflow");
     }
     layout->totalPixels = phiCount * zCount;
     return layout;
   }
 
   // Map one flat index in [0, totalPixels) to the two segmentation indices.
   // This gives every addressable pixel exactly the same selection probability.
   std::pair<std::int64_t, std::int64_t> pixelIndices(const RandomNoisePixel::PixelLayout& layout,
                                                      std::uint64_t linearIndex) {
     if (linearIndex >= layout.totalPixels) {
       throw std::out_of_range("RandomNoisePixel linear pixel index is out of range");
     }
     if (layout.rectangular) {
       const auto firstCount = static_cast<std::uint64_t>(layout.firstMax - layout.firstMin + 1);
       return {layout.firstMin + static_cast<std::int64_t>(linearIndex % firstCount),
               layout.secondMin + static_cast<std::int64_t>(linearIndex / firstCount)};
     }
 
     const auto row =
         std::upper_bound(layout.rows.begin(), layout.rows.end(), linearIndex,
                          [](std::uint64_t value, const RandomNoisePixel::PixelRow& candidate) {
                            return value < candidate.cumulativeEnd;
                          });
     if (row == layout.rows.end()) {
       throw std::out_of_range("RandomNoisePixel could not resolve compact pixel row");
     }
     const auto previousEnd = row == layout.rows.begin() ? 0 : (row - 1)->cumulativeEnd;
     return {row->firstMin + static_cast<std::int64_t>(linearIndex - previousEnd), row->secondIndex};
   }
 
 } // namespace
 
 // Build all static geometry and pixel-count metadata before event processing.
 void RandomNoisePixel::init() {
   // Step 1: reset state and validate the user-facing configuration.
   m_components.clear();
   m_layers.clear();
 
   if (!m_cfg.addNoise) {
     debug("RandomNoisePixel '{}': disabled by configuration; skipping geometry cache", name());
     return;
   }
   if (!(m_cfg.noise_rate_per_pixel_per_event >= 0.0 &&
         m_cfg.noise_rate_per_pixel_per_event <= 1.0)) {
     throw std::invalid_argument(
         "RandomNoisePixel noise_rate_per_pixel_per_event must be within [0, 1]");
   }
 
   // Step 2: CellIDPositionConverter::cellID(global) navigates shared TGeo state. JANA can
   // initialize detector factories concurrently, so serialize only cache construction.
   const std::scoped_lock geometryLock{geometryInitializationMutex()};
 
   // Step 3: obtain DD4hep services and the requested readout/segmentation.
   const auto& geo = algorithms::GeoSvc::instance();
   m_dd4hepGeo     = geo.detector();
   m_converter     = geo.cellIDPositionConverter();
   if (!m_dd4hepGeo || !m_converter) {
     throw std::runtime_error("RandomNoisePixel requires DD4hep geometry and cell-ID services");
   }
 
   m_readout = m_dd4hepGeo->readout(m_cfg.readout_name);
   if (!m_readout.isValid()) {
     throw std::invalid_argument("RandomNoisePixel invalid readout: " + m_cfg.readout_name);
   }
 
   // Step 4: collect geometry needed to identify each sensitive placement. This
   // vector owns the large TGeo transforms only during initialization.
   {
     std::vector<SensitiveComponentGeometry> componentGeometry;
     for (const auto& [name, detector] : m_dd4hepGeo->detectors()) {
       const auto sensitive = m_dd4hepGeo->sensitiveDetector(name);
       if (sensitive.isValid() && sensitive.readout().isValid() &&
           sensitive.readout().name() == m_cfg.readout_name) {
         collectDetectorComponents(detector, componentGeometry);
       }
     }
 
     // Step 5: convert initialization geometry into the compact event-time cache.
     cachePixelLayouts(componentGeometry);
   }
 
   // Step 6: the temporary transforms have now been released; construct layer totals.
   buildLayers();
   info("RandomNoisePixel '{}': cached {} sensitive components and {} layer groups for readout '{}'",
        name(), m_components.size(), m_layers.size(), m_cfg.readout_name);
 }
 
 // Enter each top-level layer of one detector system.
 void RandomNoisePixel::collectDetectorComponents(
     const dd4hep::DetElement& detector,
     std::vector<SensitiveComponentGeometry>& componentGeometry) {
   for (const auto& [_, layer] : detector.children()) {
     collectLayerComponents(detector.name(), layer, componentGeometry);
   }
 }
 
 // Collect every sensitive placement below one detector layer.
 void RandomNoisePixel::collectLayerComponents(
     const std::string& detectorName, const dd4hep::DetElement& layer,
     std::vector<SensitiveComponentGeometry>& componentGeometry) {
   for (const auto& [_, module] : layer.children()) {
     // Step 1: handle the common detector -> layer -> module hierarchy.
     const auto layerID = layerFromDetElement(layer, module);
     if (layerID) {
       auto components = findSensitiveComponents(detectorName, *layerID, module, m_cfg.readout_name);
       if (!components.empty()) {
         componentGeometry.insert(componentGeometry.end(),
                                  std::make_move_iterator(components.begin()),
                                  std::make_move_iterator(components.end()));
         continue;
       }
     }
 
     // Step 2: some geometries add one extra assembly level below the module.
     for (const auto& [__, child] : module.children()) {
       const auto childLayerID = layerID ? layerID : layerFromDetElement(layer, child);
       if (!childLayerID) {
         continue;
       }
       auto childComponents =
           findSensitiveComponents(detectorName, *childLayerID, child, m_cfg.readout_name);
       componentGeometry.insert(componentGeometry.end(),
                                std::make_move_iterator(childComponents.begin()),
                                std::make_move_iterator(childComponents.end()));
     }
   }
 }
 
 // Attach a base volume ID and a compact addressable-pixel layout to every component.
 void RandomNoisePixel::cachePixelLayouts(
     const std::vector<SensitiveComponentGeometry>& componentGeometry) {
   // Step 1: get the generic segmentation handle and the cell-ID bit-field encoder.
   const auto segmentation = m_readout.segmentation();
   const auto* decoder     = m_readout.idSpec().decoder();
   if (!segmentation.isValid() || !decoder) {
     throw std::runtime_error("RandomNoisePixel readout has no segmentation or ID decoder");
   }
 
   // Repeated Cartesian sensors share a logical TGeo volume and segmentation.
   // Store one immutable layout for all such placements to keep memory O(components).
   struct SharedLayout {
     const TGeoVolume* volume            = nullptr;
     const RawSegmentation* segmentation = nullptr;
     std::shared_ptr<const PixelLayout> layout;
   };
   std::vector<SharedLayout> sharedCartesianLayouts;
 
   m_components.reserve(componentGeometry.size());
   for (const auto& geometry : componentGeometry) {
     // Step 2: map one known interior point to a full DD4hep cell ID, then clear
     // the local pixel fields to obtain the physical sensor's base volume ID.
     const auto reference = referenceLocalPoint(geometry);
     double local[3]      = {reference[0], reference[1], reference[2]};
     if (!geometry.volume->GetShape()->Contains(local)) {
       throw std::runtime_error("RandomNoisePixel reference point is outside sensitive volume '" +
                                std::string{geometry.volume->GetName()} + "' with shape '" +
                                geometry.volume->GetShape()->ClassName() + "'");
     }
     double globalCoordinates[3];
     geometry.localToWorldTransform.LocalToMaster(local, globalCoordinates);
     const dd4hep::Position global{globalCoordinates[0], globalCoordinates[1], globalCoordinates[2]};
     const auto referenceCell = m_converter->cellID(global);
     if (!m_converter->findContext(referenceCell)) {
       throw std::runtime_error("RandomNoisePixel could not resolve a component reference cell ID");
     }
     const auto baseVolumeID = segmentation.volumeID(referenceCell);
     SensitiveComponent component;
     component.detectorName = geometry.detectorName;
     component.layer        = geometry.layer;
     component.baseVolumeID = baseVolumeID;
 
     // Step 3: select the concrete grid and reuse a Cartesian layout when possible.
     const auto& selected  = selectedSegmentation(segmentation, baseVolumeID);
     const auto findShared = [&]() -> std::shared_ptr<const PixelLayout> {
       for (const auto& cached : sharedCartesianLayouts) {
         if (cached.volume == geometry.volume && cached.segmentation == &selected) {
           return cached.layout;
         }
       }
       return {};
     };
 
     // Step 4: dispatch only to segmentation types explicitly supported by the SVT.
     if (const auto* grid =
             dynamic_cast<const dd4hep::DDSegmentation::CartesianGridXY*>(&selected)) {
       component.layout = findShared();
       if (!component.layout) {
         component.layout = makeCartesianLayout(
             geometry, GridKind::CartesianXY, grid->fieldNameX(), grid->fieldNameY(),
             grid->gridSizeX(), grid->gridSizeY(), grid->offsetX(), grid->offsetY(), *decoder);
         sharedCartesianLayouts.push_back({geometry.volume, &selected, component.layout});
       }
     } else if (const auto* grid =
                    dynamic_cast<const dd4hep::DDSegmentation::CartesianGridXZ*>(&selected)) {
       component.layout = findShared();
       if (!component.layout) {
         component.layout = makeCartesianLayout(
             geometry, GridKind::CartesianXZ, grid->fieldNameX(), grid->fieldNameZ(),
             grid->gridSizeX(), grid->gridSizeZ(), grid->offsetX(), grid->offsetZ(), *decoder);
         sharedCartesianLayouts.push_back({geometry.volume, &selected, component.layout});
       }
     } else if (const auto* grid =
                    dynamic_cast<const dd4hep::DDSegmentation::CylindricalGridPhiZ*>(&selected)) {
       component.layout = makeCylindricalLayout(geometry, *grid, *decoder);
     } else {
       throw std::runtime_error("RandomNoisePixel unsupported segmentation type: " +
                                selected.type());
     }
 
     // Step 5: cache the count, then choose representative addresses across the
     // layout. Testing the edges as well as the middle catches placement/frame
     // mistakes that a single central pixel could hide.
     component.pixelCount                             = component.layout->totalPixels;
     std::array<std::uint64_t, 5> sampleLinearIndices = {
         0, component.pixelCount / 4, component.pixelCount / 2,
         component.pixelCount / 2 + component.pixelCount / 4, component.pixelCount - 1};
     std::sort(sampleLinearIndices.begin(), sampleLinearIndices.end());
     const auto uniqueEnd = std::unique(sampleLinearIndices.begin(), sampleLinearIndices.end());
 
     // Step 6: round-trip each cell ID through DD4hep and verify that its center
     // returns inside this exact placed sensitive solid. DD4hep applies the
     // placement once when producing the global position; MasterToLocal removes
     // that same complete (possibly nested) translation and rotation here.
     for (auto sample = sampleLinearIndices.begin(); sample != uniqueEnd; ++sample) {
       const auto [sampleFirst, sampleSecond] = pixelIndices(*component.layout, *sample);
       auto sampleID                          = component.baseVolumeID;
       decoder->set(sampleID, component.layout->firstField, sampleFirst);
       decoder->set(sampleID, component.layout->secondField, sampleSecond);
       if (!m_converter->findContext(sampleID)) {
         throw std::runtime_error("RandomNoisePixel generated an invalid sample cell ID");
       }
 
       const auto sampleGlobal           = m_converter->position(sampleID);
       double sampleGlobalCoordinates[3] = {sampleGlobal.x(), sampleGlobal.y(), sampleGlobal.z()};
       double sampleLocalCoordinates[3];
       geometry.localToWorldTransform.MasterToLocal(sampleGlobalCoordinates, sampleLocalCoordinates);
       if (!containsOrTouches(*geometry.volume->GetShape(), sampleLocalCoordinates)) {
         throw std::runtime_error(
             "RandomNoisePixel sample cell center does not return inside sensitive volume '" +
             std::string{geometry.volume->GetName()} + "' in detector '" + component.detectorName +
             "' layer " + std::to_string(component.layer) + "; pixel indices=(" +
             std::to_string(sampleFirst) + ", " + std::to_string(sampleSecond) +
             "), returned local position=(" + std::to_string(sampleLocalCoordinates[0]) + ", " +
             std::to_string(sampleLocalCoordinates[1]) + ", " +
             std::to_string(sampleLocalCoordinates[2]) +
             "); check segmentation coordinates and placement transforms");
       }
     }
     m_components.push_back(std::move(component));
   }
 }
 
 // Group components by detector and layer and prepare weighted random selection.
 void RandomNoisePixel::buildLayers() {
   for (std::size_t componentIndex = 0; componentIndex < m_components.size(); ++componentIndex) {
     const auto& component = m_components[componentIndex];
     auto layer =
         std::find_if(m_layers.begin(), m_layers.end(), [&](const LayerGeometry& candidate) {
           return candidate.detectorName == component.detectorName &&
                  candidate.layer == component.layer;
         });
     if (layer == m_layers.end()) {
       m_layers.push_back({component.detectorName, component.layer, {}, {}, 0});
       layer = std::prev(m_layers.end());
     }
     if (layer->totalPixels > std::numeric_limits<std::uint64_t>::max() - component.pixelCount) {
       throw std::overflow_error("RandomNoisePixel layer pixel count overflow");
     }
     // Appending cumulative totals lets upper_bound select a component with
     // probability N_component / N_layer without visiting all components/event.
     layer->totalPixels += component.pixelCount;
     layer->componentIndices.push_back(componentIndex);
     layer->cumulativePixels.push_back(layer->totalPixels);
   }
 
   for (const auto& layer : m_layers) {
     info("RandomNoisePixel geometry: detector='{}' layer={} components={} pixels={} "
          "expected_hits={}",
          layer.detectorName, layer.layer, layer.componentIndices.size(), layer.totalPixels,
          m_cfg.noise_rate_per_pixel_per_event * static_cast<double>(layer.totalPixels));
   }
 }
 
 // Derive the event RNG seed from the run/event identity and algorithm name.
 std::uint64_t
 RandomNoisePixel::seedFromEventHeader(const edm4hep::EventHeaderCollection& headers) const {
   if (headers.empty()) {
     throw std::runtime_error("RandomNoisePixel requires a non-empty EventHeader collection");
   }
   return m_uid.getUniqueID(headers, name());
 }
 
 // Select one pixel uniformly within a sensitive component and encode its fields.
 std::uint64_t RandomNoisePixel::randomCellID(const SensitiveComponent& component,
                                              std::mt19937_64& rng) const {
   std::uniform_int_distribution<std::uint64_t> pickPixel(0, component.pixelCount - 1);
   const auto [firstIndex, secondIndex] = pixelIndices(*component.layout, pickPixel(rng));
   auto cellID                          = component.baseVolumeID;
   const auto* decoder                  = m_readout.idSpec().decoder();
   decoder->set(cellID, component.layout->firstField, firstIndex);
   decoder->set(cellID, component.layout->secondField, secondIndex);
   return cellID;
 }
 
 // Generate noise for one layer using lambda = rate_per_pixel_per_event * N_pixels.
 void RandomNoisePixel::addNoiseHitsForLayer(
     const LayerGeometry& layer, std::map<std::uint64_t, edm4eic::MutableRawTrackerHit>& hitMap,
     std::mt19937_64& rng) const {
   if (layer.totalPixels == 0 || m_cfg.noise_rate_per_pixel_per_event == 0.0) {
     return;
   }
 
   // Step 1: draw the total layer occupancy from the sparse-noise Poisson model.
   const double mean = m_cfg.noise_rate_per_pixel_per_event * static_cast<double>(layer.totalPixels);
   std::poisson_distribution<std::uint64_t> poisson(mean);
   const auto requested   = poisson(rng);
   const auto maxAttempts = std::max<std::uint64_t>(100, requested * 20);
   std::uniform_int_distribution<std::uint64_t> pickLayerPixel(0, layer.totalPixels - 1);
 
   std::uint64_t created  = 0;
   std::uint64_t attempts = 0;
   // Step 2: select a uniform layer-wide pixel. The cumulative component counts
   // first choose the sensor with probability proportional to its pixel count.
   while (created < requested && attempts < maxAttempts) {
     ++attempts;
     const auto selectedPixel     = pickLayerPixel(rng);
     const auto componentPosition = std::upper_bound(layer.cumulativePixels.begin(),
                                                     layer.cumulativePixels.end(), selectedPixel);
     const auto componentOffset =
         static_cast<std::size_t>(componentPosition - layer.cumulativePixels.begin());
     const auto& component = m_components[layer.componentIndices[componentOffset]];
     const auto cellID     = randomCellID(component, rng);
     // Step 3: electronic pixels can fire at most once in an event. At the default
     // occupancy, retries due to duplicates are extremely rare.
     if (hitMap.contains(cellID)) {
       continue;
     }
 
     // This algorithm models pixel occupancy, not the sensor pulse shape.
     edm4eic::MutableRawTrackerHit hit;
     hit.setCellID(cellID);
     hit.setCharge(1.0e6);
     hit.setTimeStamp(0);
     hitMap.emplace(cellID, hit);
     ++created;
   }
 
   if (created < requested) {
     warning("RandomNoisePixel '{}': created {}/{} requested hits for detector '{}' layer {}",
             name(), created, requested, layer.detectorName, layer.layer);
   }
 }
 
 // Event entry point: seed, sample every layer, and emit deterministic cell-ID order.
 void RandomNoisePixel::process(const Input& in, const Output& out) const {
   auto [outHits]       = out;
   const auto [headers] = in;
 
   if (!m_cfg.addNoise) {
     return;
   }
   if (!headers) {
     throw std::runtime_error("RandomNoisePixel requires an EventHeader collection");
   }
   if (m_layers.empty()) {
     throw std::runtime_error("RandomNoisePixel has no cached sensitive geometry");
   }
 
   // Step 1: make the random sequence depend on event identity, not thread scheduling.
   std::mt19937_64 rng(seedFromEventHeader(*headers));
 
   // Step 2: a sorted map removes duplicates and gives stable output ordering.
   std::map<std::uint64_t, edm4eic::MutableRawTrackerHit> noiseHits;
   for (const auto& layer : m_layers) {
     addNoiseHitsForLayer(layer, noiseHits, rng);
   }
   // Step 3: transfer the completed noise-only collection to PODIO.
   for (const auto& [_, hit] : noiseHits) {
     outHits->push_back(hit);
   }
 }
 
 } // namespace eicrecon

This page was automatically generated by the 2.3.7 LXR engine. The LXR team