Repository code AI-audit

[cuddandr]

The final markdown report from Codex for reference. We should go through it and check the validity of its findings. Several of these are on completely out-of-date code and not directly applicable.

Physics and unintended-behavior audit of develop

Audit target: develop at 3b6449466e1e8036413ad9c6750b04a68515aea3

Audit date: 2026-08-28

Executive summary

This audit found several high-confidence correctness defects in the charge,
light, and optional far-field paths. The most consequential default-path issue
is that light waveform digitization does not use the trigger time: all threshold
triggers digitize from the beginning of the simulated buffer. The previously
noted light-truth indexing concern is also confirmed and is broader than one
wrong index: the detector-response convolution tests source IDs in the input
array at the output tick, writes IDs from the output tick instead of the
convolution tick, and never tests the destination array when looking for an
existing slot.

The active-volume and pixel traversal code also disagree about partially
contained segments. A segment with one endpoint inside is retained without
being clipped, while a segment that crosses the volume with both endpoints
outside is lost. The retained, unclipped case can make get_active_pixels
write valid pixels at indices beyond an array sized only for the count of valid
pixels.

The optional far-field implementation should not be used for physics results
without fixes and validation. It loses per-segment t0, the segment mode
removes charge in infinite x/y strips because of an or condition, and voxel
sizes documented as 0.5 cm evaluate to 5 cm in the coordinate system used by
the simulation.

Priority meanings:

• P1: can substantially change simulated charge/light, trigger acceptance,
  waveform timing, or cause invalid GPU memory access.
• P2: conditional or localized correctness defect, usually affecting truth,
  edge cases, configuration variants, or high-occupancy events.
• P3: low-impact defect or latent unit/configuration hazard.

Prioritized findings

