Support pulse-height tallies in shared secondary mode

[GuySten]

Support pulse-height tallies with the shared secondary bank

Removes the restriction that forced shared_secondary_bank off whenever a
pulse-height tally was present, and implements the deferred per-history scoring
described in the TODO that guard carried.

Problem

A pulse-height tally scores once per source history, in the bin containing the
total energy that history's entire particle tree deposited in a cell. In the
default transport modes the whole tree is carried by a single Particle, so
pht_storage_ already holds the per-history total when event_death() runs.

Under the shared secondary bank each generation is transported as a fresh set of
Particle objects, load balanced across MPI ranks between generations. Three
things broke:

1. pht_storage_ is a per-Particle member reset in initialize_particle_track, so every descendant's deposition was dropped from its root's total.
2. event_death() scored unconditionally, so each fragment of one physical pulse was scored as a separate pulse, shredding the spectrum toward low energy.
3. The subtraction removing a secondary's energy from its parent was never implemented on this path. event_revive_from_secondary() skipped it on the grounds that create_secondary() had already done it; create_secondary() contained no pulse-height code. The compensating pre-add in sample_secondary_photons() still ran unconditionally.

Approach

Every particle carries a root index: the dense global index of its primary
in the source bank, in [0, n_particles). It has to survive bank sorting and
the MPI migration between generations, so it travels on SourceSite.

Each Particle still accumulates its own fragment as before. At death the
fragment is staged in a per-thread buffer keyed by root index rather than
scored. Once the generation loop drains, contributions are routed by
MPI_Alltoallv to the rank owning each root, summed, and scored once per
history.

Carrying the root without growing SourceSite

SourceSite is the element type of the fission bank and the source bank, both
sized by particles per batch, so an added field would be paid for in every run
including those that never touch the shared secondary bank. The existing
parent_id field is already free by the time the root is wanted: it and
progeny_id serve only to give a site its place in the collected bank, via
progeny_per_particle[parent_id] + progeny_id, and that key is consumed the
moment the site is placed.

The field is therefore renamed to ancestor_index and given two accessors over
the same storage:

• parent_slot() — the immediate parent's slot in the current generation's work, valid from creation until collection places the site.
• root_index() — the primary at the root of the history, valid afterwards.
  resolve_root_indices() performs the handover, once per generation, after
  collection and before the migration. A site's parent is a primary on the first
  pass, so its root comes from the phase-1 partition; on later passes the parent
  already carries its own resolved root, so the value propagates down the tree.
  Net effect on SourceSite is zero fields and zero bytes, and the MPI datatype
  is unchanged in size.

progeny_id stays: fission sites are appended concurrently, so their position
in the bank is nondeterministic and the ordinal is the only thing sort_bank()
has to work from. add_surf_source_to_bank() wrote a particle ID into one field
and a progeny count into the other, neither matching what those fields mean
anywhere else and neither read by anything, so both assignments are removed.

The phase-1 partition is computed, not stored

calculate_work() is a pure function of the primary count and the rank count —
a block partition with the remainder given to the first remainder ranks — so
the phase-1 partition can be recomputed wherever it is needed rather than
snapshotted before calculate_work() starts rewriting work_index for each
secondary generation. phase1_first_root() and phase1_owner_of_root() give
the range owned by a rank and the owner of a root in closed form. When
min_work is zero the boundary equals n_particles, so the branch that would
divide by zero is unreachable.

Reproducibility

Arrival order at a history's accumulator depends on which thread transported
each descendant and on which rank it landed, so summing fragments as they are
found would make a history's total vary in its last bits with the thread and
rank count. For most tallies that is invisible, but a pulse height is binned, so
an ulp can move a whole count across a bin edge. The regression test here is
exactly the exposed case: a delta_function(1e6) source with bin edges at
linspace(0, 1e6, 101), so every fully absorbed history lands on the top edge
and a one-ulp difference decides whether it is counted at all.

Contributions are therefore collected, not summed, and folded in track id order
at the end. A track id is the particle's global slot within its generation plus
the tracks completed in earlier generations, both global quantities, so the
order is the same for any thread or rank count, and it is unique within a batch
so the ordering needs no tie-break. This relies on the global ordering of the
secondary bank being independent of the rank count, which is the property
sort_bank() already exists to provide.

The cost is one sort per batch over the fragments that actually deposited,
which in a detector problem is a small fraction of tracks.

Changes

• particle_data.h — SourceSite::parent_id renamed to ancestor_index, with parent_slot() and root_index() accessors; ParticleData::root_index_ and accessors
• particle.cpp — propagate the root index; implement the parent-side subtraction in create_secondary(); route event_death() to staging in shared mode; drop the two dead bank-ordering writes in add_surf_source_to_bank()
• tallies/pulse_height.{h,cpp} (new) — staging buffers, owner exchange, canonical fold order, deferred scoring
• simulation.{h,cpp} — phase1_first_root(), phase1_owner_of_root(), resolve_root_indices(); hooks in both shared drivers
• tally_scoring.cpp — score_pulse_height_tally() takes the energy vector as an argument instead of reading p.pht_storage()
• initialize.cpp — remove check_pulse_height_compatibility(); mpi::source_site loses a field
• lib/core.py — _SourceSite mirrors the renamed member

Testing

New tests in tests/unit_tests/test_pulse_height.py:

• Count conservation (photon/neutron × shared/local). Fixed-source results are normalized per source particle, so the tally summed over all bins must be exactly 1.0. Per-track scoring pushes this above one; a dropped history pushes it below.
• Mode comparison, over a thin photon, thin neutron and thick photon detector. Exact comparison is not possible, since compute_particle_id() and compute_transport_seed() both branch on the shared bank and the two modes therefore sample different random number streams. Mean deposited energy and per-bin shape are compared against combined standard errors. The thick case is the sensitive one: nothing escapes and each history spawns roughly 1.7 secondaries.
  Regression references for shared_photon and shared_neutron are regenerated;
  the local_* references are unchanged.

Notes for reviewers

• SourceSite::parent_id is renamed. It is a public member of an installed header, so out-of-tree code touching it — most plausibly a compiled custom source — needs the new name. CollisionTrackSite::parent_id is a different struct and is untouched.
• Two accessors over one field give no compile-time protection. Calling root_index() on an uncollected site compiles and returns the placement key. The window is about ten lines in each driver, between collection and resolve_root_indices(). If reviewers would rather have it enforced, storing the placement key negated makes the sign a free discriminant and lets each accessor assert; that is a small follow-up.
• One allocation per secondary generation. The history-based driver used to recycle the read bank's buffer into the write bank before collection. The read bank is now needed until roots are resolved, so it is released afterwards instead.
• Memory. Staging holds one record per track that actually deposited, not per track, plus one int64 track id each; in a typical detector problem most histories never reach a pulse-height cell. The per-rank totals array is 8 * n_owned_roots * n_pulse_height_cells bytes.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)