Add high-level access to MF=8/9/10 radionuclide production data

[jon-proximafusion]

Proposal (with complete code below, developed and tested): a small interpreted layer over the existing MF=8/9/10 parsers. parse_mf8 and parse_mf9_mf10 already expose radionuclide production data as raw dicts, but every consumer has to re-implement the same join: walking the MF=8 subsections to attach excitation energies to each (ZAP, LFS) final state, and pairing the MF=9 yield with the MF=10 cross section for the same state.

Motivating consumer: OpenMC is gaining the ability to store this data verbatim on depletion chain files (shimwell/openmc#121, shimwell/openmc#122), and currently carries its own copy of this join. Activation tooling generally needs the same thing. The docstring deliberately records the classic pitfall that LFS is a level index and not an isomeric-state index, so consumers know that mapping a level to a metastable state requires comparing excitation energies against decay data.

If merged here (or upstream in paulromano/endf-python, where this would fit naturally since it is pure format interpretation with no policy), the OpenMC side can shrink its extraction module to a thin adapter.

One design note: Material.__init__ eagerly parses every section, so consumers processing thousands of files for just MF=8/9/10 may want to pre-filter by section presence before constructing Material objects. A ready branch with the commit exists locally (patch also at Jon's ~/0001-Add-high-level-access-to-MF-8-9-10-radionuclide-prod.patch); everything needed to apply it by hand is below.

New module src/endf/radionuclide_production.py

from dataclasses import dataclass
from typing import Dict, List, Optional

from .function import Tabulated1D
from .material import Material

__all__ = ['RadionuclideProduction', 'radionuclide_production']


@dataclass
class RadionuclideProduction:
    """Production data for a single final state of one reaction.

    Joins the MF=8 identification of a radioactive product (its ZAP, LFS
    level number, and excitation energy) with the energy-dependent data
    evaluated for that state: the MF=9 yield multiplicity and/or the
    MF=10 production cross section. Note that LFS is a level index of
    the product nuclide, not an isomeric-state index; matching a level
    to a metastable state requires comparing :attr:`excitation_energy`
    against decay data.

    Attributes
    ----------
    ZAP : int
        1000*Z + A of the product nuclide
    LFS : int
        Level number of the final state (0 for the ground state)
    QM : float
        Mass-difference Q value in [eV]
    QI : float
        Reaction Q value for this state in [eV]
    ELFS : float or None
        Excitation energy of the final state in [eV] from MF=8, or None
        when the evaluation has no MF=8 subsection for this state
    yields : Tabulated1D or None
        MF=9 yield multiplicity of the reaction cross section as a
        function of incident energy in [eV]
    cross_section : Tabulated1D or None
        MF=10 production cross section in [b] as a function of incident
        energy in [eV]

    """

    ZAP: int
    LFS: int
    QM: float
    QI: float
    ELFS: Optional[float] = None
    yields: Optional[Tabulated1D] = None
    cross_section: Optional[Tabulated1D] = None

    @property
    def excitation_energy(self) -> float:
        """Excitation energy of the final state in [eV], taken from the
        MF=8 ELFS value when present and otherwise from QM minus QI."""
        if self.ELFS is not None:
            return self.ELFS
        return self.QM - self.QI


def radionuclide_production(
        material: Material) -> Dict[int, List[RadionuclideProduction]]:
    """Collect radionuclide production data from MF=8/9/10.

    For every reaction that has an MF=9 or MF=10 section, the final
    states are returned in the order they appear in the evaluation, with
    the MF=9 and MF=10 data for the same (ZAP, LFS) pair merged into one
    :class:`RadionuclideProduction` and the MF=8 excitation energy
    attached when available. The tabulated functions are returned
    exactly as evaluated.

    Parameters
    ----------
    material
        Material to read production data from

    Returns
    -------
    dict
        Mapping of MT numbers to lists of :class:`RadionuclideProduction`

    """
    by_mt: Dict[int, set] = {}
    for mf, mt in material.sections:
        if mf in (9, 10):
            by_mt.setdefault(mt, set()).add(mf)

    result = {}
    for mt in sorted(by_mt):
        # MF=8 links each (ZAP, LFS) pair to an excitation energy
        elfs = {}
        if (8, mt) in material.section_data:
            for subsection in material.section_data[8, mt]['subsections']:
                key = int(subsection['ZAP']), int(subsection['LFS'])
                elfs[key] = float(subsection['ELFS'])

        states: Dict[tuple, RadionuclideProduction] = {}
        ordered = []
        for mf in sorted(by_mt[mt]):
            for level in material.section_data[mf, mt]['levels']:
                key = int(level['IZAP']), int(level['LFS'])
                state = states.get(key)
                if state is None:
                    state = RadionuclideProduction(
                        ZAP=key[0],
                        LFS=key[1],
                        QM=level['QM'],
                        QI=level['QI'],
                        ELFS=elfs.get(key),
                    )
                    states[key] = state
                    ordered.append(state)
                if mf == 9:
                    state.yields = level['Y']
                else:
                    state.cross_section = level['sigma']
        result[mt] = ordered
    return result

Export in src/endf/__init__.py

 from .product import *
+from .radionuclide_production import *
 from .reaction import *

Tests tests/test_radionuclide_production.py

All 5 pass alongside the existing test suite.

from pathlib import Path

import pytest
import endf
from endf import RadionuclideProduction, radionuclide_production


@pytest.fixture
def in115():
    # ENDF/B-VIII.1 In115 evaluation trimmed to MF=1/451, the MF=3
    # sections for MT=4/16/102, and the corresponding MF=8/9/10 sections
    filename = Path(__file__).with_name('n-049_In-115_trimmed.endf')
    return endf.Material(filename)


def test_reactions_found(in115):
    production = radionuclide_production(in115)
    assert sorted(production) == [4, 16, 102]


def test_mf9_yields(in115):
    # In115 capture gives only the In116 first metastable state, with an
    # implicit ground-state share
    (state,) = radionuclide_production(in115)[102]
    assert state.ZAP == 49116
    assert state.LFS == 1
    assert state.ELFS == pytest.approx(127269.7)
    assert state.excitation_energy == pytest.approx(127269.7)
    assert state.cross_section is None
    assert state.yields is not None
    assert state.yields(0.0253) == pytest.approx(0.79)


def test_mf10_cross_section(in115):
    # Inelastic scattering to In115_m1 is given as an MF=10 partial
    # cross section with no MF=9 counterpart
    (state,) = radionuclide_production(in115)[4]
    assert state.ZAP == 49115
    assert state.LFS == 1
    assert state.ELFS == pytest.approx(336.2e3)
    assert state.yields is None
    assert state.cross_section is not None
    assert state.QM == 0.0
    assert state.QI == pytest.approx(-336.2e3)

    (state,) = radionuclide_production(in115)[16]
    assert state.ZAP == 49114
    assert state.LFS == 1
    assert state.ELFS == pytest.approx(190268.2)


def test_excitation_energy_fallback():
    # Without an MF=8 subsection the excitation energy falls back to the
    # difference of the Q values
    state = RadionuclideProduction(ZAP=41093, LFS=1, QM=0.0, QI=-30730.0)
    assert state.ELFS is None
    assert state.excitation_energy == pytest.approx(30730.0)


def test_material_without_data():
    filename = Path(__file__).with_name('n-095_Am_244.endf')
    material = endf.Material(filename)
    assert radionuclide_production(material) == {}

Test fixture tests/n-049_In-115_trimmed.endf

A 46 kB trim of the ENDF/B-VIII.1 In115 evaluation (MF=1/451, MF=3 for MT=4/16/102, and the corresponding MF=8/9/10 sections). It covers the three interesting shapes in one file: MF=9 yields with an implicit ground-state share (capture to In116 m1, Y=0.79), MF=10-only threshold states with no ground counterpart, and excitation energy attachment from MF=8. Regenerate it from the full evaluation with:

import endf

KEEP = [(1, 451), (3, 4), (3, 16), (3, 102),
        (8, 4), (8, 16), (8, 102), (9, 102), (10, 4), (10, 16)]

src = endf.Material('n-049_In-115.endf')  # from the ENDF/B-VIII.1 neutron sublibrary
mat = src.MAT
with open('n-049_In-115.endf') as fh:
    lines = [fh.readline()]                                # TPID
last_mf = None
for mf, mt in KEEP:
    if last_mf is not None and mf != last_mf:
        lines.append(f"{' ' * 66}{mat:4d} 0  0    0\n")    # FEND
    lines.append(src.section_text[mf, mt])
    lines.append(f"{' ' * 66}{mat:4d}{mf:2d}  0    0\n")   # SEND
    last_mf = mf
lines.append(f"{' ' * 66}{mat:4d} 0  0    0\n")            # FEND
lines.append(f"{' ' * 66}{0:4d} 0  0    0\n")              # MEND
lines.append(f"{' ' * 66}{-1:4d} 0  0    0\n")             # TEND
with open('n-049_In-115_trimmed.endf', 'w') as fh:
    fh.writelines(lines)

[shimwell]

Implemented in a Rust port of this reader, in case the shape is useful.

MF=8 identifies each radioactive product and its state, MF=9 gives an energy-dependent yield and MF=10 a production cross section. The three are joined and keyed by reaction, because a consumer wants them together rather than as three parallel lookups. Turning a product's state number into a named nuclide needs decay data for the isomer excitation energies, so that is a separate input rather than something the reader guesses.

Source: crates/endf/src/radionuclide_production.rs.