Priority	File	Current likely unintended behavior	Suggested fix
P1	larndsim/light_sim.py:530-539	digitize_signal ignores both trigger_idx and the pre-trigger window, so every threshold trigger reads the same samples from tick zero rather than a window around that trigger.	Restore a trigger-relative sample index, using the already padded trigger index: sample_tick = isample * spacing / tick_size - pre_window / tick_size + trigger_idx[itrig]; add single- and multi-trigger waveform tests.
P1	larndsim/light_sim.py:355-369	Light detector-response truth uses itick where the source is jtick, consults light_sample_inc_true_track_id where it should consult the destination light_response_true_track_id, and writes the wrong source ID. Signal amplitudes are unaffected, but waveform truth is corrupted or dropped.	Cache source_id = light_sample_inc_true_track_id[idet, jtick, itrue]; match it against light_response_true_track_id[idet, itick, jtrue], then write source_id and accumulate the (jtick, itrue) photons.
P1	larndsim/light_sim.py:462-473	After the second threshold trigger, the dead-time loop slices an already shortened mask with an absolute tick. It suppresses progressively more time than digit_ticks and can miss valid later triggers.	Select from absolute candidate ticks with next_allowed_tick, or slice by the relative offset. The correction already present on the current light-trigger optimization topic branch is an appropriate model.
P1	larndsim/active_volume.py:24-46; larndsim/drifting.py:34-45	“At least partially contained” is implemented as “either endpoint strictly inside.” Crossing segments with both endpoints outside are lost; boundary-equal endpoints are lost; retained segments are not clipped or energy-apportioned; plane assignment then uses only the segment midpoint.	Use segment/AABB intersection, clip at TPC boundaries, and split/apportion dE, length, and timing when a segment crosses TPCs. Make boundary inclusivity consistent with drift.
P1	larndsim/pixels_from_track.py:139-235	max_pixels allocates by the number of in-bounds pixels, but get_active_pixels increments its output index for out-of-bounds traversal steps. A retained segment with one endpoint outside can write later in-bounds pixels beyond the allocated row.	Compact output indices only when a pixel is in bounds, or clip the segment before traversal. Add guard assertions/counters and tests with each endpoint outside each detector face.
P1	larndsim/fee.py:583-590	num_backtrack[ip] is read before checking ip < pixels_signals.shape[0]. Every launch whose pixel count is not a multiple of the block size creates out-of-range threads and an illegal/undefined read.	Move calculation of ntrks and all other ip-indexed reads inside the bounds guard, preferably with an early return.
P1 (far-field)	larndsim/far_field/signal_calculation.py:179-241	Segment-mode exclusion skips a sub-piece when either dx or dy is within the near-field radius. This removes far-field current in two infinite strips, even if the 2-D distance is large.	Use and for the intended square exclusion, or hypot(dx, dy) <= radius for a radial exclusion. The same or was changed to and in non-develop commit 3707a67.
P1 (far-field)	larndsim/far_field/signal_calculation.py:59,172; cli/simulate_pixels.py:1482-1491,1549-1577	Far-field kernels never use segment t0. Near-field pixels receive the event-wide far-field signal starting at tick zero; induction-only pixels receive the aggregate shifted by one min(t0). Tracks with different deposit times are therefore overlaid at the wrong times.	In segment mode, use t - segment['t0'] per segment and skip negative relative time. In voxel mode, retain a time dimension or avoid combining deposits with different t0. Size the output from event-wide timing, not the last near-field subbatch.
P1 (far-field)	larndsim/consts/ff_induction.py:8-14; larndsim/far_field/voxelization.py:98-113	Coordinates passed to voxelization are numeric centimetres, but 0.5 * cm evaluates to 5 in the base unit module. Voxels documented as 0.5 cm are therefore 5 cm on each axis.	Store 0.5 for values consumed in centimetres, or define and consistently apply an explicit internal unit convention. Add a test asserting the grid shape for a known 10 cm box.
P1/P2 truth	larndsim/light_sim.py:654-684; cli/simulate_pixels.py:1688-1709	i_trig advances once per file flush, not once per exported waveform. If an event has multiple threshold triggers, later light_wvfm_mc_assn.trigger_id values overlap earlier rows. A multi-event flush also labels every truth row with event_id[0].	Derive the starting trigger row from the waveform dataset length or increment by the number of exported triggers; assign event IDs with the trigger-axis index idx0.
P2	larndsim/pixels_from_track.py:271-296; larndsim/detsim.py:334-356	Duplicate neighboring pixels keep the distance from the first active pixel encountered, not the minimum distance. A pixel directly crossed later in a segment can remain labeled distance 1 or 2, making track retention direction-dependent when MAX_TRACKS_PER_PIXEL is reached. Since unretained tracks are not summed, this can affect charge as well as truth in crowded pixels.	On duplicates, update the stored radius with the minimum; alternatively generate the unique neighborhood and minimum distance deterministically.
P2	larndsim/detsim.py:256-276	The bounds check validates raw itick, but both outputs are indexed with itick + track_t0[itrk]. Negative t0 can wrap/write before the intended window, and the flattened truth index is calculated before checking the shifted tick or pixel_index.	Compute shifted_tick only after validating pixel_index, then require 0 <= shifted_tick < pixels_signals.shape[1] before calculating either destination index.
P2	cli/simulate_pixels.py:841-866	The code explicitly supports event IDs offset by multiples of MAX_EVENTS_PER_FILE, but non-spill vertex timestamps index event_times with raw uniq_ev. Globally offset IDs can index out of range or select the wrong timestamp.	Index with uniq_ev % MAX_EVENTS_PER_FILE, as the surrounding event-time logic and save path already do. Validate gaps and nonzero file offsets.
P2	larndsim/fee.py:599-617	Periodic-reset phase is sampled from 0..PERIODIC_RESET_CYCLES inclusive, but ic % PERIODIC_RESET_CYCLES can never equal the last value. About 1/(cycles+1) of pixels never reset (for example, 1/257 in FSD).	Sample int(uniform * PERIODIC_RESET_CYCLES) so the phase is in [0, cycles-1].
P2 truth	larndsim/fee.py:313-360	Segment associations are sorted by contribution, but trajectory contributions are regrouped in ascending trajectory-ID order and then truncated without re-sorting. High-contribution trajectories can be omitted from the stored top-N trajectory truth.	After summing by trajectory ID, sort the grouped trajectory fractions descending before truncation and storing.
P2	larndsim/light_sim.py:91-102	LUT time-profile bins are hard-coded to 1 ns. A LUT with another bin width produces the wrong propagation delay and profile duration. This is especially risky because LUT smearing is enabled in FSD/FSD-Cube/ND-LAr detector files.	Store/read the profile bin width in LUT metadata and pass it into the kernel; reject smearing when metadata is absent or inconsistent.
P2	larndsim/light_sim.py:103-104,120-123	Both ends of the light tick interval are strict. A photon landing exactly on a tick boundary belongs to neither adjacent tick and is dropped. Integer-ns profiles make exact boundaries plausible.	Use one half-open convention consistently, normally start_tick_time <= profile_time < end_tick_time, and test exact boundaries.
P2 truth	larndsim/light_sim.py:532-571	digitize_signal correctly resolves global optical channel idet to signal row idet_signal, but reads signal_true_photons[idet, ...] in one branch. Non-identity or reordered channel maps read another row or go out of bounds.	Use idet_signal for all signal/truth array reads; reserve idet for channel metadata only.
P2 config	larndsim/consts/light.py:111-117,171-184	Channel divisibility is checked before loading the configured OP_CHANNEL_PER_TRIG. A two-value threshold list is tiled to twice the required number of groups. Custom group sizes can therefore pass/fail the wrong validation and threshold array shape is inconsistent.	Load group size first, validate against it, and tile two thresholds by N_OP_CHANNEL // (2 * OP_CHANNEL_PER_TRIG); assert the final threshold length exactly.
P2 config	larndsim/detector_properties/fsd.yaml:29,40; fsd_cube.yaml:39,50; ndlar-module.yaml:114,124	Each file defines light_trig_mode twice. PyYAML silently keeps the last value; in fsd.yaml that changes the effective mode from 1 to 0.	Remove duplicate keys and use a loader/lint check that rejects duplicates. Document the intended mode once.
P3	larndsim/consts/units.py:117-125	gigahertz is defined as 1.e+6 * hertz, identical to megahertz rather than 1.e+9 * hertz.	Correct the multiplier and add dimensional ratio tests.

Detailed findings by folder and file

larndsim/light_sim.py

Trigger-relative waveform extraction is disabled

digitize_signal contains the correct trigger-relative expression as a
comment at line 530, but the live expression at line 531 is only a function of
isample. sim_triggers goes to considerable effort to pad the signal and to
produce padded_trigger_idx, then passes that array to a kernel that does not
use it.

Consequences:

• a beam trigger at tick zero happens to work because front padding places the
  desired window at buffer index zero;
• a threshold trigger later in the buffer reads from buffer index zero rather
  than from its trigger;
• every trigger in a multi-trigger event returns the same underlying time
  window (apart from channel choice/noise mutation).

This should be fixed before using threshold-trigger waveforms for physics.

Detector-response truth convolution uses the wrong arrays and ticks

For each convolution source (jtick, itrue), lines 366-368 instead compare and
write IDs taken from (itick, jtrue)/(itick, itrue) in the input truth
array. The destination array is never consulted for slot matching. A minimal
correct pattern is:

source_id = light_sample_inc_true_track_id[idet, jtick, itrue]
for jtrue in range(light_response_true_track_id.shape[-1]):
    destination_id = light_response_true_track_id[idet, itick, jtrue]
    if destination_id == source_id or destination_id == -1:
        light_response_true_track_id[idet, itick, jtrue] = source_id
        light_response_true_photons[idet, itick, jtrue] += (
            tick_weight * light_sample_inc_true_photons[idet, jtick, itrue]
        )
        break

This confirms the earlier concern and shows why changing only one index would
not be sufficient.

Threshold dead-time selection uses mixed coordinate systems

After the first accepted trigger, module_above_thresh is a sliced, relative
view. next_idx is converted back to an absolute tick by adding
last_trigger, but that absolute value is then used to slice the already
relative array. From the second trigger onward, the start position advances by
more than one readout window. Absolute candidate ticks plus a single
next_allowed_tick avoids this class of error.

Truth export IDs do not track waveform rows

zero_suppress_waveform_truth treats i_evt as one scalar and assigns it to
every row. Its caller supplies event_id[0]. Meanwhile, i_trig is a flush
counter even though the exported array can have more than one trigger. These
two assumptions are only valid when every flush contains exactly one event and
exactly one trigger.

Additional light concerns

• Exact tick-boundary photons are excluded at both ends of each interval.
• One truth read uses global idet instead of resolved row idet_signal.
• LUT profile spacing is fixed at 1 ns instead of coming from the LUT.
• get_nticks silently caps its caller's light buffer at 50,000 ticks in
  simulate_pixels.py:1601-1602; the cap should at least warn when it truncates
  a requested window.

larndsim/active_volume.py and larndsim/drifting.py

The selector's documented contract and its implementation differ. It does not
perform a line-box intersection; it tests only whether either endpoint is
strictly interior. It also creates unused tpc_start_mask and tpc_end_mask
arrays, which suggests the endpoint logic may have been left partially
refactored.

The downstream drift kernel assigns a plane using the segment midpoint. This
creates inconsistent cases:

1. one endpoint is inside, so the selector retains the segment;
2. its midpoint is outside, so drift leaves the default plane;
3. the batcher marks the segment as simulated in the first matching batch;
4. the main loop later filters it from charge and light response.

Clipping/splitting at the geometry boundary is preferable to adding more
tolerances because the deposited energy and segment length must remain
consistent.

larndsim/pixels_from_track.py

get_num_active_pixels increments only for valid pixels. get_active_pixels,
however, increments i on every traversal step and writes a valid pixel at
that un-compacted index. The caller allocates using the former count. This is a
direct allocation/write contract violation for partially out-of-bounds
segments.

The neighboring radius is also documented as distance from the nearest active
pixel, but duplicate suppression never updates an existing distance. For a
track whose active pixels are visited [A, B], pixel B can first be inserted
as a neighbor of A at distance 1 and remain distance 1 when B is later
visited as an active pixel. Reversing track direction can change which pixels
receive distance zero.

larndsim/detsim.py

sum_pixel_signals must validate the shifted tick, not the unshifted kernel
index. This is particularly important because input segments are only bounded
above by MAX_SEGMENT_T0; negative t0 is not rejected. Calculating
base_idx before validating pixel_index also performs an unnecessary read of
offset_backtrack[-1] for sentinel pixels.

There is a separate low-frequency edge case in tracks_current_mc: segment
direction is normalized before checking that the full segment length is
nonzero (detsim.py:158-162). Point-like deposits with positive charge can
produce divisions by zero/NaNs. A point-deposit branch should place the full
charge at the point rather than dropping or normalizing a zero vector.

larndsim/fee.py

GPU bounds guard occurs too late

The kernel block size is four in the main caller, so an arbitrary number of
pixels normally launches extra threads. Reading num_backtrack[ip] before the
guard is therefore not a theoretical case. The read should occur after an
early return for out-of-range ip.

Periodic reset leaves a subset of pixels with no resets

The random phase expression multiplies by PERIODIC_RESET_CYCLES + 1. The
modulo expression used thereafter has only PERIODIC_RESET_CYCLES possible
values. The last sampled phase is unreachable. This path is active in the FSD
and FSD-Cube detector properties.

Trajectory truth is not top-N after regrouping

Segment fractions are correctly sorted before segment truth is stored.
Trajectory truth then groups those segments using np.unique, which returns
sorted IDs, and writes the grouped values in that ID order. The subsequent
slice to ASSOCIATION_COUNT_TO_STORE is consequently not a top-contribution
slice.

For WRITE_BATCH_SIZE > 1, line 285 also creates per-message timestamp packets
using event_start_time_list[0] rather than the current row's event. Default
simulation files use a write batch of one, so this is a conditional P2 issue.

larndsim/far_field/

The far-field feature is disabled in the checked-in simulation-property files,
so these defects do not affect the default configurations. They are P1 for any
run that enables the feature.

signal_calculation.py

The segment exclusion dx <= radius or dy <= radius removes sources far from a
pixel whenever only one coordinate happens to align. An unmerged branch already
contains the direct and correction, supporting the conclusion that the
develop behavior is unintended.

Neither far-field kernel reads track t0. Voxel mode cannot recover it because
charges from all segments are spatially aggregated before signal calculation.
The main loop's single minimum-time shift for induction-only pixels cannot
represent multiple deposits at different times.

consts/ff_induction.py and voxelization.py

The rest of the detector code converts geometry into numeric centimetres; for
example, pixel pitch is loaded as yaml_mm * mm / cm. Far-field voxelization
directly subtracts those centimetre coordinates and divides by
COARSE_VOXEL_SIZE_*. In that context 0.5 * cm is numerically 5, so the
current coarse grid is ten times coarser along each dimension than documented.

cli/simulate_pixels.py

In addition to the far-field timing and light truth-counter issues above, the
non-spill vertex timestamp path fails to apply the event-ID modulus that its
own preceding comments require. This should be covered with an input whose
event IDs begin at MAX_EVENTS_PER_FILE and whose local IDs contain gaps.

The final beam-trigger export has another gap-sensitive assumption:
light_event_times is the full local event-time array, while
export_light_trig_to_hdf5 indexes it with the compact inverse of the observed
event IDs. For observed events [0, 2], event 2 receives event_times[1]
instead of event_times[2]. Pass an explicitly aligned compact time array or a
mapping keyed by event ID.

larndsim/consts/light.py and detector property YAML

OP_CHANNEL_PER_TRIG is read after divisibility checks, so the checks use the
previous/default value. The special two-threshold expansion creates twice the
expected threshold entries. It happens to be partially masked by later channel
indexing in current layouts, but it is not a sound configuration contract.

Duplicate YAML keys were found with a duplicate-rejecting loader in:

• larndsim/detector_properties/fsd.yaml
• larndsim/detector_properties/fsd_cube.yaml
• larndsim/detector_properties/ndlar-module.yaml

Only FSD's duplicate changes the value (1 followed by 0), but all should be
removed and prevented in CI.

larndsim/quenching.py

One physics-input issue needs domain confirmation rather than an immediate
code-only patch. The kernel uses the full segment dE as visible energy, while
the existing TODO notes that edep-sim's secondary/non-ionizing energy may need
to be subtracted. If dE includes SecondaryEnergyDeposit for the supported
input schema, both electrons and photons are overproduced. Confirm the edep-sim
schema/version used in production, document the contract, and add a fixture
where total and visible deposit differ.

There is also a code/test mismatch at dEdx == 0: the kernel returns zero
charge for both models because it skips the entire calculation, while
tests/testQuenching.py expects the Birks zero-density limit
BIRKS_Ab * dE / W_ION. A segment with dE > 0 and dEdx == 0 may be invalid
input, but the implementation and test should agree explicitly.

larndsim/consts/units.py

gigahertz is a straightforward factor-of-1000 typo. No use of GHZ was found
in the audited tree, so this is latent rather than an active simulation error.

Test and validation gaps

The audit attempted the checked-in test entry points against a clean archive of
develop:

• python -m pytest -q tests reported no tests ran because the files are
  named testDrifting.py, testQuenching.py, etc., which do not match pytest's
  default test_*.py convention.
• python -m unittest discover -s tests -p 'test*.py' -v found the files but
  could not import them in the audit environment because CuPy is unavailable.
• full pytest collection also treats three far-field performance scripts under
  examples/ as tests and fails at import without CuPy.
• several checked-in tests call APIs/signatures that no longer exist, including
  detsim.rho, detsim.tracks_current, and an older get_pixels signature.
• the only checked-in GitHub workflow that runs larnd-sim is a profiling job
  after a pull request is closed; there is no normal pull-request correctness
  test job.

Recommended minimum regression matrix:

1. rename/update the unit tests so pytest tests collects them;
2. run CUDA simulator-compatible kernel tests on CPU where practical and a
   small real-GPU suite for memory-safety/indexing;
3. add exact tests for light tick boundaries, one late threshold trigger, three
   threshold triggers with dead time, reordered optical channels, and more than
   one event per write batch;
4. add geometry tests for segments wholly inside, one endpoint outside, both
   endpoints outside but crossing, boundary equality, and cross-TPC splitting;
5. add far-field invariants for voxel dimensions, per-segment time translation,
   x/y exclusion symmetry, and charge conservation during voxelization;
6. reject duplicate YAML keys during configuration tests.

Audit limitations

This was a static/control-flow audit plus configuration and repository-history
cross-checking. The audit environment did not have CuPy or a CUDA device, so no
end-to-end GPU simulation was run. Priorities reflect likely result impact, not
measured event rates. The report deliberately separates direct indexing/control
flow defects from the visible-energy question in quenching.py, which requires
confirmation of the production input schema.